差分法求解 Burgers 方程（附完备MATLAB 及 Python代码）

Burgers 方程的数值解及偏差分析

弁言

Burgers 方程是一个非线性偏微分方程，在流体力学、非线性声学和交通流理论中有广泛应用。本文将通过数值方法求解带粘性的 Burgers 方程，并分析其偏差。
方程模子

Burgers 方程的情势为：
                                                    u                               t                                      +                            u                                       u                               x                                      =                            ϵ                                       u                                           x                                  x                                                 ,                                  u_t + u u_x = \epsilon u_{xx},                     ut​+uux​=ϵuxx​,
其中，                                   u                              u                  u 是待求函数，                                   x                              x                  x 是空间变量，                                   t                              t                  t 是时间变量，                                   ϵ                              \epsilon                  ϵ 是黏性系数。
初始条件和界限条件

为了求解方程，我们需要指定初始条件和界限条件。在本文中，我们选择如下初始条件：
                                         u                            (                            x                            ,                            0                            )                            =                            −                            sin                            ⁡                            (                            π                            x                            )                            ,                                  u(x, 0) = -\sin(\pi x),                     u(x,0)=−sin(πx),
界限条件设定为常数界限条件：
                                         u                            (                            −                            1                            ,                            t                            )                            =                            0                            ,                                     u                            (                            1                            ,                            t                            )                            =                            0.                                  u(-1, t) = 0, \quad u(1, t) = 0.                     u(−1,t)=0,u(1,t)=0.
数值方法

我们使用向后欧拉法进行时间离散，并使用中央差分法进行空间离散。时间步长为                                    d                         t                              dt                  dt，空间步长为                                    d                         x                              dx                  dx。
时间离散

向后欧拉法的时间离散情势为：
                                                                             u                                                   n                                        +                                        1                                                           −                                               u                                     n                                                                  d                                  t                                                 +                                       u                                           n                                  +                                  1                                                                        ∂                                               u                                                   n                                        +                                        1                                                                               ∂                                  x                                                 =                            ϵ                                                                ∂                                     2                                                           u                                                   n                                        +                                        1                                                                               ∂                                               x                                     2                                                                   \frac{u^{n+1} - u^n}{dt} + u^{n+1} \frac{\partial u^{n+1}}{\partial x} = \epsilon \frac{\partial^2 u^{n+1}}{\partial x^2}                     dtun+1−un​+un+1∂x∂un+1​=ϵ∂x2∂2un+1​
空间离散

中央差分法用于空间离散，                                             u                            x                                       u_x                  ux​ 和                                              u                                       x                               x                                                 u_{xx}                  uxx​ 的差分格式为：
                                                                ∂                                  u                                                      ∂                                  x                                                 ≈                                                                u                                                   i                                        +                                        1                                                           −                                               u                                                   i                                        −                                        1                                                                               2                                  d                                  x                                                 .                                  \frac{\partial u}{\partial x} \approx \frac{u_{i+1} - u_{i-1}}{2 dx}.                     ∂x∂u​≈2dxui+1​−ui−1​​.
                                                                             ∂                                     2                                              u                                                      ∂                                               x                                     2                                                             ≈                                                                u                                                   i                                        +                                        1                                                           −                                  2                                               u                                     i                                              +                                               u                                                   i                                        −                                        1                                                                               d                                               x                                     2                                                             .                                  \frac{\partial^2 u}{\partial x^2} \approx \frac{u_{i+1} - 2u_i + u_{i-1}}{dx^2}.                     ∂x2∂2u​≈dx2ui+1​−2ui​+ui−1​​.
差分格式

综合时间离散和空间离散，得到差分格式：
                                                                             u                                     i                                                   n                                        +                                        1                                                           −                                               u                                     i                                     n                                                                  d                                  t                                                 +                                       u                               i                                           n                                  +                                  1                                                                                     u                                                   i                                        +                                        1                                                                n                                        +                                        1                                                           −                                               u                                                   i                                        −                                        1                                                                n                                        +                                        1                                                                               2                                  d                                  x                                                 =                            ϵ                                                                u                                                   i                                        +                                        1                                                                n                                        +                                        1                                                           −                                  2                                               u                                     i                                                   n                                        +                                        1                                                           +                                               u                                                   i                                        −                                        1                                                                n                                        +                                        1                                                                               d                                               x                                     2                                                             .                                  \frac{u_i^{n+1} - u_i^n}{dt} + u_i^{n+1} \frac{u_{i+1}^{n+1} - u_{i-1}^{n+1}}{2 dx} = \epsilon \frac{u_{i+1}^{n+1} - 2u_i^{n+1} + u_{i-1}^{n+1}}{dx^2}.                     dtuin+1​−uin​​+uin+1​2dxui+1n+1​−ui−1n+1​​=ϵdx2ui+1n+1​−2uin+1​+ui−1n+1​​.
数值求解过程

为了读者方便，下面我们先给出 Matlab 差分法解 Burgers 方程的代码实现：

