用BE、FE和CN方法求解1D扩散方程的Matlab实现
使用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扩散方程求解的过程和结果差异。向前欧拉方法简单直观但稳定性条件苛刻，向后欧拉方法稳定性好但计算量相对大些，克兰克 - 尼科尔森方法则在精度和稳定性上有较好的平衡。每种方法都有其适用场景，具体使用哪种方法取决于实际问题的需求。