基础的二维微分方程的数值解法

在求解基础偏微分方程的数值方法我们研究的是一根细金属棒上的温度变化。空间中只有一个方向，所以温度可以写成 $u(x,t)$ 。

如果把金属棒换成一块薄板，温度不仅随横向位置 $x$ 变化，还会随纵向位置 $y$ 变化。这时温度场需要写成：

u(x,y,t).

从一维推广到二维并不会改变有限差分法的基本思想。我们仍然使用离散网格代替连续空间，只是原来的一排网格点变成了一张网格，每个内部节点也从两个相邻节点变成了四个相邻节点。

本文考虑一块边长为 2 的正方形薄板，并分别使用 Python 和 MATLAB 模拟热量在平面上的扩散。

二维热方程

设正方形薄板所占区域为：

\Omega=[0,2]\times[0,2].

用 $u(x,y,t)$ 表示薄板上位置 $(x,y)$ 在时刻 $t$ 的温度。没有外部热源时，二维热方程为：

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

其中 $\alpha>0$ 是材料的热扩散系数。括号中的部分也常写作拉普拉斯算子：

\nabla^2u = \frac{\partial^2u}{\partial x^2} + \frac{\partial^2u}{\partial y^2}.

因此方程也可以简写为：

u_t=\alpha\nabla^2u.

为简单起见，取热扩散系数 $\alpha=1$ 。

边界条件

假设薄板四条边始终保持为 0 度，即采用齐次 Dirichlet 边界条件：

\begin{aligned} u(0,y,t)&=0, &u(2,y,t)&=0,\\ u(x,0,t)&=0, &u(x,2,t)&=0. \end{aligned}

初始条件

设初始温度分布为：

u(x,y,0) = 3\sin\left(\frac{\pi x}{2}\right) \sin\left(\frac{\pi y}{2}\right).

薄板中心 $(1,1)$ 的初始温度最高，为 3；越靠近边缘，初始温度越低，并在四条边上恰好等于 0。

二维有限差分格式

分别在 $x$ 和 $y$ 方向划分网格：

x_i=i\Delta x,\qquad y_j=j\Delta y,

时间方向上的离散时刻为：

t_n=n\Delta t.

记

u_{i,j}^n\approx u(x_i,y_j,t_n).

也就是说，上标 $n$ 表示时间层，下标 $i,j$ 表示平面网格中的位置。

时间向前差分

在时间方向使用向前差分：

\frac{\partial u}{\partial t} \approx \frac{u_{i,j}^{n+1}-u_{i,j}^n}{\Delta t}.

空间中心差分

在两个空间方向都使用中心差分：

\frac{\partial^2u}{\partial x^2} \approx \frac{u_{i+1,j}^n-2u_{i,j}^n+u_{i-1,j}^n}{\Delta x^2},

\frac{\partial^2u}{\partial y^2} \approx \frac{u_{i,j+1}^n-2u_{i,j}^n+u_{i,j-1}^n}{\Delta y^2}.

将它们代入二维热方程，定义：

r_x=\frac{\alpha\Delta t}{\Delta x^2}, \qquad r_y=\frac{\alpha\Delta t}{\Delta y^2},

便得到显式有限差分格式：

\begin{aligned} u_{i,j}^{n+1} = {} &(1-2r_x-2r_y)u_{i,j}^n\\ &+r_x\left(u_{i-1,j}^n+u_{i+1,j}^n\right)\\ &+r_y\left(u_{i,j-1}^n+u_{i,j+1}^n\right). \end{aligned}

计算下一时刻某个节点的温度，只需要当前节点及其上、下、左、右四个相邻节点的温度。这种结构通常称为五点差分模板。

当两个方向使用相同的空间步长 $h=\Delta x=\Delta y$ 时，令

r=\frac{\alpha\Delta t}{h^2},

格式可以简化为：

u_{i,j}^{n+1} = (1-4r)u_{i,j}^n +r\left( u_{i-1,j}^n+u_{i+1,j}^n +u_{i,j-1}^n+u_{i,j+1}^n \right).

稳定性条件

二维显式格式的稳定条件为：

r_x+r_y\le\frac12.

如果 $\Delta x=\Delta y=h$ ，则 $r_x=r_y=r$ ，所以：

r\le\frac14.

这个条件也可以从加权平均的角度理解。当 $r\le 1/4$ 时，中心节点的系数 $1-4r$ 和四个相邻节点的系数 $r$ 都是非负数，而且所有系数之和为 1。新温度是旧温度的加权平均，不会凭空产生异常的极大值或极小值。

与一维问题相比，二维节点拥有更多邻居，因此相同空间步长下允许使用的时间步长更小。程序中应始终检查：

rx = alpha * dt / dx**2
ry = alpha * dt / dy**2
assert rx + ry <= 0.5

使用 Python 求解