1. % The solution surface of Burgers Equations
2. % Variables:
3. % x-space variable, t-time variable.
4. % Output:
5. % Solution surface of u_{t}+uu_x=\epsilon u_{xx}
7. % 参数设置  
8. epsilon = 0.05;  
9. dx = 5e-2; % 空间步长  
10. dt = 5e-3; % 时间步长  
11. x = -1:dx:1; % 空间网格  
12. t = 0:dt:1;  % 时间网格  
13. Nx = length(x);  
14. Nt = length(t);  
16. % 初始条件  
17. u = -sin(pi * x);  
19. % 边界条件函数（这里假设为常数边界条件）  
20. Leftbdry = @(t) 0;  
21. Rightbdry = @(t) 0;  
23. % 初始化解矩阵  
24. U = zeros(Nx, Nt);  
25. U(:, 1) = u; % 初始条件  
27. % 向后欧拉迭代求解  
28. for n = 1:Nt-1  
29.     u_n = U(:, n); % 当前时间步的解  
30.     u_np1 = u_n;   % 下一时间步的解，初始猜测为当前时间步的解  
31.    
32.     % 迭代求解隐式方程  
33.     max_iter = 100;  
34.     tol = 1e-6;  
35.     for iter = 1:max_iter  
36.         % 计算u_x和u_xx（这里使用简单的二阶中心差分，注意边界处理）  
37.         u_x = zeros(Nx, 1);
38.         u_xx = zeros(Nx, 1);
39.         
40.         % 中心差分
41.         u_x(2:end-1) = (u_np1(3:end) - u_np1(1:end-2)) / (2 * dx);
42.         u_xx(2:end-1) = (u_np1(3:end) - 2 * u_np1(2:end-1) + u_np1(1:end-2)) / (dx^2);
43.         
44.         % 边界处理
45.         u_x(1) = (u_np1(2) - u_np1(1)) / dx;
46.         u_x(end) = (u_np1(end) - u_np1(end-1)) / dx;
47.         
48.         u_xx(1) = (u_np1(2) - 2 * u_np1(1) + Leftbdry(t(n+1))) / (dx^2);
49.         u_xx(end) = (Rightbdry(t(n+1)) - 2 * u_np1(end) + u_np1(end-1)) / (dx^2);
51.         % 计算u*u_x  
52.         uu_x = u_np1 .* u_x;  
54.         % 向后欧拉方程  
55.         u_np1_new = u_n - dt * (uu_x - epsilon * u_xx);  
57.         % 检查收敛性  
58.         if norm(u_np1_new - u_np1, inf) < tol  
59.             break;  
60.         end  
61.         u_np1 = u_np1_new;  
62.     end  
64.     % 更新解矩阵  
65.     U(:, n+1) = u_np1;  
67.     % 边界条件修正（如果需要）  
68.     U(1, n+1) = Leftbdry(t(n+1));  
69.     U(end, n+1) = Rightbdry(t(n+1));  
70. end  
72. % 使用 meshgrid 生成网格数据
73. [T, X] = meshgrid(t, x);
75.   
76. % 绘图显示数值解
77. figure;
78. % 将 T, X, U 转换为列向量
79. % T_vec = T(:);
80. % X_vec = X(:);
81. % U_vec = U(:);
82. surf(T, X, U);hold on;
83. scatter3(T(:),X(:),U(:),'.' )
84. xlabel('t')  
85. ylabel('x')  
86. zlabel('u(x,t)')  
87. title('Burgers Equation Solution using Backward Euler Method')

复制代码

求解效果：

思量到 Python 用户较多，笔者也实现了 Python 版本的代码供各人参考：

