IV (2)

数値流体力学特論

COMPUTATIONAL FLUID DYNAMICS

(CFD)

IV

_{数値スキームの解析}

(2)

The Analysis of Numerical Schemes (2)

11.

_{代数方程式の反復解法}

Iterative methods for algebraic systems

岩津玲磨

Reima Iwatsu, e-mail : [email protected]

Winter Semester 2007, Graduate School, Tokyo Denki University

IV

_{数値スキームの解析}

(2)

_の構成

10.

_{発展方程式の時間積分法}

10.1 発展方程式の解析 10.2 時間積分スキームの解析 10.3 時間積分法の選択

11.

_{代数方程式の反復解法}

11.1 基本的な反復解法 11.2 過緩和法 11.3 前処理法 11.4 非線形問題 11.5 マルチ・グリッド法

11.

_{代数方程式の反復解法}

Iterative methods for algebraic systems

背景

離散化 → 格子点での変数に関する代数方程式 ⇒ CFDでは大規模

3D RANS 乱流計算 (7変数 × 100万格子点 = 700万未知変数) 陽解法 °, 陰解法, 非圧縮性流体, 陽解法でも非線形問題

代数方程式の解法 : 直接法, 反復法

Varga (1962), Dahlquist & Bjork (1974), Saad (2003) Dongarra

目的とガイドライン

• 基本的な概念の紹介

• もっともシンプルな方法への導入

11.1

_{基本的な反復解法}

11.1.1

_{直交等間隔メッシュ上のポアソン方程式}

∆u = f 0 ≤ x, y ≤ L

(1)

u = g on the bounradies

(2)

2次精度5点差分スキームで離散化

(u

_i+1,j

− 2u

_ij

+ u

_i−1,j

) + (u

_i,j+1

− 2u

_ij

+ u

_i−1,j

) = f

_ij

∆x

2

(3)

-6 } } } } }

0 1

i

M

x

0

1

j

M

y

行列表示

U = (u

₁₁

, u

₂₁

, · · · , u

_M,1

, u

_1,2

, u

_2,2

, · · · , u

_1,j

, u

_2,j

, · · · , u

_M,j

, · · · ) (4)

SU = F ∆x

2

+ G ≡ −Q

(5)

     −4 1 · · · 1 1 −4 1 · · · · · · 1 · · · 1 −4 1 · · · 1 · · · · · · ·      ×      u₁₁ ... u_i1 ...      = ∆x 2      f₁₁ ... f_i1 ...      −      g₁₀ + g₀₁ ... g_i0 ...      (6) (i) S は対角優位 (ii) S は対称 (したがって非特異, 正定値) 線形代数の表示 N = M2, I = i + (j − 1)M N

X

J=1

s

_IJ

u

_J

= −Q

_I

I = 1, · · · , N

(7)

11.1.2

_{ポイント・ヤコビ法}

/

_{ポイント・ガウスザイデル法}

ポイント・ヤコビ法 初期近似 点毎順番に値を修正 (スイープ) 近似値 un_ij の修正 un+1_ij (n 反復回数)

u

n+1_ij

=

1

4

(u

n

i+1,j

+ u

ni−1,j

+ u

ni,j+1

+ u

ni,j−1

) +

1

4

q

n ij

(8)

u

n+1_i

=

1

s

_ij



−q

n i

−

N

X

j=1,j6=i

s

_ij

u

n_j





₍₉₎

4

} } } }

4 : Unknown at iteration n + 1

• : Known from iteration n

行列表示 D: 対角行列, E: 下三角行列, F : 上三角行列

S = D + E + F

(10)

DU

n+1

= −Q

n

− (E + F )U

n

(11)

(D + E + F )U = −Q

(12)

∆-フォーム 残差 Rn

R

n

≡ R(U

n

) = SU

n

+ Q

n

= (D + E + F )U

n

+ Q

n

(13)

D∆U

