热传导方程的差分格式.docVIP

热传导方程的差分格式 第 PAGE 2页  热传导方程的差分格式 ——上机实习报告 二零一四年五月 热传导方程的差分格式 第 PAGE 1页 一维抛物方程的初边值问题 分别用向前差分格式、向后差分格式、六点对称格式，求解下列问题： MACROBUTTON MTEditEquationSection2 SEQ MTEqn \r \h SEQ MTSec \r 1 \h SEQ MTChap \r 3 \h MACROBUTTON MTEditEquationSection2 SEQ MTEqn \r \h SEQ MTSec \r 3 \h SEQ MTChap \h 在时刻的数值解，并与解析解。 1差分格式形式 设空间步长, 时间步长，，网比. 向前差分格式 向前差分格式，即 （1） ， 其中，以表示网比。（1）式可改写成如下： 此格式为显格式。 其矩阵表达式如下： 向后差分格式 向后差分格式，即 （2） 其中（2）式可改写成 此种差分格式被称为隐格式。 其矩阵表达式如下： 六点对称格式 六点差分格式： （3） 将（3）式改写成 其矩阵表达式如下: 2利用MATLAB求解问题的过程 对每种差分格式依次取，，，，用MATLAB求解并图形比较数值解与精确解，用表格列出不同剖分时的误差。 （1）向前差分格式 （2）向后差分格式 （3）六点对称格式 3 方法总结及分析 六点对称方法克服了向前差分格式与向后差分格式方法误差偏大的缺点，更接近准确数值 本文向前差分格式，向后差分格式及六点差分格式都是使用三对角系数矩阵，计算简单。根据matlab作，特别明显的是，时，图像解析解波动特别大，与真解差距很大。 附件程序 （1）向前差分格式 源程序： function[e]=forward(h,t,T) N=1/h; x=zeros(1,N+1); for i=1:N x(i+1)=i*h; end u=sin(pi.*x); r=t/(h*h); for j=1:T/t; u=for_climb(u,r,t); end plot(x,u,o) hold u_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-2 b(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); r=dt/dx/dx; r2=1+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 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; 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,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;