使用BE（向后欧拉），FE（向前欧拉），C N方法求解1d扩散方程 Matlab

在数值求解偏微分方程领域，1D扩散方程是一个经典的研究对象。而向前欧拉（FE）、向后欧拉（BE）和克兰克 - 尼科尔森（CN）方法是常用的求解手段。今天咱们就用Matlab来实现这三种方法对1D扩散方程的求解。

1D扩散方程

1D扩散方程的一般形式为：

\[

\frac{\partial u}{\partial t} = D \frac{\partial^2 u}{\partial x^2}

\]

其中 \( u(x, t) \) 是在位置 \( x \) 和时间 \( t \) 的扩散量， \( D \) 是扩散系数。

向前欧拉（FE）方法

向前欧拉方法是一种显式的时间推进方法。在空间和时间上进行离散化后，1D扩散方程的向前欧拉格式为：

\[

ui^{n + 1} = ui^n + \Delta t D \frac{u{i + 1}^n - 2ui^n + u_{i - 1}^n}{\Delta x^2}

\]

这里 \( u_i^n \) 表示在时间步 \( n \) 和空间点 \( i \) 的值，\( \Delta t \) 是时间步长，\( \Delta x \) 是空间步长。

下面是Matlab代码实现：

% 参数设置
L = 1; % 区域长度
T = 0.1; % 总时间
nx = 101; % 空间节点数
nt = 1000; % 时间步数
D = 1; % 扩散系数

dx = L / (nx - 1);
dt = T / nt;
x = linspace(0, L, nx);
t = linspace(0, T, nt + 1);

u = zeros(nx, nt + 1);
u(:, 1) = exp(-100 * (x - 0.5).^2); % 初始条件

for n = 1:nt
    for i = 2:nx - 1
        u(i, n + 1) = u(i, n) + D * dt / dx^2 * (u(i + 1, n) - 2 * u(i, n) + u(i - 1, n));
    end
    % 边界条件
    u(1, n + 1) = 0;
    u(nx, n + 1) = 0;
end

% 绘图
figure;
for n = 1:nt + 1
    plot(x, u(:, n));
    hold on;
end
xlabel('x');
ylabel('u(x, t)');
title('向前欧拉方法求解1D扩散方程');
hold off;

代码分析：首先设定了区域长度、总时间、空间和时间节点数以及扩散系数。通过 linspace 函数生成空间和时间向量。初始化 \( u \) 矩阵，并设置初始条件为一个高斯分布。在双重循环中，按照向前欧拉格式更新 \( u \) 的值，并在每次更新后应用边界条件（这里设为0）。最后通过循环绘图展示不同时间步的解。

向后欧拉（BE）方法

向后欧拉方法是一种隐式的时间推进方法。其格式为：

使用BE（向后欧拉），FE（向前欧拉），C N方法求解1d扩散方程 Matlab

\[

\frac{ui^{n + 1} - ui^n}{\Delta t} = D \frac{u{i + 1}^{n + 1} - 2ui^{n + 1} + u_{i - 1}^{n + 1}}{\Delta x^2}

\]

整理后可以写成矩阵形式 \( A \mathbf{u}^{n + 1} = \mathbf{b} \)，其中 \( A \) 是系数矩阵，\( \mathbf{u}^{n + 1} \) 是未知向量，\( \mathbf{b} \) 是已知向量。

Matlab代码如下：

% 参数设置同FE方法
L = 1;
T = 0.1;
nx = 101;
nt = 1000;
D = 1;

dx = L / (nx - 1);
dt = T / nt;
x = linspace(0, L, nx);
t = linspace(0, T, nt + 1);

u = zeros(nx, nt + 1);
u(:, 1) = exp(-100 * (x - 0.5).^2);

% 构建系数矩阵A
A = zeros(nx, nx);
A(1, 1) = 1;
A(nx, nx) = 1;
for i = 2:nx - 1
    A(i, i - 1) = -D * dt / dx^2;
    A(i, i) = 1 + 2 * D * dt / dx^2;
    A(i, i + 1) = -D * dt / dx^2;
end

for n = 1:nt
    b = u(:, n);
    b(1) = 0; % 边界条件
    b(nx) = 0;
    u(:, n + 1) = A \ b;
end

% 绘图
figure;
for n = 1:nt + 1
    plot(x, u(:, n));
    hold on;
end
xlabel('x');
ylabel('u(x, t)');
title('向后欧拉方法求解1D扩散方程');
hold off;

代码分析：参数设置部分和FE方法类似。重点在于构建系数矩阵 \( A \)，按照向后欧拉格式的系数填充矩阵。在时间推进循环中，每次根据前一时间步的 \( u \) 值构建向量 \( b \)，并应用边界条件，然后通过矩阵除法求解 \( u^{n + 1} \)。绘图部分和FE方法类似，展示不同时间步的解。

克兰克 - 尼科尔森（CN）方法

克兰克 - 尼科尔森方法是一种隐式的、二阶精度的时间推进方法。其格式为：

\[

\frac{ui^{n + 1} - ui^n}{\Delta t} = \frac{D}{2} \left( \frac{u{i + 1}^{n + 1} - 2ui^{n + 1} + u{i - 1}^{n + 1}}{\Delta x^2} + \frac{u{i + 1}^n - 2ui^n + u{i - 1}^n}{\Delta x^2} \right)

\]

同样可以写成矩阵形式求解。

Matlab代码：

% 参数设置同前
L = 1;
T = 0.1;
nx = 101;
nt = 1000;
D = 1;

dx = L / (nx - 1);
dt = T / nt;
x = linspace(0, L, nx);
t = linspace(0, T, nt + 1);

u = zeros(nx, nt + 1);
u(:, 1) = exp(-100 * (x - 0.5).^2);

% 构建系数矩阵A
A = zeros(nx, nx);
A(1, 1) = 1;
A(nx, nx) = 1;
for i = 2:nx - 1
    A(i, i - 1) = -D * dt / (2 * dx^2);
    A(i, i) = 1 + D * dt / dx^2;
    A(i, i + 1) = -D * dt / (2 * dx^2);
end

% 构建矩阵B
B = zeros(nx, nx);
B(1, 1) = 1;
B(nx, nx) = 1;
for i = 2:nx - 1
    B(i, i - 1) = D * dt / (2 * dx^2);
    B(i, i) = 1 - D * dt / dx^2;
    B(i, i + 1) = D * dt / (2 * dx^2);
end

for n = 1:nt
    b = B * u(:, n);
    b(1) = 0; % 边界条件
    b(nx) = 0;
    u(:, n + 1) = A \ b;
end

% 绘图
figure;
for n = 1:nt + 1
    plot(x, u(:, n));
    hold on;
end
xlabel('x');
ylabel('u(x, t)');
title('克兰克 - 尼科尔森方法求解1D扩散方程');
hold off;

代码分析：参数设置依旧相同。这里构建了两个矩阵 \( A \) 和 \( B \)，分别对应克兰克 - 尼科尔森格式中的相关系数。在时间推进循环中，通过矩阵乘法得到向量 \( b \) 并应用边界条件，再通过矩阵除法求解 \( u^{n + 1} \)。绘图部分展示不同时间步的解。

通过这三种方法的Matlab实现，我们可以直观地看到它们对1D扩散方程求解的过程和结果差异。向前欧拉方法简单直观但稳定性条件苛刻，向后欧拉方法稳定性好但计算量相对大些，克兰克 - 尼科尔森方法则在精度和稳定性上有较好的平衡。每种方法都有其适用场景，具体使用哪种方法取决于实际问题的需求。