n

= −R

n

(14)

∆U

n

≡ U

n+1

− U

n

(15)

ポイント・ガウスザイデル法

u

n+1_ij

=

1

4

(u

n+1

i+1,j

+ u

ni−1,j

+ u

ni,j+1

+ u

n+1i,j−1

) +

1

4

q

n ij

(17a)

u

n+1_I

=

1

s

_II



−q

n I

−

I−1

_X

J=1

s

_IJ

u

n+1_J

−

N

X

J=I+1

s

_IJ

u

n_J





_(17b)

(D + E)U

n+1

= −Q

n

− F U

n

(18)

(D + E)∆U

n

= −R

n

(19)

4

m m } }

4 : Unknown at iteration n + 1

• : Known from iteration n + 1

° : Known from iteration n

11.1.3

_{反復スキームの収束解析}

反復回数 n → ∞ のときに誤差 en → 0 のとき, 反復法は収束するという

e

n

= U

n

− ¯

U

U :

¯

厳密解

(20)

行列 S の任意の分割

S = P + A

(21)

反復法

P U

n+1

= −Q

n

− AU

n

(22)

P ∆U

n

= −R

n

(23)

P : 収束条件行列 ( (22)式が簡単に解けるように選ぶ) 直接法と比較すると

SU

n+1

= −Q

n

or S∆U

n

= −R

n

(24)

厳密解がみたす(25)式を(22)式からひいて誤差に関する式を導く

S ¯

U = (P + A) ¯

U = −Q

(25)

誤差についての式

P e

n+1

= −Ae

n

(26)

e

n+1

= −(P

−1

A)e

n

= −(P

−1

A)

n

e

1

= (1 − P

−1

S)

n

e

1

(27)

もし n を増加させるときに ||en+1|| が0に近づくためには行列 P−1A が収束行列 増幅行列 (反復行列) G

G = 1 − P

−1

S

(28)

スペクトル半径

ρ(G) ≤ 1

(29)

固有値

|λ

_J

(G)| ≤ 1 for all J

(30)

行列 S を他の行列 P にとりかえて解きやすくする 「解きやすい」とはどんなことなのか？ ˆP−1 _{を計算するのに必要な演算数が O(N) より大きくない} ˆP はできるだけ S に似ていること

残差 R と誤差 e の関係式

R

n

= Se

n

(31)

誤差が0に近づくときに残差も0に近づく 実際の計算では誤差が直接わからない (厳密解がわかっていない！) 残差のノルム ||R|| を収束への目安とする ||R|| がマシンゼロでも Un が U¯ と同じ程度正確なわけではない SU が2次精度空間離散化ならば ∆x = 10−2 で U に 10−4 程度の誤差 反復法

U

n+1

= (1 − P

−1

S)U

n

− P

−1

Q = GU

n

− P

−1

Q

(32)

残差と誤差の関係式

e

n+1

= Ge

n

= G

n+1

e

0

(33)

反復回数 n を擬似時間インデックスとみなせば時間発展方程式

ヤコビ法とガウスザイデル法の収束条件

ヤコビ法

G = −D

−1

(E + F ) = 1 − D

−1

S

ガウスザイデル法

G = −(D + E)

−1

F = 1 − (D + E)

−1

S

ヤコビ法,ガウスザイデル法 : S が優位対角ならば収束する S が対称行列ならば ヤコビ法: S または −S, (2D − S) または (S − 2D) が正定値ならば収束 ガウスザイデル法 : S または −S, (2D − S) が正定値ならば収束

収束率の評価

反復の回数とともにどれだけ誤差が小さくなるのだろうか？ n 回の反復による平均誤差減少率

¯

∆e

_n

≡

µ

||e

n

||

||e

0

||

¶

_1/n

(34)

残差と誤差の関係式(33)より

¯

∆e

_n

≤ ||G

n

||

1/n

(35)

