二维抛物方程的初边值问题的数值解法（附Matlab程序，六点对称格式，Du Fort-Frankel 格式，ADI格式，LOD格式）

二维抛物方程的初边值问题的数值解法有很多种，下面介绍几种比较常用的方法。

1. 六点对称格式

对于二维抛物方程$u_t = a(u_{xx} + u_{yy}) + f(x,y,t)$，可以采用六点对称格式进行数值求解。该格式的离散方程为：

$$
\frac{u_{i,j}^{n+1}-u_{i,j}^n}{\Delta t} = a\frac{1}{h^2}(u_{i+1,j}^n-4u_{i,j}^n+u_{i-1,j}^n+u_{i,j+1}^n+u_{i,j-1}^n)+f_{i,j}^n
$$

其中，$i,j$为空间网格点的下标，$n$为时间层的下标，$\Delta t$为时间步长，$h$为空间步长。

2. Du Fort-Frankel 格式

Du Fort-Frankel 格式是隐式差分格式，可以有效避免数值解的不稳定性。离散方程为：

$$
\frac{u_{i,j}^{n+1}-u_{i,j}^{n-1}}{2\Delta t} = a\frac{1}{h^2}(u_{i+1,j}^n-2u_{i,j}^n+u_{i-1,j}^n+u_{i,j+1}^n-2u_{i,j}^n+u_{i,j-1}^n)+f_{i,j}^n
$$

3. ADI格式

ADI格式是一种交替方向显式差分格式，可以降低计算复杂度。具体地，先沿$x$方向进行隐式差分，再沿$y$方向进行隐式差分，离散方程为：

$$
\begin{aligned}
\frac{u_{i,j}^{n+1/2}-u_{i,j}^n}{\Delta t/2} &= a\frac{1}{h^2}(u_{i+1,j}^{n+1/2}-2u_{i,j}^{n+1/2}+u_{i-1,j}^{n+1/2})+f_{i,j}^n \\
\frac{u_{i,j}^{n+1}-u_{i,j}^{n+1/2}}{\Delta t/2} &= a\frac{1}{h^2}(u_{i,j+1}^{n+1}-2u_{i,j}^{n+1}+u_{i,j-1}^{n+1})+f_{i,j}^{n+1/2}
\end{aligned}
$$

4. LOD格式

LOD格式是一种基于有限元方法的数值解法，可以适用于非规则网格。具体地，将二维区域分解成若干个三角形，利用Galerkin方法得到离散方程。LOD格式的主要优点是可以处理非规则区域和非线性问题，但计算量较大。

下面给出Matlab程序，以六点对称格式为例：

%% 二维抛物方程的数值解法
% 六点对称格式
clear all;

a = 1; % 热传导系数
Lx = 1; % 区域长度
Ly = 1; % 区域宽度
T = 1; % 求解时间
Nx = 50; % x方向网格数
Ny = 50; % y方向网格数
Nt = 1000; % 时间层数
dx = Lx/Nx; % x方向网格步长
dy = Ly/Ny; % y方向网格步长
dt = T/Nt; % 时间步长

% 初始化
u = zeros(Nx+1,Ny+1,Nt+1);
x = linspace(0,Lx,Nx+1);
y = linspace(0,Ly,Ny+1);
t = linspace(0,T,Nt+1);
f = zeros(Nx+1,Ny+1,Nt+1);
u(:,:,1) = sin(pi*x)'*sin(pi*y); % 初始条件
u(:,1,:) = 0; % 边界条件
u(:,Ny+1,:) = 0;
u(1,:,:) = 0;
u(Nx+1,:,:) = 0;

% 求解
for n = 1:Nt
    for i = 2:Nx
        for j = 2:Ny
            u(i,j,n+1) = u(i,j,n) + a*dt/dx^2*(u(i+1,j,n)-4*u(i,j,n)+u(i-1,j,n)+u(i,j+1,n)+u(i,j-1,n)) + dt*f(i,j,n);
        end
    end
end

% 可视化
for n = 1:Nt+1
    surf(x,y,u(:,:,n));
    axis([0 Lx 0 Ly -1 1]);
    title(['t = ',num2str(t(n))]);
    drawnow;
end

其他格式的程序类似，只需要修改离散方程即可。