T1(i,j)=M1*(T(i-1,j)+T(i+1,j))+M2*(T(i,j-1)+T(i,j+1))+(1-2*M1-2*M2)*T(i,j); end end
for i=2:8 for j=2:8
T(i,j)=T1(i,j); end end if(k==5) l=l+1;
Tc(l)=T(5,5); k=0; end
end
i=1:9;j=1:9;
[x,y]=meshgrid(i); figure(1);
subplot(1,2,1);
mesh(x,y,T(i,j))% 画出10min后的温度场 axis tight;
xlabel('x','FontSize',14);ylabel('y','FontSize',14);zlabel('T/℃','FontSize',14) title('1min后二维温度场模拟图','FontSize',18) subplot(1,2,2);
[C,H]=contour(x,y,T(i,j)); clabel(C,H);axis square;
xlabel('x','FontSize',14);ylabel('y','FontSize',14); title('1min后模拟等温线图','FontSize',18) figure(2); xx=1:600;
plot(xx,Tc,'k-','linewidth',2)
xlabel('时间/s','FontSize