Varga (1962)による証明 log s を平均収束率とよぶ

s ≡ lim

n→∞

||G

n

_||

1/n

_{= ρ(G)}

₍₃₆₎

n 回の反復で誤差は1桁 (10倍) 小さくなる

µ

1

10

¶

_n

≤ s = ρ(G)

(37)

n ≥ −1/ log ρ(G)

(38)

スペクトル半径は ρ(G) < 1

11.1.4

_{反復法の固有値解析}

増幅行列 G の周波数応答 (P−1 がどのように S の固有値に近づくか)

U

n

=

X

J=1

(λ

_J

)

n

U

₀

V

(J)

−

X

J=1

Q

_J

Ω

_J

[1 − (λ

J

)

n

_]V

(J)

₍₃₉₎

V (J) : S の固有モード, U0 : 初期成分, ΩJ: Sの固有値, S, G交換,同じ固有関数仮定 過渡状態 (λJ: Gの固有値) 収束

lim

n→∞

U

n

₌

X

J=1

Q

_J

Ω

_J

V

(J)

₍₄₀₎

解 U = −S−1Q の固有値展開 ¯ U = X J=1 UJV (J) (41) S ¯U = X J=1 Ω_JU_JV (J) = −X J=1 Q_JV (J) (42) U_J = −QJ Ω_J (43)

n 回の反復ののちに誤差はどれだけ小さくなるか (e0J: 初期の誤差)

e

n

=

X

J=1

(λ

_J

)

n

e

_0J

V

(J)

(44)

残差も固有関数 V (J) で展開できる (誤差と残差の関係式を使うと)

R

n

=

X

J=1

(λ

_J

)

n

e

_0J

SV

(J)

=

X

J=1

(λ

_J

)

n

Ω

_J

e

_0J

V

(J)

(45)

固有ベクトルの直交を仮定すると残差のL2ノルムは

||R

n

||

_L2

=

X

J=1

[(λ

_J

)

n

Ω

_J

e

_0J

]

2

(46)

反復回数 n を増加させるとともに最も小さな(λJ)がより早く減衰する もし最も大きい(λJ)が1に近ければ漸近的な収束はかなり遅い -6

n

||R||

a

b

G の固有値とSの固有値の関係

λ

_J

= λ(Ω

_J

)

は一般には求めにくい ヤコビ法の場合

G

_J

= 1 − D

−1

S

(47)

D は対角行列なので

λ

_J

= 1 −

1

d

_J

Ω

J

(S)

(48)

11.1.5

_{反復法のフーリエ解析}

周期境界条件ではフーリエ・モードがS, Gの固有ベクトル Sの固有値

Ω = −4(sin

2

φ

_x

/2 + sin

2

φ

_y

/2)

φ

_x

= lπ/M φ

_y

= mπ/M J = l+(m−1)M l, m = 1, · · · , M (49)

GJ の固有ベクトル

λ(G

_J

) = 1 − (sin

2

φ

_x

/2 + sin

2

φ

_y

/2) =

1

2

(cos φ

x

+ cos φ

y

)

(50)

