对流matlab程序,热对流代码MATLAB

PI=3.1415192756;

h=0.025;%网格间距

T=0.001;%时间间隔

t=10;%时间

M=1/h;

N=t/T;

u=zeros(M+1,M+1);

v=zeros(M,M);

d=zeros(M);

q=zeros(M);

w=zeros(M);

E=zeros(N);

r=T/(h*h);

a=-r/2;

b=1+r;

c=a;

%边界值上的初始值

n=1;

for i=1:M+1

u(1,i)=20+80*(i-1)*h;

end

for i=1:M+1

u(M+1,i)=20+80*((i-1)*h-exp(-0.5*PI*PI*(n-1)*T)*sin(PI*(i-1)*h/2));

end

for i=1:M+1

u(i,1)=20;

end

for i=1:M+1

u(i,M+1)=20+80*(1-exp(-0.5*PI*PI*(n-1)*T)*sin(PI*(i-1)*h/2));

end

for n=1:N

%x方向

for j=2:M

v(1,j)=0.5*(u(1,j)+u(1,j))-0.25*r*((u(1,j-1)-2*u(1,j)+