T=zeros(10,10);T1=zeros(10,10);T0=zeros(10,10);
t=0;k=0;
for x=1:10
for y=1:10
T(x,y)=1.05-y/10;
T1(x,y)=1.05-y/10;
T0(x,y)=1.05-y/10;
end
end
for n=1:700%时层
for x=1:10
for y=1:10
if x==1&&y==1
T1(x,y)=(2+T0(x,y+1)+T0(x+1,y)+10*T0(x,y)+0.1)/14;
elseif x==1&&y==10
T1(x,y)=(T0(x,y-1)+T0(x+1,y)+10*T0(x,y)+0.1)/14;
elseif x==10&&y==1
T1(x,y)=(2+T0(x,y+1)+T0(x-1,y)+10*T0(x,y)+0.1)/14;
elseif x==10&&y==10
T1(x,y)=(T0(x,y-1)+T0(x-1,y)+10*T0(x,y)+0.1)/14;
elseif x==1&&y~=1&&y~=10
T1(x,y)=(T0(x+1,y)+T0(x,y+1)+T0(x,y-1)+10*T0(x,y)+0.1)/13;
elseif x==10&&y~=1&&y~=10
T1(x,y)=(T0(x-1,y)+T0(x,y+1)+T0(x,y-1)+10*T0(x,y)+0.1)/13;
elseif y==1&&x~=1&&x~=10
T1(x,y)=(2+T0(x-1,y)+T0(x,y+1)+T0(x+1,y)+10*T0(x,y))/15;
elseif y==10&&x~=1&&x~=10
T1(x,y)=(T0(x-1,y)+T0(x,y-1)+T0(x+1,y)+10*T0(x,y))/15;
elseif x~=1&&x~=10&&y~=1&&y~=10
T1(x,y)=(T0(x+1,y)+T0(x-1,y)+T0(x,y+1)+T0(x,y-1)+10*T0(x,y))/14;
end
end
end
k=k+1;
E=0;
for x=1:10
for y=1:10
if (abs(T1(x,y)-T0(x,y))/T0(x,y))>E
E=(abs(T1(x,y)-T0(x,y))/T0(x,y));
end
end
end
if E<0.000025%判断时层达到稳态
评论0
最新资源