低い周波数(φx, φy ∼ 0)では収束率は低い(λ ' 1) スペクトル半径 (Ω-スペクトルのうちで最も低い周波数に対応)

ρ(G

_J

) = cos πM

(51)

収束率は ∆x = ∆y でメッシュ数が大きいときに

ρ(G

_J

) ' 1 −

π

2

2M

2

= 1 −

π

2

2N

= 1 − O(∆x

2

₎

₍₅₂₎

ガウス・ザイデル法

ノイマンの安定性解析を適用すると

(4 − e

−iφx

_{− e}

−iφy

_)λ

GS

= e

iφx

+ e

iφy

(53)

ガウスザイデル法の増幅行列の固有値はヤコビ法の固有値の2乗に等しい (Young, 1971)

λ(G

_GS

) =

1

4

(cos φ

x

+ cos φ

y

)

2

_{= λ(G}

J

)

2

(54)

スペクトル半径

ρ(G

_GS

) = ρ(G

_J

)

2

= cos

2

π/M

(55)

M が大きいときに

ρ(G

_GS

) ' 1 −

π

2

M

2

= 1 −

π

2

N

(56)

ガウスザイデル法はヤコビ法の2倍速く収束する 誤差を1桁減衰させるのに必要な反復回数は N/(0.43π2) もっといろいろな場合について Varga (1962)

10.2

_過緩和法

収束率を上げたい

修正量 ∆Un = Un+1 − Un をより速くメッシュを伝播させる

過緩和法のもとになっているアイディア (Frankel & Young, 1950) 過緩和法はもとになる反復法によって得られる値Un+1¯ を使って次の回の値を

U

n+1

= ω ¯

U

n+1

+ (1 − ω)U

n

(57)

で求める ω は過緩和係数

新しい修正量を

∆U

n

= U

n+1

− U

n

with ∆U

n

= ω ¯

∆U

n

(58)

で求める 最大の収束率をもつように ω を調節すればかなりの加速が可能

10.2.1

_{ヤコビ過緩和法}

ポアソン方程式

u

n+1_ij

=

ω

4

(u

n

i+1,j

+ u

ni−1,j

+ u

ni,j+1

+ u

ni,j−1

+ q

ijn

) + (1 − ω)u

nij

(59)

行列表示 (演算子表示)

DU

n+1

= −ωQ

n

−ω(E +F )U

n

+(1−ω)DU

n

= −ω(SU

n

+Q

n

)+DU

n

(60)

D∆U

n

= −ωR

n

(61)

増幅行列

G

_J

(ω) = (1 − ω)I + ωG

_J

= 1 − ωD

−1

S

(62)

収束するためには

ρ(G

_J

(ω)) ≤ |1 − ω| + ωρ(G

_J

) < 1

(63)

収束条件

0 < ω <

2

1 + ρ(G

_J

)

(64)

ヤコビ法の固有値とヤコビ過緩和法の固有値の関係は

λ(G

_J

(ω)) = (1 − ω) + ωλ(G

_J

)

(65)

与えられた周波数の範囲内で減衰が最大になるように最適なωの値を選べる 

ω = 1/(1 − λ(G

_J

))

ととれば λ(GJ) = 0 対応する周波数成分は寄与しない 高波数を選択的に減衰させたいとき λ ∼ −1, ω < 1 に選ぶ 低波数を選択的に減衰させたいとき λ ≤ 1, ω ∼ ωmax に選ぶ

ω の最適値 ωopt ⇒ 全周波数域での最適値の平均

−1 + ω

_opt

(1 − λ

_min

) = 1 − ω

_opt

(1 − λ

_max

)

(66)

λ_min に対するグラフとλmax に対するグラフの交点

ω

_opt

=

2

2 − (λ

_min

+ λ

_max

)

(67)

ラプラス演算子に対しては(50)式より

λ

_max

= −λ

_min

= cos π/M

(68)

最適値 ωopt = 1 ヤコビ法が最適 -6 ω |λ(G_J(ω))| 1

1−λ_min ωopt 1−λ1max

A A A A A A A A A A AA¢¢ ¢¢ ¢¢ ¢¢ ¢¢ ¢¢ HH HH HH HH HH HH HH HH HH HH HH_H©©©©©© ©©©© ©© Unstable

10.2.2

_{ガウス・ザイデル過緩和法}

: SOR (Successive

Overrelax-ation)

_法

ポアソン方程式 (ラプラス演算子)

¯

u

n+1_ij

=

ω

4

(u

n

i+1,j

+ u

n+1i−1,j

+ u

ni,j+1

+ u

n+1i,j−1

+ q

ijn

) + (1 − ω)u

nij

u

n+1_ij

= ω ¯

u

n+1_ij

+ (1 − ω)u

n_ij

(69)

演算子形式

D ¯

∆U

n

= −R

n

− E∆U

n

(70)

∆-フォーム (反復演算子 (D + ωE))

(D + ωE)∆U

n

= −ωR

n

(71)

増幅行列

緩和パラメータの条件

三角行列の積 行列積は対角要素の積

det G

_SOR

(ω) = det(I + ωD

−1

E)

−1

det[(1 − ω)I − ωD

−1

F ]

= I det[(1 − ω)I − ωD

−1

F ] = (1 − ω)

N

(73)

一方で行列積は固有値の積

ρ(G

_SOR

(ω))

N

≥ | det G

_SOR

(ω)| = (1 − ω)

N

(74)

スペクトル半径が1以下になるのは

|(1 − ω)| ≤ ρ(G

_SOR

(ω)) < 1

(75)

0 < ω < 2

(76)

のとき 対角優位行列のときSOR法が収束する条件 (Young, 1971)

0 < ω <

2

1 + ρ(G

_GS

)

(77)

対称行列の場合

S =

Ã

D

₁

F

F

T

D

₂

!

(78)

D₁, D₂: ブロック対角行列, λ(GSOR(ω)) = λ(ω) と書く  i) SOR法の増幅行列の固有値

λ(ω) = 1 − ω + ωλ

1/2

(ω)λ(G

_J

)

(79)

ii) ω = 1 の場合にはガウスザイデル法とヤコビ法の固有値の関係式

