1、热传导方程的差分格式上机实习报告二零一四年五月一维抛物方程的初边值问题分别用向前差分格式、向后差分格式、六点对称格式，求解以下问题:在r=0.05,0.1和0.2时刻的数值解，并与解析解(x,1)=eSin(G)进行比较。1差分格式形式设空间步长力=1N,时间步长ZO,T=Mr,网比r=*.(1)向前差分格式向前差分格式，即ft = f(xi) 其中，=12一,-1次=1,2,-,加一1.以厂=22表示网比。(1)式可改写成如下:此格式为显格式。其矩阵表达式如下:(2)向后差分格式向后差分格式，即.*+1Hk.+17fi+,-4-j+,Uj_UjiUAZUj+My-1-二-a+%其中/=1,2
2、N-1M=1,2,/-1.(2)式可改写成此种差分格式被称为隐格式。其矩阵表达式如下：(3)六点对称格式六点差分格式:a ukj -2wy+1 +wl1IT将(3)式改写成其矩阵表达式如下:1 + r-r2r/21-r2利用MATLAB求解问题的过程对每种差分格式依次取N=40.,r=l1600,r=l3200,%1/6400,用MATLAB求解并图形比拟数值解与精确解，用表格列出不同剖分时的Z3误差。(1)向前差分格式(2)向后差分格式(3)六点对称格式3方法总结及分析六点对称方法克服了向前差分格式与向后差分格式方法误差偏大的缺点，更接近准确数值本文向前差分格式，向后差分格式及六点差分格式
3、都是使用三对角系数矩阵，计算简单。根据mat根b作,特别明显的是，/=0.05:=l1600:时，图像解析解波动特别大，与真解差距很大。附件程序(1)向前差分格式源程序：functione=forward(h,t,T)N=lh;x=zeros(1,N+1);fori=l:Nx(i+l)=i*h;u=sin(pi.*x);r=t(h*h);forj=l: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)1
4、Iength(a);b=zeros(1,L);fori=l:L-2b(i+l)=r*a(i+2)+(l-2*r)*a(i+l)+r*a(i);end(2)向后差分格式源程序：functione=back(dx,dt,T)M=ldx;N=Tdt;ul=zeros(1,M+l);x=1:M-1*dx;ul(2:M)=sin(pi.*x);r=dtdxdx;r2=l+2*r;fori=l:M-1A(i,i)=r2;ifilA(i-l,i)=-r;A(i,i-l)=-r;endendu=zeros(N+l,M+l);u(N+l,:)=ul;fork=l:Nb=u(N+2-k,2:M);u(N+l-
5、k,2:M)=inv（八）*b,;enduT=u(1,:);xl=0,x,1;plot(xl,uT,o)holdu_xt=exp(-pi*pi*T)*sin(pi.*xl);plot(xl,u_xt,r)e=u_xt-uT;六点对称格式：源程序：functione=six(dx,dt,T)M=ldx;N=Tdt;ul=zeros(1,M+l);x=1:M-1*dx;ul(2:M)=sin(pi*x);r=dtdxdx;r2=2+2*r;r3=2-2*r;fori=l:M-1A(i,i)=r2;ifilA(i-l,i)=-r;A(i,i-l)=-r;endendfori=l:M-1B(izi)=r3;ifilB(i-l,i)=r;B(i,i-l)=r;endendu=zeros(N+l,M+l);u(N+l,:)=ul;fork=l:Nb=B*(u(N+2-k,2:M),;u(N+l-kz2:M)=inv（八）*b;enduT=u(1,:);xl=0,x,1;plot(xl,uT,o)holdu_xt=exp(-pi*pi*T)*sin(pi.*xl);plot(xl,u_xt,r)e=uxt-uT;
