编写MATLAB代码，实现二维平面热传导中绝热壁面条件下，壁面温度的计算

假设热传导方程为：

$\frac{\partial T}{\partial t} = \alpha (\frac{\partial^2 T}{\partial x^2} + \frac{\partial^2 T}{\partial y^2})$

其中 $T$ 为温度，$t$ 为时间，$x$ 和 $y$ 为空间坐标，$\alpha$ 为热扩散系数。

在绝热壁面条件下，壁面温度不发生变化，因此可以将边界温度设置为常数。

下面是 MATLAB 代码实现：

% 设置参数
Lx = 1; % x方向长度
Ly = 1; % y方向长度
Nx = 50; % x方向网格数
Ny = 50; % y方向网格数
alpha = 1; % 热扩散系数
T0 = 0; % 初始温度
Tb = 100; % 边界温度

% 初始化温度分布
T = T0 * ones(Ny,Nx);
T(:,1) = Tb; % x=0边界
T(:,end) = Tb; % x=Lx边界
T(1,:) = Tb; % y=0边界
T(end,:) = Tb; % y=Ly边界

% 设置时间步长和总时间
dt = 0.001;
t = 0:dt:1;

% 迭代计算温度分布
for i = 1:length(t)
    % 计算中心部分的二阶偏导数
    d2Tdx2 = diff(T, 2, 2) / (Lx/Nx)^2;
    d2Tdy2 = diff(T, 2, 1) / (Ly/Ny)^2;
    % 补全边缘部分的二阶偏导数
    d2Tdx2 = [zeros(Ny,1), d2Tdx2, zeros(Ny,1)];
    d2Tdy2 = [zeros(1,Nx+2); d2Tdy2; zeros(1,Nx+2)];
    % 计算下一个时间步的温度分布
    T = T + alpha * dt * (d2Tdx2 + d2Tdy2);
    % 保持边界温度不变
    T(:,1) = Tb;
    T(:,end) = Tb;
    T(1,:) = Tb;
    T(end,:) = Tb;
end

% 绘制温度分布图
x = linspace(0,Lx,Nx+2);
y = linspace(0,Ly,Ny+2);
[X,Y] = meshgrid(x,y);
figure;
surf(X,Y,T);
xlabel('x');
ylabel('y');
zlabel('Temperature');

在运行代码后，会得到一个二维平面的温度分布图。其中红色代表高温，蓝色代表低温。在边缘处可以看到温度始终保持为100度，符合绝热壁面条件。