隐式差分方程课件

这是一个六点差分格式（如图2.7所示），称为加权六点差分格式。

m-1,n+1 m,n+1 m+1,n+1
图2.7：
1
m-1,n m,n m+1,n
显然，当 0时，加权六点格式为古典显示格式； 1 当 2 时，加权六点格式为Crank-Nicolson隐式格式； 当 1 时，加权六点格式为古典隐式格式。 加权六点格式亦可直接由差商代替导数得到
0.413 637 929 568
0.171 096 336 778
0.000 683
0.000 564
-0.000 000 001 712
-0.000 000 001 417
Hale Waihona Puke Baidu
0.629 273 956 459
0.012 108 818 740
0.000 194
0.000 100
-0.000 000 000 485
表2.2
tn t1 t2 t4 t8 t16 t 80 t160 t 320 t 640 t 800
Sr
0.994 497 915 630 0.489 026 104 192 0.978 172 634 773 0.956 821 703 419 0.915 507 772 134 0.643 146 895 793
-0.000 000 000 257
2.3.3 加权六点隐式格式
前面，我们已经推导了热传导方程（2.26）的古典显示格式，古 典显示格式及Crank-Nicolson格式等。实际上，它们都可以作为本 节推导的加权六点隐式格式的特殊情形。 n 1 2 n u exp( kD ) u m x m 由 2 n 1 2 n exp(kDx )um exp(( 1 )kDx )um ,1 1 得 1 2 2 4 1 2 n 1 2 2 2 4 n [ 1 kD k D ] u [ 1 ( 1 ) kD ( 1 ) k D ] u x x m x x m 即 2 2 1 2 2 D 2 x x 两边去掉高于二阶导数的项，且用 h 代替 ，则得差分格式
与显式差分格式不同，隐式差分格式中 包括了（n+1）时间层上二个或二个以上结 n 1 n1 n 1 点处的未知值（例如 U m1 ,U m ,U m1 ）， 使用隐式差分格式和使用显式差分格式求 解完全不同。相对而言，使用隐式差分格 式求解，每时间层包含有较多的计算工作 量。从后面对差分格式的稳定性分析可知， 隐式格式的优点在于，其稳定性要求对步 长比的限制大为放宽，而这正是我们所期 望的。
2 x
n m
（2.44） 由于格式（2.44）中包括六个结点，故也可 称为六点格式（如图2.6所示）。
m-1,n+1 m,n+1 m+1,n+1
1 1 n 1 n 1 n 1 n n n (1 r )U m r (U m U ) ( 1 r ) U r ( U U 1 m 1 m m 1 m 1 ) 2 2
n
n 1 )U m ] 1 2
1 2 x [(1 2
1 6ra
1 n 2 m
（2.49.1）
格式（2.49.1）具有截断误差阶 (h k ) ，可写 成更方便的形式
4 2
1 n 2 2 n 2 1 1 n 2 2 n1 1 n 2 2 n 2 1 1 n 2 2 n [1 am x (am ) ram x ]U m [1 am x (am ) ram x ]U m 12 2 12 2
n 1 n (1 r x2 )U m [1 (1 )r x2 ]U m
n 1 n n n 1 n 1 n ( 1 2 r ) U r ( 1 )( U U ) r ( U U ) [ 1 2 r ( 1 )] U m m 1 m 1 m 1 m 1 m ,0 1 （2.46） 或者
2.3.1 古典隐式格式
现在对热传导方程 n1 2 n 推导其最简单的隐式差分逼近——古典隐式格式。由 um exp(kDx )um exp( kD )u u 故
2 x n 1 m n m 2 (1 kDx
u 2u t x 2
1 2 4 n 1 n k Dx )u m um 2
ECN
0.000011 0.000 022 0.000 040 0.000 079 0.000 151 0.000 531
ED
0.000 000 000 026 -0.000 000 000 051 -0.000 000 000 101 -0.000 000 000 198 -0.000 000 000 379 -0.000 000 000 331
（2.45） 它称为Douglas差分格式，具有截断误差阶 (k 2 h4 )。
1 1 2 n 1 1 1 2 n [1 (r ) x ]U m [1 (r ) x ]U m 2 6 2 6
例2.1 解初边值问题
u 2u t x 2 u t 0 sin x u (0, t ) u ( , t ) 0
t k n 1 n 1 n 1 2u n 1 um 1 2u m u m 1 ( 2 )m x h2
n 1 m n 1 m n m
n 1 Um
代入微分方程，得到格式（2.41）。
2.3.2 Crank-Nicolson隐式格式
exp( 1 1 n 1 n kL)um exp( kL)um 2 2
2 2 4 1 n 2 m 2 2 n 1 n x2 (um um )(
u 2u a ( x, t ) 2 t x
)m
2

h2 (
)m
2
(h 4 k 2 )
（2.48）
令
代入式（2.48），则 因此得差分方程
1 ra
1 n 2 m n 1 n (U m Um )
式中左边如果仅保留二阶导数项，且以 ，则得差分 k 2 n 1 n ( 1 ) U U x m m 格式 h2 n 1 n 1 n 1 n 或者 rU m （2.41） 1 (1 2r )U m rU m 1 U m 格式用图2.5表示，其截断误差阶为 (k 2 h2 ) ,与古典差分格式相同。
n 1 n Um Um 1 2 n 1 1 2 n 2 x U m (1 ) 2 x U m k h h
2.3.4 系数依赖于x,t的一维热传导方程的一 个隐式格式的推导
考虑方程 （2.47） 的差分逼近。 1 已知 x 2 sinh( hDx ) 4 2 h 2 2 2 2 由其Taylor展开式可得 x h Dx Dx (h 6 ) 12 据此，可得 1 2u n 1 1 4u n 1
{
0 x ;0 t T 0 x 0t T
应用（1） Crank-Nicolson差分格式，（2） Douglas差分格式解上述问题。对每一种情况， 令 h 20, r 1 20 （r的这个值对Douglas格式 有最小的截断误差），由初值条件和边值条件通 n U 过上述二个格式的每一个逐层求出 m 的值。一般 而言，当由第n层去求第（n+1）层的解时，二个 格式的每一个都需解一线性代数方程组，其系数 是三对角阵，可用追赶法求解（见2.4）。已知上 述定解问题的理论解，记为 Sr ， 有 u e t sin x 记 SCN , S D 分别为用高速数字计算机解出的CrankNicolson格式的解，而 ECN Sr SCN , ED Sr SD 分别 x 表示它们对精确解的误差，在 ，时间层n 2 上， t n nk 。 它们的值由表2.2给出。
Crank-Nicolson隐式差分格式是解热传导方程（2.26） 的常用的差分格式，为了推导它，由式（2.24），有 由 得
2 L Dx
[1
1 1 1 1 1 1 2 2 2 n 1 2 2 2 n kDx ( kDx ) ]um [1 kDx ( kDx ) ]u m 2 2 2 2 2 2
m-1,n+1 m,n+1 m+1,n+1
1 2 2 x h 2 替代 Dx
图2.5：
m,n
为了求得第（n+1）时间层上的 的值， 必须通过解线性代数方程组。这是一个隐 式差分格式，必须联合其初边值条件求解。 格式（2.41）通常称为古典隐式格式。 2 D , D 我们也可以通过直接用差分算子代替 x x u u u 的方法，即 ( )
1
1
1
1
1
1
（2.49.2）
这是一个隐式差分格式（如图2.8所示）。
m-1,n+1 m,n+1 m+1,n+1
图2.8
m-1,n m,n m+1,n
x 12 x 1 1 u h 1 u n 2 （ ） ( ) m （h 4 k 2） 2 a t 12 x a t 1 2 1 u n 2 1 2 1 1 n1 n ( ) m 2 x [ 1 (um um ) (h 2 k 2 ) 2 n k x a t h am 2 2h
2 x
（2.42） 1 两边仅保留前二项，用 h 代替 D ，则得差分格式 1 1 (1 r )U (1 r )U （2.43） 2 2 这是一个隐式差分格式，称为Crank-Nicolson差分格式， 截断误差阶为 (k 2 h2 ) ，也可写为
2
2 x
2 x
n 1 m
代入微分方程（2.26），得到Crank-Nicolson格式。
基于如同Crank-Nicolson格式一样的六个网 格结点可获得另一精度较高的差分格式，如在 2 D 前式（2.42）中仅保留直到 x 的项，即有
(1
由式(2.19.3),可令 则可得 代入上式，则有如下差分格式：
1 1 2 n 1 2 n kDx )um (1 kDx )um 2 2 1 2 1 2 n 2 n Dx um 2 x (1 x )um h 12 1 2 1 2 1 2 Dx 2 x (1 x ) h 12
1 2 n 1 n x (um um ) 2 2h
1
am 2 k
n
1
n 1 n (um um )
1 1 2 1 n 1 n x [ 1 (um um )] (h 4 k 2 ) 2 n 12r h am 2
n )U m ]
1 2 x [(1 2
1 6ram
图2.6
m-1,n m,n m+1,n
n 1 n n 1 n 1 n 1 n n n um um 2u n 2 1 u m 也可将( u )nm 2 um 1 2um um 1 1 2um um 1 ( 2 )m [ ] 2 2 t k t 2 h h 1 1