Crank Nicolson求二维热传导方程

二维热传导方程可以写成如下形式：

<img src="https://latex.codecogs.com/svg.image?\begin{aligned}&space;\frac{\partial&space;u}{\partial&space;t}&space;=&space;\alpha&space;\left(\frac{\partial^2&space;u}{\partial&space;x^2}&space;&plus;&space;\frac{\partial^2&space;u}{\partial&space;y^2}\right)\\&space;u(x,y,0)&space;=&space;f(x,y)\\&space;u(a,y,t)&space;=&space;g_1(y,t)\\&space;u(b,y,t)&space;=&space;g_2(y,t)\\&space;u(x,c,t)&space;=&space;g_3(x,t)\\&space;u(x,d,t)&space;=&space;g_4(x,t)&space;\end{aligned}" title="\begin{aligned} \frac{\partial u}{\partial t} = \alpha \left(\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2}\right)\\ u(x,y,0) = f(x,y)\\ u(a,y,t) = g_1(y,t)\\ u(b,y,t) = g_2(y,t)\\ u(x,c,t) = g_3(x,t)\\ u(x,d,t) = g_4(x,t) \end{aligned}" />

其中，<img src="https://latex.codecogs.com/svg.image?\alpha" title="\alpha" /> 为热扩散系数，<img src="https://latex.codecogs.com/svg.image?f(x,y)" title="f(x,y)" /> 为初始温度分布，<img src="https://latex.codecogs.com/svg.image?g_1(y,t),g_2(y,t),g_3(x,t),g_4(x,t)" title="g_1(y,t),g_2(y,t),g_3(x,t),g_4(x,t)" /> 分别为四边界的热传导条件。

Crank-Nicolson方法是一种求解偏微分方程的有限差分方法，其中采用了向前差分和向后差分的平均值来计算时间步长上的值。下面是 Crank-Nicolson 方法的 MATLAB 程序示例：

% 定义问题参数
a = 0; b = 1; c = 0; d = 1;
alpha = 0.1;
dx = 0.1; dy = 0.1;
dt = 0.01;
x = a:dx:b;
y = c:dy:d;
Nx = length(x);
Ny = length(y);

% 定义初始温度分布
u = zeros(Nx,Ny);
for i=1:Nx
    for j=1:Ny
        u(i,j) = sin(pi*x(i))*sin(pi*y(j));
    end
end

% 定义边界条件
g1 = zeros(Ny,1);
g2 = zeros(Ny,1);
g3 = zeros(Nx,1);
g4 = zeros(Nx,1);
for j=1:Ny
    g1(j) = sin(pi*y(j));
    g2(j) = sin(pi*y(j))*exp(-pi);
end
for i=1:Nx
    g3(i) = 0;
    g4(i) = 0;
end

% 定义矩阵
r = alpha*dt/(2*dx*dx);
s = alpha*dt/(2*dy*dy);
A = -r*diag(ones(Nx-1,1),-1) - r*diag(ones(Nx-1,1),1) + (1+2*r)*eye(Nx);
B = -s*diag(ones(Ny-1,1),-1) - s*diag(ones(Ny-1,1),1) + (1+2*s)*eye(Ny);

% 迭代求解
for k=1:100
    u_old = u;
    for j=1:Ny
        f = (1-2*s)*u(:,j) + s*[g1(j); zeros(Nx-2,1); g2(j)] + s*dt/2*[g3; zeros(Nx-2,1); g4] + dt/2*[zeros(1,Nx); u(2:Nx-1,1:Ny-1)]*B;
        u(:,j) = (A+B)\f;
    end
    if max(max(abs(u_old-u))) < 1e-6
        break;
    end
end

% 显示结果
surf(x,y,u');
xlabel('x');
ylabel('y');
zlabel('u');

在这个程序中，我们首先定义了问题的参数，初始温度分布和边界条件。然后，我们定义了矩阵 A 和 B，用于计算迭代中每个时间步长的解。最后，我们使用双重循环迭代求解二维热传导方程，并使用 surf 函数绘制了最终的温度分布图像。