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 0def 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 = 100tol = 1e-6for 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]) / dxu_x[-1] = (u_np1[-1] - u_np1[-2]) / dxu_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_xuu_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:breaku_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 为整个图像添加总标题，并显示图像。
数值结果如下：
效果不错！
本专栏致力于普及各种偏微分方程的不同数值求解方法，所有文章包含全部可运行代码。欢迎大家支持、关注！
作者 ：计算小屋
个人主页 ： 计算小屋的主页