λ(ω = 1) ≡ λ(G

_GS

) = λ

2

(G

_J

)

(80)

iii) 過緩和係数の最適値

ω

_opt

=

2

1 +

p

1 − ρ

2

(G

_J

)

(81)

ラプラス演算子に対する点ヤコビ法のスペクトル半径

ρ(G

_J

) = cos π/M

(83)

過緩和係数の最適値

ω

_opt

=

2

1 + sin π/M

' 2(1 −

π

M

+

π

2

M

2

)

(84)

最適係数におけるスペクトル半径

ρ(G

_SOR

(ω

_opt

)) =

1 − sin π/M

1 + sin π/M

' 1 −

2π

M

+ O(

1

M

2

)

(85)

誤差のノルムが1桁減少するのに必要な反復回数は 2.3M/π ∼ √N 最適SOR法は誤差をk桁減少させるのに kN√N 演算を要する 最適値の推定

ρ(G) = lim

n→∞

||e

n+1

||

||e

n

||

(86)

n が大きなとき誤差は最も大きな固有値に支配される (79)式より ρ(GJ) を推定 (81)式より ω の最適値

10.2.3 SSOR (Symmetric Successive Overrelaxation)

SORはスイープの方向がいつも同じ バイアスによって誤差が蓄積しないか？ 両方向への交互スイープ 予測値 ∆U¯ n+1はSORスイープ

(D + ωE) ¯

∆U

n+1/2

= −ωR

n

(87)

反対方向へのスイープが次に続く

(D + ωF ) ¯

∆U

n

= −ωR( ¯

U

n+1/2

)

(88)

ここで

¯

U

n+1/2

= U

n

+ ¯

∆U

n+1/2

U

n+1

= ¯

U

n+1/2

+ ¯

∆U

n

= U

n

+ ¯

∆U

n+1/2

+ ¯

∆U

n

(90)

SSOR法はSORと同じ条件で収束 最適緩和係数は

ω

_opt

'

2

1 +

p

2(1 − ρ

2

(G

_J

))

(91)

10.2.4 SLOR (Successive Line Overrelaxation)

修正量 ∆Un が伝播される方法を変えて収束を加速 ラプラス演算子に対するSLOR (VLOR)

¯

u

n+1_ij

=

ω

4

(u

n

i+1,j

+ u

n+1i−1,j

+ ¯

u

ni,j+1

+ ¯

u

ni,j−1

+

1

4

q

n ij

)

