有限差分法

一、直观理解：用网格上的差商代替导数

有限差分法（Finite Difference Method, FDM）把连续的时空区域切成网格，用邻近网格点的函数值差商来近似偏导数，从而把 PDE 变成每个网格点上的代数方程或递推式。

以温度 $u(x,t)$ 为例：不再追求处处可微的解析表达式，而是计算

$$u_i^n \approx u(x_i,t_n),\qquad x_i=i\Delta x,\ t_n=n\Delta t.$$

空间二阶导数常用中心差：

$$u_{xx}(x_i,t_n)\approx\frac{u_{i+1}^n-2u_i^n+u_{i-1}^n}{\Delta x^2}.$$

时间一阶导数用前差或后差，就得到显式或隐式时间推进。这像用“相邻气象站温差”估计温度梯度，站点越密，估计通常越准——但时间步若与空间步搭配不当，预报也可能数值炸掉。

有限差分是规则区域上最容易上手的 PDE 数值法，也是理解 CFL、稳定性与收敛阶的最佳入口。

二、适用场景与优缺点

类型	说明
适用场景	矩形/立方体等规则区域；热传导、扩散、波动、对流—扩散的教学与快速原型；需要完全掌控格式与边界离散；一维/二维结构网格问题。
核心优点	推导直观、编程直接；精度阶可通过 Taylor 分析清楚给出；显式格式实现极快；便于做网格加密与稳定性实验。
局限与缺点	复杂几何与局部加密不自然；显式格式受稳定性限制；高维时网格点数膨胀；守恒性、单调性需挑选合适格式，不是自动保证。

使用前应确认：

1. 方程类型（抛物/双曲/椭圆）与定解条件完整；
2. 区域是否可用结构网格拟合；
3. 显式还是隐式，以及稳定性条件；
4. 有解析解、稳态或守恒量用于验证。

三、核心原理：从差商到时间推进

1. 常见差商

$$\begin{aligned} u_x &\approx \frac{u_{i+1}-u_i}{\Delta x} \quad\text{（前差，$O(\Delta x)$）},\\ u_x &\approx \frac{u_{i+1}-u_{i-1}}{2\Delta x} \quad\text{（中心，$O(\Delta x^2)$）},\\ u_{xx} &\approx \frac{u_{i+1}-2u_i+u_{i-1}}{\Delta x^2} \quad\text{（中心，$O(\Delta x^2)$）}. \end{aligned}$$

2. 一维热方程显式格式

方程 $u_t=\alpha u_{xx}$。FTCS（Forward Time Central Space）：

$$\frac{u_i^{n+1}-u_i^n}{\Delta t} =\alpha\frac{u_{i+1}^n-2u_i^n+u_{i-1}^n}{\Delta x^2}.$$

令 $r=\alpha\Delta t/\Delta x^2$，则

$$u_i^{n+1}=u_i^n+r\big(u_{i+1}^n-2u_i^n+u_{i-1}^n\big).$$

稳定性（Fourier 分析）要求

$$r\le \frac12.$$

不满足时，数值解通常出现指数增长的网格振荡。

3. 隐式与 Crank–Nicolson

• 全隐式（BTCS）：右端用第 $n+1$ 层，无条件稳定（线性常系数），但每步解三对角组；
• Crank–Nicolson：显式与隐式平均，时间二阶，常用于热方程。

4. 双曲型与 CFL

对对流方程 $u_t+c u_x=0$，显式格式通常要求 CFL 条件

$$\nu=\frac{c\Delta t}{\Delta x}\le 1$$

（具体上限依赖格式）。CFL 的直观含义：数值依赖域必须覆盖物理依赖域。

5. Lax 等价定理（概念）

适定线性问题中，一致性 + 稳定性 $\Rightarrow$ 收敛性。因此只加密网格却不稳定，解不会收敛到真解。

四、关键决策

1. 格式选择
   扩散问题：显式（小规模）/ CN 或隐式（更稳健）；强对流：注意迎风与振荡。

2. 步长比
   显式热方程必须监控 $r=\alpha\Delta t/\Delta x^2$；不要只加密空间而忘记缩小时间步。

3. 边界离散
   Dirichlet 直接赋值；Neumann 常用虚拟点或单侧差，阶数应与内点匹配。

4. 收敛验证
   固定 $r$ 或按理论约束同时加密，比较某时刻 $L^2$/$L^\infty$ 误差。

5. 何时换方法
   复杂区域 → 有限元；强调通量守恒与间断 → 有限体积。

五、MATLAB 示例

%% 一维热方程 u_t = alpha u_xx
% x in [0,1], u(0,t)=u(1,t)=0, u(x,0)=sin(pi x)
% 真解 u=exp(-alpha pi^2 t)*sin(pi x)