首先导入依赖：

import matplotlib.pyplot as plt
import numpy as np

设置空间和时间网格

alpha = 1.0

x0, xL = 0.0, 2.0
y0, yL = 0.0, 2.0

Nx = 40
Ny = 40

x = np.linspace(x0, xL, Nx + 1)
y = np.linspace(y0, yL, Ny + 1)

dx = (xL - x0) / Nx
dy = (yL - y0) / Ny

dt = 0.0005
steps = 1000
t = np.arange(steps + 1) * dt

rx = alpha * dt / dx**2
ry = alpha * dt / dy**2

assert rx + ry <= 0.5, (
    f"Unstable: rx + ry = {rx + ry:.4f} > 0.5"
)

print(f"dx = {dx:.4f}, dy = {dy:.4f}")
print(f"dt = {dt:.6f}, rx + ry = {rx + ry:.4f}")

这里 $\Delta x=\Delta y=0.05$ ，所以 $r_x=r_y=0.2$ ，满足：

r_x+r_y=0.4\le0.5.

设置初始条件

meshgrid 可以将两个一维坐标数组组合成二维坐标网格：

X, Y = np.meshgrid(x, y, indexing="ij")

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

u = initial_profile(X, Y)

数组 u[i, j] 保存网格点 $(x_i,y_j)$ 的温度。

为了绘制不同时刻的结果，可以只保存需要观察的时间切片，而不必保留每一步的完整二维数组：

snapshot_steps = {0, 100, 400, 1000}
snapshots = {0: u.copy()}

时间迭代

只更新内部节点，边界节点始终保持为 0：

for n in range(steps):
    u_new = u.copy()

    u_new[1:-1, 1:-1] = (
        (1 - 2 * rx - 2 * ry) * u[1:-1, 1:-1]
        + rx * (u[2:, 1:-1] + u[:-2, 1:-1])
        + ry * (u[1:-1, 2:] + u[1:-1, :-2])
    )

    u_new[0, :] = 0
    u_new[-1, :] = 0
    u_new[:, 0] = 0
    u_new[:, -1] = 0

    u = u_new

    if n + 1 in snapshot_steps:
        snapshots[n + 1] = u.copy()

NumPy 的切片运算同时更新所有内部节点，代码与五点差分公式一一对应，也避免了两层 Python 循环。

可视化温度场

使用等高填色图展示四个时刻的温度分布：

fig, axes = plt.subplots(2, 2, figsize=(10, 8))

for ax, step in zip(axes.flat, sorted(snapshot_steps)):
    contour = ax.contourf(
        X,
        Y,
        snapshots[step],
        levels=30,
        cmap="hot",
        vmin=0,
        vmax=3,
    )
    ax.set_title(f"t = {step * dt:.2f}")
    ax.set_xlabel("x")
    ax.set_ylabel("y")
    ax.set_aspect("equal")

fig.colorbar(contour, ax=axes, label="u(x, y, t)")
plt.show()

随着时间增加，温度分布仍然关于薄板中心对称，但最高温度不断下降，热量逐渐从四条低温边界散失。

也可以使用三维曲面展示最终时刻：

fig = plt.figure(figsize=(8, 6))
ax = fig.add_subplot(projection="3d")
ax.plot_surface(X, Y, u, cmap="hot")

ax.set_title(f"Solution at t = {steps * dt:.2f}")
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_zlabel("u")

plt.show()

使用解析解验证

这组初始条件不仅形状简单，而且拥有解析解：

u(x,y,t) = 3e^{-\frac{\pi^2}{2}t} \sin\left(\frac{\pi x}{2}\right) \sin\left(\frac{\pi y}{2}\right).

二维情况下， $x$ 和 $y$ 两个方向的扩散作用相加，因此指数衰减率是一维问题的两倍。温度场的空间形状保持不变，振幅则随时间指数衰减。

可以在最终时刻计算最大绝对误差：

final_time = steps * dt

u_exact = (
    3
    * np.exp(-(np.pi**2 / 2) * final_time)
    * np.sin(np.pi * X / 2)
    * np.sin(np.pi * Y / 2)
)

max_error = np.max(np.abs(u - u_exact))
print(f"Maximum error at t={final_time:.2f}: {max_error:.6e}")

还可以绘制误差场：

plt.figure(figsize=(6, 5))
plt.contourf(X, Y, np.abs(u - u_exact), levels=30, cmap="viridis")
plt.colorbar(label="absolute error")
plt.xlabel("x")
plt.ylabel("y")
plt.title(f"Absolute error at t = {final_time:.2f}")
plt.axis("equal")
plt.show()

显式 FTCS 格式在时间方向具有一阶精度，在两个空间方向具有二阶精度，其截断误差可以写为：

O(\Delta t)+O(\Delta x^2)+O(\Delta y^2).

