一维热传导方程的差分格式

《微分方程数值解》

课程论文

学生姓名1：许慧卿学号：20144329 学生姓名2：向裕学号：20144327学生姓名3：邱文林学号：20144349学生姓名4：高俊学号：20144305学生姓名5：赵禹恒学号：20144359学生姓名6：刘志刚学号： 20144346 学院：理学院

专业：14级信息与计算科学

指导教师：陈红斌

2017年6 月25日

《偏微分方程数值解》课程论文

《一维热传导方程的差分格式》论文

一、《微分方程数值解》课程论文的格式

1）引言：介绍研究问题的意义和现状

2）格式：给出数值格式

3）截断误差：给出数值格式的截断误差

4）数值例子：按所给数值格式给出数值例子

5）参考文献：论文所涉及的文献和教材

二、《微分方程数值解》课程论文的评分标准

1）文献综述：10分；

2）课题研究方案可行性：10分；

3）数值格式：20分；

4）数值格式的算法、流程图：10分；

5）数值格式的程序：10分；

6）论文撰写的条理性和完整性：10分；

7）论文工作量的大小及课题的难度：10分；

8）课程设计态度：10分；

9）独立性和创新性：10分。

评阅人：

- 2 -

一维热传导方程的差分格式

1 引言

考虑如下一维非齐次热传导方程Dirichlet 初边值问题

22(,),u u

a f x t t x

??=+?? ,c x d << 0,t T <≤ (1.1) (,0)(),u x x ?= ,c x d ≤≤ (1.2)

(,)(),u c t t α= (,)(),u d t t β= 0t T <≤ (1.3)

的有限差分方法, 其中a 为正常数,(,),(),(),

()f x t x t t ?αβ为已知常数, ()(0),c ?α=

()(0).d ?β= 称(1.2)为初值条件, (1.3)为边值条件.

本文将给出(1.1) (1.3)的向前Euler 格式, 向后Euler 格式和Crank Nicolson -格式, 并给出其截断误差和数值例子. 经对比发现, Crank Nicolson -格式误差最小, 向前

Euler 格式次之, 向后Euler 格式误差最大.

2 差分格式的建立

2.1 向前Euler 格式

将区间[,]c d 作M 等分, 将[]0,T 作N 等分, 并记 ()/h d c M =-, /T N τ=,

j x c jh =+,0j M ≤≤, k t k τ=,0k N ≤≤. 分别称h 和τ为空间步长和时间步长.用

两组平行直线

j x x =, 0j M ≤≤, k t t =, 0k N ≤≤

将Ω分割成矩形网格.记{}

|0h j x j M Ω=≤≤, {}|0k t k N τΩ=≤≤, h h ττΩ=Ω?Ω. 称()

,j k x t 为结点[1]

.

定义h τΩ上的网格函数 {}|0,0k j U j M k N Ω=≤≤≤≤, 其中()

,k

j j k U u x t =.

在结点()

,j k x t 处考虑方程(1.1),有

()()()22

,,,,j k j k j k u x t u x t a

f x t t

x

??=+?? 11,j M ≤≤- 11k N ≤≤-. (2.1)

将()1,j k u x t +以结点()

,j k x t 为中心关于t 运用泰勒级数展开, 有

()()()()2

21,,,,().2!

j k j k j k j k u x t u x t u x t u x t o τττ+'''=++

+

整理有

()()

()()212

,,,,().2j k j k j k j k u x t u x t u x t u x t o t t

τττ

+-??=++?? (2.2) 再将()1,j k x t -, ()1,j k x t +分别以结点()

,j k x t 为中心关于x 运用泰勒级数展开, 有

()()()()()

()2

3

1,,,=,,(2!

3!

j k j k j k j k j k u x t h u x t h u x t u x t u x t h -'''''-'+-+

+

（-）

）

()()

4

(4)4,()4!

j k u x t h o h -+

+, ()()()()()2

3

1,,,=,,2!3!

j k j k j k j k j k u x t h u x t h u x t u x t u x t h +''''''++

+

()(4)4

4,().4!

j k u x t h o h +

+

由上述两式可得

()()()()()242112

224

,2,,,,=()12j k j k j k j k j k u x t u x t u x t u x t x t h o h h x x

-+-+??++??. (2.3) 将(2.2), (2.3)两式代入(2.1)中, 得

()()

()()()

()1112

,-u ,,2,,,j k j k j k j k j k k j k j u x t x t u x t u x t u x t a

f x t R h

τ

+-+-+=++. (2.4)

其中()()4222242

,,()122j k

j k k j

u x t x t ah R o h x t ττ??=-+++?? 为方程(2.1)的截断误差. 舍去截断误差, 用k

j u 代替()

,j k u x t , 得到如下差分方程

111

2

2,k k k k k

j j

j j j k j u u u u u a

f h τ

+-+--+=+ 11,j M ≤≤- 0 1.k N ≤≤- (2.5)

结合初边值条件, 可得如下差分格式

111

2

2,k k k k k

j j

j j j k j u u u u u a

f h τ

+-+--+=+ 11,j M ≤≤- 11,k N ≤≤- (2.6)

0(),j j u x ?= 0,j M ≤≤ (2.7)

0(),k j u t α= (),k M k u t β= 1.k N ≤≤ (2.8)

2.2 向后Euler 差分格式

在结点()

1,j k x t +处考虑方程(1.1), 有

()()()21112

,,,,j k j k j k u x t u x t a

f x t t

x +++??=+?? 11,j M ≤≤- 1.k N ≤≤ (2.9)

将(

),j k u x t 以()

1,j k x t +为中心关于t 运用泰勒级数展开, 有

()()()()2

1211,(),,,()().2!

j k j k j k j k u x t u x t u x t u x t o τττ+++''-'=+-+

+

将上式整理得

()()

()()

21112

,,,,().2j k j k j k j k u x t u x t u x t u x t o t t

τττ

+++-??=-+?? (2.10) 再将()11,j k u x t -+, ()11,j k u x t ++分别以()

1,j k x t +为中心关于x 运用泰勒级数展开, 有

()()()()()

()2

3

