热传导方程的差分格式

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