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

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

引言

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

方程模型

Burgers 方程的形式为：
$$u_t+u{u}_x=ϵu_{xx},$$
其中，$u$ 是待求函数，$x$ 是空间变量，$t$ 是时间变量，$ϵ$ 是黏性系数。

初始条件和边界条件

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

数值方法

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

时间离散

向后欧拉法的时间离散形式为：
$$\frac{u^{n+1}-u^n}{dt}+u^{n+1}\frac{\partial u^{n+1}}{\partial x}=ϵ\frac{{\partial }^2{u}^{n+1}}{\partial x^2}$$

空间离散

中心差分法用于空间离散，$u_x$ 和 $u_{xx}$ 的差分格式为：
$$\frac{\partial u}{\partial x}\approx \frac{u_{i+1}-u_{i-1}}{2dx}.$$
$$\frac{{\partial }^{2}u}{\partial x^2}\approx \frac{u_{i+1}-2u_i+u_{i-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}}{2dx}=ϵ\frac{u_{i+1}^{n+1}-2u_i^{n+1}+u_{i-1}^{n+1}}{d{x}^2}.$$

数值求解过程

为了读者方便，下面我们先给出 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 (\pi x)$.
• 边界条件：定义常数边界条件函数 Leftbdry 和 Rightbdry.
• 初始化解矩阵：创建 U 矩阵，并设置初始条件。
• 向后欧拉迭代求解：
  对每个时间步，初始猜测下一时间步的解。
  使用迭代方法求解隐式方程，计算 u_x 和 u_xx，并处理边界条件。
  计算 uu_x 和更新 u_np1.
  检查收敛性，如果满足收敛条件则跳出循环。
• 更新解矩阵：将 u_np1 的值存储到解矩阵中，并处理边界条件。
• 生成网格数据：使用 meshgrid 生成 T 和 X.
• 绘图：
  创建一个包含四个子图的 3D 图像，每个子图展示从不同视角观察的解。设置各个子图的坐标轴标签、标题和视角。
  使用 plt.suptitle 为整个图像添加总标题，并显示图像。

数值结果如下：

效果不错！

---

本专栏致力于普及各种偏微分方程的不同数值求解方法，所有文章包含全部可运行代码。欢迎大家支持、关注！

作者 ：计算小屋
个人主页 ： 计算小屋的主页