椭圆型方程的差分格式

考虑Laplace方程 ∂2u + ∂2u = 0 (x, y)∈ Ω
∂x 2 ∂y 2
(3.3)
设Ω为正方形区域，0<x<1，0<y<1，求方程(3.3)
满足边值条件
的解。
u(x, y) = f (x, y) (x, y)∈ ∂Ω
(3.4)
因此Laplace方程的五点差分格式为
( ) 1
h2
U l +1,m + U l −1,m + U l ,m +1 + U l ,m −1 − 4U l ,m
Ul−1,m +Ul+1,m +Ul,m+1 +Ul,m−1 −4Ul,m =0 l, m = 1,2,…, M −1
( ) Ul,1 +Ul+1,0 2 +Ul−1,0 2 − 2 + hql Ul,0 = hf1,l l =1,2, , M −1
( ) U1,m +U0,m+1 2 +U0,m−1 2 − 2 + hpm U0,m = hf0,m m = 1,2, , M −1
( ) ( ) U1,0
2
+U0,1 2
−
2 + hp0 + hq0
2 U0,0 = h
f0,0 + f1,0
2
令
[ ] U = U0,0,U1,0,
,UM−1,0;U0,1,
,UM−1,1;
,U0,M−1,
,U T M−1,M−1
下面把所有差分方程写成矩阵形式，于是U满足方程
AU=hg
(3.25)
⎨ ⎪
∂u
= g (x , y )
⎪⎩ ∂ n ∂ Ω
解 令h=1/2，应用图3.2中结点次序，则方程(3.18)为
⎡ 4 − 2 0 − 2 0 0 0 0 0 ⎤⎡U1 ⎤
⎢⎢−1 4 −1 0 − 2 0
0
0
0
⎥ ⎥
⎢⎢U
2
⎥ ⎥
⎡2g1 ⎤
⎢ ⎢
g2
⎥ ⎥
⎢0 −2 4 ⎢⎢−1 0 0
0 0 −2 0 0 4 − 2 0 −1 0
] ; g0,M −1,0, ,0, gM ,M −1;2g0,M , g1,M , , gM −1,M ,2gM ,M T
gl,m = g(lh,mh)
例 3.1 在单位正方形区域Ω上解Laplace方程Neumann 问题
⎧ ∂ 2u ⎪⎪ ∂ x 2
+
∂ 2u ∂y 2
=
0
0 ≤ x, y ≤ 1
( ) 2U1,m +U0,m+1 +U0,m−1 − 4 + 2hpm U0,m = 2hf0,m m = 1,2, , M −1
相似地对于y=0上的结点(lh,0)，我们有 ( ) 2Ul,1 +Ul+1,0 +Ul−1,0 − 4 + 2hql Ul,0 = 2hf1,l l =1,2, , M −1
第三章 椭圆型方程的差分方法
• 3.1 正方形区域中的Laplace方程Dirichlet边值问 题的差分模拟
• 3.2 Neumann边值问题的差分模拟 • 3.3 混合边值条件 • 3.4 非矩形区域 • 3.5 极坐标形式的差分格式 • 3.6 矩形区域上的Poisson方程的五点差分逼近
的敛速分析
(3.1)
其中，系数a(x,y)，b(x,y)，c(x,y)满足
b2 −ac< 0 (x, y)∈Ω
(3.2)
对应方程(3.1)的定解问题有下面三类：
第一边值问题，或称Drichlet问题
⎧ 方程(3.1) (x, y)∈ Ω ⎩⎨u = f (x, y) (x, y)∈ ∂Ω
第二边值问题，或称Neumann问题
(3.10)
∂u 表示函数u沿着边界的外法线方向导数。在正方
形的∂n ∂四Ω 个顶点上法向没有定义，事实上，g(x,y)在那里
将不连续，以后将取平均值作为不连续点上的值的定
义。
Neumann边值问题(3.10)的解存在，仅当
∫ g(x, y)dl＝0 ∂Ω
且除了一个任意常数外，解唯一。因为容易看到，如 果u(x,y)是式(3.10)的解，于是， u(x,y)+C(C是一个任 意常数)也是其解。为了唯一性，需要规定u(x,y)在区域 中某一点上的值。
其中矩阵A为M 2 阶对称方阵。
⎡E0 K
⎤
⎢ ⎢
K
E1
K
⎥ ⎥
⎢ A=⎢
K E2 K
⎥ ⎥
⎢
⎥
⎢ ⎢
K EM −2
K⎥ ⎥
⎢⎣
K EM −1 ⎥⎦
⎡−(2+ pm) 1
⎤
⎢ ⎢
1
ห้องสมุดไป่ตู้
−4 1
⎥ ⎥
Em
=
⎢ ⎢
⎢
⎥ m =1,2, M −1
1
−4
1
⎥ ⎥
⎢⎣
1 −4⎥⎦
E0,E1, ,EM−1,K 是M阶方阵。其中
=0
它具有截断误差:
1 12
h 2 ⎜⎜⎝⎛
∂ 4u ∂x 4
+
∂ 4u ∂y 4
⎟⎟⎠⎞ l ,m
+
我们引进记号◇，有
◇Ul,m
=
1 h2
(Ul+1,m
+Ul−1,m
+Ul,m+1
+
Ul,m−1
−
4Ul,m
)
因此差分方程(3.6)即◇Ul,m =0。
（3.6） (3.7)
如图3.1所示
在区域Ω的每一内部结点(l,m)上
⎪⎩⎪⎨⎧∂∂un方=程g((3x.,1y))
(x, y)∈ Ω (x, y)∈ ∂Ω
第三边值问题，或称Robin问题
其中
⎪⎧方程(3.1) ⎪⎩⎨α (x, y)u +
β
(x,
y)
∂u ∂n
=
γ
(x,
y)
(x, y)∈ Ω (x, y)∈ ∂Ω
α(x, y), β (x, y) > 0
3.1 正方形区域中的Laplace方程Dirichlet边 值问题的差分模拟
⎤
⎢⎢−1 4 −1
⎥ ⎥
B=⎢
⎥
⎢ ⎢
−1 4 −1⎥⎥
⎢⎣
−1 4 ⎥⎦
3.2 Neumann边值问题的差分模拟
现在我们考虑Laplace方程Neumann边值问题，即
⎪⎪⎩⎪⎪⎨⎧∂∂∂∂unx2u2|∂+Ω=∂∂y2gu2(x=, y0) (x, y)∈Ω;Ω={(x, y)| 0< x <1,0< y <1}
=
0
l , m ＝ 1，2 ，… , M − 1
在x=0上的导数边值条件的差分模拟为
( ) 1
2h
U −1,m −U1,m
= g0,m
这里 g0,m = g (0, mh )。
m = 1,2,
M −1
（3.13）
在五点差分格式(3.12)中令l=0，于是有
( ) δ
2 x
+
δ
2 y
U0,m
=
0
即
U−1,m =2U0,m −U1,m −U0,m−1 +2U0,m −U0,m+1
,U T M−1,M−1
单位正方形中的内部结点上的 (M −1)2 个线性方程
(3.8)写成矩阵形式为
AU=K
（3.9）
其中，A是 (M −1)2 阶方阵
⎡ B −I
⎤
⎢⎢− I B − I
⎥ ⎥
A=⎢
⎥
⎢ ⎢
−I
B
−
I
⎥ ⎥
⎢⎣
− I B ⎥⎦
I 是(M-1)阶单位方阵；B是(M-1)阶方阵。
⎡ 4 −1
⎡4 −2
⎤
⎢ ⎢
−
1
4
−1
⎥ ⎥
B=⎢
⎥
⎢ ⎢
− 1 4 − 1⎥⎥
⎢⎣
− 2 4 ⎥⎦
方程组(3.18)中的向量U和g由以下给出：
[ ] U=U0,0,U1,0, ,UM,0;U0,1, ,UM,1; ;U0,M, ,UM,M T [ g = 2 g0,0 , g1,0 , , g M −1,0 ,2 g M ,0 ; g0,1,0, ,0, g M ,1;
4Ul,M − 2Ul,M−1 −Ul−1,M −Ul+1,M = 2hgl,M l =1, , M −1
在四个顶点上，有
4U0,0 −2U1,0 −2U0,1 = 4hg0,0
4U0,M −2U1,M −2U0,M−1 = 4hg0,M 4UM,0 −2UM,1 −2UM−1,0 =4hgM,0 4U M ,M − 2U M −1,M − 2U M ,M −1 = 4hgM ,M
∂u ∂x
−
p(y)u
=
f0(y)
的差分模拟和五点差分公式，即
U−1,m +U1,m +U0,m+1 +U0,m−1 −4U0,m =0
( ) ( ) U1,m −U −1,m 2h − pmU0,m = f0,m pm = p(mh); f0,m = f0 (mh) （3.22）
消去U −1,m ，得
⎢⎡− ⎢
⎢⎣⎡1+
1 2
h(
p0
+
q0
)⎥⎦⎤
1 2
⎤ ⎥ ⎥
⎢
1
⎢
2
− (2 + hq1 )
1 2
⎥ ⎥
E0 = ⎢
⎥
⎢ ⎢ ⎢
1 2
− (2 + hqM−2 )
0 0
⎥ ⎥ ⎥
⎢U ⎢⎢U
3 4
⎥ ⎥ ⎥
⎢⎢⎢2gg43
⎥ ⎥ ⎥
⎢0 ⎢
−1 0
−1
4
−1 0
−1
0
⎥ ⎥
⎢U ⎢
5
⎥ ⎥
=
2h⎢ ⎢
0
⎥ ⎥
⎢ 0 0 −1 0 − 2 4 0 0 −1⎥⎢U 6 ⎥ ⎢ g6 ⎥
⎢ ⎢
0
0
0 −2 0
0
4
−2
0
⎥ ⎥
⎢⎢U
7
⎥ ⎥
⎢⎢2
g
7
⎥ ⎥
⎢ 0 0 0 0 − 2 0 −1 4 −1⎥⎢U8 ⎥ ⎢⎣ 0 0 0 0 0 − 2 0 − 2 4 ⎥⎦⎢⎣U 9 ⎥⎦
其中，ql = q(lh), f1,l = f1(lh) 。
(3.23)
在原点(0,0)上，两边值条件相遇，则
U−1,0 +U1,0 +U0,−1 +U0,1 − 4U0,0 = 0
U1,0 −U−1,0 −2hp0U0,0 = 2hf0,0
U 0,1 − U 0,−1 − 2hq0U 0,0 = 2hf1,0
Neumann边值问题的差分模拟
先在区域Ω中给定一个正方形网格区域，步长为
h，Mh=1,于是必须确定解的结点为（M +1)2个，结点上
的差分方程的解为
U l,m (0 ≤ l, m ≤ M )
（3.12）
在内点上，Laplace方程由差分方程(3.6)代替：
( ) 1
h2
δ
2 x
+
δ
2 y
U l,m
消去 和 ,则 U−1,0
U 0,−1
( ) ( ) 2U1,0 + 2U0,1 − 4 + 2hp0 + 2hq0 U0,0 = 2h f0,0 + f1,0
(3.24)
且对l=0和m=0上成立的方程(3.22)，(3.23)用1/2乘 之，对l=m=0上的方程(3.24)用1/4乘之。
这样在整个计算区域及相应边界网格点上建立了差 分方程：
⎧x = 1;0 ≤ y ≤ 1
⎨ ⎩
y
=
1;0
≤
x
≤
1
p, q, f0 , f1, g 是给定的函数。
在求解区域Ω内由逼近Laplace方程的五点差分公式
( ) 1
h2
δx2
+δ
2 y
Ul,m
=
0，l,
m
=
1,2,…,
M
−1
给出函数u在结点(lh,mh)的近似值Ulm所满足的差分方程。
对于在x=0上的结点(0,mh)，应用边值条件
(l = 1, , M −1; m = 1, , M −1。)
建立差分方程，由此在区域Ω内部 (M −1)2 个点上建立(M −1)2 个方程。
◇Ul,m =0 l,m=1, ,M−1
（3.8）
定义向量
[ ] U = U1,1,U2,1,
,UM−1,1;U1,2,U2,2,
,UM−1,2;U1,M−1,U2,M−1,
代入式(3.13)，则
4U0,m −2U1,m −U0,m−1 −U0,m+1 =2hg0,m m=1, ,M−1
（3.14）
同理，在x=1,y=0,y=1时分别有
4UM,m −2UM−1,m −UM,m−1 −UM,m+1 = 2hgM,m m=1, ,M −1 4Ul,0 − 2Ul,1 −Ul−1,0 −Ul+1,0 = 2hgl,0 l = 1, , M −1
• 3.7 一般二阶线性椭圆型方程差分逼近及其性质 研究
设 Ω 是平面中的具有边界的一个有界区域，本章 考虑如下椭圆型方程的差分解法：
a (x, y ) ∂ 2u
∂x 2
+
2b(x, y ) ∂ 2u
∂x∂y
+ c(x, y ) ∂ 2u
∂y 2
=
d
⎜⎜⎝⎛
x,
y,
u,
∂u ∂x
,
∂u ∂y
⎟⎟⎠⎞
（3.15） （3.16） （3.17）
由此，正方形 0 ≤ x, y ≤ 1区域的 (M +1)2 个结点上差分 方程解Ulm 满足线性方程组
AU=2hg
(3.18)
这里A是(M +1)2阶方阵
⎡ B − 2I
⎤
⎢ ⎢
−
I
B
−I
⎥ ⎥
A=⎢
⎥
⎢ ⎢
−I
B
−
I
⎥ ⎥
⎢⎣
− 2 I B ⎥⎦
I是(M+1)阶单位方阵；B 是如下(M+1)的阶方阵：
⎪ ⎪
∂x
2
+
∂2u ∂y 2
=
0
(x, y)∈ Ω;Ω = {(x, y) | 0 < x, y < 1}
⎪∂u
⎪ ⎨ ⎪
∂x ∂u
− −
p( y )u q(x)u
⎪ ∂y
= =
f0(y) f1(x)
x = 0;0 ≤ y ≤ 1 y = 0;0 ≤ x ≤ 1
（3.21）
⎪ ⎪
u
=
g
(x,
y)
⎩
⎢ g8 ⎥ ⎢⎣2g9 ⎥⎦
(3.19)
或简写成 AU=2hg。 显然A是一奇异矩阵。
3.3 混合边值问题
在xy平面的某区域Ω中，未知函数u满足Laplace方 程，将边界 ∂Ω分成若干弧段，要求u在每一弧段上满足 不同类型的边界条件。讨论此类定解问题的差分模拟。
例如，求解如下定解问题：
⎧ ∂2u