热传导方程的差分格式

1、热传导方程的差分格式 上机实习报告二零一四年五月热传导方程的差分格式 第13页一维抛物方程的初边值问题分别用向前差分格式、向后差分格式、六点对称格式，求解下列问题：在时刻的数值解，并与解析解。1差分格式形式 设空间步长, 时间步长，网比.（1） 向前差分格式向前差分格式，即 （1） ，其中，以表示网比。（1）式可改写成如下：此格式为显格式。其矩阵表达式如下：（2） 向后差分格式向后差分格式，即（2）其中（2）式可改写成此种差分格式被称为隐格式。其矩阵表达式如下：（3） 六点对称格式六点差分格式： （3）将（3）式改写成其矩阵表达式如下:2利用MATLAB求解问题的过程 对每种差分格式依次取，用

2、MATLAB求解并图形比较数值解与精确解，用表格列出不同剖分时的误差。（1）向前差分格式（2）向后差分格式（3）六点对称格式3 方法总结及分析六点对称方法克服了向前差分格式与向后差分格式方法误差偏大的缺点，更接近准确数值本文向前差分格式，向后差分格式及六点差分格式都是使用三对角系数矩阵，计算简单。根据matlab作，特别明显的是，时，图像解析解波动特别大，与真解差距很大。附件程序（1）向前差分格式源程序：functione=forward(h,t,T)N=1/h;x=zeros(1,N+1);for i=1:Nx(i+1)=i*h;endu=sin(pi.*x);r=t/(h*h);for j

3、=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;function b=for_climb(a,r,t)L=length(a);b=zeros(1,L);for i=1:L-2b(i+1)=r*a(i+2)+(1-2*r)*a(i+1)+r*a(i);end（2）向后差分格式源程序：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)

4、;r=dt/dx/dx;r2=1+2*r;for i=1:M-1A(i,i)=r2;if i1 A(i-1,i)=-r;A(i,i-1)=-r;end end u=zeros(N+1,M+1);u(N+1,:)=u1;for k=1:N b=u(N+2-k,2:M); u(N+1-k,2:M)=inv(A)*b;end uT=u(1,:); x1=0,x,1; plot(x1,uT,o)hold u_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

5、;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;for i=1:M-1 A(i,i)=r2; if i1 A(i-1,i)=-r; A(i,i-1)=-r;end end for i=1:M-1 B(i,i)=r3;if i1 B(i-1,i)=r;B(i,i-1)=r; end end u=zeros(N+1,M+1);u(N+1,:)=u1;for k=1:N b=B*(u(N+2-k,2:M); u(N+1-k,2:M)=inv(A)*b;end uT=u(1,:); x1=0,