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

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

《热传导方程的差分格式.docx》由会员分享，可在线阅读，更多相关《热传导方程的差分格式.docx（16页珍藏版）》请在冰豆网上搜索。

热传导方程的差分格式

热传导方程的差分格式——上机实习报告

二零一四年五月

一维抛物方程的初边值问题

分别用向前差分格式、向后差分格式、六点对称格式，求解下列问题：

在

时刻的数值解，并与解析解

。

1差分格式形式

设空间步长

时间步长

，

，网比

.

（1）向前差分格式

向前差分格式，即

（1）

，

其中，

以

表示网比。

（1）式可改写成如下：

此格式为显格式。

其矩阵表达式如下：

（2）向后差分格式

向后差分格式，即

（2）

其中

（2）式可改写成

此种差分格式被称为隐格式。

其矩阵表达式如下：

（3）六点对称格式

六点差分格式：

（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）;

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'）

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）;

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）;

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-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））';

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;