偏微分课程课件9_椭圆型方程的有限差分方法(I)

lg
10 9
lg 16 9
...
...
17
...
计算机实现
• 五点差分格式的系数矩阵 • 算法 • 程序
返回 18
• 五点差分格式的系数矩阵
对possion方程第一边值问题建立五点差分格式
1 h2
[ui , j1
ui 1, j
4uij
ui 1, j
ui ,
j1 ]
f ij
1 r2
2u
2
f (r, ),
(3.17)
并对得到的方程在
D
(r, )
r
r,
j
1 2
j
1
2
上积分
j1
2 j1
2
r
r
r
u r
1 r
2u
2
dr
d
r
j
1 2
rf
(
r
,
)d
dr
j1
2
左
第五章 椭圆型方程的差分方法
(一) Poisson方程 （二）差分格式的性质 （三）边界条件的处理 （四）变系数方程 （五）双调和方程 （六）特征值问题
1
(一) Poisson方程
2u 2u u x2 y2 f ( x, y), ( x, y) D
D {( x, y) | 0 x a,0 y b}, D {( x, y) | x 0,a,0 y b, y 0,b,0 x a}
d
hn d n
u( x h, y) (1 h dx
L
n!
dxn
L
)u( x, y)
hd e dx u( x, y)
ex 1 x x2 x3 L xn L
1 2 3!
n!
u1 =e u0 , u2 =e u0 , u3 =e u0 , u4 =e u0 , u5 =e u0等
f ( xi , y j ) fi, j
hu( xi , y j ) u( xi , y j )
u( xi h, y j ) 2u( xi , y j ) u( xi h, y j ) h2
+
u( xi
,
yj
k)
2u( xi , k2
yj
)
u( xi
,
yj
k)
2u
2u
x 2
u( xi ,
1
3 4
u1 4u2 u4 6 7 2
u3 u4
5
u1 4u3 u4 2 3 1
u1 u2
6
u2 4u4 u3 4 5
联立, 有线性代数方程组
0
8 7 1 x
4
1
1
0
1 4 0 1
1 0 4 1
0
1
1
4
u1 u2 u3 u4
1
6
2 4
8 7 3 5
f1,1
j2 j 1
j 1时
j 2 ...
4 1 0 L 1 0 L 0
I个
I个 ...
(1,1)对应H的第一行 12
分析系数矩阵H i 1,L , I; j 1,L , J;
1 h2 [ui, j1 ui1, j 4ui, j ui1, j ui, j1 ] fi, j
对于第二个结点(2,1),
x,
y)
log
(1
x2
)
y2
( x, y) D
D {( x, y) | 0 x, y 1}
解：
ui1, j
2ui, j h2
ui1, j
ui ,
j1
2ui , k2
j
ui, j1
fi, j
i 1,L , I; j 1,L , J;
改写为代数方程, 仅考虑 h k
1 h2 [ui, j1 ui1, j 4ui, j ui1, j ui, j1 ] fi, j
4
4！
1
2
2
3
3！
4
4！
1
2
2
3
3！
4
4！)
O(h6 ) S2 u5 u6 u7 u8 O(h6 )
S1
=4u0
+h2u0
+
1 12
h4
(
2
2D4
)u0
+O(h6
)
S2
=4u0
+2h2
u0
+
1 6
h4
(
2
+4D4
)u0
+O(h6
)
消去D4u0得
u0
4S1
S2 6h2
20u0
1 12
yj)
y 2
u( xi ,
yj
)
1 12
h2
4 x4
k2
4 y4
u( xi ,
yj
)
O(h4
k4)
所以五点差分格式是相容的.
24
2.九点差分格式 令 =h ,=h ,D2 = 2 ,
x y xy
625 30 1
那么有 2 + 2 =h2, =h2 D2
748
4 + 4 =( 2 + 2 )2 2 2 2 =h4 (2 2D4 )
fij
其中fij =f (xi ,y j ),
则 T( xi , yj ) Lh[u]i, j [Lu]i, j O(h2 k2 )
ຫໍສະໝຸດ Baidu uij
ui1, j 2uij ui1, j h2
ui , j1
2uij k2
ui , j1
fij
(i, j)
由于上差分方程中 只出现u在(i,j)及其四个 临点上的值, 故称为五 点差分格式
h2 2 u0
O(h4 )
由u0 f 2u0 f
4S1
S2
20u0
6h2
f
1 h4f 2
四阶精度
27
3 极坐标形式的差分格式
若求解域是圆环、环形域或扇形域，
则采用极坐标是方便的，此时Poisson方程形如
r, u
1
(r
u )
r r r
1 r2
2u
2
f (r, ),
(3.17)
r x2 y2 , tan y .
任意点 (ri , j )(i 0), 用中心差商公式
1 u
[ r
r
(r
r
)](ri , j )
r , u
1 r
r
(r
u) r
1 r2
2u
2
1 ri
r
i
1
u(ri
1
,
j
)
(r
i
1
2
2
r
i
1
)u(ri
,
j
)
2
(r )2
r
i
1
u(
ri
1
,
j
)
2
1 [r2
2u
]2 (ri , j )
1 ri2
i 1,L , I; j 1,L , J
将双下标数列ui,j按自然顺序排成单下标uk,那么
1 h2 [ukI uk1 4uk uk1 ukI ] fk
k 1,L , IJ
19
这样有代数方程组Hu=g, 其中
1 h2 [ukI
uk1 4uk uk1 ukI ]
fk
k 1,L , IJ
20
这样有代数方程组Hu=g, 其中
B I
H
1 h2
I
B O
I O
O
I B I
I B
4 1
1
4
1
H是稀疏矩阵, 每一
行最多5个数不为0
B
OOO
H是对角占优矩阵
1 4 1
1 4 H是对称正定
21
课堂练习
1.试求出解 Poisson 方程 u f ( x, y) 的五点差分格式
取h k 1 3
y
1
3 4
对于第一个结点,
2
u3 u4
5
1
u2
8
h2
u3
4u1
0
1
0
即 4u1 u2 u3 1 8
u1 u2
6
8 7 1 x
4u1 u2 u3 1 8 u( x, y) lg[(1 x)2 y2], ( x, y) D
y 同理, 对于另外三个结点, 分别有
对于前I个结点, 系数矩阵H为
j 1
j2
4 1
1
B
1
4
1
OOO
I1 O
L L
1 4 1
1
1 4
1
五点差分格式
i 1,L , I; j 1,L , J;
1 h2 [ui, j1 ui1, j 4uij ui1, j ui, j1 ] fij
等价于
1 h2
Huh
g
B I
I
u(ri , j1 ) 2u(ri , j ) u(ri , j1 ) ( )2
r, u
1 (r u) r r r
1 r2
2u
2
f (r, ),
(3.17)
代到(3.17),就i 1得到逼近它的差分方程
1
r u i 1 i1, j 2
(r
i
1
2
r
i
1
)uij
2
r
i
1
ui
1,
j
2
ri
(r )2
B
I
H
OOO
I B I I B
B,I为I维方阵, H为I*J维方阵
2u 2u
例：
x
2
y2
0,
( x, y) D [0,1][0,1]
u( x, y) lg[(1 x)2 y2 ], ( x, y) D
1
h2 [ui, j1 ui1, j 4uij ui1, j ui, j1 ] fij
x
{0 r , 0 2 }
方程(3.17)的系数于r 0奇异，因此只当r 0时
有意义。为了定出有意义的解，需补充u于r 0有
界的条件, 可设u满足
u lim r 0. r0 r
(3.18)
r, 分别取等步长 r和 .令
ri (i 0.5)r, i 0,1, 2,L ,
j ( j 1) , j 0,1,L , J 1, h 2 / J
1
于是差分方程为: h2 Huh g
11
分析系数矩阵H
i 1,L , I; j 1,L , J;
1 h2 [ui, j1 ui1, j 4ui, j ui1, j ui, j1 ] fi, j
未知节点 对于第一个结点(1,1),
1 h2
[u1,0
u0,1
4u1,1
u2,1
u1,2
]
u1 =e u0 , u2 =e u0 , u3 =e u0 , u4 =e u0 , u5 =e u0
其中 =h ,=h ,D2 = 2 ,
x y xy
625
S1 u1 u2 u3 u4
30 1
e u0 e u0 e u0 e u0
748
(1
2
2
3
3！
4
4！
1
2
2
3
3！
1 ri2
ui1, j
2ui, j ui, j1
( )2
fi, j .
下面用有限体积法获得在点 (r0 , j ) 的 差分方程 ri (i 0.5)r, i 0,1, 2,L ,
j ( j 1) , j 0,1,L , J 1, h 2 / J
用r乘方程
r, u
1 (r u) r r r
8
Possion方程第一边值问题的差分逼近为
huij =fij , (xi，y j ) Dh
uij
=
ij
,
(xi，y j ) Dh
其中uij u(xi , y j ) = (xi , y j )=ij , (x, y) D.
例：五差分格式求解
2u 2u
x
2
y2
0
(x, y) D
u(
内点 : I J
a
x 返回 5
1.五点差分格式
利用Taylor级数展开有
1 h2
u( xi
h,
yj
)
2u( xi
,
yj
)
u( xi
h,
yj
)
2u x 2
ij
h2 24
4
x
4
u(1 ,
yj
)
4 x4
u(2 ,
yj
)
xi1 1 ,2 xi1
1 k2
u( xi ,
yj
k)
2u( xi ,
2
区域以及边界离散 h a , k b I 1 J 1
内点为
Dh
( xi
,
yj
):
xi yj
ih, i 1 i jk, j 1
I; jJ
.
边界点为
( xi , y j ) : xi ih, y j jk;
Dh
i 0, I 1, j 0,1,L , J , J 1
1 h2 [u2,0 u1,1 4u2,1 u3,1 u2,2 ] f2,1
j2 j 1
4 1 0 L 1 0 L 0
1
4 1 L
j 1
01L j2
0
...
I个
I个
...
(2,1)节点对应H的第二行
分析系数矩阵H i 1,L , I; j 1,L , J;
1 h2 [ui, j1 ui1, j 4ui, j ui1, j ui, j1 ] fij
j 0, J 1, i 0,1, L , I , I 1
u u( x, y), ( x, y) [0, a][0, b]
y b
内点
总:(I 2)(J 2)
边界点
a
x 返回 4
u u( x, y), ( x, y) [0, a][0, b]
y
b
L uIJ
uI1 uI2 L u2I u1 u2 L uI
先仅下标 j 小的未知量放前面, 相同的 j 的再按 i
从小到大顺利排列
i 1,L , J; j 1,L , I;
按自然顺序排列网点(i,j)
j 1, i 1,L , I; j 2, i 1,L , I;
L j J , i 1,L , I;
定义向量
uh u1,1,L , uI ,1, u1,2 ,L , uI ,2 ,L u1,J ,L , uI ,J
1 h2 [ui, j1 ui1, j 4uij ui1, j ui, j1 ] fij
的局部截断误差（相容性）。
实用性分析
对于双曲与抛物方程：
微分方程
相 容 性
差分格式 稳定性
真解u=u(x,y)
收 敛 性
真解u=un
椭圆型方程收敛性分析区别于双曲与抛物方程
23
相容性
Ri, j u Lh[u]i, j [Lu]i, j
yj)
u( xi ,
yj
k )
2u y2
ij
h2 24
4
y
4
u( xi
,1 )
4 x4
u( xi
,2
)
y j1 1 ,2 y j1
Lu u
2u x 2
2u y2
f (x, y)
Lhuij
huij
ui1, j
2uij h2
ui1, j
ui ,
j1
2uij k2
ui, j1