椭圆方程迭代法介绍

合集下载
  1. 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。

第三章 椭圆型问题的差分法
§3-1 流体力学中的椭圆型问题
·无旋流场中 速度势02=∇ϕ (Laplace Eq.) ·二维不可压定常流动,利用涡-流函数表示:
⎪⎩⎪
⎨⎧-=∇∇=∇⋅+∂∂方程Poisson V t
ζ
ϕζυζζ22
·不可压分离流问题中,扰动压力场:02='∇p
·定常的N -S 方程求解问题
·在网格自动生成中,求解椭圆型方程的网格生成方法 由于椭圆型方程的数学性质:求解域内部任何一点的解函数依赖于所有边界上的边界条件,因此从数值计算方法来看,就不能从一部分边界起步进行推进计算到另外的边界,这与发展方程的求解方法有很大的差别,椭圆型方程的数值求解方法,只能是在整个流场中进行迭代计算来求解。

§3-2 椭圆型问题的迭代法求解
(一)迭代法的基本概念
例:方程ζψ=∇2 ( Poisson 方程) 二维
ζψψ=∂∂+
∂∂2
22
2y
x
差分离散
j i j i j i j i j
i j i j i y x ,2
1
,1,1,2
,1,,1)
(2)
(2ζψψψψψψ=∆+-+
∆+---+-+ ………………..(*)
写成矩阵形式代数方程组为: B A =ψ (1)
其中 ⎥⎥⎥
⎥⎥
⎥⎦

⎢⎢⎢⎢

⎢⎣
⎡=
A ⎥⎥⎥⎥⎥⎥⎥⎦

⎢⎢⎢⎢⎢⎢⎢⎣⎡=ψJ I j i ,,2,11,1ψψψψ 一般地,对于线性方程组有B A =ψ,欲求未知函数ψ的解矢量 若A 为非奇异矩阵,即:∃1
-A ,则B A 1
-=ψ
由于A 是个阶数甚大的矩阵(非三对角),直接求解,或利用Gauss 消去法求逆矩
阵,计算量及所需计算机的内存都将十分巨大,所以在实际计算中不希望采用直接法求解。

迭代法的基本思想是:定义一个序列)
(k ψ,当∞→k 时,B A k 1
)(-→ψ,从而得到
方程(1)的解。

迭代法设法给出[]
)()2()1()(,,,,r k k k k k B A F ---ψψψ=ψ 的迭代关系。

(通常为计算方便,迭代法采取1=r ,使之简单)
[]
)1()(,,-ψ=ψk k k B A F
若k F (即迭代关系式)与迭代步k 无关,则称为平稳迭代; 若k F 是)
1(-ψ
k 的线性函数关系,则称为线性迭代。

例如最简单的线性迭代关系可设为:)()1()()
(k k k k V H +ψ=ψ- (2)
若迭代是有效的,则)()
(k k V H
+ψ=ψ
即)(1)
(1
k k V B A H
B A +=-- (3)
B M B A H I V k M k k k )(1)()()
()(=-=∴-
即 B M H k k k )()()
(+ψ=ψ
且 1)
()(--=A H I M
k
或 I MA H =+
● 研究迭代的收敛性:
引入误差:B A E
k k 1)()
(--ψ=
而由(2)-(3)得:)(1)1()(1)
(B A H B A k k k ----ψ=-ψ
即 )1()()
(-⋅=k k k E H E
或有递推关系式:)
0()2()1()()
(E H H H HE H E H E k k k k k

⋅==⋅=⋅=-- 由于)
0(E
是初始解与精确解的误差,应是一个有界的任意函数,故迭代矩阵H 应具有:
当∞→k 时,0lim =⋅∞
→Z H H H k
k
,Z 为任意的有界向量函数。

可以证明:(参阅 “偏微分方程的有限差分方法”P239)
● 对于任意的向量Z ,
0)1()1()
(→-Z H H H k k 的充分必要条件是H 的所有的特征值i λ的
绝对值(即谱半径)都小于1。

推论 当k 很大时,)(~)()
1(H E E
k k ρ+ i i
H λρm a x )(=
所以若1~ρ,则迭代法的收敛速率很慢。

二、几种迭代法介绍
1. Jacobi 迭代 (简单点迭代)
由方程 B A =ψ
将矩阵分解为:A=L+D+U
L :主对角线以下的元素 ij a (i>j 时等于A ,其余为零) D : 主对角线元素
U : 主对角线以上的元素 ij a (i<j 时等于A ,其余为零)
B U D L =ψ++)(
B U D L k k k =ψ+ψ+ψ--)1()()1(
B D U L D k k 1)1(1)()(---+ψ+-=ψ∴
)(1D L D H +-=∴-, 1-=D M , B D V 1-=
H ,M 可以验证满足迭代有效性条件,即I MA H =+
2、Gauss-Seidel 点迭代
类似1 但是 B U D L k k k =ψ
+ψ+ψ-)
1()()( B D L UW D L k k 1)1(1)()()(---+++-=ψ∴
在实际计算中
L 中(i>j )只要遵循已有新值时,用新值,没有新值时
用旧值,即为G-S 。

*往返扫描的Gauss-Seidel 迭代,即
step1: B U D L k k k =ψ
+ψ+ψ-)
1()()( step2: B U D L k k k =ψ+ψ
+ψ++)1()
1()(
3、SOR (逐点松弛迭代)
step1. 用G -S 迭代法求中间值,即
B U D L k k k =ψ+ψ+ψ-)1()()( …………………………………(a)
step2. )1((*))
()1(-ψ-+ψ=ψ
k k ωω …. …………………………….(b)
消去中间结果(*)
ψ 即 )1((*))()1()(-ψ-+ψ=ψ⇒⋅k k D D D D b ωω 将(a)代入 )1()1()()()1(][--ψ-+ψ-ψ-=ψk k k k D U L B D ωω
B
U D D U B L D k k k k ωωωωωωω+ψ
--=ψ-+ψ-=ψ+---)
1()
1()1()(])1[()1()(
B L D U D L D k k ωωωωω1)1(1)()(])1[()(---++ψ--+=ψ∴
其中ω为松弛因此, 10<<ω为亚松弛,
21<<ω时为超松弛。

4、线迭代和线松弛迭代
将A 分解为U L D A ++=
但 ⎥⎥⎥⎦
⎤⎢⎢⎢⎣⎡= D 保留主对角元素在D 中,L ,U 则仍为余下元素的上三角与下三角矩阵
则B U D L k k k =ψ+ψ+ψ--)1()()1(→导出线迭代
而B
U D L k k k =ψ+ψ+ψ-)1()()(→导出线D -S 迭代
*往返扫描的G -S 线迭代(线松弛迭代)
而⎩⎨⎧ψ
-+ψ=ψ=ψ+ψ+ψ--)
1((*))()1((*))()1(k k k k D D B U D L ωω→导出松弛迭代
三、迭代法的收敛性及松弛因子的选择
1. 迭代法收敛的几个充分条件
对于方程 B A =ψ
① 若矩阵A 满足强对角优势条件,则Jacobi 迭代和G -S 迭代均收敛
② 若矩阵A 满足对角占优条件,且矩阵A 为不可约矩阵,则Jacobi 迭代和G -S 迭代均收敛
③ 若矩阵A 是对称正定矩阵,则G -S 迭代收敛。

④ 若),2,1(0
N i a ii =≠且有
12
)(12
≤∑≠=ii N
i j j ij a a 则Jacobi ,G –S 迭代收敛
⑤ 若对于B A =ψ Jacobi 迭代收敛的,则21<<ω的松弛迭代也总是收敛的。

只证明①(余略)
A 为强对角占优,及11>∑≠=N
i
j j ij
ii
a a
将A 的每一行元素均用该行的主对角线元素去除,可得到主对角元素为1,且不改变将对角占优的性质,B A '=ψ',然后将 A '分解为 U D L A ++=',
且有:
[]
1<+∑ij ij
U L
对于Jacobi 迭代,B U L B D U L D k k k +ψ+-='+ψ+-=ψ----)1(1)1(1)()()(
即)(U L G +-=
利用矩阵的特征值分布定理(Gerschgorin 圆盘定理),可知A '的所有特征值均在单位圆内,证毕!
2、对于Poisson 方程Jacobi 迭代矩阵的特征分析
结论:1、Jacobi 迭代矩阵的特征值为:(参考苏煜诚,吴启光,偏微分方程数值解) 1
,,2,11
,,2,1cos 21cos 21,-=-=+=
l r m s l r
m s r s ππλ
x 方向总网格数为s+1(0,1,2,…,s),0,s 为边界 y 方向总网格数为l+1(0,1,2,…,l),0,l 为边界
2、逐点松弛迭代法中迭代矩阵的特征值(G -S 或SOR )
由})1{()()(1
U D L D G ωωωω--+=- 设)(ωG 的特征值为η
则结论为:0)1(,,,=-+-ωηωληs r s r s r (1) 3、SOR 方法中松弛因子ω的最优化迭代
⎥⎥⎦

⎢⎢⎣⎡
)
(,,,)(max Min )(ωρωωηωηs r s r opt s r Max =谱半径
由0)1(=-+-ωηωλη,使1<η的充要条件是:
11<-ω(20<<ω) )1(1-+<ωωλ
结论:2*
2
*
2
*112)
11(2λ
λ
λω-+=
--=
opt (2)
r s r
s ,,*max λλ= Jacobi 迭代矩阵的谱半径
并有2*
2
*
11111)(λ
λωωη-+--=
-=opt opt
但是由于实际过程中)]([*J G ρλ,未知,所以opt ω不能预先知晓。

4、优选松弛因子ω的两个近似方法
方法1:利用的关系仍为上面之(1)(2)两式
步骤① 取200<<ω用SOR 迭代计算若干步,然后用下面的计算近似
的)(0ωη
∑∑-+--=
ij
k j
i k j
i ij k j i k j i u
u
u u
)1(,)(,)
(,)
1(,0)(ωη
②以0ω,)(0ωη代入(1)式求η
ωωηλ⋅+-=2
02
02
)0()1(][B
③根据(2)
2)1(112B
o p t λω-+=
即为)
1(opt ω的第一次近似值
可以类似求出)
1(ω
,)
2(ω,…,直至)
(k ω
,)
1(+k ω
之差小于ε为止。

方法2、令 )
1()(,)(m a x --=k ij
k ij j
i k u u δ 由1=ω开始,近似认为)2()
1()
1()(---==k k k k δ
δδδλ
直至a k k k k ∆
---=)2()1()1()(~δ
δδδ
取2
112a
opt -+=
ω
5、几种主要的迭代算法的收敛速度比较
设而为问题求解域为 ππ≤≤≤≤y x 00 n
h π= 内点共有2
)1(-n 个
a. Jacobi 迭代:
))cos()(cos(2
1
sh rh rs +=λ
21,1,,2
1
1cosh max h s r s r J -====∴λλρ
所以收敛速度 2
12
1ln h R J ≈-=ρ
b. G -S 迭代
由 0)1(=-+-ωηωλ
η
当1=ω时,即为G -S ,2λη=
所以 2
2m a x m a x J
rs S G ρληρ===- 2ln 2ln h R J S G S G ≈-=-=--ρρ
c .SOR h G opt 21sinh
1sinh
11111)]([2
*
2
*-=+-=
-+--=
λλωρ
h R 2ln =-=∴ρ
§3-3 定常问题的迭代法求解与(伪不定常)时间推进法
计算的一致性讨论
一、 概述
例1,定常方程x
u
x u a 22∂∂=∂∂ν (*)
采用Jacobi 迭代,差分格式用中心差分
2
)
(1
)1()(1)
(1
)(122x u u u x
u u a
k j k j k j k j k j ∆+-=∆--++-+υ
)(2
1)(4)(1)(1)(1)(1)1(k j k j k j k j k j u u u u x a u -+-++++-∆-
=υ )2(21)(4)
(1)()(1)(1)(1)()1(k j k j k j k j k j k j k j u u u u u x a u u -+-+++-+-∆-=-υ
2
)
(22)(1)(12)
()1(22)(4x u x x u u x a u u k j
x k j k j k j k j ∆∆+
∆-∆-=∆-∆-++δυτ
τ
其中τ∆可以视为虚拟的时间步长
)(2)(2)(22
222
2)
(x O x
u x x O x u x a O u
k j
∆+∂∂∆+∆+∂∂∆-=∆+∂∂∆υττ
τ
2
22222x u x x u x a u ∂∂∆∆=∂∂∆∆+∂∂ττυτ 即从方程222
2x
u
x u a u x ∂∂=∂∂+∂∂∆∆υττυ 出发的FTCS 格式,与从方程(*)出发的Jacobia 迭代得到相同的差分。

或稳定条件:1222
≤∆=∆∆∆∆=∆∆'='υ
ττυx a x x a x t a c
2
122
2
2=∆∆∆∆=∆∆'=x x x t s ττυ 12'2≤s c ⇒ 12≤∆υ
x a 稳定条件即取: 2≤∆=
∆υx
a R x e
例2
ςψψ
-∇=∂∂2t
k ij
k
ij
y k
ij
x k ij
k ij y
x
t
ςψδψδψψ-∆+
∆=
∆-+2
22
21 为简单起见 令∆=∆=∆y x
]4[2,1,!,,!,!21k
ij k j i k j i k j i k j i k j i k ij k ij t ςψψψψψψψ∆--+++∆
∆+
=-+-++ 若取 41=s 即4
1
2=∆∆t
][4
1]
4[4121,!,,!,!2,1,!,,!,!1k ij k j i k j i k j i k j i k
ij k j i k j i k j i k j i k j i k ij k ij ςψψψψςψψψψψψψ∆-+++=∆--++++=-+-+-+-++ 将k ,k+1视为相邻两个迭代步的解,则上式是原方程ςψ=∇2
的Jacobi 迭代,即
t
∂∂ψ
可视为一个虚拟时间项(时间相关!)。

例3 0=∆ϕ
采用线松弛迭代求解椭圆方程(y 方向是隐式求解,x 方向是G -S 迭代) step1:
0222
(*)
1
,(*),(*)1,2
)1(,1(*),)1(,1=∆+-+
∆+--+---+y x j i j i j i n j
i j i n j i ϕϕϕϕϕϕ (1)
step2:
)
1(,(*),)(,)1(--+=n j
i j i n j i ϕωωϕϕ (2)
由(2)式 )1(,)(,(*)
,)1
1(1
--
+=
n j i n j i j i ϕω
ϕωϕ 代入到(1)
[][]
0)11(11)1
1(22
)1(,)(,22)
(,1)1(,)(,)1(,1=⎭⎬⎫⎩⎨⎧-+∆+
∆+-
--
----+
II
n j i y n j i y I
n j
i n j i n j i n j i y x ϕδωϕδωϕϕω
ϕω
ϕ[]
))(12(1)()(21)11(12121)1(,)(,2)1(,1)(,1)1(,1)(,12)1(22)(222)
(,2)1(,2-----++----∆+---∆=⎪⎭⎪⎬⎫⎪⎩⎪⎨⎧⎪⎪⎭⎫ ⎝⎛∂∂-+⎪⎪⎭⎫ ⎝⎛∂∂+∆+n j i n j i n j i n j i n j i n j i n n n j i x n j i x x x y y x ϕϕω
ϕϕϕϕϕωϕωϕδϕδ 右边:
j
i j
i n j
i n j
i n j
i n j i n j i n j i n j i n j i n j i t
x
x
t x t x x t
t x x x
,2
,21,21,11,1)1(,)(,2)1(,1)(,1)1(,1)(,12)
12
()
12(2)()12()()(2∂∂-∆∆+
∂∂∂∆∆=
∂∂-∆∆+∆∂∂-∂∂⋅∆∆=∆--∆∆+⎥⎥⎦
⎤⎢⎢⎣⎡∆--∆-∆∆=----+-----++ϕ
ω
τϕτϕ
ω
τ
ϕϕττϕϕωττϕϕτϕϕτ
所以原方程相容于:222222)12(y
x t x x t x ∂∂+∂∂=∂∂-∆∆+∂∂∂∆∆ϕ
ϕϕωτϕτ 或
][212222222y
x t x x t x t ∂∂+∂∂-∆∆=∂∂∂-∆+∂∂ϕϕωϕωωϕ 注意到粘性系数应为正,021
>-ω
2<ω
二、 计算流体力学中的伪时间推进方法
定常、不可压时,连续方程 0=⋅∇V
c o n s t
=ρ 改用
02=⋅∇+∂∂V c t p
+∇=∇⋅+∂∂p V V t V
ρ1保留此项。

相关文档
最新文档