热传导方程的差分格式

热传导方程的差分格式——上机实习报告二零一四年五月热传导方程的差分格式第1页一维抛物方程的初边值问题分别用向前差分格式、向后差分格式、六点对称格式，求解下列问题：22,01,uuaxtx(,0)sin,01uxxx(0,)(1,)0,0ututt在0.05,0.10.2t和时刻的数值解，并与解析解2(,)sin()tuxtex进行比较。1差分格式形式设空间步长1/hN,时间步长0，TM，网比2/rh.（1）向前差分格式向前差分格式，即ikjkjkjkjkjfhuuuauu21112（1）)(iixff，0),(00kNkjjjuuxu其中，.1,,2,1,1,,2,1MkNj以2/har表示网比。（1）式可改写成如下：jkjkjkjkjfruurruu111)21(此格式为显格式。其矩阵表达式如下：111121112121212121jNjNjjjNjNjjuuuuuuuurrrrrrrrr（2）向后差分格式热传导方程的差分格式第2页向后差分格式，即,221111k11jkjkjjkjkjfhuuuauu（2）,0),(N00kkjjjuuxu其中.1,,2,1,1,,2,1MkNj（2）式可改写成jkjkjkjkjfuruurru11111)21(此种差分格式被称为隐格式。其矩阵表达式如下：jNjNjjjNjNjjuuuuuuuurrrrrrrrr121111121121212121（3）六点对称格式六点差分格式：jkjkjkjkjkjkjkjkjfhuuuhuuuauu]22[22112111111（3）.0),(00kNkjjjuuxu将（3）式改写成jkjkjkjkjkjkjfurururururur11111112)1(22)1(2其矩阵表达式如下:jNjNjjjNjNjjuuuurrrrrrrrruuuurrrrrrrrr1211111211212/2/12/12/2/1212/12/12/2/12利用MATLAB求解问题的过程对每种差分格式依次取40.N，=1/1600，=1/3200，=1/6400，用MATLAB求解并图形比较数值解与精确解，用表格列出不同剖分时的2L误差。（1）向前差分格式热传导方程的差分格式第3页热传导方程的差分格式第4页热传导方程的差分格式第5页热传导方程的差分格式第6页（2）向后差分格式热传导方程的差分格式第7页热传导方程的差分格式第8页热传导方程的差分格式第9页（3）六点对称格式热传导方程的差分格式第10页热传导方程的差分格式第11页热传导方程的差分格式第12页3方法总结及分析六点对称方法克服了向前差分格式与向后差分格式方法误差偏大的缺点，更接近准确数值本文向前差分格式，向后差分格式及六点差分格式都是使用三对角系数矩阵，计算简单。根据matlab作，特别明显的是，:05.0t:1600/1时，图像解析解波动特别大，与真解差距很大。附件程序（1）向前差分格式源程序：function[e]=forward(h,t,T)N=1/h;x=zeros(1,N+1);fori=1:Nx(i+1)=i*h;endu=sin(pi.*x);r=t/(h*h);forj=1:T/t;u=for_climb(u,r,t);endplot(x,u,'o')holdu_xt=exp(-pi*pi*T)*sin(pi.*x);plot(x,u_xt,'r')e=u_xt-u;functionb=for_climb(a,r,t)L=length(a);b=zeros(1,L);fori=1:L-2b(i+1)=r*a(i+2)+(1-2*r)*a(i+1)+r*a(i);end（2）向后差分格式热传导方程的差分格式第13页源程序：function[e]=back(dx,dt,T)M=1/dx;N=T/dt;u1=zeros(1,M+1);x=[1:M-1]*dx;u1([2:M])=sin(pi.*x);r=dt/dx/dx;r2=1+2*r;fori=1:M-1A(i,i)=r2;ifi1A(i-1,i)=-r;A(i,i-1)=-r;endendu=zeros(N+1,M+1);u(N+1,:)=u1;fork=1:Nb=u(N+2-k,2:M);u(N+1-k,2:M)=inv(A)*b';enduT=u(1,:);x1=[0,x,1];plot(x1,uT,'o')holdu_xt=exp(-pi*pi*T)*sin(pi.*x1);plot(x1,u_xt,'r')e=u_xt-uT;（3）六点对称格式：源程序：function[e]=six(dx,dt,T)M=1/dx;N=T/dt;u1=zeros(1,M+1);x=[1:M-1]*dx;u1([2:M])=sin(pi*x);r=dt/dx/dx;r2=2+2*r;r3=2-2*r;fori=1:M-1A(i,i)=r2;ifi1热传导方程的差分格式第14页A(i-1,i)=-r;A(i,i-1)=-r;endendfori=1:M-1B(i,i)=r3;ifi1B(i-1,i)=r;B(i,i-1)=r;endendu=zeros(N+1,M+1);u(N+1,:)=u1;fork=1:Nb=B*(u(N+2-k,2:M))';u(N+1-k,2:M)=inv(A)*b;enduT=u(1,:);x1=[0,x,1];plot(x1,uT,'o')holdu_xt=exp(-pi*pi*T)*sin(pi.*x1);plot(x1,u_xt,'r')e=u_xt-uT;