時間差分スキーム(4)

GFDワークゼミ 時間差分スキーム(4) 1

時間差分スキーム (4)

アダムス - バッシュフォース (Adams-Bashforth) スキームの安定性と 位相

ここでは2次のアダムス-バッシュフォーススキームについて振動方程式に当ては めた場合の安定性と位相比について述べる. 2 次のアダムスバッシュフォースス キームは

Uⁿ⁺¹ =Uⁿ+ ∆t (3

2fⁿ−1 2fⁿ⁻¹

)

と表せる. 今,f =iωU とすると, Uⁿ⁺¹ =Uⁿ+ ∆t

(3

2fⁿ−1 2fⁿ⁻¹

)

=Uⁿ+iω∆t (3

2Uⁿ− 1 2Uⁿ⁻¹

)

. (1)

増幅係数λを用いてUⁿ⁺¹とUⁿをUⁿ⁻¹で表すと, Uⁿ =λUⁿ⁻¹,

Uⁿ⁺¹ =λUⁿ =λ²Uⁿ⁻¹. これらを(1)式に代入すると

λ²Uⁿ⁻¹ =λUⁿ⁻¹+iω∆t (3

2λUⁿ⁻¹−1 2Uⁿ⁻¹

)

= (

1 + 3

2(iω∆t) )

λUⁿ⁻¹− 1

2(iω∆t)Uⁿ⁻¹ となる. ここでp≡ω∆tとすると

λ²− (

1 +i3 2p

)

λ+i1

2p= 0. (2)

2014˙0911˙kawahara.tex 2014/09/11(川原 健史)

GFDワークゼミ 時間差分スキーム(4) 2

(2)式からアダムス-バッシュフォーススキームもリープフロッグスキームと同様に 2つのλをもつ. (2)式をλについて解くと,

λ₁ = 1 2

[ 1 +i3

2p+

√ 1− 9

4p²+ip ]

,

λ₂ = 1 2

[ 1 +i3

2p−

√ 1−9

4p²+ip ]

.

(3)

p → 0のとき, λ₁ → 1, λ₂ → 0である. したがって, λ₁ に対応するモードを物理 モード, λ₂ に対応するモードを計算モードと解釈できる. pが十分小さいとき, 計 算モードは減衰する. これはアダムス-バッシュフォーススキームの利点である. そ こで, |p| < 1のときの λ₁ とλ₂ の振る舞いを(3)式を用いて調べる. (3)式 の根号 の部分をテイラー展開し,地道に計算すると,

λ₁ = 1 +ip− 1

2p²+i1 4p³− 1

8p⁴+· · · , λ₂ =i1

2p+1

2p²−i1 4p³+1

8p⁴+· · ·. 実部と虚部に分けて表すと,

λ₁ = (

1− 1 2p²− 1

8p⁴− · · · )

+i (

p+1

4p³+· · · )

, λ₂ =

(1 2p²+ 1

8p⁴+· · · )

+i (1

2p− 1

4p³− · · · )

2014˙0911˙kawahara.tex 2014/09/11(川原 健史)

GFDワークゼミ 時間差分スキーム(4) 3

となる. ここで,増幅率を計算すると,

|λ₁|=√ λ₁λ^∗₁

∼

√{(

1− 1 2p²− 1

8p⁴ )

+i (

p+ 1 4p³

)} {(

1− 1 2p²− 1

8p⁴ )

−i (

p+1 4p³

)}

=

√(

1− 1 2p²− 1

8p⁴ )2

+ (

p+1 4p³

)2

=

√ 1 + 1

2p⁴+o(p⁶)

∼

√ 1 + 1

2p⁴

∼1 + 1

4p² >1,

|λ₂|=√ λ₂λ^∗₂

∼

√(

i1 2p+1

2p² ) (

−i1 2p+ 1

2p² )

=

√1 4p⁴ +1

4p²

= p 2

√1 +p²

∼ p 2

( 1 + 1

2p² )

<1

となる. つまりアダムス-バッシュフォーススキームの物理モードは不安定であり, 計算モードは安定で減衰する.

また,リープフロッグスキームの時と同様にアダムスバッシュフォームスキームの 物理モードの位相の振る舞いを調べると,

2014˙0911˙kawahara.tex 2014/09/11(川原 健史)

GFDワークゼミ 時間差分スキーム(4) 4

θ₁ p = 1

parctanλ_1,im λ_1,re

∼ 1

parctan p+¹₄p³ 1−¹₂p²− ¹₈p⁴

= 1 p

{ p+¹₄p³

1− ¹₂p²−¹₈p⁴ − 1 3

( p+ ¹₄p³ 1− ¹₂p²− ¹₈p⁴

)3

+1 5

( p+¹₄p³ 1−¹₂p² −¹₈p⁴

)5

+· · · }

∼ 1 + ¹₄p²

1− ¹₂p²− ¹₈p⁴ − 1 3p

( p+¹₄p³ 1−¹₂p²− ¹₈p⁴

)3

= (

1 + 1 4p²

) ( 1 + 1

2p² +o(p⁴) )

− 1 3p

( p+ 1

4p³ )3(

1 + 1

2p² +o(p⁴) )3

= 1 + 1 2p²+ 1

4p²− 1

3p²+o(p³)

= 1 + 5

12p²+o(p³)>1

となり, 真の解よりも位相が早く進む. ここで, 4行目から5行目の変形には,以下 のテーラー展開の式を用いた.

1

1− ¹₂p²− ¹₈p⁴ = 1 + (1

2p²+ 1 8p⁴

) + 1

2 (1

2p²+1 8p⁴

)2

+· · ·

= 1 + 1

2p²+o(p⁴)

2014˙0911˙kawahara.tex 2014/09/11(川原 健史)