MATLAB 代码

下面给出对应的完整 MATLAB 实现：

% 二维热方程的显式有限差分解法
%
% 方程：
%   du/dt = alpha * (d^2u/dx^2 + d^2u/dy^2)
%
% 稳定条件：
%   rx + ry <= 1/2

clear;
clc;
close all;

%% 参数设置

alpha = 1.0;

x0 = 0.0;
xL = 2.0;
y0 = 0.0;
yL = 2.0;

Nx = 40;
Ny = 40;

x = linspace(x0, xL, Nx + 1);
y = linspace(y0, yL, Ny + 1);

dx = (xL - x0) / Nx;
dy = (yL - y0) / Ny;

dt = 0.0005;
steps = 1000;

rx = alpha * dt / dx^2;
ry = alpha * dt / dy^2;

assert(rx + ry <= 0.5, ...
    'Unstable: rx + ry = %.4f > 0.5', rx + ry);

fprintf('dx = %.4f, dy = %.4f\n', dx, dy);
fprintf('dt = %.6f, rx + ry = %.4f\n', dt, rx + ry);

%% 初始条件

[X, Y] = ndgrid(x, y);

u = 3 .* sin(pi .* X ./ 2) .* sin(pi .* Y ./ 2);

snapshotSteps = [0, 100, 400, 1000];
snapshots = zeros(Nx + 1, Ny + 1, length(snapshotSteps));
snapshots(:, :, 1) = u;

%% 时间迭代

for n = 1:steps
    uNew = u;

    uNew(2:end-1, 2:end-1) = ...
        (1 - 2 * rx - 2 * ry) .* u(2:end-1, 2:end-1) ...
        + rx .* (u(3:end, 2:end-1) + u(1:end-2, 2:end-1)) ...
        + ry .* (u(2:end-1, 3:end) + u(2:end-1, 1:end-2));

    % 四条边保持为 0
    uNew(1, :) = 0;
    uNew(end, :) = 0;
    uNew(:, 1) = 0;
    uNew(:, end) = 0;

    u = uNew;

    index = find(snapshotSteps == n, 1);
    if ~isempty(index)
        snapshots(:, :, index) = u;
    end
end

%% 绘制温度场

figure;
tiledlayout(2, 2);

for k = 1:length(snapshotSteps)
    nexttile;
    contourf(X, Y, snapshots(:, :, k), 30, 'LineColor', 'none');
    axis equal tight;
    clim([0, 3]);
    colorbar;
    xlabel('x');
    ylabel('y');
    title(sprintf('t = %.2f', snapshotSteps(k) * dt));
end

colormap(hot);

%% 与解析解比较

finalTime = steps * dt;

uExact = 3 .* exp(-(pi^2 / 2) * finalTime) ...
    .* sin(pi .* X ./ 2) .* sin(pi .* Y ./ 2);

maxError = max(abs(u - uExact), [], 'all');
fprintf('Maximum error at t=%.2f: %.6e\n', finalTime, maxError);

figure;
surf(X, Y, u);
shading interp;
colormap(hot);
colorbar;
xlabel('x');
ylabel('y');
zlabel('u(x,y,t)');
title(sprintf('Numerical solution at t = %.2f', finalTime));
view(45, 30);

MATLAB 可视化的结果

由于现在多了一个维度，总共有 4 个维度， $(x,y,u,t)$ 无法在三维上完成完整的可视化。可以通过设定不同的时间 t 感受温度变化：

不同时间的温度变化

从一维到二维发生了什么

从算法上看，一维与二维热方程的主要区别可以概括为：

项目	一维	二维
空间区域	线段	平面区域
解的表示	$u(x,t)$	$u(x,y,t)$
内部节点的邻居	左、右	上、下、左、右
差分模板	三点	五点
等距网格稳定条件	$r\le 1/2$	$r\le 1/4$
单个时间层的数据量	$O(N)$	$O(N^2)$

二维问题不仅需要更多网格点，而且稳定条件更加严格。当空间步长缩小一半时，为保持稳定，时间步长最多只能取原来的四分之一；与此同时，平面上的网格点数量大约变成原来的四倍。因此，高分辨率二维模拟的计算量会迅速增加。

总结

二维热方程仍然遵循与一维问题相同的离散思路：

将连续平面划分成有限个网格点；
使用向前差分近似时间导数；
使用中心差分近似两个方向的二阶空间导数；
根据五点差分模板逐层推进温度场；
检查 $r_x+r_y\le1/2$ ，确保显式格式稳定；
使用解析解和误差场验证程序。

显式有限差分法直观、容易实现，非常适合用来理解二维扩散过程。但它的时间步长受到稳定条件限制。对于更细的网格或更长时间的模拟，可以进一步考虑后向欧拉法、Crank–Nicolson 方法，或者交替方向隐式方法（ADI）。