111111,,,=,,(2!

3!

j k j k j k j k j k u x t h u x t h u x t u x t u x t h ++-+++'''''-'+-+

+

（-）

）

()()

4

(4)14,(),4!

j k u x t h o h +-+

+ ()()()()()2

3

111111,,,=,,2!3!

j k j k j k j k j k u x t h u x t h u x t u x t u x t h ++++++''''''++

+

()(4)4

14,().4!

j k u x t h o h ++

+

由上述两式可得

()()()()()

24211111112224

,2,,,,=()12j k j k j k j k j k u x t u x t u x t u x t x t h o h h x x

-++++++-+??++??. (2.11) 将(2.10), (2.11)两式代入(2.9)中, 得

()()

()()()

1111112

,,,2,,j k j k j k j k j k u x t u x t u x t u x t u x t a

h τ

+-++++--+=

()1,.k j k j f x t R +++ (2.12)

其中()()

24

21122

4

,,()2

12

i k i k k

j

x t u x t ah R o h t x ττ++??=-

-

++?? 为方程(2.9)的截断误差.

舍去截断误差, 用k

j u 代替()

,j k u x t ,得到如下差分方程

1111

11

12

2,k k k k k j j

j j j k j u u u u u a

f h

τ

++++-++--+=+ 11,j M ≤≤- 1.k N ≤≤ (2.13)

结合初边值条件, 可得如下差分格式

1

111

11

12

2,k k k k k j j

j j j k j u u u u u a

f h

τ

++++-++--+=+ 11,j M ≤≤- 1,k N ≤≤ (2.14)

0(),j j u x ?= 0,j M ≤≤ (2.15)

0(),k k u t α= (),k M k u t β= 1.k N ≤≤ (2.16)

2.3 Crank Nicolson -差分格式 在结点()

1/2,j k x t +处考虑方程(1.1), 有

()()()21/21/21/22

,,,,j k j k j k u x t u x t a

f x t t

x

+++??=+?? 11,j M ≤≤- 0 1.k N ≤≤- (2.17)

将()1,j k u x t +, (

),j k u x t 以()

1/2,j k x t +为中心关于t 运用泰勒级数展开, 有

()()()()2

1/211/21/2,,,,22!2j k j k j k j k u x t u x t u x t u x t τ

τ++++''??'=++ ???

()()3

1/23

,.3!2j k u x t o ττ+'''??++ ???

()()()()2

1/21/21/2,,,,22!2j k j k j k j k u x t u x t u x t u x t ττ+++''?

???'=+-+- ? ?????

()()3

1/23

,.3!2j k u x t o ττ+'''??+-+ ???

将上述两式整理得

()()

()()3211/21/233

,,,,().24j k j k j k j k u x t u x t u x t u x t o t t

τττ

+++-??=++?? (2.18) 再将()1,j k u x t -, ()1,j k u x t +分别以()

,j k x t 为中心关于x 运用泰勒级数展开, 有

()()()()()

2

3

1,(,),=,,(2!

3!

j k j k j k j k j k u x t h u x t h u x t u x t u x t h -''-''''+-+

+

（-）

）

()()

4

(4)4,()4!j k u x t h o h -+

+, ()()()()()()2

3

(4)4

41,,,,=,,().2!

3!

4!

j k j k j k j k j k j k u x t h u x t h u x t h u x t u x t u x t h o h +''''''++

++

+

由上述两式得

()()()()()242112

224

,2,,,,=().12j k j k j k j k j k u x t u x t u x t u x t x t h o h h x x

-+-+??++?? (2.19) 同理, 将()11,j k u x t -+, ()11,j k u x t ++分别以()

1,j k x t +为中心关于x 泰勒级数展开, 整理得

()()()()()

24211111112224

,2,,,,=().12j k j k j k j k j k u x t u x t u x t u x t x t h o h h x x

-++++++-+??++?? (2.20) 此时，分别将()1,j k u x t +, ()

,j k u x t 以()

1/2,j k u x t +为中心关于t 泰勒级数展开，有

()()()()()22342

11/21/21/22

22222

,,,,1=,22!2j k j k j k j k u x t u x t u x t u x t o x x x t x t τττ++++??????+++ ?????????

()()()2231/21/2222

,,,=2j k j k j k u x t u x t u x t x x x t τ++?????

+- ???????

()()42

1/2222,1.2!2j k u x t o x t ττ+???+-+ ?????

利用上述两式得

()()()()()22242

1/211/22

2

2222

,,,,11.222j k j k j k j k u x t u x t u x t u x t o x x x x t ττ+++?????????? ? ?=+-+ ? ? ???????????

? (2.21)

利用(2.19), (2.20)两式, 整理有

()()()()()

221112

2

2

,,,2,,j k j k j k j k j k u x t u x t u x t u x t u x t x

x

h

+-+??-++

=

??

()()()

()4

21111124

,2,,,12j k j k j k j k u x t u x t u x t x t h h x -++++-+?+

-

? ()()4

2124,.12j k x t h o h x

+?-+? (2.22) 结合(2.21),(2.22)两式, 整理得

()()()()

21/2112

2

,,2,,2j k j k j k j k u x t u x t u x t u x t x

h

+-+?-+=

?

()()()

()42

111111/22

22,2,,,1222j k j k j k j k u x t u x t u x t u x t h x t τ-+++++??

-+??? ?+

- ? ??????

? ()()()44

2212

44

,,1.21212j k j k x t x t h h o h x x +???? ?-++ ?????

(2.23) 将(2.18), (2.23)式代入(2.17)得

()()

()()()

1112

,,,2,,2j k j k j k j k j k u x t u x t u x t u x t u x t a

h τ

+-+--+=

()()()

()111112

,2,,,.2j k j k j k k j k j u x t u x t u x t a

f x t R h -++++-++++ (2.24)

其中

()()()()4443222211/21/244223

,,,,21212224j k j k j k j k k j x t x t u x t u x t a

h h R x x x t t ττ+++???????? ?=-+++ ? ?????????? 32()o h τ++为方程(2.17)的截断误差.

舍去截断误差, 用k

j u 代替(,)j k u x t , 可得如下差分方程

1

111

11

11

1/22

2

22,22k k k k k

k k k j j

j j j j j j k j u u u u u u u u a

a

f h h τ

++++-+-++--+-+=++

11,j M ≤≤- 0 1.k N ≤≤- (2.24)

结合初边值条件, 可得如下差分格式

1

111

11

11

1/22

2

22,22k k k k k

k k k j j

j j j j j j k j u u u u u u u u a

a

f h h τ

++++-+-++--+-+=++

11,j M ≤≤- 01,k N ≤≤- (2.25)

0(),j j u x ?= 0,j M ≤≤ (2.26)

0(),k k u t α= (),k M k u t β= 1.k N ≤≤ (2.27)

3 差分格式的求解

3.1 向前Euler 格式 记2

a r h τ

= , 称r 为步长比. 差分格式(2.6) (2.8)中(2.6)可改写为

111(12)k k k k k

j j j j j u r u r u r u f τ+-+=?+-+?+? ,

11j M ≤≤- , 01k N ≤≤- .

(3.1)

将(3.1)写成如下形式

11k k k U A U F r F τ+=?+?+? .

其中

+1111121=,,,T

k k k k M U u u u +++-???? , 121=,,,T

k k k k

M U u u u -???

? , 121,,,T k k k k

M F f f f -??=?? , 10

(1)*1,0,0,,0,T

k k M M F u u -??=?? , (1)*(1)

121212M M r r r r r A r r r r ---??

?- ?

?= ?

?

?-??

. 3.2 向后Euler 格式 记2

a r h τ

= , 称r 为步长比. 差分格式(2.14) (2.16)中(2.14)可改写为

1111

11(12)k k k k k j j j j j

u r u r u r u f τ++++-+-=?-++?+? ,

11j M ≤≤- , 1k N ≤≤.

(3.2)

将(3.2)写成如下形式

111k k k A U U F r F τ++?=--?-? .

其中

+1111121=,,,T k k k k M U u u u +++-???? , 121=,,,T

k k k k

M U u u u -???

? , 1111

121,,,T k k k k M F f f f ++++-??=?? , 1110(1)*1

,0,0,,0,T k k M M F u u ++-??=??

,

(1)*(1)

(12)(12)(12)M M r r r r r A r r r r ---+??

?

-+ ?

?= ?

?

?-+??

. 3.3 Crank Nicolson -格式 记2

a r h τ

= , 称r 为步长比. 差分格式(2.25) (2.27)中(2.25)可改写为

1111/2

1111-[/2(1)/2]/2(1)/2k k k k k k k j j j j j j j

r u r u r u r u r u r u f τ++++-+-+?-++?=?+-+?+? ,

11j M ≤≤- , 01k N ≤≤- .

(3.3)

将(3.3)写成如下形式

11/22112-(k k k A U A U F r F F τ++?=?+?++）

. 其中

+1111121=,,,T

k k k k M U u u u +++-???? , 121=,,,T

k k k k

M U u u u -???

? , 10(1)*1,0,0,,0,T k k M M F u u -??=?? , 11

20(1)*1,0,0,,0,T

k k M M F u u ++-??=?? ,

1/2

1/21/21/2

1

21,,,T

k k k k M F

f f f ++++-??=?? , 1(1)*(1)

1/2/21/2/2/2/21M M r r r r r A r r r r ---??

?- ?

?= ?

?

?-??

, 2(1)*(1)

(1)/2/2(1)/2/2/2/2(1)M M r r r r r A r r r r ---+?? ?-+ ?

?= ?

?

?-+??

. 4 数值例子

在本文中, 将应用向前Euler 格式(2.6) (2.8), 向后Euler 格式(2.14) (2.16)和

Crank Nicolson -格式(2.25) (2.27)分别来求解如下算例.

算例 现有如下定解问题

222,x t

u u e t x

+??=+?? 01,x << 01,t <≤ (,0),x u x e = 01,x ≤≤

2(0,),t u t e = 12(1,),t u t e += 0 1.t <≤

上述定解问题的精确解为2(,)x t

u x t e +=.

4.1 向前Euler 格式

在表4.1.1中给出了取步长1/10h =和1/200τ=(步长比1/2r =)时计算得到的部分数值结果. 表4.1.2给出了取步长1/10h =和1/100τ=(步长比1r =)时计算得到的部分数值结果, 随着计算层数的增加, 误差越来越大, 数值结果是发散的. 经步长比r 的多次取值发现, 当1/2r ≤时, 其数值结果是收敛的, 当1/2r >时, 其数值结果是发散的.表4.1.3给出了1/2r =时, 取不同步长, 数值解的最大误差

111(,)max |(,)|.k j k j j M k N

E h u x t u τ∞≤≤-≤≤=-

从表4.1.3可以看出当空间步长缩小到原来的1/2, 时间步长缩小到原来的1/4时, 最大误差约缩小到原来的1/4.

图4.1.1给出了1t =时的数值解曲线(1/10,1/200)h τ==与精确解曲线. 图4.1.2给出了1t =时不同步长数值解的误差曲线图(1/2)r =, 由图4.1.2可知, 当步长比r 保持不变, 随着时间与空间方向上的加密, 误差越来越小. 图4.1.3给出了不同步长的误差曲面图.

表4.1.1 算例在部分结点处数值解、精确解和绝对误差(1/10,1/200h τ==)

k

(,)x t

数值解 精确解 |精确解-数值解|

20

(0.5,0.1)

1.802810

1.803988

1.178e-3

40 (0.5,0.2) 2.201438 2.203396 1.959e-3 60 (0.5,0.3) 2.688651 2.691234 2.583e-3 80 (0.5,0.4) 3.283856 3.287081 3.225e-3 100 (0.5,0.5) 4.010886 4.014850 3.965e-3 120 (0.5,0.6) 4.898897 4.903749 4.852e-3 140 (0.5,0.7) 5.983523 5.989452 5.929e-3 160 (0.5,0.8) 7.308290 7.315534 7.243e-3 180 (0.5,0.9) 8.926366 8.935213 8.847e-3 200 (0.5,1.0)

10.902687

10.913494

1.081e-3

表4.1.2 算例在部分结点处数值解、精确解和绝对误差(1/10,1/100h τ==)

k

(,)x t

数值解 精确解 |精确解-数值解|

1 (0.5,0.01) 1.521674 1.52196

2 2.879e-4 2 (0.5,0.02) 1.55212

3 1.552707 5.845e-

4 3 (0.5,0.03) 1.583184 1.584074 8.901e-4 4 (0.5,0.04) 1.614870 1.616074 1.205e-3

5 (0.5,0.05) 1.64738

6 1.648721 1.336e-3 6 (0.5,0.06) 1.679785 1.682028 2.242e-3

7 (0.5,0.07) 1.716057 1.716007 5.051e-5

8 (0.5,0.08) 1.741446 1.750673 9.227e-3

9 (0.5,0.09) 1.807089 1.786038 2.105e-2 10 (0.5,0.10) 1,744143 1.822119 7.798e-2 11 (0.5,0.11) 2.092546 1.858928 2.336e-1 12 (0.5,0.12) 1.164923 1.896481 7.316e-1 13 (0.5,0.13) 4.148393 1.934792 2.213e+0 14 (0.5,0.14) -4.709415 1.973878 6.684e+0 15 (0.5,0.15) 21.998093 2.013753 1.998e+1 16 (0.5,0.16) -57.422399 2.054433 5.948e+1 17

(0.5,0.17)

178.294623

2.095936

1.762e+2

18 (0.5,0.18) -518.143503 2.138276 5.203e+2

表4.1.3 算例取不同步长时数值解的最大误差(1/2r =)

h

τ

(,)E h τ∞

(2,4)/(,)E h E h ττ∞∞

1/10 1/200 1.178e-2 * 1/20 1/800 2.973e-3 3.964 1/40 1/3200 7.431e-4 4.000 1/80

1/12800

1.858e-4

4.000

图4.1.1 算例取1t =时的数值解曲线(1/10,1/200)h τ==与精确解曲线

图4.1.2 算例取1t =时不同步长数值解的误差曲线(1/2)r =

图4.1.3 算例取不同步长数值解的误差曲面(1/2)r =

4.2 向后Euler 格式

在表 4.2.1和表 4.2.2中分别给出了1/10,1/200h τ==(步长比1/2r =)和

1/10,h =1/100τ=(步长比1r =)时计算得到的部分数值结果, 发现两表的数值结果都

是收敛的, 经步长比r 的多次取值发现, 无论r 取任何值, 其数值结果都是收敛的. 表4.2.3 给出了1/2r =时, 取不同步长, 数值解的最大误差

111(,)max |(,)|.k j k j j M k N

E h u x t u τ∞≤≤-≤≤=-

从表4.2.3可以看出当空间步长缩小到原来的1/2, 时间步长缩小到原来的1/4时, 最大误差缩小到原来的1/4. 表4.2.4给出了1r =时的类似结论.

图4.2.1给出了1t =时的数值解曲线(1/10,1/200)h τ==与精确解曲线. 图4.2.2给出了1t =时, 取不同步长数值解的误差曲线图(1/2)r =, 由图4.2.2可知, 当步长比r 保持不变, 随着时间与空间方向上的加密, 误差越来越小. 图4.2.3和图4.2.4分别当给出了步长比1/2r =和1r =下不同步长的误差曲面图.

表4.2.1 算例在部分结点处数值解、精确解和绝对误差(1/10,1/200h τ==)

k

(,)x t

数值解 精确解 |精确解-数值解|

20 (0.5,0.1) 1.805343 1.803988 1.355e-3 40 (0.5,0.2) 2.205674 2.203396 2.278e-3 60 (0.5,0.3) 2.694258 2.691234 3.023e-3 80 (0.5,0.4) 3.290867 3.287081 3.785e-3 100 (0.5,0.5) 4.019509 4.014850 4.659e-3 120 (0.5,0.6) 4.909453 4.903749 5.704e-3 140 (0.5,0.7) 5.996425 5.989452 6.973e-3 160 (0.5,0.8) 7.324052 7.315534 8.518e-3 180 (0.5,0.9) 10.926203 8.935213 1.041e-2 200 (0.5,1.0)

10.926203

10.913494

1.271e-2

表4.2.2 算例在部分结点处数值解、精确解和绝对误差(1/10,1/100h τ==)

k

(,)x t

数值解 精确解 |精确解-数值解|

10 (0.5,0.1) 1.788500 1.786038 2.461e-3 20 (0.5,0.2) 2.185742 2.181472 4.269e-3 30 (0.5,0.3) 2.670171 2.664456 5.715e-3 40 (0.5,0.4) 3.261551 3.254374 7.177e-3 50 (0.5,0.5) 3.983745 3.974902 8.843e-3 60 (0.5,0.6) 4.865788 4.854956 1.083e-2 70 (0.5,0.7) 5.943098 5.929856 1.324e-2 80 (0.5,0.8) 7.258921 7.242743 1.618e-2 90 (0.5,0.9) 8.866068 8.846306 1.976e-2 100 (0.5,1.0)

10.829041

10.804903

2.414e-2

表4.2.3 算例取不同步长时数值解的最大误差(1/2r =)

h

τ

(,)E h τ∞

(2,4)/(,)E h E h ττ∞∞

1/10 1/200 1.386e-2 * 1/20 1/800 3.509e-3 3.950 1/40 1/3200 8.779e-4 3.997 1/80

1/12800

2.195e-4

3.999

表4.2.4 算例取不同步长时数值解的最大误差(1r =)

h

τ

(,)E h τ∞

(2,4)/(,)E h E h ττ∞∞

1/10 1/100 2.659e-2 * 1/20 1/400 6.744e-3 3.943 1/40

1/1600

1.688e-3

3.955

1/80

1/6400 4.222e-4 3.999

图4.2.1 算例取1t =时的数值解曲线(1/10,1/200)h τ== 与精确解曲线

图4.2.2 算例取1t = 时不同步长数值解的误差曲线(1/2)r =

图 4.2.3 算例取不同步长数值解的误差曲面(1/2)r =

图 4.2.4 算例取不同步长数值解的误差曲面(1)r =

4.3 Crank Nicolson -格式

在表4.3.1中给出了取步长1/10,1/10h τ==时计算得到的部分数值结果, 表4.3.2给出了取步长1/100,1/100h τ==时计算得到的部分数值结果, 发现两表的数值结果都是收敛的, 经步长比r 的多次取值发现, 无论r 取任何值, 其数值结果都是收敛的. 表4.3.3给出了取不同步长时, 所得数值解的最大误差

111(,)max |(,)|.k j k j j M k N

E h u x t u τ∞≤≤-≤≤=-

从表可以看出当空间步长缩小到原来的1/2, 时间步长缩小到原来的1/2时, 最大误差缩小

到原来的1/4.

图4.3.1给出了1t =时的数值解曲线(1/10,1/10)h τ==与精确解曲线. 图4.3.2给出了1t =时, 取不同步长数值解的误差曲线图, 由图 4.3.2可知, 当步长比r 保持不变, 随着时间与空间方向上的加密, 误差越来越小. 图4.3.3给出了取不同步长时所得数值解的误差曲面图. 图4.3.4给出了向前Euler , 向后Euler 和Crank Nicolson -三种方法数值解的误差对比图(1/10,1/200)h τ==, 由图4.3.4分析可得, 当时间方向t 与空间方向

x 都取相同步长时, Crank Nicolson -方法误差最小, 精度最高, 向前Euler 方法次之,

向后Euler 方法所得误差最大, 精度最低.

表4.3.1 算例在部分结点处数值解、精确解和绝对误差(1/10,1/10h τ==)

k

(,)x t

数值解 精确解 |精确解-数值解|

1 (0.5,0.1) 1.823114 1.822119 9.949e-4

2 (0.5,0.2) 2.227196 2.225541 1.656e-

3 3 (0.5,0.3) 2.72041

4 2.718282 2.132e-3 4 (0.5,0.4) 3.322776 3.320117 2.659e-3

5 (0.5,0.5) 4.058460 4.055200 3.260e-3

6 (0.5,0.6) 4.957021 4.953032 3.988e-3

7 (0.5,0.7) 6.054520 6.049647 4.873e-3

8 (0.5,0.8) 7.395008 7.389056 5.952e-3

9 (0.5,0.9) 9.032284 9.025013 7.271e-3 10

(0.5,1.0)

11.032056

11.023176

8.880e-3

表4.3.2 算例在部分结点处数值解、精确解和绝对误差(1/100,1/100h τ==)

k

(,)x t

数值解 精确解 |精确解-数值解|

10 (0.5,0.1) 1.993726 1.993715 1.081e-5 20 (0.5,0.2) 2.435147 2.435130 1.749e-5 30 (0.5,0.3) 2.974297 2.974297 2.296e-5 40 (0.5,0.4) 3.632815 3.632786 2.864e-5 50 (0.5,0.5) 4.437131 4.437096 3.521e-5 60 (0.5,0.6) 5.919524 5.419481 4.308e-5 70 (0.5,0.7) 6.619421 6.619369 5.266e-5 80 (0.5,0.8) 8.084979 8.084915 6.433e-5 90 (0.5,0.9) 9.875016 9.874938 7.857e-5 100 (0.5,1.0)

12.061372

12.061276

9.597e-5

表4.3.3 算例取不同步长时数值解得最大误差

h

τ

(,)E h τ∞

(2,2)/(,)E h E h ττ∞∞

1/10 1/10 9.588e-3 * 1/20 1/20 2.429e-3 3.9457 1/40 1/40 6.078e-4 3.9961 1/80 1/80 1.520e-4 3.9990 1/160 1/160 3.799e-5 3.9998 1/320 1/320 9.499e-6 3.9999 1/640

1/640

2.375e-6

4.0000

图4.3.1 算例取1t =时数值解曲线(1/10,1/10)h τ==与精确解曲线

t 时不同步长数值解的误差曲线图4.3.2 算例取1

图4.3.3 算例取不同步长数值解的误差曲面

一维热传导方程的差分格式

《微分方程数值解》 课程论文 学生姓名1：许慧卿学号：20144329 学生姓名2：向裕学号：20144327学生姓名3：邱文林学号：20144349学生姓名4：高俊学号：20144305学生姓名5：赵禹恒学号：20144359学生姓名6：刘志刚学号： 20144346 学院：理学院 专业：14级信息与计算科学 指导教师：陈红斌 2017年6 月25日

《偏微分方程数值解》课程论文 《一维热传导方程的差分格式》论文 一、《微分方程数值解》课程论文的格式 1）引言：介绍研究问题的意义和现状 2）格式：给出数值格式 3）截断误差：给出数值格式的截断误差 4）数值例子：按所给数值格式给出数值例子 5）参考文献：论文所涉及的文献和教材 二、《微分方程数值解》课程论文的评分标准 1）文献综述：10分； 2）课题研究方案可行性：10分； 3）数值格式：20分； 4）数值格式的算法、流程图：10分； 5）数值格式的程序：10分； 6）论文撰写的条理性和完整性：10分； 7）论文工作量的大小及课题的难度：10分； 8）课程设计态度：10分； 9）独立性和创新性：10分。 评阅人： - 2 -

一维热传导方程的差分格式 1 引言 考虑如下一维非齐次热传导方程Dirichlet 初边值问题 22(,),u u a f x t t x ??=+?? ,c x d << 0,t T <≤ (1.1) (,0)(),u x x ?= ,c x d ≤≤ (1.2) (,)(),u c t t α= (,)(),u d t t β= 0t T <≤ (1.3) 的有限差分方法, 其中a 为正常数,(,),(),(), ()f x t x t t ?αβ为已知常数, ()(0),c ?α= ()(0).d ?β= 称(1.2)为初值条件, (1.3)为边值条件. 本文将给出(1.1) (1.3)的向前Euler 格式, 向后Euler 格式和Crank Nicolson -格式, 并给出其截断误差和数值例子. 经对比发现, Crank Nicolson -格式误差最小, 向前 Euler 格式次之, 向后Euler 格式误差最大. 2 差分格式的建立 2.1 向前Euler 格式 将区间[,]c d 作M 等分, 将[]0,T 作N 等分, 并记 ()/h d c M =-, /T N τ=, j x c jh =+,0j M ≤≤, k t k τ=,0k N ≤≤. 分别称h 和τ为空间步长和时间步长.用 两组平行直线 j x x =, 0j M ≤≤, k t t =, 0k N ≤≤ 将Ω分割成矩形网格.记{} |0h j x j M Ω=≤≤, {}|0k t k N τΩ=≤≤, h h ττΩ=Ω?Ω. 称() ,j k x t 为结点[1] . 定义h τΩ上的网格函数 {}|0,0k j U j M k N Ω=≤≤≤≤, 其中() ,k j j k U u x t =. 在结点() ,j k x t 处考虑方程(1.1),有

一维热传导方程

一维热传导方程 一． 问题介绍 考虑一维热传导方程： （1） ,0),(22T t x f x u a t u ≤<+??=?? 其中a 是正常数，)(x f 是给定的连续函数。按照定解条件的不同给法，可将方程（1）的定解问题分为两类： 第一类、初值问题（也称Cauthy 问题）：求具有所需次数偏微商的函数),(t x u ，满足方程（1）（∞<<∞-x ）和初始条件： （2） ),()0,(x x u ?= ∞<<∞-x 第二类、初边值问题（也称混合问题）：求具有所需次数偏微商的函数),(t x u ，满足方程（1）（l x <<0）和初始条件： （3） ),()0,(x x u ?= l x <<0 及边值条件 (4) .0),(),0(==t l u t u T t ≤≤0 假定)(x ?在相应区域光滑，并且在l x ,0=满足相容条件，使上述问题有唯一充分光滑的解。 二． 区域剖分 考虑边值问题（1），（4）的差分逼近。去空间步长N l h /=和时间步长M T /=τ，其中N,M 都是正整数。用两族平行直线： 将矩形域}0;0{T t l x G ≤≤≤≤=分割成矩形网格，网格节点为),(k j t x 。以h G 表示网格内点集合，即位于开矩形G 的网点集合；h G 表示所有位于闭矩形G 的网点集合；h Γ=h G --h G 是网格界点集合。 三． 离散格式 第k+1层值通过第k 层值明显表示出来，无需求解线性代数方程组，这样的格式称为显格式。 第k+1层值不能通过第k 层值明显表示出来，而由线性代数方程组确定，这样的格式称为隐格式。 1. 向前差分格式 （5） ,221 11j k j k j k j k j k j f h u u u a u u ++-=--++τ

【文献综述】热传导方程差分格式的收敛性和稳定性

文献综述 信息与计算科学 热传导方程差分格式的收敛性和稳定性在实际研究物理问题过程中, 往往能给出问题相应的数学表达式, 但是由于实际物理问题的复杂性, 它的解却一般不容易求出. 由此计算物理应运而生, 计算物理是以计算机为工具, 应用数学的方法解决物理问题的一门应用性学科, 是物理、数学和计算机三者结合的交叉性学科. 它产生于二战期间美国对核武器的研究, 伴随着计算机的发展而发展. 计算物理的目的不仅仅是计算, 而是要通过计算来解释和发现新的物理规律. 这一点它与传统的实验物理和理论物理并无差别, 所不同的只是使用的工具和方法. 计算物理早已与实验物理和理论物理形成三足鼎立之势, 甚至有人提出它将成为现代物理大厦的“栋梁”. 在一个物理问题中一个数值解往往比一个式子更直观, 更有价值. 在实际求解方程时, 除了一些特殊的情况下可以方便地求得其精确解外, 在一般情况下, 当方程或定解条件具有比较复杂的形式, 或求解区域具有比较复杂的形状时, 往往求不到, 或不易求到其精确解. 这就需要我们去寻找方程的近似解, 特别是数值近似解, 简称数值解. 这里主要研究的是热传导方程. 有限差分法是微分方程和积分微分方程数值解的方法. 其基本思想是把连续的定解区域用有限个离散点构成的网格来代替, 这些离散点称作网格的节点;把连续定解区域上的连续变量的函数用在网格上定义的离散变量函数来近似;把原方程和定解条件中的微商用差商来近似, 积分用积分和来近似, 于是原微分方程和定解条件就近似地代之以代数方程组, 即有限差分方程组, 解此方程组就可以得到原问题在离散点上的近似解. 然后再利用插值方法便可以从离散解得到定解问题在整个区域上的近似解. 热传导的差分法是求解热传导方程的重要方法之一. 对于差分格式的的求解, 我们首先要关注差分格式的收敛性和稳定性. 对于一个微分方程建立的各种差分格式, 为了有实用意义, 一个基本要求是它们能够任意逼近微分方程, 即相容性要求. 一个差分格式是否有用, 就要看差分方程的精确解能否任意逼近微分方程的解, 即收敛性的概念. 此外, 还有一个重要的概念必须考虑, 即差分格式的稳定性. 因为差分格式的计

热传导方程向后差分格式的MATLAB程序

向后差分格式MATLAB编程： c lear;clc; format short e a=input('请输入系数a的值'); l=input('请输入长度l的值'); M=input('请输入将区间[0,1]等分的个数M '); ot=input('请输入时间增量ot的值'); n=input('请输入运行次数n的值'); ox=1/M; x0=zeros(M+1,1) for ii=1:M x0(ii+1)=ii*ox; end u=sin(pi*x0/l); r=a*ot/(ox)^2; for ii=1:n %数据的输入 B=zeros(M-1,1); A=zeros(M-2,1); C=zeros(M-2,1); S=zeros(M-1,1); for ii=1:M-2 B(ii)=1+2*r;A(ii)=-r;C(ii)=-r; S(ii)=u(ii+1,1); end B(M-1,1)=1+2*r;S(M-1,1)=u(M,1);u(1,2)=0;u(M+1,2)=0; S(1,1)=S(1,1)+r*u(1,2);S(M-1,1)=S(M-1,1)+r*u(M+1,2); %追赶法 S(1)=S(1)/B(1);T=B(1);k=2; while k~=M B(k-1)=C(k-1)/T; T=B(k)-A(k-1)*B(k-1); S(k)=(S(k)-A(k-1)*S(k-1))/T; k=k+1 end k=1; while k~=M-1 S(M-1-k)=S(M-1-k)-B(M-1-k)*S(M-k); k=k+1; end u(2:M,2)=S; u(:,1)=u(:,2); end %计算精确解 for x=0:M

一维导热方程 有限差分法 matlab实现

第五次作业（前三题写在作业纸上） 一、用有限差分方法求解一维非定常热传导方程，初始条件和边界条件见说明.pdf 文件，热扩散系数α=const ， 22T T t x α??=?? 1. 用Tylaor 展开法推导出FTCS 格式的差分方程 2. 讨论该方程的相容性和稳定性，并说明稳定性要求对求解差分方程的影响。 3. 说明该方程的类型和定解条件，如何在程序中实现这些定解条件。 4. 编写M 文件求解上述方程，并用适当的文字对程序做出说明。（部分由网络搜索得到，添加，修改后得到。） function rechuandaopde %以下所用数据，除了t 的范围我根据题目要求取到了20000，其余均从pdf 中得来 a=0.00001;%a 的取值 xspan=[0 1];%x 的取值范围 tspan=[0 20000];%t 的取值范围 ngrid=[100 10];%分割的份数，前面的是t 轴的，后面的是x 轴的 f=@(x)0;%初值 g1=@(t)100;%边界条件一 g2=@(t)100;%边界条件二 [T,x,t]=pdesolution(a,f,g1,g2,xspan,tspan,ngrid);%计算所调用的函数 [x,t]=meshgrid(x,t); mesh(x,t,T);%画图，并且把坐标轴名称改为x ，t ，T xlabel('x') ylabel('t') zlabel('T') T%输出温度矩阵 dt=tspan(2)/ngrid(1);%t 步长 h3000=3000/dt;

h9000=9000/dt; h15000=15000/dt;%3000,9000,15000下，温度分别在T矩阵的哪些行T3000=T(h3000,:) T9000=T(h9000,:) T15000=T(h15000,:)%输出三个时间下的温度分布 %不再对三个时间下的温度-长度曲线画图，其图像就是三维图的截面 %稳定性讨论,傅里叶级数法 dx=xspan(2)/ngrid(2);%x步长 sta=4*a*dt/(dx^2)*(sin(pi/2))^2; if sta>0,sta<2 fprintf('\n%s\n','有稳定性') else fprintf('\n%s\n','没有稳定性') error end %真实值计算 [xe,te,Te]=truesolution(a,f,g1,g2,xspan,tspan,ngrid); [xe,te]=meshgrid(xe,te); mesh(xe,te,Te);%画图，并且把坐标轴名称改为xe，te，Te xlabel('xe') ylabel('te') zlabel('Te') Te%输出温度矩阵 %误差计算 jmax=1/dx+1;%网格点数 [rms]=wuchajisuan(T,Te,jmax) rms%输出误差

热传导方程向前差分格式的MATLAB程序

向前差分格式MATLAB编程： c lear;clc; format short e a=input('请输入系数a的值'); l=input('请输入长度l的值'); M=input('请输入将区间[0,1]等分的个数M '); ot=input('请输入时间增量ot的值'); n=input('请输入运行次数n的值'); ox=1/M; x0=zeros(M+1,1) for ii=1:M x0(ii+1)=ii*ox; end u=sin(pi*x0/l); r=a*ot/(ox)^2; for ii=1:n %数据的输入 B=zeros(M-1,1); A=zeros(M-2,1); C=zeros(M-2,1); S=zeros(M-1,1); for ii=1:M-2 B(ii)=1+2*r;A(ii)=-r;C(ii)=-r; S(ii)=u(ii+1,1); end B(M-1,1)=1+2*r;S(M-1,1)=u(M,1);u(1,2)=0;u(M+1,2)=0; S(1,1)=S(1,1)+r*u(1,2);S(M-1,1)=S(M-1,1)+r*u(M+1,2); %追赶法 S(1)=S(1)/B(1);T=B(1);k=2; while k~=M B(k-1)=C(k-1)/T; T=B(k)-A(k-1)*B(k-1); S(k)=(S(k)-A(k-1)*S(k-1))/T; k=k+1 end k=1; while k~=M-1 S(M-1-k)=S(M-1-k)-B(M-1-k)*S(M-k); k=k+1; end D=(1-2*r)*eye(M-1); temp=r*linspace(1,1,M-2); D=D+diag(temp,1)+diag(temp,-1); S=D*S

一维热传导方程

一维热传导方程 一． 问题介绍 考虑一维热传导方程： （1） ,0),(22 T t x f x u a t u ≤<+??=?? 其中a 是正常数，)(x f 是给定的连续函数。按照定解条件的不同给法，可将方程（1）的定解问题分为两类： 第一类、初值问题（也称Cauthy 问题）：求具有所需次数偏微商的函数),(t x u ，满足方 程（1）（∞<<∞-x ）和初始条件： （2） ),()0,(x x u ?= ∞<<∞-x 第二类、初边值问题（也称混合问题）：求具有所需次数偏微商的函数),(t x u ，满足方 程（1）（l x <<0）和初始条件： （3） ),()0,(x x u ?= l x <<0 及边值条件 (4) .0),(),0(==t l u t u T t ≤≤0 假定)(x ?在相应区域光滑，并且在l x ,0=满足相容条件，使上述问题有唯一充分光滑 的解。 二． 区域剖分 考虑边值问题（1），（4）的差分逼近。去空间步长N l h /=和时间步长M T /=τ，其中N,M 都是正整数。用两族平行直线： ),,1,0(N j jh x x j === ),,1,0(M k k t t k ===τ 将矩形域}0;0{T t l x G ≤≤≤≤=分割成矩形网格，网格节点为),(k j t x 。以h G 表示网格内点集合，即位于开矩形G 的网点集合；h G 表示所有位于闭矩形G 的网点集合； h Γ=h G --h G 是网格界点集合。 三． 离散格式 第k+1层值通过第k 层值明显表示出来，无需求解线性代数方程组，这样的格式称为

一维热传导方程

一维热传导方程Last revision on 21 December 2020

一维热传导方程 一． 问题介绍 考虑一维热传导方程： （1） ,0),(22T t x f x u a t u ≤<+??=?? 其中a 是正常数，)(x f 是给定的连续函数。按照定解条件的不同给法，可将方程（1）的定解问题分为两类： 第一类、初值问题（也称Cauthy 问题）：求具有所需次数偏微商的函数),(t x u ，满足方程（1）（∞<<∞-x ）和初始条件： （2） ),()0,(x x u ?= ∞<<∞-x 第二类、初边值问题（也称混合问题）：求具有所需次数偏微商的函数),(t x u ，满足方程（1）（l x <<0）和初始条件： （3） ),()0,(x x u ?= l x <<0 及边值条件 (4) .0),(),0(==t l u t u T t ≤≤0 假定)(x ?在相应区域光滑，并且在l x ,0=满足相容条件，使上述问题有唯一充分光滑的解。 二． 区域剖分 考虑边值问题（1），（4）的差分逼近。去空间步长N l h /=和时间步长M T /=τ，其中N,M 都是正整数。用两族平行直线： 将矩形域}0;0{T t l x G ≤≤≤≤=分割成矩形网格，网格节点为),(k j t x 。以h G 表示网格内点集合，即位于开矩形G 的网点集合；h G 表示所有位于闭矩形G 的网点集合；Γ=G --G 是网格界点集合。

三． 离散格式 第k+1层值通过第k 层值明显表示出来，无需求解线性代数方程组，这样的格式称为显格式。 第k+1层值不能通过第k 层值明显表示出来，而由线性代数方程组确定，这样的格式称为隐格式。 1. 向前差分格式 （5） ,22111j k j k j k j k j k j f h u u u a u u ++-=--++τ )(j j x f f =， )(0 j j j x u ??==， 00==k N k u u ， 其中j = 1,2,…，N-1，k = 1,2,…,M-1。以2/h a r τ=表示网比。则方程（5）可以改写为： 易知向前差分格式是显格式。 2. 向后差分格式 （6） ,11111)21(j k j k j k j k j f u ru u u ru τ+=-++-+-+++ )(0 j j j x u ??==， 00==k N k u u ， 其中j = 1,2,…，N-1，k = 1,2,…,M-1，易知向前差分格式是显格式。 3. 六点对称格式（Grank-Nicolson 格式） 将向前差分格式和向后差分格式作算术平均，即得到六点对称格式： （7） 111112)1(2+-+++-++-k j k j k j u r u r u r =j k j k j k j f u r u r u r τ++-+-+112 )1(2 利用0j u 和边值便可逐层求到k j u 。六点对称格式是隐格式，由第k 层计算第k+1层时需解线性代数方程组（因系数矩阵严格对角占优，方程组可唯一求解）。

本科毕业设计--求解热传导方程的高精度隐式差分格式

新疆大学毕业论文(设计) 题目:求解热传导方程的高精度隐式差分格式所属院系：数学与系统科学学院 专业：信息与计算科学

声明 本人郑重声明该毕业论文（设计）是本人在开依沙尔老师指导下独立完成的，本人拥有自主知识产权，没有抄袭、剽窃他人成果，由此造成的知识产权纠纷由本人负责。 声明人（签名）： 年月日 亚库甫江.买买提同学在指导老师的指导下，按照任务书的内容，独立完成了该毕业论文（设计），指导教师已经详细审阅该毕业论文（设计）。 指导教师（签名）： 年月日

新疆大学 毕业论文（设计）任务书 班级：信计07-2 姓名：亚库甫江.买买提论文（设计）题目：求解热传导方程的高精度隐式差分格式 专题：毕业设计 论文（设计）来源：教师自拟 要求完成的内容：学习和掌握一维热传导方程已有的各种差分 格式的基础上，扩散方程对空间变量应用紧 致格式离散，对时间变量应用梯形方法，构 造热传导方程的精度为() 24 τ+数值格式， O h 讨论格式的稳定性，最后数值例子来验证。发题日期：2012 年12月25日完成日期：2012 年5月28 日实习实训单位：数学学院地点：数学学院 论文页数：19页；图纸张数：4 指导教师：开依沙尔老师 教研室主任 院长（系主任）

摘要 本文首先对热传导方程经典差分格式进行复习和讨论，然后热传导方程对空间变量四阶紧致格式进行离散,时间变量保持不变，把一维热传导方程转化为常微分方程组的初值问题, 再利用梯形方法构造热传导方程方程的时间二阶空间四阶精度的一种差分格式，并稳定性进行分析，数值结果与Crank-Nicholson 格式进行比较，数值结果表明, 该方法是有效求解热传导方程的数值计算. 关键词: 热传导方程,高精度紧致格式; 梯形方法；两层隐格式; Crank-Nicolson格式 ABSTRACT This paper first study on some classical finite difference for the heat conduction equation, secondely secondely we apply compact finite difference approximation of fourth order for discretizing spatial derivatives but leave the time variable Continuous. This approach results in a system of ODEs, which can then be used trapezodial formula derived fourth order in space and second order in time unconditionally stable implicit scheme .the stability and local truncation error of the obtained method are analysied. Numerical experiments shows that this method Useful, efficient method for solving diffusion equation Keywords: Heat conduction eqution;Higher- oder compact scheme; Trapezodial formula ;Two- level implict scheme; Crank- Nicolson scheme

有限差分和Matlabpde求解一维稳态传热问题.(优选)

有限差分和pde 函数求解一维定态热传导方程 分别用有限差分方法和pde 函数求解一维定态热传导方程，初始条件和边界条件，热扩散系数α=0.00001， 22 T T t x α??=?? （1） 求解过程： 1. 用Tylaor 展开法推导出FTCS 格式的差分方程 首先对T 进行泰勒展开得到如下两式子： 2 3 1231 2 3 ... 232! 3! 2 3 ... 232!3!n n n n n j j j j j n n n n n j j j j j t t T T t x x T T x T T T t t t T T T x x x ++??=+?+ + +??=+?+ + +????????? ? ? ?????? ???? ????????? ? ? ?????? ?? ?? 上述两个方程变换得： ()11223 23...23n n n n n n n j j j j j j j T T T T T t T t T o t t t t t t ++--???? ???????= --=+? ? ? ???????????? （2） 223 123...23n n n n n j j j j j T T T x T x T x x x x --???? ???????= -- ? ? ??????????? ()1232422 342222...3!4!n n n n n n j j j j j j T T T T x T x T x x x x x x +-?? ????????????=--- ? ? ? ??????????????? ()()2112 22 22-n n n j j j T T T T o x x x +--+???=+? ????? （3） 将上述式子（2）（3）代入（1）得：

向前差分格式求解二维热传导方程

用向前差分格式求解二维热传导方程function varargout=liu(varargin) T=1;a=1;h=1/30;dt=1/150; [X,T,Z]=chfenmethed(h,dt,a,T); mesh(X,T,Z(:,:,3)); shading flat; % xlabel('X','FontSize',14); % ylabel('t','FontSize',14); % zlabel('error','FontSize',14); % title('误差图'); function [X,Y,Z]=chfenmethed(h,dt,a,T); %求解下问题 %u_t-a*(u_xx+u_yy)=f(x,y,t) 0

n=length(t); r=a*dt/h^2; [X,Y]=meshgrid(x,y); Z=zeros(m,m,n); U=zeros(m,m,n); for i=1:m for j=1:m U(i,j,1)=d(x(i),y(j)); end end for j=2:n for k=1:m U(1,k,j)=g0(y(k),t(j)); U(m,k,j)=g1(y(k),t(j)); U(k,1,j)=h0(x(k),t(j)); U(k,m,j)=h1(x(k),t(j)); end end for k=2:n for i=2:m-1 for j=2:m-1

热传导方程地差分格式

一维抛物方程的初边值问题 分别用向前差分格式、向后差分格式、六点对称格式，求解下列问题： 22,01,u u a x t x ??=< 在0.05,0.10.2t =和时刻的数值解，并与解析解2 (,)sin()t u x t e x ππ-=进行比较。 1差分格式形式 设空间步长1/h N =, 时间步长0τ>，T M τ=，网比2/r h τ=. （1）向前差分格式 该问题是第二类初边值问题（混合问题），我们要求出所需次数的偏微商的函数 (,)u x t ，满足方程22,01,u u a x t x ??=<。 已知sin x π在相应区域光滑，并且在0,x l =与边值相容，使问题有唯一充分光滑的 解。 取空间步长1/h N =，和时间步长/T M τ=，其中,N M 都是正整数。用两族平行直 线 (0,1,,) j x x jh j N ===L 和 (0,1,,) k t t k k M τ===L 将矩形域 {01,0}G x t =≤≤≥分割成矩形网络，网络格节点为(,)j k x t 。以h G 表示网格内点集合， 即位于矩形G 的网点集合；h G 表示闭矩形G 的网格集合；h h G G -=Γh 是网格界点的集合。 向前差分格式，即 i k j k j k j k j k j f h u u u a u u ++-=--++2 1 112τ （1）

综合实验2 热传导方程的有限差分数值模拟

微分方程数值解实验报告 专业信息与计算科学班级信计101 姓名学号 协作队员实验日期2013 年 4 月18 日星期四成绩评定教师签名批改日期 题目 一、问题提出

二、 模型建立 三、 求解方法 使用古典显格式：)2(111n m n m n m n m n m U U U U U -+++-+=τ 其中22/h k a =γ

(k 和h 分别为时间与空间方向的步长，取k=0.005，h=0.1使得2/1/2≤h k ) 有12=a L=1，细杆各处的初始温度为)sin(x π，两端截面上的温度为0。 Matlab 程序如下： clc; k=0.005; h=0.1; r=k/h^2; t=0:k:0.1; n=length(t); x=0:h:1; Un=sin(pi*x); un=[]; for i=1:n u=[]; for p=1:11 u1=exp(-pi^2*t(i))*sin(pi*x(p)); u=[u u1]; end un1=[]; for j=2:10 Un1=r*Un(j-1)+(1-2*r)*Un(j)+r*Un(j+1); un1=[un1 Un1]; end e=abs(u-Un); un=[un;u;Un;e]; Un=[0 un1 0]; end un 四、 输出结果

五、结果分析 在同一个时间下，细杆内的温度分布为：细杆内中间的温度最高，往两边逐渐下降到0，并且温度值关于x=0.5这条直线对称。 在不同的时间下，细杆内的温度分布为：随着时间的增加，在同一点的温度逐渐减少。 模拟值与真实值之间的误差不超过1％。

一维热传导方程的前向 、紧差分格式

中南林业科技大学 本科课程论文 学院：理学院 专业年级：09信息与计算科学一班 课程：偏微分方程数值解法 论文题目：一维热传导方程的前向Euler和紧差分格式指导教师：陈红斌 2012年7月

学生姓名：唐黎学号： 20093936分工：程序编写，数值例子 学生姓名：何雄飞学号：20093925分工：格式建立，资料收集 学生姓名：汪霄学号：20093938分工：文档编辑，资料整理 学生姓名：毛博伟学号：20093931分工：公式编辑，查找资料 学生姓名：倪新东学号：20093932分工：数据分析，查找资料 学生姓名：何凯明学号：20093924分工：数据分析，查找资料

目录 1引言 (1) 2物理背景 (1) 3网格剖分 (2) 4．1.1向前Euler格式建立 (2) 4.1.2差分格式的求解 (4) 4.1.3收敛性与稳定性 (4) 4.1.4 数值例子 (7) 4.2.1紧差分格式建立 (10) 4.2.2差分格式求解 (12) 4.2.3数值例子 (13) 总结 (17) 参考文献 (18) 附录 (19)

1 引言 本文考虑的一维非齐次热传导方程的定解问题： 22(,),0,0,u u a f x t x l t T t x ??-=<<<≤?? (,0)(),0,u x x x l φ=≤≤ (0,)(), (1,)(), 0.u t t u t t t T αβ==<≤ 其中a 为正常数，(,),(),(),()f x t x t t ?αβ为已知函数，(0)(0),(1)(0).?α?β== 目前常用的求解热传导方程的差分格式有前向Euler 差分格式、向后Euler 差分格式、Crank-Nicolson 格式、Richardson 格式[1,2,3]．本文将给出前向Euler 格式和紧差分格式，并给出其截断误差和数值例子． 2 物理背景 热传导是由于物体内部温度分布不均匀,热量要从物体内温度较高的点流向温度较低的点处.以函数(),,,u x y z t 表示物体在t 时刻,(),M M x y =处的温度,并假设 (),,u x y z 关于,,x y z 具有二阶连续偏导数,关于t 具有一阶连续偏导数.() ,,k k x y z =是物体在(),,M x y z 处的热传导系数,取正值.设物体的比热容为(),,c c x y z =,密度为 (),,x y z ρ.根据Fourier 热传导定律,热量守恒定律以及Gauss 公式得 ,u u u u c kx k k t x x y y z z ρ ????????????? =++ ? ? ???????????? ?? 如果物体是均匀的,此时,k c 以及ρ均为常数.令2 k a c ρ = ,上式方程化为 2222 2222,t u u u u a a u x y z ?? ???=++=? ?????? 若考虑物体内有热源,其热源密度函数为(),,F F x y z =,则有热源的热传导方程为 ()2,,,,t u a u f x y z t =?+ 其中F f c ρ = .

最简差分格式

最简差分格式 朱勇军 20083738 考虑问题： 考虑一维热传导方程： ()x f x u a t u +??=??22， T ≤

差分格式： 向前差分格式： 向量形式： j k j k j k j k j k j f h u u u a u u ++-=--++2 1 11 2τ () ()0 ,00 =====k N k j j j j j u u x u x f f ?? 其中1,,2,1,1,2,1-=-=M k N j ，以2h a r τ=表示网比，则该格式可变为：()j k j k j k j k j f ru u r ru u τ++-+=-++111 21 矩阵形式： ?????????? ??????+++??????????????????????????????? ?----=????????????????----+-+-++k N N N k k N k N k k k N k N k k ru f f f ru f u u u u r r r r r r r r r u u u u 1220112211112121121002100210021ττττ 向后差分格式： 向量形式： j k j k j k j k j k j f h u u u a u u ++-=-+-++++2 1 1 1 1 1 1 2τ () ()0 ,00 =====k N k j j j j j u u x u x f f ?? 其中1,,2,1,1,2,1-=-=M k N j ，以2 h a r τ=表示网比，则该格式可 变为：()j k j k j k j k j f u ru u r ru τ+=-++-+-+++1 1 1 1 1 21 矩阵形式： ??????? ?????????+++????????????????=????????????????????????????????+--+-+--+----+-+-++k N N N k k N k N k k k N k N k k ru f f f ru f u u u u u u u u r r r r r r r r r 122011********* 1121002100210021ττττ

热传导方程差分格式

热传导方程的差分格式第2页 一维抛物方程的初边值问题 分别用向前差分格式、向后差分格式、六点对称格式，求解下列问题： .:u ；：2u a 2，0 ：：： x ：： 1, .:t ；x u(x,0) =sin 二x, 0 ：： x :: 1 u(O,t) =u(1,t) =0, t 0 2 在t =0.05,0.1和0.2时刻的数值解，并与解析解u(x,t)=ef t s in (二x)进行比较 1差分格式形式 2 设空间步长h =1/ N ,时间步长? ? 0, T =M ?，网比r = ? / h . (1) 向前差分格式 该问题是第二类初边值问题(混合问题)，我们要求出所需次数的偏微商的函数 Eu c2u u(x,t)，满足方程—a—^, 0 ：：：x ：：：1,和初始条件u(x,0)= sin x , 0 x ：：： 1抚ex 及边值条件u(0,t) =u(1,t) =0, t 0。 已知sin二x在相应区域光滑，并且在x =0,丨与边值相容，使问题有唯一充分光滑的解。 取空间步长h =1/ N，和时间步长-T /M，其中N,M都是正整数。用两族平行 直线x= j X =( j h0Hj 1 ，和tlNt k =X(k=0,1,|||,M) 将矩形域 G二｛0 — x — 1,t — 0｝分割成矩形网络，网络格节点为(X j,t k)。以G h表示网格内点集合， 即位于矩形G的网点集合；G h表示闭矩形G的网格集合；丨h=G h-G h是网格界点的集合。 向前差分格式，即 k 1 k k k k u, -u, u, 1 -2比■ u, 4 - j二a 亠2亠t ( 1) ￡h

热传导方程的差分格式

热传导方程的差分格式——上机实习报告 二零一四年五月

热传导方程的差分格式 第1页 一维抛物方程的初边值问题 分别用向前差分格式、向后差分格式、六点对称格式，求解下列问题： 22,01,u u a x t x ??=< 在0.05,0.10.2t =和时刻的数值解，并与解析解2(,)sin()t u x t e x ππ-=进行比较。 1差分格式形式 设空间步长1/h N =, 时间步长0τ >，T M τ=，网比2/r h τ=. （1）向前差分格式 向前差分格式，即 i k j k j k j k j k j f h u u u a u u ++-=--++21112τ （1） )(i i x f f =， 0),(00 ====k N k j j j u u x u ?? 其中，.1,,2,1,1,,2,1-=-=M k N j 以2 /h a r τ=表示网比。（1）式可改写成如下： j k j k j k j k j f ru u r ru u τ++-+=-++111)21( 此格式为显格式。 其矩阵表达式如下： ??????? ? ??=???????? ?????????? ??----++-++-111121112121212121j N j N j j j N j N j j u u u u u u u u r r r r r r r r r

热传导方程的差分格式 第2页 （2）向后差分格式 向后差分格式，即 ,2211 11k 11j k j k j j k j k j f h u u u a u u ++-=-+-++++τ （2） ,0),(N 00 ====k k j j j u u x u ?? 其中.1,,2,1,1,,2,1-=-=M k N j （2）式可改写成 j k j k j k j k j f u ru u r ru τ+=-++-+-+++11111)21( 此种差分格式被称为隐格式。 其矩阵表达式如下： ??????? ? ??=???????? ?????????? ??+--+-+--+-++-++j N j N j j j N j N j j u u u u u u u u r r r r r r r r r 121111121121212121 （3）六点对称格式 六点差分格式： j k j k j k j k j k j k j k j k j f h u u u h u u u a u u ++-++-=--++-++++]22[22112111111τ （3） .0),(00====k N k j j j u u x u ?? 将（3）式改写成 j k j k j k j k j k j k j f u r u r u r u r u r u r τ++-+=-++--++-+++11111112 )1(22)1(2 其矩阵表达式如下: ???????? ?????????? ??----=???????? ?????????? ??+--+-+--+-++-++j N j N j j j N j N j j u u u u r r r r r r r r r u u u u r r r r r r r r r 1211111211212/2/12/12/2/1212/12/12/2/1 2利用MATLAB 求解问题的过程 对每种差分格式依次取40.N =，=1/1600τ，=1/3200τ，=1/6400τ，用MATLAB

一维热传导方程的前向 、紧差分格式

页眉内容 中南林业科技大学 本科课程论文学院：理学院 专业年级：09信息与计算科学一班 课程：偏微分方程数值解法 论文题目：一维热传导方程的前向Euler和紧差分格式指导教师：陈红斌 2012年7月 学生姓名：唐黎学号： 分工：程序编写，数值例子 学生姓名：何雄飞学号： 分工：格式建立，资料收集 学生姓名：汪霄学号： 分工：文档编辑，资料整理 学生姓名：毛博伟学号： 分工：公式编辑，查找资料 学生姓名：倪新东学号： 分工：数据分析，查找资料 学生姓名：何凯明学号：

页眉内容 分工：数据分析，查找资料 目录 1引言 (1) 2物理背景 (1) 3网格剖分 (2) 4．1.1向前Euler格式建立 (2) (4) 4.1.4 数值例子 (7) (10) (12) (13) 总结 (17) 参考文献 (18) 附录 (19)

页眉内容 1 引言 本文考虑的一维非齐次热传导方程的定解问题： 其中a 为正常数，(,),(),(),()f x t x t t ?αβ为已知函数，(0)(0),(1)(0).?α?β== 目前常用的求解热传导方程的差分格式有前向Euler 差分格式、向后Euler 差分格式、Crank-Nicolson 格式、Richardson 格式[1,2,3]．本文将给出前向Euler 格式和紧差分格式，并给出其截断误差和数值例子． 2 物理背景 热传导是由于物体内部温度分布不均匀,热量要从物体内温度较高的点流向温度较低的点处.以函数(),,,u x y z t 表示物体在t 时刻,(),M M x y =处的温度,并假设(),,u x y z 关于,,x y z 具有二阶连续偏导数,关于t 具有一阶连续偏导数.(),,k k x y z =是物体在(),,M x y z 处的热传导系数,取正值.设物体的比热容为(),,c c x y z =,密度为(),,x y z ρ.根据Fourier 热传导定律,热量守恒定律以及Gauss 公式得 如果物体是均匀的,此时,k c 以及ρ均为常数.令2k a c ρ =,上式方程化为 若考虑物体内有热源,其热源密度函数为(),,F F x y z =,则有热源的热传导方程为 其中F f c ρ =. 3 网格剖分 取空间步长N l h /=和时间步长M T /=τ，其中M N ,都是正整数．用两族平行直线),1,0(N j jh x j Λ==和),,1,0(M k k t k Λ==τ将矩形域}0,0{T t l x G ≤≤≤≤=分割成矩形网格，网格节点为),(k j t x ．记),(k j k j t x u u =．以h G 表示网格内点集合，即 位于开矩形G 的网点集合；h G 表示所有位于闭矩形的网点集合；h h h G G -=Γ是网格界点集合. 引进如下记号: