一维热传导方程的差分格式
《微分方程数值解》
课程论文
学生姓名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 初边值问题
∂u ∂2u
=a 2+f (x , t ), c
u (x ,0) =ϕ(x ), c ≤x ≤d , (1.2)
u (c , t ) =α(t ), u (d , t ) =β(t ), 0
的有限差分方法, 其中a 为正常数, f (x , t ), ϕ(x ), α(t ),
β(t ) 为已知常数, ϕ(c ) =α(0),
ϕ(d ) =β(0). 称(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 ,
x j =c +jh , 0≤j ≤M , t k =k τ, 0≤k ≤N . 分别称h 和τ为空间步长和时间步长. 用
两组平行直线
x =x j , 0≤j ≤M , t =t k , 0≤k ≤N
将Ω分割成矩形网格. 记Ωh =x j |0≤j ≤M , Ωτ={t k |0≤k ≤N }, Ωh τ=Ωh ⨯Ωτ. 称x j , t k 为结点.
k k
定义Ωh τ上的网格函数 Ω=U j |0≤j ≤M ,0≤k ≤N , 其中U j =u x j , t k .
{}
()
[1]
{}()
在结点x j , t k 处考虑方程(1.1),有
()
∂u (x j , t k )=a
∂2u (x j , t k )+f (x j , t k ), 1≤j ≤M -1, 1≤k ≤N -1. (2.1)
∂t
∂x
2
将u (x j , t k +1)以结点(x j , t k )
为中心关于t 运用泰勒级数展开, 有
u (x u ''(x j , t k )τ2
j , t k +1)=u (x j , t k )+u '(x j , t k )τ+
2!
+o (τ2).
整理有
u (x 2j , t k +1)-u (x j , t k )
=∂u (x j , t k )τ∂u (x j , t k )τ
∂t +2∂t
2
+o (τ). (2.2) 再将(x j -1, t k ), (x j +1, t k )分别以结点(x j , t k )
为中心关于x 运用泰勒级数展开, 有
2
u (x j , t k （)-h ）
j -1, t k )=u (x j , t k )+u '(x j , t k )(-h )
j , t k )(-h ）+
u ''(x 2!
+
u '''(x 3
3!
4
+
u (4)(x j , t k )(-h )
4!
+o (h 4) , (x u (x u ''(x j , t k )h 2
u '''(x j , t k )h 3
u j +1, t k )=j , t k )+u '(x j , t k )h +
2! +
3!
+
u (4)(x j , t k )h 4
4!
+o (h 4).
由上述两式可得
u (x j -1, t k )-2u (x j , t k )+u (x j +1, t k )∂2u (x j , t k h 2=)∂x 2+h 2∂4(x j , t k )12∂x
4
+o (h 2
) . (2.3) 将(2.2), (2.3)两式代入(2.1)中, 得
u (x j , t k +1)-u (x j , t k )
k )+u (x j +1, t k )
τ
=a
u (x j -1, t k )-2u (x j , t h
2
+f (x j , t k )+R k j . (2.4)
2其中R k ah ∂4u (x j , t k
)τ2∂2(x j , t k )j
=-12∂x 4+2∂t 2
+o (τ+h 2) 为方程(2.1)的截断误差. 舍去截断误差, 用u k
j 代替u (x j , t k )
, 得到如下差分方程
u k +1-u k j j
=a
u k k k
j -1-2u j +u j +1
τ
h 2
+f k j , 1≤j ≤M -1, 0≤k ≤N -1. 结合初边值条件, 可得如下差分格式
(2.5)
+1u k -u k j j
τ
=a
k k
u k j -1-2u j +u j +1
h 2
+f j k , 1≤j ≤M -1, 1≤k ≤N -1, (2.6)
u 0j =ϕ(x j ), 0≤j ≤M , (2.7)
k k u 0=α(t j ), u M =β(t k ), 1≤k ≤N . (2.8)
2.2 向后Euler 差分格式
在结点x j , t k +1处考虑方程(1.1), 有
()
∂u (x j , t k +1)∂t
=a
∂2u (x j , t k +1)∂x 2
+f (x j , t k +1), 1≤j ≤M -1, 1≤k ≤N . (2.9)
将u x j , t k 以x j , t k +1为中心关于t 运用泰勒级数展开, 有
(
)()
u (x j , t k )=u (x j , t k +1)+u '(x j , t k +1)(-τ) +
将上式整理得
u ''(x j , t k +1)(-τ) 2
2!
+o (τ2).
u (x j , t k +1)-u (x j , t k )
τ
∂u (x j , t k +1)τ∂2u (x j , t k +1)
=-+o (τ). (2.10) 2
∂t 2∂t
再将u x j -1, t k +1, u x j +1, t k +1分别以x j , t k +1为中心关于x 运用泰勒级数展开, 有
()()()
u (x j -1, t k +1)=u (x j , t k +1)+u '(x j , t k +1)(-h ）+
u ''(x j , t k +1)(-h )
2!
4
2
+
3
u '''(x j , t k +1（)-h ）
3!
+
u (4)(x j , t k +1)(-h )
4!
+o (h 4), u ''(x j , t k +1)h 2
2! +o (h 4).
u '''(x j , t k +1)h 3
3!
u (x j +1, t k +1)=u (x j , t k +1)+u '(x j , t k +1)h +
u (4)(x j , t k +1)h 4
4!
+
+
由上述两式可得
u (x j -1, t k +1)-2u (x j , t k +1)+u (x j +1, t k +1)∂2u (x j , t k +1)h 2∂4(x j , t k +1)
=++o (h 2) . (2.11) 224
h ∂x 12∂x
将(2.10), (2.11)两式代入(2.9)中, 得
τ
=a
h 2
+f (x j , t k +1)+R k j . (2.12)
其中R =-
k
j
24
τ∂(x i , t k +1)ah 2∂u (x i , t k +1)
2
∂t 2
-
12
∂x 4
+o (τ+h 2) 为方程(2.9)的截断误差.
k
舍去截断误差, 用u j 代替u x j , t k , 得到如下差分方程
+1u k -u k j j
+1k +1+1
u k +u k j -1-2u j j +1
()
τ
=a
h
2
+f j k +1, 1≤j ≤M -1, 1≤k ≤N . (2.13)
结合初边值条件, 可得如下差分格式
+1
u k -u k j j
τ
=a
+1k +1+1
u k +u k j -1-2u j j +1
h
2
+f j k +1, 1≤j ≤M -1, 1≤k ≤N , (2.14)
u 0j =ϕ(x j ), 0≤j ≤M , (2.15)
k k u 0=α(t k ), u M =β(t k ), 1≤k ≤N . (2.16)
2.3 Crank -Nicolson 差分格式 在结点x j , t k +1/2处考虑方程(1.1), 有
()
∂u (x j , t k +1/2)∂t
=a
∂2u (x j , t k +1/2)∂x
2
+f (x j , t k +1/2), 1≤j ≤M -1, 0≤k ≤N -1. (2.17)
将u x j , t k +1, u x j , t k 以x j , t k +1/2为中心关于t 运用泰勒级数展开, 有
()(
)()
u ''(x j , t k +1/2)⎛τ⎫2
u (x j , t k +1)=u (x j , t k +1/2)+u '(x j , t k +1/2)+ ⎪
22! ⎝2⎭
τ
u '''(x j , t k +1/2)⎛τ⎫3
3
+ ⎪+o (τ).
3! ⎝2⎭
2
τ⎫u ''(x j , t k +1/2)⎛τ⎫⎛
u (x j , t k )=u (x j , t k +1/2)+u '(x j , t k +1/2) -⎪+ -⎪
22! ⎝⎭⎝2⎭
u '''(x j , t k +1/2)⎛τ⎫3
3
+ -⎪+o (τ).
3! ⎝2⎭
将上述两式整理得
τ
=++o (τ3). (2.18) 3
∂t 24∂t
再将u x j -1, t k , u x j +1, t k 分别以x j , t k 为中心关于x 运用泰勒级数展开, 有
()()()
u (x j -1, t k )=u (x j , t k )+u '(x j , t k )(-h ）+
u ''(x j , t k )(-h )
2!
4
2
+
3
u '''(x j , t k （) -h ）
3!
+
u (4)(x j , t k )(-h )
4! u ''(x j , t k )h 2
2!
+o (h 4) , u '''(x j , t k )h 3
3!
u (4)(x j , t k )h 4
4!
u (x j +1, t k )=u (x j , t k )+u '(x j , t k )h +
由上述两式得
++
+o (h 4).
u (x j -1, t k )-2u (x j , t k )+u (x j +1, t k )∂2u (x j , t k )h 2∂4(x j , t k )2
=++o (h ). (2.19) 224
h ∂x 12∂x
同理, 将u x j -1, t k +1, u x j +1, t k +1分别以x j , t k +1为中心关于x 泰勒级数展开, 整理得
()()()
u (x j -1, t k +1)-2u (x j , t k +1)+u (x j +1, t k +1)∂2u (x j , t k +1)h 2∂4(x j , t k +1)
=++o (h 2). (2.20) 224
h ∂x 12∂x
此时，分别将u x j , t k +1, u x j , t k 以u x j , t k +1/2为中心关于t 泰勒级数展开，有
()()
()
∂2u (x j , t k +1)∂2u (x j , t k +1/2)∂3u (x j , t k +1/2)τ1∂4u (x j , t k +1/2)⎛τ⎫2
2
=+++o τ, () ⎪22222
∂x ∂x ∂x ∂t 22! ∂x ∂t ⎝2⎭
∂2u (x j , t k )∂2u (x j , t k +1/2)∂3u (x j , t k +1/2)⎛τ⎫
=+ -⎪ 222
∂x ∂x ∂x ∂t ⎝2⎭
42
1∂u (x j , t k +1/2)⎛τ⎫2+ -⎪+o (τ). 222! ∂x ∂t ⎝2⎭
利用上述两式得
∂2u (x j , t k +1/2)∂x 2
2242
1⎛∂u (x j , t k +1)∂u (x j , t k )⎫1⎛∂u (x j , t k +1/2)⎛τ⎫⎫2
⎪ ⎪= +-+o τ. () ⎪2222
⎪ ⎪2∂x ∂x 2∂x ∂t ⎝2⎭⎭⎝⎭⎝
(2.21)
利用(2.19), (2.20)两式, 整理有
∂2u (x j , t k +1)∂x
2
+
∂2u (x j , t k )∂x
2
=
u (x j -1, t k )-2u (x j , t k )+u (x j +1, t k )
h
2
+
u (x j -1, t k +1)-2u (x j , t k +1)+u (x j +1, t k +1)
h 24
h 2∂(x j , t k ) -
12∂x 4
4
h 2∂(x j , t k +1)2-+o h . (2.22) ()412∂x
结合(2.21),(2.22)两式, 整理得
∂2u (x j , t k +1/2)∂x
+
2
=
u (x j -1, t k )-2u (x j , t k )+u (x j +1, t k )
2h
2
u (x j -1, t k +1)-2u (x j , t k +1)+u (x j +1, t k +1)
2h 2
42
1⎛∂u (x j , t k +1/2)⎛τ⎫⎫
- ⎪⎪⎪2 ∂x 2∂t 22⎝⎭⎭⎝
44
⎫1⎛h 2∂(x j , t k )h 2∂(x j , t k +1)2
- ++o (h )⎪. (2.23) 44
⎪212∂x 12∂x ⎝⎭
将(2.18), (2.23)式代入(2.17)得
u (x j , t k +1)-u (x j , t k )
τ
+a
其中
=a
u (x j -1, t k )-2u (x j , t k )+u (x j +1, t k )
2h 2
u (x j -1, t k +1)-2u (x j , t k +1)+u (x j +1, t k +1)
2h 2
+f (x j , t k )+R k j . (2.24)
⎛h 2∂4(x j , t k )h 2∂4(x j , t k +1)∂4u (x j , t k +1/2)⎛τ⎫2⎫τ2∂3u (x j , t k +1/2)a
⎪+R k ++ j =- ⎪44223
⎪2 12∂x 12∂x ∂x ∂t 224∂t ⎝⎭⎭⎝+o (τ3+h 2) 为方程(2.17)的截断误差.
舍去截断误差, 用u j 代替u (x j , t k ) , 可得如下差分方程
+1
u k -u k j j
k
τ
=a
k k
u k j -1-2u j +u j +1
2h 2
+a
+1k +1+1
u k +u k j -1-2u j j +1
2h 2
+f j k +1/2,
1≤j ≤M -1, 0≤k ≤N -1. (2.24)
结合初边值条件, 可得如下差分格式
+1
u k -u k j j
τ
=a
k k
u k j -1-2u j +u j +1
2h 2
+a
+1k +1+1
u k +u k j -1-2u j j +1
2h 2
+f j k +1/2,
1≤j ≤M -1, 0≤k ≤N -1, (2.25)
u 0j =ϕ(x j ), 0≤j ≤M , (2.26)
k k u 0=α(t k ), u M =β(t k ), 1≤k ≤N . (2.27)
3 差分格式的求解
3.1 向前Euler 格式 记r =a τ
h 2
, 称r 为步长比. 差分格式(2.6) (2.8)中(2.6)可改写为
u k +1j =r ⋅u k ) u k k k
j -1+(1-2r j +r ⋅u j +1+τ⋅f j ,
1≤j ≤M -1 , 0≤k ≤N -1 .
将(3.1)写成如下形式
U k +1=A ⋅U k +τ⋅F k +r ⋅F 1 .
其中
U k +1=⎡u k +1, u k +1, , u k +1⎤T
U k =⎡u k , u k k
⎤T
⎣12M -1⎦, ⎣
12, , u M -1⎦, F k =⎡f k , f k , , f k
⎤T F =⎡u k , 0, k ⎤T
⎣12M -1⎦, 1⎣0
0, , 0, u M ⎦(M -1)*1, ⎛ 1-2r r ⎫
r 1-2r r ⎪A = ⎪
r
⎪. r ⎪
⎪
⎝r 1-2r ⎪⎭(M -1)*(M -1)
3.2 向后Euler 格式 记r =a τ
h 2
, 称r 为步长比. 差分格式(2.14) (2.16)中(2.14)可改写为
-u k j =r ⋅u k +1k +1k +1k +1
j -1-(1+2r ) u j +r ⋅u j +1+τ⋅f j
,
1≤j ≤M -1 , 1≤k ≤N .
将(3.2)写成如下形式
A ⋅U k +1=-U k -τ⋅F k +1-r ⋅F 1 .
其中
(3.1)
(3.2)
k +1k +1k +1k k k k
⎤⎡⎤U k +1=⎡u , u , , u U =u , u , , u , 12M -112M -1⎦, ⎣⎦⎣
T T
k +1k +1k +1k +1k +1
⎡⎤⎤F k +1=⎡f , f , , f F =u , 0, 0, , 0, u , 12M -110M ⎣⎦⎣⎦
T T (M -1)*1
,
r ⎛-(1+2r ) ⎫
⎪
r -(1+2r ) r ⎪
⎪A = . r
⎪
r ⎪
r -(1+2r ) ⎪⎝⎭(M -1)*(M -1)
3.3 Crank -Nicolson 格式 记r =a τ
, 称r 为步长比. 差分格式(2.25) (2.27)中(2.25)可改写为
h 2
+1k +1k +1k k k k +1/2
, -[r /2⋅u k -(1+r ) u +r /2⋅u ]=r /2⋅u +(1-r ) u +r /2⋅u +τ⋅f j -1j j +1j -1j j +1j
将(3.3)写成如下形式
1≤j ≤M -1 , 0≤k ≤N -1 .
(3.3)
. -A 2⋅U k +1=A 1⋅U k +τ⋅F k +1/2+r (F 1+F 2）
其中
k +1k +1k +1k k k k
⎤⎡⎤U k +1=⎡u , u , , u U =u , u , , u , 12M -112M -1⎦, ⎣⎦⎣
T
T
k k k +1k +1
F 1=⎡⎣u 0, 0, 0, , 0, u M ⎤⎦(M -1)*1, F 2=⎡⎣u 0, 0, 0, , 0, u M ⎤⎦(M -1)*1,
T T
F
k +1/2
k +1/2k +1/2
⎤=⎡, f 2k +1/2, , f M -1⎦, ⎣f 1
T
⎛1-r r /2⎫
⎪r /21-r r /2 ⎪
⎪A 1= , r /2
⎪
r /2 ⎪
r /21-r ⎪⎝⎭(M -1)*(M -1)
r /2⎛-(1+r ) ⎫ ⎪r /2-(1+r ) r /2 ⎪
⎪A 2= . r /2
⎪
r /2 ⎪
r /2-(1+r ) ⎪⎝⎭(M -1)*(M -1)
4 数值例子
在本文中, 将应用向前Euler 格式(2.6) (2.8), 向后Euler 格式(2.14) (2.16)和
Crank -Nicolson 格式(2.25) (2.27)分别来求解如下算例.
算例 现有如下定解问题
∂u ∂2u x +2t
=2+e , 0
u (x ,0) =e x , 0≤x ≤1,
u (0,t ) =e 2t , u (1,t ) =e 1+2t , 0
上述定解问题的精确解为u (x , t ) =e 4.1 向前Euler 格式
在表4.1.1中给出了取步长h =1/10和τ=1/200(步长比r =1/2) 时计算得到的部分数值结果. 表4.1.2给出了取步长h =1/10和τ=1/100(步长比r =1) 时计算得到的部分数值结果, 随着计算层数的增加, 误差越来越大, 数值结果是发散的. 经步长比r 的多次取值发现, 当r ≤1/2时, 其数值结果是收敛的, 当r >1/2时, 其数值结果是发散的. 表4.1.3给出了r =1/2时, 取不同步长, 数值解的最大误差
x +2t
.
E ∞(h , τ) =max |u (x j , t k ) -u j k |.
1≤j ≤M -11≤k ≤N
从表4.1.3可以看出当空间步长缩小到原来的1/2, 时间步长缩小到原来的1/4时, 最大误差约缩小到原来的1/4.
图4.1.1给出了t =1时的数值解曲线(h =1/10,τ=1/200) 与精确解曲线. 图4.1.2给出了t =1时不同步长数值解的误差曲线图(r =1/2) , 由图4.1.2可知, 当步长比r 保持不变, 随着时间与空间方向上的加密, 误差越来越小. 图4.1.3给出了不同步长的误差曲面图.
表4.1.1 算例在部分结点处数值解、精确解和绝对误差(h =1/10,τ=1/200)
k
20
(x , t )
(0.5,0.1)
数值解 1.802810
精确解 1.803988
|精确解-数值解|
1.178e-3
60 80 100 120 140 160 180 200
(0.5,0.3) (0.5,0.4) (0.5,0.5) (0.5,0.6) (0.5,0.7) (0.5,0.8) (0.5,0.9) (0.5,1.0)
2.688651 3.283856 4.010886 4.898897 5.983523 7.308290 8.926366 10.902687
2.691234 3.287081 4.014850 4.903749 5.989452 7.315534 8.935213 10.913494
2.583e-3 3.225e-3 3.965e-3 4.852e-3 5.929e-3 7.243e-3 8.847e-3 1.081e-3
表4.1.2 算例在部分结点处数值解、精确解和绝对误差(h =1/10,τ=1/100)
k
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
(x , t )
(0.5,0.01) (0.5,0.02) (0.5,0.03) (0.5,0.04) (0.5,0.05) (0.5,0.06) (0.5,0.07) (0.5,0.08) (0.5,0.09) (0.5,0.10) (0.5,0.11) (0.5,0.12) (0.5,0.13) (0.5,0.14) (0.5,0.15) (0.5,0.16) (0.5,0.17)
数值解 1.521674 1.552123 1.583184 1.614870 1.647386 1.679785 1.716057 1.741446 1.807089 1,744143 2.092546 1.164923 4.148393 -4.709415 21.998093 -57.422399 178.294623
精确解 1.521962 1.552707 1.584074 1.616074 1.648721 1.682028 1.716007 1.750673 1.786038 1.822119 1.858928 1.896481 1.934792 1.973878 2.013753 2.054433 2.095936
|精确解-数值解|
2.879e-4 5.845e-4 8.901e-4 1.205e-3 1.336e-3 2.242e-3 5.051e-5 9.227e-3 2.105e-2 7.798e-2 2.336e-1 7.316e-1 2.213e+0 6.684e+0 1.998e+1 5.948e+1 1.762e+2
表4.1.3 算例取不同步长时数值解的最大误差(r =1/2)
h
1/10 1/20 1/40 1/80
τ
1/200 1/800 1/3200 1/12800
E ∞(h , τ)
1.178e-2 2.973e-3 7.431e-4 1.858e-4
E ∞(2h ,4τ) /E ∞(h , τ)
* 3.964 4.000 4.000
图4.1.1 算例取t =1时的数值解曲线(h =1/10,τ=1/200) 与精确解曲线
图4.1.2 算例取t =1时不同步长数值解的误差曲线(r =1/2)
图4.1.3 算例取不同步长数值解的误差曲面(r =1/2)
4.2 向后Euler 格式
在表4.2.1和表4.2.2中分别给出了h =1/10,τ=1/200(步长比r =1/2) 和
h =1/10, τ=1/100(步长比r =1) 时计算得到的部分数值结果, 发现两表的数值结果都
是收敛的, 经步长比r 的多次取值发现, 无论r 取任何值, 其数值结果都是收敛的. 表4.2.3 给出了r =1/2时, 取不同步长, 数值解的最大误差
E ∞(h , τ) =max |u (x j , t k ) -u j k |.
1≤j ≤M -11≤k ≤N
从表4.2.3可以看出当空间步长缩小到原来的1/2, 时间步长缩小到原来的1/4时, 最大误差缩小到原来的1/4. 表4.2.4给出了r =1时的类似结论.
图4.2.1给出了t =1时的数值解曲线(h =1/10,τ=1/200) 与精确解曲线. 图4.2.2给出了t =1时, 取不同步长数值解的误差曲线图(r =1/2) , 由图4.2.2可知, 当步长比r 保持不变, 随着时间与空间方向上的加密, 误差越来越小. 图4.2.3和图4.2.4分别当给出了步长比r =1/2和r =1下不同步长的误差曲面图.
表4.2.1 算例在部分结点处数值解、精确解和绝对误差(h =1/10,τ=1/200)
k
20 40 60 80 100 120 140 160 180 200
(x , t )
(0.5,0.1) (0.5,0.2) (0.5,0.3) (0.5,0.4) (0.5,0.5) (0.5,0.6) (0.5,0.7) (0.5,0.8) (0.5,0.9) (0.5,1.0)
数值解 1.805343 2.205674 2.694258 3.290867 4.019509 4.909453 5.996425 7.324052 10.926203 10.926203
精确解 1.803988 2.203396 2.691234 3.287081 4.014850 4.903749 5.989452 7.315534 8.935213 10.913494
|精确解-数值解|
1.355e-3 2.278e-3 3.023e-3 3.785e-3 4.659e-3 5.704e-3 6.973e-3 8.518e-3 1.041e-2 1.271e-2
表4.2.2 算例在部分结点处数值解、精确解和绝对误差(h =1/10,τ=1/100)
k
10 20 30 40 50 60 70 80 90 100
(x , t )
(0.5,0.1) (0.5,0.2) (0.5,0.3) (0.5,0.4) (0.5,0.5) (0.5,0.6) (0.5,0.7) (0.5,0.8) (0.5,0.9) (0.5,1.0)
数值解 1.788500 2.185742 2.670171 3.261551 3.983745 4.865788 5.943098 7.258921 8.866068 10.829041
精确解 1.786038 2.181472 2.664456 3.254374 3.974902 4.854956 5.929856 7.242743 8.846306 10.804903
|精确解-数值解|
2.461e-3 4.269e-3 5.715e-3 7.177e-3 8.843e-3 1.083e-2 1.324e-2 1.618e-2 1.976e-2 2.414e-2
表4.2.3 算例取不同步长时数值解的最大误差(r =1/2)
h
1/10 1/20 1/40 1/80
τ
1/200 1/800 1/3200 1/12800
E ∞(h , τ)
1.386e-2 3.509e-3 8.779e-4 2.195e-4
E ∞(2h ,4τ) /E ∞(h , τ)
* 3.950 3.997 3.999
表4.2.4 算例取不同步长时数值解的最大误差(r =1)
h
1/10 1/20 1/40
τ
1/100 1/400 1/1600
E ∞(h , τ)
2.659e-2 6.744e-3 1.688e-3
E ∞(2h ,4τ) /E ∞(h , τ)
* 3.943 3.955
1/80
1/6400 4.222e-4 3.999
图4.2.1 算例取t =1时的数值解曲线(h =1/10,τ=1/200) 与精确解曲线
图4.2.2 算例取t =1 时不同步长数值解的误差曲线(r =1/2)
图 4.2.3 算例取不同步长数值解的误差曲面(r =
1/2)
图 4.2.4 算例取不同步长数值解的误差曲面(r =1)
4.3 Crank -Nicolson 格式
在表4.3.1中给出了取步长h =1/10,τ=1/10时计算得到的部分数值结果, 表4.3.2给出了取步长h =1/100,τ=1/100时计算得到的部分数值结果, 发现两表的数值结果都是收敛的, 经步长比r 的多次取值发现, 无论r 取任何值, 其数值结果都是收敛的. 表4.3.3给出了取不同步长时, 所得数值解的最大误差
E ∞(h , τ) =max |u (x j , t k ) -u j k |.
1≤j ≤M -11≤k ≤N
从表可以看出当空间步长缩小到原来的1/2, 时间步长缩小到原来的1/2时, 最大误差缩小
到原来的1/4.
图4.3.1给出了t =1时的数值解曲线(h =1/10,τ=1/10)与精确解曲线. 图4.3.2给出了t =1时, 取不同步长数值解的误差曲线图, 由图4.3.2可知, 当步长比r 保持不变, 随着时间与空间方向上的加密, 误差越来越小. 图4.3.3给出了取不同步长时所得数值解的误差曲面图. 图4.3.4给出了向前Euler , 向后Euler 和Crank -Nicolson 三种方法数值解的误差对比图(h =1/10,τ=1/200) , 由图4.3.4分析可得, 当时间方向t 与空间方向
x 都取相同步长时, Crank -Nicolson 方法误差最小, 精度最高, 向前Euler 方法次之,
向后Euler 方法所得误差最大, 精度最低.
表4.3.1 算例在部分结点处数值解、精确解和绝对误差(h =1/10,τ=1/10)
k
1 2 3 4 5 6 7 8 9 10
(x , t )
(0.5,0.1) (0.5,0.2) (0.5,0.3) (0.5,0.4) (0.5,0.5) (0.5,0.6) (0.5,0.7) (0.5,0.8) (0.5,0.9) (0.5,1.0)
数值解 1.823114 2.227196 2.720414 3.322776 4.058460 4.957021 6.054520 7.395008 9.032284 11.032056
精确解 1.822119 2.225541 2.718282 3.320117 4.055200 4.953032 6.049647 7.389056 9.025013 11.023176
|精确解-数值解|
9.949e-4 1.656e-3 2.132e-3 2.659e-3 3.260e-3 3.988e-3 4.873e-3 5.952e-3 7.271e-3 8.880e-3
表4.3.2 算例在部分结点处数值解、精确解和绝对误差(h =1/100,τ=1/100)
k
10 20 30 40 50 60 70 80 90 100
(x , t )
(0.5,0.1) (0.5,0.2) (0.5,0.3) (0.5,0.4) (0.5,0.5) (0.5,0.6) (0.5,0.7) (0.5,0.8) (0.5,0.9) (0.5,1.0)
数值解 1.993726 2.435147 2.974297 3.632815 4.437131 5.919524 6.619421 8.084979 9.875016 12.061372
精确解 1.993715 2.435130 2.974297 3.632786 4.437096 5.419481 6.619369 8.084915 9.874938 12.061276
|精确解-数值解|
1.081e-5 1.749e-5 2.296e-5 2.864e-5 3.521e-5 4.308e-5 5.266e-5 6.433e-5 7.857e-5 9.597e-5
表4.3.3 算例取不同步长时数值解得最大误差
h
1/10 1/20 1/40 1/80 1/160 1/320 1/640
τ
1/10 1/20 1/40 1/80 1/160 1/320 1/640
E ∞(h , τ)
9.588e-3 2.429e-3 6.078e-4 1.520e-4 3.799e-5 9.499e-6 2.375e-6
E ∞(2h ,2τ) /E ∞(h , τ)
* 3.9457 3.9961 3.9990 3.9998 3.9999 4.0000
图4.3.1 算例取t =1时数值解曲线(h =1/10,τ=1/10)与精确解曲线
图4.3.2 算例取t 1时不同步长数值解的误差曲线
图4.3.3 算例取不同步长数值解的误差曲面
图4.3.4 算例取t =1时三种方法数值解的误差对比曲线(h =1/10,τ=1/200)
5 总结
在本文中, 给出了向前Euler , 向后Euler 和Crank -Nicolson 三种差分格式, 分别对抛物型方程进行求解, 其中向前Euler 格式是条件稳定的, 后两者是无条件稳定的. 经对比发现, Crank -Nicolson 格式的精度最高, 误差最小, 向前Euler 格式次之, 向后Euler 格式精度最低, 误差最大.
◆ 参考文献
[1] 孙志忠. 偏微分方程数值解法[M]. 第二版. 北京: 科学出版社, 2005.
[2] 李荣华. 偏微分方程数值解法[M]. 第二版. 北京: 高等教育出版社, 2005.
[3] 刘卫国. MATLAB程序设计教程[M]. 第二版. 北京: 中国水利水电出版社, 2010. ◆ 程序流程图
1 向前Euler 法
2 向后Euler 法
3 Crank Nicolson 法
程序代码
附录A
(向前Euler 法）
建立函数文件ForwardEuler.m ，如下：
% 求解一维非齐次热传导方程 ut-uxx=e^(x+2t), 0
% u(x,0)=e^x, 0
% u(0,t)=e^(2t), u(1,t)=e^(1+2t), 0
% 该问题的定解为 u(x,t)=e^(x+2t).
% 其中取a=1, h1表示时间t 方向上的步长, h2表示空间x 方向上的步长.
function [uxy,u,eps]=ForwardEuler(h1,h2,M,N,T)
% 求解一维非齐次热传导方程 ut-uxx=f(x,t),(a>0)
% 其中M,N 分别为x 与t 方向上的网格数。
% uxy为精确解, uij为数值解, eps为数值解与精确解的绝对误差.
u=zeros(N+1,M+1); % 初始化矩阵u, 存储数值解的值.
x=0:h2:1; % h2=(1-0)/M
t=0:h1:T; % h1=(T-0)/N
% 给边值条件赋值 u(0,t)=e^(2t), u(1,t)=e^(1+2t).
for k=1:N+1
u(k,1)=exp(2.*t(k));
u(k,M+1)=exp(1+2.*t(k));
end
% 给初值条件赋值 u(x,0)=e^x
for j=1:M+1
u(1,j)=exp(x(j));
end
f=zeros(N+1,M+1); % 初始化源函数f(x,t).
% 给源函数(Source Function)赋值, f(x,t)=e^(x+2t).
for k=1:N+1
for j=1:M+1
f(k,j)=exp(x(j)+2.*t(k));
end
end
% r=(a*h1)/(h2^2)为步长比, 这里取a=1,即 r=h1/(h2^2);
% 当 r≤1/2 时, 数值解是越来越接近精确解的,
% 当 r > 1/2 时, 两者的误差会越来越大.
r=h1/(h2^2);
for k=1:N
for j=2:M
u(k+1,j)=r.*( u(k,j-1)+u(k,j+1) ) + (1-2*r).*u(k,j) + h1.*f(k,j); end
end
uxy=zeros(N+1,M+1); %给精确解赋初值.
eps=zeros(N+1,M+1); %给数值解与精确解的绝对误差赋初值.
for k=1:N+1
for j=1:M+1
uxy(k,j)=exp(x(j)+2.*t(k)); %给精确解赋值.
eps(k,j)=abs( u(k,j)-uxy(k,j) ); %数值解与精确解的绝对误差. end
end
附录B
（向后Euler 法）
建立函数文件BackwordEuler.m ，如下：
% 求解一维非齐次热传导方程 ut-uxx=e^(x+2t), 0
% u(x,0)=e^x, 0
% u(0,t)=e^(2t), u(1,t)=e^(1+2t), 0
% 该问题的定解为 u(x,t)=e^(x+2t).
% 其中取a=1, h表示时间t 方向上的步长, c表示空间x 方向上的步长.
function [U,ux,eps]=BackwordEuler(M,N,h,c,x,t)
% 求解一维非齐次热传导方程 ut-uxx=f(x,t),(a>0)
% 其中M,N 分别为x 与t 方向上的网格数。
% ux为精确解, U为数值解, eps为数值解与精确解的绝对误差.
U=zeros((M+1)*(N+1),1);
r=c./h^2;
aa1=zeros(1,M-2); %给矩阵A 对角线序号为-1的位置赋值 for i=1:M-2
aa1(i)=r;
end
aa2=zeros(1,M-1); %给矩阵A 主对角线上元素赋值. for i=1:M-1
aa2(i)=-1-2*r;
end
aa3=zeros(1,M-2); %给矩阵A 对角线序号为1的位置赋值 for i=1:M-2
aa3(i)=r;
end
a=zeros(M-1,M-1);
A=diag(aa1,-1)+diag(aa2)+diag(aa3,1); %构建矩阵A.
k=1;%记录行数
count=1;%记录总行数
for i=1:M+1 %给u （x,0）赋值
U(count,1)=exp(x(i));
count=count+1;
end
U0=zeros(M-1,1);
for i=1:M-1 %U向量初始值
U0(i,1)=exp(x(i+1));
end
f1=zeros(M-1,1);
f=zeros(M-1,1);
while k~=N+1
f1(1,1)=exp(2*t(k+1));
f1(M-1,1)=exp(1+2*t(k+1));
F1=r*f1;
for j=1:M-1
f(j,1)=exp(x(j+1)+2*t(k+1));
end
F=c*f;
U1=inv(A)*(-U0-F-F1);
U(count,1)=exp(2*t(k+1)); %给最左边界赋值
count=count+1;
for i=1:M-1
U(count,1)=U1(i,1);
count=count+1;
end
U(count,1)=exp(1+2*t(k+1)); %给最右边界赋值
count=count+1;
k=k+1; %进去下一行
U0=U1;
end
U=reshape(U',M+1,N+1);
p=1; %求精确解
for i=1:N+1
for j=1:M+1
ux(p,1)=exp(x(j)+2*t(i));
p=p+1;
end
end
ux=reshape(ux',M+1,N+1);
eps=abs(ux-U); %精确解与数值解的绝对误差
附录C
(Crank Nicolson 法）
建立函数文件GrankNicolson.m ，如下：
% 求解一维非齐次热传导方程 ut-uxx=e^(x+2t), 0
% u(x,0)=e^x, 0
% u(0,t)=e^(2t), u(1,t)=e^(1+2t), 0
% 该问题的定解为 u(x,t)=e^(x+2t).
% 其中取a=1, h表示时间t 方向上的步长, c表示空间x 方向上的步长.
function [U,ux,eps]=GrankNicolson(M,N,h,c,x,t)
% 求解一维非齐次热传导方程 ut-uxx=f(x,t),(a>0)
% 其中M,N 分别为x 与t 方向上的网格数。
% ux为精确解, U为数值解, eps为数值解与精确解的绝对误差.
U=zeros((M-1)*(N-1),1);
r=c./h^2;
aa1=zeros(1,M-2); %给矩阵A2对角线序号为-1的位置赋值 for i=1:M-2
aa1(i)=-r./2;
end
aa2=zeros(1,M-1); %给矩阵A2主对角线上元素赋值. for i=1:M-1
aa2(i)=1+r;
end
aa3=zeros(1,M-2); %给矩阵A2对角线序号为1的位置赋值 for i=1:M-2
aa3(i)=-r./2;
end
A2=diag(aa1,-1)+diag(aa2)+diag(aa3,1); %构造矩阵A2.
a1=zeros(1,M-2); %给矩阵A1对角线序号为-1的位置赋值 for i=1:M-2
a1(i)=r./2;
end
a2=zeros(1,M-1); %给矩阵A1主对角线上元素赋值. for i=1:M-1
a2(i)=1-r;
end
a3=zeros(1,M-2); %给矩阵A1对角线序号为1的位置赋值 for i=1:M-2
a3(i)=r./2;
end
A1=diag(a1,-1)+diag(a2)+diag(a3,1); %构造矩阵A1.
k=1;%记录行数
count=1;%记录总行数
for i=1:M+1 %给u （x,0）赋值
U(count,1)=exp(x(i));
count=count+1;
end
U0=zeros(M-1,1);
for i=1:M-1 %U向量初始值
U0(i,1)=exp(x(i+1));
end
f1=zeros(M-1,1); %F(k=0)
for j=1:M-1
f1(j,1)=exp(x(j+1)+t(k));
end
f2=zeros(M-1,1); %F(k=1)
for j=1:M-1
f2(j,1)=exp(x(j+1)+2*t(k+1)); end
F1=zeros(M-1,1); %F1(k=0)
F1(1,1)=1;
F1(M-1,1)=exp(1);
F2=zeros(M-1,1); %F2(k=1)
while k~=N+1
F=(f1+f2)/2;
F2(1,1)=exp(2*t(k+1));
F2(M-1,1)=exp(1+2*t(k+1));
U1=inv(A2)*(A1*U0+c*F+r./2*(F1+F2));
U(count,1)=exp(2*t(k+1)); %给最左边界赋值 count=count+1;
for i=1:M-1
U(count,1)=U1(i,1);
count=count+1;
end
U(count,1)=exp(1+2*t(k+1)); %给最右边界赋值 count=count+1;
if k~=N
F1(1,1)=exp(2*t(k+1));
F1(M-1,1)=exp(1+2*t(k+1));
for j=1:M-1
f1(j,1)=exp(x(j+1)+2*t(k+1));
end
for j=1:M-1
f2(j,1)=exp(x(j+1)+2*t(k+2));
end
end
k=k+1;
U0=U1;
end
U=reshape(U',M+1,N+1);
p=1; %求精确解
for i=1:N+1
for j=1:M+1
ux(p,1)=exp(x(j)+2*t(i)); p=p+1;
end
end
ux=reshape(ux',M+1,N+1);
eps=abs(ux-U); %精确解与数值解的绝对误差