matlab 利用crank-nicolson格式求解抛物型方程∂u/∂t = u_xx+u_yy,0<=x<=1,0<=y<=1,u(x,y,0)=sin(pi*x)sin(pi*y),u|∂ G = 0，给出示例代码

以下是利用Crank-Nicolson格式求解抛物型方程的Matlab代码示例：

% 定义参数和网格
Nx = 50;    % x方向网格数
Ny = 50;    % y方向网格数
Lx = 1;     % x方向长度
Ly = 1;     % y方向长度
dx = Lx/Nx; % x方向步长
dy = Ly/Ny; % y方向步长
dt = 0.01;  % 时间步长
T = 1;      % 最终时间
x = linspace(0, Lx, Nx+1);
y = linspace(0, Ly, Ny+1);
[X, Y] = meshgrid(x, y);

% 初始条件
u0 = sin(pi*X).*sin(pi*Y);

% 定义系数矩阵
r = dt/(2*dx^2);
s = dt/(2*dy^2);
a = -r*ones(Nx-1,1);
b = 1+2*r+2*s;
c = -r*ones(Nx-1,1);
d = -s*ones(Ny-1,1);
e = 1+2*r+2*s;
f = -s*ones(Ny-1,1);

% 初始化解
u = u0;

% 迭代求解
for t = dt:dt:T
    % 定义矩阵A和B
    A = spdiags([a b c], [-1 0 1], Nx-1, Nx-1);
    B = spdiags([d e f], [-1 0 1], Ny-1, Ny-1);
    
    % 定义右端项
    f = zeros(Nx-1, Ny-1);
    for j = 1:Ny-1
        f(:,j) = u(2:Nx,j+1)+r*(u(3:Nx+1,j+1)-2*u(2:Nx,j+1)+u(1:Nx-1,j+1));
    end
    
    % 利用Thomas算法求解A*u_n+1=A*u_n+f/2
    for j = 1:Ny-1
        u(2:Nx,j+1) = thomas(A, u(2:Nx,j+1)+f(:,j)/2);
    end
    
    % 利用Thomas算法求解B*u_n+1=B*u_n+f/2
    for i = 1:Nx-1
        u(i+1,2:Ny) = thomas(B, u(i+1,2:Ny)'+f(i,:)'/2)';
    end
end

% 可视化结果
surf(X, Y, u);
xlabel('x');
ylabel('y');
zlabel('u');

在上述代码中，我们利用了Thomas算法来求解三对角矩阵，可以通过以下函数实现：

function x = thomas(A, b)
% Thomas算法求解三对角线性方程组Ax=b
% A为三对角矩阵，b为列向量

n = length(b);
a = diag(A, -1);
d = diag(A);
c = diag(A, 1);

% 前向消元
for k = 2:n
    factor = a(k-1)/d(k-1);
    d(k) = d(k)-factor*c(k-1);
    b(k) = b(k)-factor*b(k-1);
end

% 回代求解
x(n) = b(n)/d(n);
for k = n-1:-1:1
    x(k) = (b(k)-c(k)*x(k+1))/d(k);
end

希望能对你有所帮助！