用有限差分方法中的ADI格式求解偏微分方程，并用matlab写代码

ADI（Alternating Direction Implicit）方法是一种常用的求解偏微分方程的方法，它将二阶偏微分方程分解为两个一维方程，分别用Crank-Nicolson格式和向后Euler格式求解。下面我们介绍如何用ADI方法求解二维热传导方程：

$$\frac{\partial u}{\partial t}=k(\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2})$$

我们将以上方程用Crank-Nicolson格式离散化：

$$\frac{u_{i,j}^{n+1}-u_{i,j}^n}{\Delta t}=\frac{k}{2}(\frac{u_{i-1,j}^{n+1}-2u_{i,j}^{n+1}+u_{i+1,j}^{n+1}}{(\Delta x)^2}+\frac{u_{i,j-1}^n-2u_{i,j}^n+u_{i,j+1}^n}{(\Delta y)^2}+\frac{u_{i-1,j}^n-2u_{i,j}^n+u_{i+1,j}^n}{(\Delta x)^2}+\frac{u_{i,j-1}^{n+1}-2u_{i,j}^{n+1}+u_{i,j+1}^{n+1}}{(\Delta y)^2})$$

将上式进行整理得：

$$-au_{i-1,j}^{n+1}+(2a+2b+1)u_{i,j}^{n+1}-bu_{i+1,j}^{n+1}=cu_{i,j-1}^n+(2c+2d+1)u_{i,j}^n-du_{i,j+1}^n$$

其中，$a=\frac{k\Delta t}{2(\Delta x)^2}$，$b=\frac{k\Delta t}{2(\Delta x)^2}$，$c=\frac{k\Delta t}{2(\Delta y)^2}$，$d=\frac{k\Delta t}{2(\Delta y)^2}$，$1+2a+2b=2c+2d+1=1$。

我们可以通过将时间步长$\Delta t$和空间步长$\Delta x$、$\Delta y$的比值保持不变，来提高数值解的精度。因此，我们将$\Delta x=\Delta y=h$，$\Delta t=\frac{1}{2}(h^2)/k$。

根据以上推导，我们可以得到ADI格式的迭代公式：

$$u_{i,j}^{*}=\frac{1}{2}(au_{i-1,j}^{n+1}+(2a+2b+1)u_{i,j}^{n}+bu_{i+1,j}^{n})$$

$$u_{i,j}^{n+1}=\frac{1}{2}(cu_{i,j-1}^{n}+(2c+2d+1)u_{i,j}^{*}+du_{i,j+1}^{n})$$

其中，上述公式中的$*$表示该点的值是通过Crank-Nicolson格式得到的中间值。

最后，我们可以通过以下的Matlab代码来实现ADI格式的求解过程：

% 定义参数
Lx = 1; % x轴方向长度
Ly = 1; % y轴方向长度
nx = 41; % x轴方向网格数
ny = 41; % y轴方向网格数
dx = Lx / (nx - 1); % x轴方向网格间距
dy = Ly / (ny - 1); % y轴方向网格间距
x = 0:dx:Lx; % x轴网格
y = 0:dy:Ly; % y轴网格
k = 1; % 热传导系数
dt = 0.5 * (dx^2) / k; % 时间步长
nt = 1200; % 时间步数
u = zeros(nx, ny, nt); % 温度分布矩阵

% 边界条件
u(1, :, :) = 100; % 左边界
u(nx, :, :) = 100; % 右边界
u(:, 1, :) = 0; % 下边界
u(:, ny, :) = 0; % 上边界

% 初始化
u(:, :, 1) = 0; % 初始温度分布

% 迭代求解
for n = 1:nt
    % x方向求解
    for j = 2:ny-1
        % 构造三对角矩阵
        a = -0.5 * k * dt / (dx^2);
        b = -0.5 * k * dt / (dx^2);
        c = 1 + k * dt / (dx^2);
        A = diag(a*ones(nx-2,1),-1) + diag(c*ones(nx-1,1),0) + diag(b*ones(nx-2,1),1);
        % 更新中间值
        for i = 2:nx-1
            B = [a*u(i-1,j,n+1); (c+2*b)*u(i,j,n) + a*u(i-1,j,n+1) + b*u(i+1,j,n); b*u(i+1,j,n+1)];
            u(i,j,n+1) = A\B;
        end
    end
    
    % y方向求解
    for i = 2:nx-1
        % 构造三对角矩阵
        a = -0.5 * k * dt / (dy^2);
        b = -0.5 * k * dt / (dy^2);
        c = 1 + k * dt / (dy^2);
        A = diag(a*ones(ny-2,1),-1) + diag(c*ones(ny-1,1),0) + diag(b*ones(ny-2,1),1);
        % 更新温度分布
        for j = 2:ny-1
            B = [a*u(i,j-1,n+1); (c+2*b)*u(i,j,n+1) + a*u(i,j-1,n+1) + b*u(i,j+1,n); b*u(i,j+1,n+1)];
            u(i,j,n+1) = A\B;
        end
    end
end

% 绘制温度分布图像
[X, Y] = meshgrid(x, y);
for n = 1:nt
    surf(X, Y, u(:, :, n));
    zlim([0 100]);
    drawnow;
end