1. import numpy as np
2. import matplotlib.pyplot as plt
4. # 参数设置
5. epsilon = 0.05
6. dx = 5e-2  # 空间步长
7. dt = 5e-3  # 时间步长
8. x = np.arange(-1, 1 + dx, dx)  # 空间网格
9. t = np.arange(0, 1 + dt, dt)  # 时间网格
10. Nx = len(x)
11. Nt = len(t)
13. # 初始条件
14. u = -np.sin(np.pi * x)
16. # 边界条件函数（这里假设为常数边界条件）
17. def Leftbdry(t):
18.     return 0
20. def Rightbdry(t):
21.     return 0
23. # 初始化解矩阵
24. U = np.zeros((Nx, Nt))
25. U[:, 0] = u  # 初始条件
27. # 向后欧拉迭代求解
28. for n in range(Nt - 1):
29.     u_n = U[:, n]  # 当前时间步的解
30.     u_np1 = u_n  # 下一时间步的解，初始猜测为当前时间步的解
31.    
32.     # 迭代求解隐式方程
33.     max_iter = 100
34.     tol = 1e-6
35.     for iter in range(max_iter):
36.         # 计算 u_x 和 u_xx（这里使用简单的二阶中心差分，注意边界处理）
37.         u_x = np.zeros(Nx)
38.         u_xx = np.zeros(Nx)
39.         
40.         # 中心差分
41.         u_x[1:-1] = (u_np1[2:] - u_np1[:-2]) / (2 * dx)
42.         u_xx[1:-1] = (u_np1[2:] - 2 * u_np1[1:-1] + u_np1[:-2]) / (dx**2)
43.         
44.         # 边界处理
45.         u_x[0] = (u_np1[1] - u_np1[0]) / dx
46.         u_x[-1] = (u_np1[-1] - u_np1[-2]) / dx
47.         
48.         u_xx[0] = (u_np1[1] - 2 * u_np1[0] + Leftbdry(t[n+1])) / (dx**2)
49.         u_xx[-1] = (Rightbdry(t[n+1]) - 2 * u_np1[-1] + u_np1[-2]) / (dx**2)
51.         # 计算 u * u_x
52.         uu_x = u_np1 * u_x
54.         # 向后欧拉方程
55.         u_np1_new = u_n - dt * (uu_x - epsilon * u_xx)
57.         # 检查收敛性
58.         if np.linalg.norm(u_np1_new - u_np1, np.inf) < tol:
59.             break
60.         u_np1 = u_np1_new
62.     # 更新解矩阵
63.     U[:, n+1] = u_np1
65.     # 边界条件修正（如果需要）
66.     U[0, n+1] = Leftbdry(t[n+1])
67.     U[-1, n+1] = Rightbdry(t[n+1])
69. # 使用 meshgrid 生成网格数据
70. T, X = np.meshgrid(t, x)
72. # 绘图显示数值解
73. fig = plt.figure(figsize=(12, 12))
75. # 第一个视角
76. ax1 = fig.add_subplot(221, projection='3d')
77. ax1.plot_surface(T, X, U, cmap='viridis', edgecolor='none')
78. ax1.set_title('View from angle 1')
79. ax1.set_xlabel('t')
80. ax1.set_ylabel('x')
81. ax1.set_zlabel('u(x,t)')
82. ax1.view_init(elev=30, azim=45)
84. # 第二个视角
85. ax2 = fig.add_subplot(222, projection='3d')
86. ax2.plot_surface(T, X, U, cmap='viridis', edgecolor='none')
87. ax2.set_title('View from angle 2')
88. ax2.set_xlabel('t')
89. ax2.set_ylabel('x')
90. ax2.set_zlabel('u(x,t)')
91. ax2.view_init(elev=30, azim=135)
93. # 第三个视角
94. ax3 = fig.add_subplot(223, projection='3d')
95. ax3.plot_surface(T, X, U, cmap='viridis', edgecolor='none')
96. ax3.set_title('View from angle 3')
97. ax3.set_xlabel('t')
98. ax3.set_ylabel('x')
99. ax3.set_zlabel('u(x,t)')
100. ax3.view_init(elev=60, azim=45)
102. # 第四个视角
103. ax4 = fig.add_subplot(224, projection='3d')
104. ax4.plot_surface(T, X, U, cmap='viridis', edgecolor='none')
105. ax4.set_title('View from angle 4')
106. ax4.set_xlabel('t')
107. ax4.set_ylabel('x')
108. ax4.set_zlabel('u(x,t)')
109. ax4.view_init(elev=60, azim=135)
111. plt.suptitle('Burgers Equation Solution from Multiple Angles')
112. plt.show()

复制代码

详细阐明：

• 参数设置：界说 epsilon、dx、dt、x 和 t.
  
• 初始条件：设置初始条件为                                         u                            (                            x                            ,                            0                            )                            =                            −                            sin                            ⁡                            (                            π                            x                            )                                  u(x, 0) = -\sin(\pi x)                     u(x,0)=−sin(πx).
  
• 界限条件：界说常数界限条件函数 Leftbdry 和 Rightbdry.
  
• 初始化解矩阵：创建 U 矩阵，并设置初始条件。
  
• 向后欧拉迭代求解：
  对每个时间步，初始推测下一时间步的解。
  使用迭代方法求解隐式方程，盘算 u_x 和 u_xx，并处置惩罚界限条件。
  盘算 uu_x 和更新 u_np1.
  查抄收敛性，如果满足收敛条件则跳出循环。
  
• 更新解矩阵：将 u_np1 的值存储到解矩阵中，并处置惩罚界限条件。
  
• 天生网格数据：使用 meshgrid 天生 T 和 X.
  
• 绘图：
  创建一个包罗四个子图的 3D 图像，每个子图展示从不同视角观察的解。设置各个子图的坐标轴标签、标题和视角。
  使用 plt.suptitle 为整个图像添加总标题，并表现图像。
  

数值效果如下：

效果不错！

---

本专栏致力于普及各种偏微分方程的不同数值求解方法，全部文章包罗全部可运行代码。欢迎各人支持、关注！
作者 ：盘算小屋
个人主页 ： 盘算小屋的主页

免责声明：如果侵犯了您的权益，请联系站长，我们会及时删除侵权内容，谢谢合作！更多信息从访问主页：qidao123.com:ToB企服之家，中国第一个企服评测及商务社交产业平台。