苏佳园：二维热传导方程有限差分法的MATLAB实现

附 录

追赶法程序

function x=chase(A,f) n=rank(A); b=zeros(n,1); a=zeros(n-1,1); c=zeros(n-1,1); for i=1:n-1 b(i,1)=A(i,i); a(i,1)=A(i+1,1); c(i,1)=A(i,i+1); end

b(n,1)=A(i,i); for i=2:n

b(i,1)=f(i,1)-(a(i-1,1)/b(i-1,1))*c(i-1,1); f(i,1)=f(i,1)-(a(i-1,1)/b(i-1,1))*f(i-1,1); end

x(n)=f(n,1)-b(n,1); for i=(n-1):-1:1

x(i)=(f(i,1)-c(i,1)*x(i+1))/b(i,1); end

26

求解热传导方程的程序

function [u,x,y,t]=sjy(a,D,T,u_xy0,u_xyt,Mx,My,N) %a为方程系数

%D为在x轴和y轴上的边界值,D(1)<=x<=D(2),D(3)<=y<=D(4) %T为时间上限 %u_xy0为t=0时的初值 %u_xyt为边界取值函数 %Mx为x轴的等分段数 %My为y轴的等分段数 %N为时间轴t的等分段数

ox=(D(2)-D(1))/Mx;x=D(1)+[0:Mx]*ox; oy=(D(4)-D(3))/My;y=D(3)+[0:My]*oy; ot=T/N;t=[0:N]*ot; %初始化u for j=1:Mx+1 for i=1:My+1

u(i,j)=u_xy0(x(j),y(i)); end end

rx=a*ot/(ox*ox) ry=a*ot/(oy*oy) rx1=1+2*rx;rx2=1-2*rx; ry1=1+2*ry;ry2=1-2*ry; for j=1:Mx-1 A(j,j)=ry1; if j>1

A(j-1,j)=-ry; A(j,j-1)=-ry; end end

%A为y方向隐式时的系数矩阵

27

苏佳园：二维热传导方程有限差分法的MATLAB实现

for i=1:My-1 B(i,i)=rx1; if i>1

B(i-1,i)=-rx; B(i,i-1)=-rx; end end

%B为X方向隐式时的系数矩阵 for k=1:N

u_1=u; t=k*ot; for i=1:My+1

u(i,1)=feval(u_xyt,x(1),y(i),t); u(i,Mx+1)=feval(u_xyt,x(Mx+1),y(i),t); end

for j=1:Mx+1

u(1,j)=feval(u_xyt,x(j),y(1),t); u(My+1,j)=feval(u_xyt,x(j),y(My+1),t); end

if mod(k,2)==0 for i=2:My jj=2:Mx;

bx=[ry*u(i,1),zeros(1,Mx-3),ry*u(i,My+1)]+rx*(u_1(i-1,jj)+u_1(i+1,jj))+rx2*u_1(i,jj);

u(i,jj)=linsolve(A,bx');%bx为y方向隐式时的常数项 end else

for j=2:Mx ii=2:My;

by=[rx*u(1,j);zeros(My-3,1);rx*u(Mx+1,j)]+ry*(u_1(ii,j-1)+u_1(ii,j+1))+ry2*u_1(i

i,j);

u(ii,j)=linsolve(B,by);%by为x方向隐式时的常数项 end

28

end end

29