緩和法を定義する行列 前処理行列とみなせる 反復法の一般形より P の選び方を考える
P ∆Un
τ = −ω(SUn+Qn) = −ωRn Un+1 = Un −τ ωP−1(SUn+Qn)
これは演算子 B = τ P−1S に対するリチャードソン法と考えられる 反復行列 G = 1 − τ P−1S の固有値は関係式
λ(G) = 1 − ωλ(τ P −1S) (110)
最適緩和係数, G のスペクトル半径は
ωopt = 2
λmin(τ P −1S) + λmax(τ P −1S) ρ(G) = κ(τ P−1s) − 1
κ(τ P−1s) + 1 (111)
しかしリチャードソン法の収束率は低い 非常に大きい条件数
実際問題での条件数 = 105 ∼ 106 も稀でない メッシュ間隔の大きい変化
↓
さまざまな方法が考案 狙い : τ P−1S の固有値の大きさが揃うように P を選ぶ
不完全コレスキー
(Choleski)分解
Meijerink & Van der Vorst
Un+1 = Un − τ ωP −1(SUn + Qn) (112) Strongly Implicit Procedure (SIP)
Stone, Schneider & Zedan
λ(G) = 1 − ωλ(τ P −1S) (113) CG (Conjugate Gradient)
法
Reid, Concus, Kershaw
ωopt = 2
λmin(τ P−1S) + λmax(τ P−1S) (114)
GMRES (Generalized Minimum Residual) 法 Saad & Schultz
ρ(G) = κ(τ P −1s) − 1
κ(τ P −1s) + 1 (115)
10.4
非線形問題
非線形問題
S(U) = −Q (116)
ニュートンの線形化
S(Un+1) = S(Un + ∆U) = S(Un) +
µ∂S
∂U
¶
∆U = −Q (115)
ヤコビアン (Jacobian)
J(U) ≡
µ∂S
∂U
¶
(116)
反復法 (ニュートン法)
J(U)∆Un = −Rn (117)
直接法で解けば厳密解 U¯ が得られる筈
U¯ = Un − J(U)−1Rn (118)
誤差 en = Un − U¯ の方程式
J(U)en = −Rn (119)
式(117) J(U)∆Un = −Rn を条件演算子 P/τ で表わされる反復法で解く
P
τ ∆Un = −Rn (120)
誤差 en は演算子 G で増幅される
en+1 = Un+1−U¯ = en+∆U = en−τ P−1Rn = (1−τ P−1J)en ≡ Gen (121)
非線形問題では条件演算子 P はヤコビアンを近似 ニュートン法は
G = 1 − τ P −1J (122)
のスペクトル半径が1以下ならば収束する
10.5
マルチグリッド法
Brandt (1972,1977,1982), Thomas et al. (2003), Briggs (2000)
反復法 漸近的に遅い収束 低い周波数の誤差 = 長い波長の誤差
誤差の漸近的な振る舞い = 増幅行列Gの1に近い固有値 = 空間離散化演算子Sの 絶対値が最も小さい固有値
ラプラス演算子の|ΩJ|min 低波数成分φx = φy = π/M
↓
空間離散化Sの低周波数誤差成分が反復法では最も減衰しにくい 一方で高波数成分は数回の反復で効果的に減衰
それぞれのフーリエ・モードは(λJ)nの割合で減衰 (Gの固有値)
もしλJ = 0.5ならn回の反復後に2nだけ減少 λJ = 0.999 → 101 に23000回の反復
-6
x
Error after smoothing
10.5.1 平滑化 (smoothing) スムージング・ファクター
µ = max
π/2≤φ≤π |λ(G)| (123)
ヤコビ法: φx = φy = π で −1 ヤコビ過緩和法 高周波数領域で λ(GJ) = −1,+1/2
µ(GJ(ω)) = max[|1 − 2ω|, |1 − ω/2|] (124)
緩和係数 ω = 4/5 で µ = 0.6
線ヤコビ法 0.6
ガウス・ザイデル法 0.5 線ガウス・ザイデル法 0.447 レッド・ブラック点緩和法 0.25 ゼブラ線緩和法 0.25 交替ゼブラ法 0.048
10.5.2 CGC (Coarse Grid Correction) 法
h: 細かいメッシュ, H: 粗いメッシュ, ex. H = 2h 細かいメッシュ上の問題
ShUh = −Qh (125)
反復法で解く (残差形式, residual form)
P ∆Uh = −Rh (126)
1. 細かいメッシュから粗いメッシュへの制約 (restriction)
RH = IhHRh (127)
2. 粗いメッシュでの問題を解く
SHUH = −QH (128)
SH∆UH = −RH (129)
3. 修正量 ∆UH を粗いメッシュから細かいメッシュへ延長 (prolongation)する
∆UH = IHh ∆UH (130)
10.5.3 二重格子法 Two-grid Iteration Method
ふたつのグリッドh, H 上のマルチグリッド法
1. 細かいグリッド上の解Uhnに平滑演算子S1でn1回の緩和スイープをおこなう 2. 粗いグリッドによる修正によりUhn+1 = Uhn + ∆Uhを得る
3. 細かいグリッド上の解Uhn+1に平滑演算子S2でn2回の緩和スイープをおこなう
}
} HH }
HHHH
HHHHHj ©©©©©©©©©©©*
n1
S1 S2 n2
coarse grid correstion
収束特性
漸近収束率がグリッドサイズh (したがってグリッド点数 N)によらない 全計算量はN に線形に変化する (実質的な計算量はスイープ部分∝ N)
10.5.4 マルチグリッド法
FAS (Full Approximation Scheme) 与えられた問題
Sh(Uh) = −Qh (131)
Sh(Uh + ∆Uh) − Sh(Uh) = −Rh (132)
制約演算子 IhH, 延長演算子 IˆHh
RH = IhHRh (133)
UH = ˆIhHUH (134)
粗いグリッド上の問題
SH(UH + ∆UH) − SH(UH) = −RH (135)
代数方程式の反復解法 まとめ 反復解法についての短い導入
• それぞれの反復法 (緩和法とも呼ばれる) は収束行列 P で特長づけられる
• 反復法は増幅行列 G = 1 − P−1S の最も小さい固有値を効果的に減衰させる
• G の小さい固有値は空間離散化行列 S の高い振動数に対応する 振動をフーリエ・モードで表すと
• 前処理法には多くの種類がある
• 現在知られている方法のなかではマルチグリッド法が他の方法と比べてはるかに 効果的で一般性がある
メッシュ数に依存しない収束率
計算コスト O(N) という理想的状態