热传导方程差分格式.docx

热传导方程差分格式.docx

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

热传导方程差分格式

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

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

.:

u；：

2u

a2，0：

：

：

x：

：

1,

.:

t；x

u（x,0）=sin二x,0：

：

x:

:

1

u（O,t）=u（1,t）=0,t0

2

在t=0.05,0.1和0.2时刻的数值解，并与解析解u（x,t）=eftsin（二x）进行比较

1差分格式形式

2

设空间步长h=1/N,时间步长••0,T=M•，网比r=•/h.

（1）向前差分格式

该问题是第二类初边值问题（混合问题），我们要求出所需次数的偏微商的函数

Euc2u

u（x,t），满足方程—a—^,0：

：

：

x：

：

：

1,和初始条件u（x,0）=sinx,0x：

：

：

1

抚ex

及边值条件u（0,t）=u（1,t）=0,t0。

已知sin二x在相应区域光滑，并且在x=0,丨与边值相容，使问题有唯一充分光滑的解。

取空间步长h=1/N，和时间步长-T/M，其中N,M都是正整数。

用两族平行

直线x=jX=（jh0Hj1，和tlNtk=X（k=0,1,|||,M）将矩形域

G二｛0—x—1,t—0｝分割成矩形网络，网络格节点为（Xj,tk）。

以Gh表示网格内点集合，

即位于矩形G的网点集合；Gh表示闭矩形G的网格集合；丨h=Gh-Gh是网格界点的集合。

向前差分格式，即

k1kkkk

u,-u,u,1-2比■u,4

-j二a亠2亠t

（1）

£h

fi=f（x）,

0kk

Ujj=（Xj）,Uo=Un=0

其中，j=1,2，…，N_1,k=1,2,…，M一1.以r=a./h2表示网比。

（1）式可改写成如下：

u：

+=ru+（1-2r）u：

+ru：

」j

此格式为显格式。

其矩阵表达式如下:

1-2r

r

（j、

U1

厂j*、

U1

r

1-2r

+

u2

*

u2*

I-

r

1-2rr

j

UN」

j*

UNJ

r1-2r丿

j

j*

（2）向后差分格式

（3）

向后差分格式，即

u0

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

其矩阵表达式如下：

勺+2r

—r

\

fj八

U1

fj、

U1

—r

1+2r

+

r

U2

9

j

U2

■

-r1+2r-r

j+

UN4

j

UN4

-r1+2r>

j+

j

（4）六点对称格式

六点差分格式:

0

Uj

=（Xj）,U0=Un=0.

将（3）式改写成

rk.1一、kirk1rk_、krk

_2出1（1r）Uj―刁山/=?

山i（1_r）Uj2uj-

其矩阵表达式如下：

广1+r-r/2、

-r/21+r

+

j出

U2

■

'1-rr/2、

r/21-r

+

5\u2

*

—r/21+r-r/2

uNA

r/21-r

r/2

uN_i

、、_r1+2r/

uj+凹丿

r/2

1-2r/

2利用MATLAB求解问题的过程

对每种差分格式依次取N=40.,=1/1600,=1/3200，.=1/6400，用MATLAB

求解并图形比较数值解与精确解，用表格列出不同剖分时的L2误差。

向前差分格式：

t705:

=1/1600:

*数值解

——精确解

-1

00,10.2030.40.60.60.7D.日0.91

=1/3200：

-

、1

*数值解

精确解

/

■

/

-//

/

L/

/

-4

/

/

-/

/

#

fl111

li

\

\-\

\

\\\

06

05

0.4

03

02

01

°00.10.20.3Q.40.6G.60.70.00.91

.=1/6400:

*数值解

06

精确解

05

0.4

02

n£I，$*nlin■I中

00.10.2030.40.60.60.70,00.91

t=0.1:

=1/3600：

x10

2.5|

2-

*数值解

精确解

15-

-1--

+*

-1.5-*-

*

显-**-

■

9片|iIa||||i

''00.10.2030.40.50.60.70,00.91

.=1/3200:

*数值解

036

精确解-

0.25

0.15

0.05

n£H*nlin■I]

00.10.20.30.40.60.60.70,00.91

■=1/6400:

*数值解

精确解一

-/

Z*

V

*

///

X

\

1

j////

#

/

/

/

/

■/

11|]

h

\

\

\

%

\1

\・\-

\

11

0.36

03

0.25

02

0.15

005

°00.10.2030.40.50.60.70,00.9

t=0.2

.=1/1600：

=1/3200：

D.14

0.12

0,08

006

004

002

Q]IIHiIIIII

00.10.20.30.40.50.60.70,00.91

.=1/6400:

0.14

0.12

*数值解

精确解

0,08

006

004

002

n£h，i*nknij

00.10.20.30.40.60.60.70,00.91

向后差分格式:

t二0.05:

;U32OO

m6400

06

*数值解

精确解

05

0.4

Q]IIHiIIII]

00.10.20.30.40.50.60.70,00.91

t=0.1:

.=1/1600:

*数值解

0.36

精确解

0.25

0.15

0.05

Q]IIHiIIII]

=1/3200

0.36

*数值解

精确解

0.25

0.15

0.05

Q]IIHiliIIII]

00.10.20.30.40.50.60.70,00.91

.=1/6400

*数值解

036

精确解

0.25

02

0.15

005

t=0.2

=1/1600:

0.12

*数值解

精确解

0,08

006

004

002