alpha = 0.1;
nx = 51;
x = linspace(0, 1, nx);
dx = x(2) - x(1);
rTarget = 0.4;                      % < 0.5
dt = rTarget * dx^2 / alpha;
r = alpha * dt / dx^2;
tEnd = 0.5;
nt = ceil(tEnd / dt);
dt = tEnd / nt;                     % 落到准确终点
r = alpha * dt / dx^2;

u = sin(pi * x);
Uhist = zeros(nt+1, nx);
Uhist(1,:) = u;

for n = 1:nt
    un = u;
    u(2:end-1) = un(2:end-1) + r * (un(3:end) - 2*un(2:end-1) + un(1:end-2));
    u(1) = 0; u(end) = 0;
    Uhist(n+1,:) = u;
end

t = linspace(0, tEnd, nt+1);
uTrue = exp(-alpha * pi^2 * tEnd) * sin(pi * x);
err = max(abs(u - uTrue));
fprintf("r=%.3f, max error at t=%.2f: %.3e\n", r, tEnd, err);

figure;
subplot(1,2,1);
plot(x, u, "o-", x, uTrue, "-"); grid on;
legend("FDM 显式", "解析解"); title("t=tEnd 截面");
subplot(1,2,2);
imagesc(x, t, Uhist); axis xy; colorbar;
xlabel("x"); ylabel("t"); title("u(x,t) 时空图");

%% 网格收敛：加密空间并保持 r 固定
nxs = [21, 41, 81, 161];
errs = zeros(size(nxs));
for k = 1:numel(nxs)
    nxk = nxs(k);
    xk = linspace(0,1,nxk);
    dxk = xk(2)-xk(1);
    dtk = 0.4 * dxk^2 / alpha;
    ntk = ceil(tEnd/dtk); dtk = tEnd/ntk;
    rk = alpha*dtk/dxk^2;
    uk = sin(pi*xk);
    for n = 1:ntk
        un = uk;
        uk(2:end-1) = un(2:end-1) + rk*(un(3:end)-2*un(2:end-1)+un(1:end-2));
        uk(1)=0; uk(end)=0;
    end
    ut = exp(-alpha*pi^2*tEnd)*sin(pi*xk);
    errs(k) = sqrt(dxk * sum((uk-ut).^2));  % 近似 L2
end
disp(table(nxs(:), errs(:), [NaN, errs(1:end-1)./errs(2:end)]', ...
    'VariableNames', {'nx', 'L2err', 'ratio'}));
% 空间二阶时，dx 减半 ratio 约 4

%% 不稳定演示：r > 0.5（短时间即可看到振荡）
% 仅作教学，勿用于正式结果

隐式三对角（Thomas 思路可用矩阵左除）：

% BTCS 一步： (I - r L) u^{n+1} = u^n
% 对内点构造三对角
nx = 51; x = linspace(0,1,nx); dx = x(2)-x(1);
alpha = 0.1; dt = 0.01; r = alpha*dt/dx^2;
e = ones(nx-2,1);
A = spdiags([-r*e, (1+2*r)*e, -r*e], -1:1, nx-2, nx-2);
u = sin(pi*x)';
for n = 1:50
    rhs = u(2:end-1);
    u(2:end-1) = A \ rhs;
    u(1)=0; u(end)=0;
end
plot(x, u, "o-"); grid on; title("隐式 BTCS");

六、如何解释结果

有限差分数值实验的完整叙述应包含：

1. PDE、初边值与物理参数；
2. 网格 $\Delta x,\Delta t$、格式名称与理论阶；
3. 稳定性参数 $r$ 或 CFL 数及其是否满足判据；
4. 与解析解/参考解的误差及加密比；
5. 边界处理方式；
6. 结果的物理解读（扩散变平、峰值衰减率等）；
7. 已知局限（网格取向、数值耗散/色散）。

若误差不随加密下降，先查稳定性和边界实现，再查代码索引。

七、常见误区

1. 只加密空间、时间步不变

显式热格式中 $\Delta x$ 减半会使 $r$ 变为 4 倍，极易跨过稳定阈值。

2. 把数值振荡解释成物理不稳

先检查 $r$、CFL、边界与初值是否兼容。

3. Neumann 边界降阶

内点二阶、边界一阶，整体精度可能被边界拖累。

4. 二维仍用过期显式步长公式

维数增加后稳定条件更严（系数变化），不能照搬一维 $r\le 1/2$ 的同一数值而不推导。

5. 未验证就进入复杂非线性项

先在线性热方程上证明代码收敛，再加反应项或对流项。

八、小结

有限差分法的核心思想可以概括为：

在结构网格上用差商替换导数，把 PDE 化为可递推或可求解的代数系统，并通过稳定性与一致性保证收敛。

它是规则区域 PDE 的第一选择，也是理解所有更高级离散方法的基础。使用时必须同时报告网格、格式、稳定性约束和收敛证据；只有曲线“看起来光滑”远远不够。