u

n+1_ij

= ω ¯

u

n+1_ij

+ (1 − ω)u

n_ij

(92)

4∆U

_ijn

− ∆U

_i,j−1n

− ∆U

_i,j+1n

− ω∆U

_i−1,jn

= ωR

n_ij

(93)

SLORは同じ値のωoptに対して √ 2倍速く収束する 4 4 4 } m 4 : Unknown at iteration n + 1

• : Known from iteration n

レッド

-

ブラック点緩和法

赤色の点 : I = i + (j − 1)M が偶数, 黒色の点 : I が奇数 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 • • • • • • • • • • • • • • • 4 : Red point • : Black point

ゼブラ線緩和法

4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 • • • • • • • • • • • • • • • • • •

10.3

_前処理法

反復法を残差形 (∆-form)に書く 反復回数 n を時間ステップ数とみなす

P

∆U

n

τ

= −ω(SU

n

_{+ Q}

n

_{) = −ωR}

n

₍₉₄₎

τ を時間刻みとした時間発展微分方程式の陽的オイラー法による時間積分

P

dU

dt

= −ω(SU + Q)

(95)

収束行列 (前処理行列) P は増幅行列の固有値がみな1以下

G = 1 − ωτ P

−1

S

(96)

ρ(G) ≤ 1

である限り任意にとれる もうひとつの条件として (−ωτ P−1)S の固有値の実部が非正であること (安定性) P = (ωτ S), G = 0 ならば1回の反復で厳密解が求まる ⇒ P/ωτ が S に近いほどよい反復法 ただし同時に P は解きやすいこと 

10.3.1

_{リチャードソン法}

P

∆U

n

τ

= −ω(SU

n

_{+ Q}

n

_{) = −ωR}

n P/τ = −1 である場合 (−はSの固有値が負の実部をもつため)

U

n+1

= U

n

+ ω(SU

n

+ Q

n

) = (1 + ωS)U

n

+ ωQ

n

≡ G

_R

U

n

+ ωQ

n

(97)

増幅演算子(反復演算子) GR = 1 + ωS GR の固有値と S の固有値 ΩJ の関係

λ

_J

= 1 + ωΩ

_J

(98)

過緩和係数の値に対する安定性の条件

0 < ω <

2

|Ω

_J

|

_max

=

2

ρ(S)

(99)

-6

ω

|λ|

1

Ω_max

ω

opt Ω_min1

A A A A A A A A A A AA¢¢ ¢¢ ¢¢ ¢¢ ¢¢ ¢¢ HH HH HH HH HH HH HH HH HH HH HH_H©©©©©© ©©©© ©©

Unstable

ω

_opt

=

2

|Ω

_J

|

_max

+ |Ω

_J

|

_min

(100)

対応する増幅演算子(反復演算子)のスペクトル半径は

ρ(G

_R

) = 1 −

2|Ω

J

|

max

|Ω

_J

|

_max

+ |Ω

_J

|

_min

=

κ(S) − 1

κ(S) + 1

(101)

行列 S の条件数

κ(S) =

|Ω

J

|

max

|Ω

_J

|

_min

(102)

ラプラス演算子では

|Ω

_J

|

_max

= 8 |Ω

_J

|

_min

= 8 sin

2

π

2M

'

2π

2

条件数

κ(S) =

1

sin

2 π_2M

'

4N

π

2

(104)

M が大きいときになりたつ 50 × 50 点のメッシュで κ(S) ∼ 1000 収束率

ρ(G

_R

) = 1 −

2π

2

N

= 1 −

2

κ(S)

(105)

は点ヤコビ法と類似の(遅い)収束 反復毎に異なる値の ω を 1/|ΩJ|max と 1/|ΩJ|min の間から選ぶ ⇒ non-stationary 緩和法 cf. stationary

10.3.2 ADI (Alternating Direction Implicit)

