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+12dxui+1n+1−ui−1n+1=ϵdx2ui+1n+1−2uin+1+ui−1n+1.
数值求解过程
为了读者方便，下面我们先给出 Matlab 差分法解 Burgers 方程的代码实现：
% The solution surface of Burgers Equations
% Variables:
% x-space variable, t-time variable.
% Output:
% Solution surface of u_{t}+uu_x=\epsilon u_{xx}
% 参数设置
epsilon = 0.05;
dx = 5e-2; % 空间步长
dt = 5e-3; % 时间步长
x = -1:dx:1; % 空间网格
t = 0:dt:1; % 时间网格
Nx = length(x);
Nt = length(t);
% 初始条件
u = -sin(pi * x);
% 边界条件函数（这里假设为常数边界条件）
Leftbdry = @(t) 0;
Rightbdry = @(t) 0;
% 初始化解矩阵
U = zeros(Nx, Nt);
U(:, 1) = u; % 初始条件
% 向后欧拉迭代求解
for n = 1:Nt-1
u_n = U(:, n); % 当前时间步的解
u_np1 = u_n; % 下一时间步的解，初始猜测为当前时间步的解
% 迭代求解隐式方程
max_iter = 100;
tol = 1e-6;
for iter = 1:max_iter
% 计算u_x和u_xx（这里使用简单的二阶中心差分，注意边界处理）
u_x = zeros(Nx, 1);
u_xx = zeros(Nx, 1);
% 中心差分
u_x(2:end-1) = (u_np1(3:end) - u_np1(1:end-2)) / (2 * dx);
u_xx(2:end-1) = (u_np1(3:end) - 2 * u_np1(2:end-1) + u_np1(1:end-2)) / (dx^2);
% 边界处理
u_x(1) = (u_np1(2) - u_np1(1)) / dx;
u_x(end) = (u_np1(end) - u_np1(end-1)) / dx;
u_xx(1) = (u_np1(2) - 2 * u_np1(1) + Leftbdry(t(n+1))) / (dx^2);
u_xx(end) = (Rightbdry(t(n+1)) - 2 * u_np1(end) + u_np1(end-1)) / (dx^2);
% 计算u*u_x
uu_x = u_np1 .* u_x;
% 向后欧拉方程
u_np1_new = u_n - dt * (uu_x - epsilon * u_xx);
% 检查收敛性
if norm(u_np1_new - u_np1, inf) < tol
break;
end
u_np1 = u_np1_new;
end
% 更新解矩阵
U(:, n+1) = u_np1;
% 边界条件修正（如果需要）
U(1, n+1) = Leftbdry(t(n+1));
U(end, n+1) = Rightbdry(t(n+1));
end
% 使用 meshgrid 生成网格数据
[T, X] = meshgrid(t, x);
% 绘图显示数值解
figure;
% 将 T, X, U 转换为列向量
% T_vec = T(:);
% X_vec = X(:);
% U_vec = U(:);
surf(T, X, U);hold on;
scatter3(T(:),X(:),U(:),'.' )
xlabel('t')
ylabel('x')
zlabel('u(x,t)')
title('Burgers Equation Solution using Backward Euler Method')
求解结果：
考虑到 Python 用户较多，笔者也实现了 Python 版本的代码供大家参考：
import numpy as np
import matplotlib.pyplot as plt
# 参数设置
epsilon = 0.05
dx = 5e-2 # 空间步长
dt = 5e-3 # 时间步长
x = np.arange(-1, 1 + dx, dx) # 空间网格
t = np.arange(0, 1 + dt, dt) # 时间网格
Nx = len(x)
Nt = len(t)
# 初始条件
u = -np.sin(np.pi * x)
# 边界条件函数（这里假设为常数边界条件）
def Leftbdry(t):
return 0
def Rightbdry(t):
return 0
# 初始化解矩阵
U = np.zeros((Nx, Nt))
U[:, 0] = u # 初始条件
# 向后欧拉迭代求解
for n in range(Nt - 1):
u_n = U[:, n] # 当前时间步的解
u_np1 = u_n # 下一时间步的解，初始猜测为当前时间步的解
# 迭代求解隐式方程
max_iter = 100
tol = 1e-6
for iter in range(max_iter):
# 计算 u_x 和 u_xx（这里使用简单的二阶中心差分，注意边界处理）
u_x = np.zeros(Nx)
u_xx = np.zeros(Nx)
# 中心差分
u_x[1:-1] = (u_np1[2:] - u_np1[:-2]) / (2 * dx)
u_xx[1:-1] = (u_np1[2:] - 2 * u_np1[1:-1] + u_np1[:-2]) / (dx**2)
# 边界处理
u_x[0] = (u_np1[1] - u_np1[0]) / dx
u_x[-1] = (u_np1[-1] - u_np1[-2]) / dx
u_xx[0] = (u_np1[1] - 2 * u_np1[0] + Leftbdry(t[n+1])) / (dx**2)
u_xx[-1] = (Rightbdry(t[n+1]) - 2 * u_np1[-1] + u_np1[-2]) / (dx**2)
# 计算 u * u_x
uu_x = u_np1 * u_x
# 向后欧拉方程
u_np1_new = u_n - dt * (uu_x - epsilon * u_xx)
# 检查收敛性
if np.linalg.norm(u_np1_new - u_np1, np.inf) < tol:
break
u_np1 = u_np1_new
# 更新解矩阵
U[:, n+1] = u_np1
# 边界条件修正（如果需要）
U[0, n+1] = Leftbdry(t[n+1])
U[-1, n+1] = Rightbdry(t[n+1])
# 使用 meshgrid 生成网格数据
T, X = np.meshgrid(t, x)
# 绘图显示数值解
fig = plt.figure(figsize=(12, 12))
# 第一个视角
ax1 = fig.add_subplot(221, projection='3d')
ax1.plot_surface(T, X, U, cmap='viridis', edgecolor='none')
ax1.set_title('View from angle 1')
ax1.set_xlabel('t')
ax1.set_ylabel('x')
ax1.set_zlabel('u(x,t)')
ax1.view_init(elev=30, azim=45)
# 第二个视角
ax2 = fig.add_subplot(222, projection='3d')
ax2.plot_surface(T, X, U, cmap='viridis', edgecolor='none')
ax2.set_title('View from angle 2')
ax2.set_xlabel('t')
ax2.set_ylabel('x')
ax2.set_zlabel('u(x,t)')
ax2.view_init(elev=30, azim=135)
# 第三个视角
ax3 = fig.add_subplot(223, projection='3d')
ax3.plot_surface(T, X, U, cmap='viridis', edgecolor='none')
ax3.set_title('View from angle 3')
ax3.set_xlabel('t')
ax3.set_ylabel('x')
ax3.set_zlabel('u(x,t)')
ax3.view_init(elev=60, azim=45)
# 第四个视角
ax4 = fig.add_subplot(224, projection='3d')
ax4.plot_surface(T, X, U, cmap='viridis', edgecolor='none')
ax4.set_title('View from angle 4')
ax4.set_xlabel('t')
ax4.set_ylabel('x')
ax4.set_zlabel('u(x,t)')
ax4.view_init(elev=60, azim=135)
plt.suptitle('Burgers Equation Solution from Multiple Angles')
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 为整个图像添加总标题，并显示图像。
数值结果如下：
效果不错！
本专栏致力于普及各种偏微分方程的不同数值求解方法，所有文章包含全部可运行代码。欢迎大家支持、关注！
作者 ：计算小屋
个人主页 ： 计算小屋的主页