Q]IIHiliIIIII

00.10.2030.40.50.60.70,00.91

.=1/3200

D.1J

0.12

0,08

006

004

002

n£h，$*nknij

■=1/6400

六点差分格式：

t=0.05:

.=1/1600:

Q]IIHiIIII]

00.10.20.30.40.60.60.70,00.91

=1/3200

07

06

05

0.4

03

02

1

0.2030.40.50.60.70,00.91

.=1/6400

06

05

0.4

*数值解

精确眸

02

t=0.1:

=1/1600:

0.36

0.25

*数值解

精确解

0.15

0.05

Q]IIHiliIIII]

00.10.20.30.40.50.60.70,00.91

.=1/3200

036

0.25

*数值解

精确解

02

0.15

0.05

n£【I，$*nlini]

00.10.20.30.40.60.60.70,00.91

■=1/6400

0.36

0.25

*数值解

精确解-

0.15

0.05

Q]IIHiliIIII]

00.10.20.30.40.50.60.70,00.91

t=0.2

.=1/1600:

0.12

0,08

006

004

002

Q]IIHiliIIIII

=1/3200

.=1/6400

0.14

0.12

*数值解

精确解

0,08

006

004

002

n£I，$*nRn■I]

00.10.20.30.40.60.60.70,00.91

L2误差:

AL

Tun

向前差分格式

现21

+

8e

1

0匸

2

0匸

向后差分格式

8

2

兀3

5

O2

1

0匸

2

0匸

六点差分格式

4

兀2

3

03

O

07

I诃

1

0匸

o

2a

MO

O

6

4

1

o

O9

2

0匸

9

09

8

8

3

O1

1o

5

3方法总结及分析

本文向前差分格式，向后差分格式及六点差分格式都是使用三对角系数矩阵，计算简单。

根据matlab作，特别明显的是，t=0.05®=1/1600:

时，图像解析解波动特别大，与真

1

解差距很大。

这是因为r=1一一差分格式不稳定。

根据误差比较，解本题时向后差分格式

2

误差最小。

衡量一个差分格式是否经济实用，有点多方面因素决定，应考虑不同的情况决定。

附件程序

向前差分格式：

function[u,err]=xq（t,t1）

%t是时刻值；

%t1是时间步长；N=40;

h=1/N;x=[h:

h:

（1-h）]';

r=t1/（hA2）;

u（:

1）=sin（pi*x）;%t=0时刻

R=zeros（N-1,N-1）;

fori=2:

N-2

R（i,i）=1-2*r;

R（i,i+1）=r;R（i,i-1）=r;

end

R（1,1）=1-2*r;R（1,2）=r;

R（N-1,N-1）=1-2*r;R（N-1,N-2）=r;k=t/t1;

fori=1:

k

u（:

i+1）=R*u（:

i）;endu=[0;u（:

i+1）;0];fori=1:

kforj=1:

N+1up（j,i）=exp（-pi*pi*t）*sin（pi*（j-1）*h）;endendx=0:

h:

1;

'r--'）;

plot（x,u,'b.',x,up（:

k）,

legend（'数值解','精确解'）;err=norm（u-up（:

k））;

end

向后差分格式：

function[u,err]=xh（t,t1）N=40;

h=1/N;

x=[h:

h:

（1-h）]';

r=t1/（hA2）;

u（:

1）=sin（pi*x）;%t=0时刻

R=zeros（N-1,N-1）;

fori=2:

N-2

R（i,i）=1+2*r;

R（i,i+1）=-r;

R（i,i-1）=-r；

end

R（1,1）=1+2*r;

R（1,2）=-r;

R（N-1,N-1）=1+2*r;

R（N-1,N-2）=-r;

k=t/t1;

fori=1:

k

u（:

i+1）=inv（R）*u（:

i）;

end

u=[0;u（:

i+1）;0];

fori=1:

k

forj=1:

N+1

up（j,i）=exp（-pi*pi*t）*sin（pi*（j-1）*h）;

end

end

x=0:

h:

1;

plot（x,u,b',x,up（:

k）,'r--'）;

legend（'数值解’，‘精确解'）;

err=norm（u-up（:

k））;

end

六点差分格式：

function[u,err]=ld（t,t1）

N=40;

h=1/N;

x=[h:

h:

（1-h）]';

r=t1/（hA2）;

u（:

1）=sin（pi*x）;%t=0时刻

R=zeros（N-1,N-1）;

fori=2:

N-2

R（i,i）=1+r;

R（i,i+1）=-r/2;

R（i,i-1）=-r/2;

end

R（1,1）=1+r;

R（1,2）=-r/2;

R（N-1,N-1）=1+2*r;

R（N-1,N-2）=-r;

fori=2:

N-2

R1（i,i）=1-r;

R1（i,i+1）=r/2;

R1（i,i-1）=r/2;

end

R1（1,1）=1-r;

R1（1,2）=r/2;

R1（N-1,N-1）=1-2*r;

R1（N-1,N-2）=r/2;

k=t/t1;

fori=1:

k

u（:

i+1）=inv（R）*R1*u（:

i）;

end

u=[0;u（:

i+1）;0];

fori=1:

k

forj=1:

N+1

up（j,i）=exp（-pi*pi*t）*sin（pi*（j-1）*h）;

end

end

x=0:

h:

1;

'r--'）;

plot（x,u,'b.',x,up（:

k）,legend（'数值解','精确解'）err=norm（u-up（:

k））;

end