_法

条件行列 P としてADI因数分解演算子をとる

(1 − τ S

_x

)(1 − τ S

_y

)(1 − τ S

_z

) = τ ω(SU

n

+ Q

n

)

(106)

順番に解く

(1 − τ S

_x

) ¯

∆U = τ ω(SU

n

+ Q

n

)

(1 − τ S

_y

)¯¯∆U = ¯

∆U

(1 − τ S

_z

)∆U =¯¯∆U

(107)

ω, τ の最適値に対してADIは N log N のオーダーで収束する (しかし一般の問題 に対して最適化が難しい) 過緩和係数の最適値 ωopt ' 2 τ の値のガイドラインは

τ

_opt

=

p

1

|Ω

_J

|

_min

|Ω

_J

|

_max

(108)

ラプラス演算子の場合

λ(G

_ADI

) =

(1 − 4τ

1

sin

2

_φ

x

/2)(1 − 4τ

1

sin

2

φ

y

/2)

(1 + 4τ

₁

sin

2

φ

_x

/2)(1 + 4τ

₁

sin

2

φ

_y

/2)

(109)

10.3.3

_{そのほかの前処理法と緩和法}

緩和法を定義する行列 前処理行列とみなせる 反復法の一般形より P の選び方を考える

P

∆U

n

τ

= −ω(SU

n

_{+ Q}

n

_{) = −ωR}

n

_U

n+1

_{= U}

n

_{− τ ωP}

−1

_(SU

n

_{+ Q}

n

₎

これは演算子 B = τ P−1S に対するリチャードソン法と考えられる 反復行列 G = 1 − τ P−1S の固有値は関係式

λ(G) = 1 − ωλ(τ P

−1

S)

(110)

最適緩和係数, G のスペクトル半径は

ω

_opt

=

2

λ

_min

(τ P

−1

S) + λ

_max

(τ P

−1

S)

ρ(G) =

κ(τ P

−1

s) − 1

κ(τ P

−1

s) + 1

(111)

しかしリチャードソン法の収束率は低い 非常に大きい条件数 実際問題での条件数 = 105 ∼ 106 も稀でない メッシュ間隔の大きい変化 ↓ さまざまな方法が考案 狙い : τ P−1S の固有値の大きさが揃うように P を選ぶ

不完全コレスキー

(Choleski)

分解

Meijerink & Van der Vorst

U

n+1

= U

n

− τ ωP

−1

(SU

n

+ Q

n

)

(112)

Strongly Implicit Procedure (SIP)

Stone, Schneider & Zedan

λ(G) = 1 − ωλ(τ P

−1

S)

(113)

CG (Conjugate Gradient)

_法

Reid, Concus, Kershaw

ω

_opt

=

2

λ

_min

(τ P

−1

S) + λ

_max

(τ P

−1

S)

(114)

GMRES (Generalized Minimum Residual) 法

Saad & Schultz

ρ(G) =

κ(τ P

−1

_{s) − 1}

10.4

_{非線形問題}

非線形問題

S(U) = −Q

(116)

ニュートンの線形化

S(U

n+1

) = S(U

n

+ ∆U) = S(U

n

) +

µ

∂S

∂U

¶

∆U = −Q

(115)

ヤコビアン (Jacobian)

J(U) ≡

µ

∂S

∂U

¶

(116)

反復法 (ニュートン法)

J(U)∆U

n

= −R

n

(117)

直接法で解けば厳密解 U¯ が得られる筈

¯

U = U

n

− J(U)

−1

R

n

(118)

誤差 en = Un − ¯U の方程式

J(U)e

n

= −R

n

(119)

式(117) J(U)∆Un = −Rn を条件演算子 P/τ で表わされる反復法で解く

P

τ

∆U

n

_{= −R}

n

₍₁₂₀₎

誤差 en は演算子 G で増幅される

e

n+1

= U

n+1

− ¯

U = e

n

+∆U = e

n

−τ P

−1

R

n

= (1−τ P

−1

J)e

n

≡ Ge

n

(121)

非線形問題では条件演算子 P はヤコビアンを近似 ニュートン法は

G = 1 − τ P

−1

J

(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 → ₁₀1 に23000回の反復 -6 x

10.5.1 平滑化 (smoothing) スムージング・ファクター

µ =

max

π/2≤φ≤π

|λ(G)|

(123)

ヤコビ法: φx = φy = π で −1 ヤコビ過緩和法 高周波数領域で λ(GJ) = −1, +1/2

µ(G

_J

(ω)) = 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 細かいメッシュ上の問題

S

_h

U

_h

= −Q

_h

(125)

反復法で解く (残差形式, residual form)

P ∆U

_h

= −R

_h

(126)

1. 細かいメッシュから粗いメッシュへの制約 (restriction)

R

_H

= I

_hH

R

_h

(127)

2. 粗いメッシュでの問題を解く

S

_H

U

_H

= −Q

_H

(128)

S

_H

∆U

_H

= −R

_H

(129)

3. 修正量 ∆UH を粗いメッシュから細かいメッシュへ延長 (prolongation)する

∆U

_H

= I

_Hh

∆U

_H

(130)

10.5.3 二重格子法 Two-grid Iteration Method ふたつのグリッドh, H 上のマルチグリッド法 1. 細かいグリッド上の解U_hnに平滑演算子S1でn1回の緩和スイープをおこなう 2. 粗いグリッドによる修正によりU_hn+1 = U_hn + ∆Uhを得る 3. 細かいグリッド上の解U_hn+1に平滑演算子S2でn2回の緩和スイープをおこなう } } } HH HH HH HH HHHj _©©©© ©©©© ©©©* n1 S1 S2 n2

