求解基础偏微分方程的数值方法

偏微分方程的求解非常困难，无论是解析还是数值方法，然而在多维度、连续变化且相互耦合的复杂现象（如流体流动、电磁波传播、金融市场波动等）的建模中又需要偏微分方程。

偏微分方程听起来可能有些抽象，但它研究的其实是一些非常具体的问题：

一根金属棒被加热后，温度怎样随时间变化？
水面上的波纹怎样传播？
空气中的污染物怎样扩散？
一块薄膜受力后怎样振动？
电磁场怎样在空间中传播？

这些问题的共同特点是：我们关心的量不仅随时间变化，也随空间位置变化。

一维热方程

这算是偏微分方程最经典（简单）的例子了。想象一根长度为 2 米的细金属棒。

我们用 $x$ 表示金属棒上的位置： $0\le x\le 2$ 。

用 $t$ 表示时间，用函数 $u(x,t)$ 表示位置 $x$ 在时刻 $t$ 的温度。

如果 $u(1,0)=3$ ，则表示在初始时刻 $t=0$ ，金属棒中点 $x=1$ 的温度为 3。

一维热方程：

\frac{\partial u}{\partial t} = \alpha\frac{\partial^2u}{\partial x^2} + f(x, t)

为简单起见，我们取 $f(t,x)=0$ ，即系统不被加热或冷却。方程简化为：

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

假设热扩散系数 $\alpha = 1$ ，进一步简化为：

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

设置边界条件：

u(0,t)=u(2,t)=0

意味着无论何时，两端的温度都是 0。

设定初始温度曲线函数为：

u(0, x) = 3sin(\frac{\pi}{2}x), (0\le x\le 2)

有限差分法（FEM）

有限差分法核心思想是用有限的差分代替无限的微分，把微积分问题变成代数方程组。这个方法是数值分析和科学计算中最经典、最直观，也是历史最悠久的偏微分方程（PDE）数值求解方法。

程序涉及的依赖：

import matplotlib.pyplot as plt
import numpy as np
from mpl_toolkits import mplot3d
from scipy import sparse

设定条件

创建系统物理约束条件的变量，包括金属棒的长度范围和扩散系数 α 的值：

alpha = 1
x0 = 0
xL = 2

划分区间

使用 N + 1 个点将 x 的范围划分为 N 个等间隔的区间（相当于对空间的划分）：

N = 10
x = np.linspace(x0, xL, N + 1)
h = (xL - x0) / N

划分时间

我们需要设置时间方向上的步长。在这里设置了时间步长 k 和步数（在隐含意义上，我们假设从时间 0 开始）。

k = 0.01
steps = 100
t = np.array([i * k for i in range(steps + 1)])

稳定性

由于这里用的是显式有限差分格式，也叫 FTCS：

Forward Time：时间向前差分
Centered Space：空间中心差分

稳定条件为 $r \le 0.5$ ，其中

r = \frac{ak}{h^2}

所以有：

r = alpha * k / h ** 2
assert r < 0.5, f"Must have r < 0.5, currently r = {r}"

如果是不稳定的，微小的扰动就会像雪崩一样迅速放大，最终导致系统崩溃或计算结果炸掉，比如输出为无穷大或 NaN。

为什么 $r=\dfrac{\alpha k}{h^2}$

一阶差分：

\Delta f(x_0) = f(x_0 + \Delta x) - f(x_0)

二阶差分：

\begin{align*} \Delta^2 f(x_0) &= \Delta f(x_0 + \Delta x) - \Delta f(x_0) \\ &= (f(x_0 + 2 \Delta x) - f(x_0 + \Delta x)) - (f(x_0 + \Delta x) - f(x_0)) \\ &= f(x_0 + 2 \Delta x) - 2f(x_0 + \Delta x) + f(x_0) \end{align*}

对时间导数使用向前差分：

\frac{\partial u}{\partial t} \approx \frac{u_j^{n+1}-u_j^n}{k},

对空间二阶导数使用中心差分：

\frac{\partial^2u}{\partial x^2} \approx \frac{u_{j+1}^n-2u_j^n+u_{j-1}^n}{h^2}.