coarse grid correstion

収束特性

漸近収束率がグリッドサイズh (したがってグリッド点数 N)によらない 全計算量はN に線形に変化する (実質的な計算量はスイープ部分∝ N)

10.5.4 マルチグリッド法

FAS (Full Approximation Scheme)

与えられた問題

S

_h

(U

_h

) = −Q

_h

(131)

S

_h

(U

_h

+ ∆U

_h

) − S

_h

(U

_h

) = −R

_h

(132)

制約演算子 I_hH, 延長演算子 Iˆ_Hh

R

_H

= I

_hH

R

_h

(133)

U

_H

= ˆ

I

_hH

U

_H

(134)

粗いグリッド上の問題

S

_H

(U

_H

+ ∆U

_H

) − S

_H

(U

_H

) = −R

_H

(135)

代数方程式の反復解法 まとめ 反復解法についての短い導入 • それぞれの反復法 (緩和法とも呼ばれる) は収束行列 P で特長づけられる • 反復法は増幅行列 G = 1 − P−1S の最も小さい固有値を効果的に減衰させる • G の小さい固有値は空間離散化行列 S の高い振動数に対応する 振動をフーリエ・モードで表すと • 前処理法には多くの種類がある • 現在知られている方法のなかではマルチグリッド法が他の方法と比べてはるかに 効果的で一般性がある メッシュ数に依存しない収束率 計算コスト O(N) という理想的状態

さらに進んだ学習のために勧めること n簡単なベンチマーク問題を解くプログラムを自分でつくってみること ˆ非粘性流れ 1次元の衝撃波管 円柱まわりの縮まないポテンシャル流れ ˆ粘性流れ クエット流 平板をすぎる境界層流れ 2次元キャビティ流れ n進んだ読みもの ˆ教科書 数値流体力学シリーズ(全6冊), 河村哲也著「流体解析I」 藤井孝蔵著「流体力学の数値解法」, Toro パタンカー ˆ論文

Chorin (1967,1968), Harlow & Welch (1965) Patankar & Spalding (1972)