代入热方程：

\frac{u_j^{n+1}-u_j^n}{k} =\alpha\frac{u_{j+1}^n-2u_j^n+u_{j-1}^n}{h^2}

两边乘以 $k$ ：

u_j^{n+1}-u_j^n = \frac{\alpha k}{h^2} \left( u_{j+1}^n-2u_j^n+u_{j-1}^n \right)

于是定义

r=\frac{\alpha k}{h^2}

公式就可以简写为差分形式：

u_j^{n+1} = r u_{j-1}^n + (1-2r)u_j^n + r u_{j+1}^n

保存有限差分方案的系数

可以构建一个矩阵来保存有限差分方案的系数。

为此，我们使用 scipy.sparse 模块中的 diags 来创建一个稀疏的对角矩阵：

diag = [1, *(1 - 2*r for _ in range(N - 1)), 1]
abv_diag = [0, *(r for _ in range(N - 1))]
blw_diag = [*(r for _ in range(N - 1)), 0]

A = sparse.diags([blw_diag, diag, abv_diag], (-1, 0, 1),
                 shape=(N + 1, N + 1), dtype=np.float64, format="csr")

空白矩阵保存解

u = np.zeros((steps + 1, N + 1), dtype=np.float64)

初始温度曲线

def initial_profile(x):
    return 3 * np.sin(np.pi * x / 2)

u[0, :] = initial_profile(x)

求解

循环执行每一步，通过将矩阵与前一行相乘计算矩阵的下一行：

for i in range(steps):
    u[i + 1, :] = A @ u[i, :]

可视化

X, T = np.meshgrid(x, t)
fig = plt.figure()
ax = fig.add_subplot(projection="3d")
ax.plot_surface(T, X, u, cmap="gray")

ax.set_title("Solution of the heat equation")
ax.set_xlabel("t")
ax.set_ylabel("x")
ax.set_zlabel("u")

MATLAB 代码

% 偏微分方程的数值求解
% 一维热方程
%
% 方程：
%     du/dt = alpha * d^2u/dx^2
%
% 使用显式有限差分格式 FTCS：
%     u_j^(n+1) = r*u_(j-1)^n + (1-2r)*u_j^n + r*u_(j+1)^n

clear;
clc;
close all;

%% 参数设置

alpha = 1;

x0 = 0;
xL = 2;

%% 空间差分

N = 10;                          % 空间区间数量
x = linspace(x0, xL, N + 1);     % N+1 个网格点
h = (xL - x0) / N;               % 空间步长

%% 时间步长

k = 0.01;                        % 时间步长
steps = 100;                     % 时间迭代次数
t = (0:steps) * k;               % 时间网格

%% 稳定性参数

r = alpha * k / h^2;

assert(r < 0.5, ...
    sprintf('Must have r < 0.5, currently r = %.6f', r));

fprintf('h = %.4f\n', h);
fprintf('k = %.4f\n', k);
fprintf('r = %.4f\n', r);

%% 构造有限差分系数矩阵

% 初始化稀疏矩阵
A = sparse(N + 1, N + 1);

% 左边界保持不变
A(1, 1) = 1;

% 内部网格点
for j = 2:N
    A(j, j - 1) = r;
    A(j, j)     = 1 - 2 * r;
    A(j, j + 1) = r;
end

% 右边界保持不变
A(N + 1, N + 1) = 1;

%% 初始化解矩阵

% 行表示时间，列表示空间
u = zeros(steps + 1, N + 1);

%% 初始条件

initial_profile = @(x) 3 * sin(pi * x / 2);

u(1, :) = initial_profile(x);

%% 时间迭代

for n = 1:steps
    % MATLAB 中列向量与矩阵相乘，因此需要进行转置
    u(n + 1, :) = (A * u(n, :)')';
end

%% 绘制三维曲面

[X, T] = meshgrid(x, t);

figure;

surf(T, X, u);

colormap(gray);
shading interp;

title('Solution of the heat equation');
xlabel('t');
ylabel('x');
zlabel('u(x,t)');

grid on;
view(45, 30);

MATLAB 可视化的结果