GFD ワークノート
ルンゲクッタ型公式の安定性
1
ルンゲクッタ型公式の安定性
ここでは, ルンゲクッタ型公式の, 振動方程式及び一次元移流方程式においての安
定性について議論する.
1 段 1 位のルンゲクッタ型公式の安定性
1 段 1 位のルンゲクッタ型公式とは, オイラー法のことである. 時間差分スキーム
においてオイラー法を用いた場合, 差分式は以下の式で表される.
un+1 = un + ∆tf (un , tn )
(1)
振動方程式
以下の様な振動方程式を考える.
du(t)
= iωu(t).
dt
(2)
ここで, ω は定数である. 式 (2) にオイラー法を適用すると, 以下のような式になる.
un+1 = un + iω∆tun .
増幅係数 λ =
un+1
un
(3)
を導入すると,
2015 0402-kawahara.tex
2015/04/02(川原 健史)
GFD ワークノート
ルンゲクッタ型公式の安定性
λ = 1 + iω∆t.
2
(4)
|λ| の大きさを求めると,
|λ| =
√
1 + (ω∆t)2 > 1
(5)
よって, |λ| > 1 となるため, 改良オイラー法及びホイン法は振動方程式に対して不
安定である.
一次元移流方程式
以下のような一次元移流方程式を考える.
∂u(x, t)
∂u(x, t)
+c
= 0,
∂t
∂x
c = const.
(6)
また, 領域は x = [0, L] の有限領域であり, 境界条件は周期境界とする. ここで, 空
間差分には 2 次の中心差分, 時間差分にはオイラー法を用いて差分化する. このと
き, 空間方向に領域を J 分割し, 周期境界より xJ = x0 とする このとき, 一次元移
流方程式は以下のようになる.
un+1
= unj − c∆t
j
unj+1 − unj−1
2∆x
(7)
ここで, j は空間方向の添字である. u を複素数に拡張し, U を振幅, k を波数として
∞
∑
Ukn e−ikj∆x を代入すると,
unj =
k=−∞
∞
∑
Ukn+1 e−ikj∆x =
k=−∞
∞
∑
k=−∞
Ukn e−ikj∆x −
∞
)
c∆t ∑ n ( −ik(j+1)∆x
Uk e
− e−ik(j−1)∆x .
2∆x k=−∞
(8)
2015 0402-kawahara.tex
2015/04/02(川原 健史)
GFD ワークノート
ルンゲクッタ型公式の安定性
3
′
両辺に eik j∆x をかけ, j について総和を取ると,
J
∞
∑
∑
′
Ukn+1 e−ikj∆x eik j∆x
j=1 k=−∞
=
J
∞
∑
∑
′
Ukn e−ikj∆x eik j∆x
j=1 k=−∞
J
∞
)
c∆t ∑ ∑ n ( −ik(j+1)∆x ik′ j∆x
′
Uk e
e
− e−ik(j−1)∆x eik j∆x ,
−
2∆x j=1 k=−∞
∞
∑
Ukn+1
J
∑
j=1
k=−∞
=
′
ei(k −k)j∆x
∞
∑
Ukn
J
∑
′
ei(k −k)j∆x
j=1
k=−∞
)
∑( ′
c∆t
i(k j−k(j+1))∆x
i(k′ j−k(j−1))∆x
n
e
−e
U
−
2∆x k=−∞ k j=1
∞
∑
=
∞
∑
Ukn
k=−∞
J
∑
J
′
ei(k −k)j∆x
j=1
(
) ∑ i(k′ −k)j∆x
c∆t
−
e
.
Ukn e−ik∆x − eik∆x
2∆x k=−∞
j=1
∞
∑
J
(9)
三角関数の直交性より,
Ukn+1
′
ここで, p =
c∆t
∆x
)
c∆t n ( −ik′ ∆x
ik′ ∆x
=
−
Uk′ e
−e
)
( 2∆x
c∆t
sin k′ ∆x
= Ukn′ 1 + i
∆x
Ukn′
sin k′ ∆x とし, 増幅係数 λ =
Ukn+1
′
Ukn′
を導入すると,
λ = 1 + ip
√
|λ| = 1 + p2 ≥ 1
2015 0402-kawahara.tex
(10)
(11)
2015/04/02(川原 健史)
GFD ワークノート
ルンゲクッタ型公式の安定性
4
よって, |λ| ≥ 1 となるため, オイラー法は一次元移流方程式に対して, k∆x = nπ(n
は整数) となる波数においては中立, それ以外の波数においては不安定である.
2 段 2 位のルンゲクッタ型公式の安定性
2 段 2 位のルンゲクッタ型公式のうち, ここでは改良オイラー法とホイン法の 2 通
りの数値解法について議論する. 時間差分スキームにおいて改良オイラー法を用い
た場合, 差分式は以下の式で表される.
k1n = ∆tf (un , tn ),
1
1
k2n = ∆tf (un + k1n , tn + ∆t),
2
2
n+1
n
n
u
= u + k2 .
(12)
また, ホイン法を用いた場合には, 以下の式で表される.
k1n = ∆tf (un , tn ),
k2n = ∆tf (un + k1n , tn + ∆t),
1
un+1 = un + (k1n + k2n ).
2
(13)
振動方程式
以下の様な振動方程式を考える.
du(t)
= iωu(t).
dt
(14)
式 (14) に式 (12) 及び式 (13) を適用すると, 両式共に以下のような式になる.
1
un+1 = un + iω∆tun − (ω∆t)2 un .
2
2015 0402-kawahara.tex
(15)
2015/04/02(川原 健史)
GFD ワークノート
増幅係数 λ =
un+1
un
ルンゲクッタ型公式の安定性
5
を導入すると,
1
λ = 1 + iω∆t − (ω∆t)2
2
1 2
= 1 + ip − p .
2
(16)
ここで, p = ω∆t とした. ここで, |λ| の大きさを求める.
√
1
(1 − p2 )2 + p2
2
√
1
= 1 + p4 > 1
4
|λ| =
(17)
よって, |λ| > 1 となるため, 改良オイラー法及びホイン法は振動方程式に対して不
安定である.
一次元移流方程式
以下のような一次元移流方程式を考える.
∂u(x, t)
∂u(x, t)
+c
= 0,
∂t
∂x
c = const.
(18)
また, 領域は x = [0, L] の有限領域であり, 境界条件は周期境界とする. ここで, 空
間差分には 2 次の中心差分, 時間差分に改良オイラー法を用いて差分化する. この
とき, 空間方向に領域を J 分割し, 周期境界より xJ = x0 とする このとき, 一次元
移流方程式は以下のようになる.
unj+1 − unj−1
,
2∆x
n
n
unj+1 + 12 k1,j+1
− unj−1 − 12 k1,j−1
= −c∆t
2∆x
n
n
= uj + k2,j .
n
k1,j
= −c∆t
n
k2,j
un+1
j
2015 0402-kawahara.tex
(19)
2015/04/02(川原 健史)
GFD ワークノート
ルンゲクッタ型公式の安定性
6
また, 時間差分にホイン法を用いた場合は, 以下のようになる.
n
k1,j
n
k2,j
un+1
j
unj+1 − unj−1
= −c∆t
,
2∆x
n
n
unj+1 + k1,j+1
− unj−1 − k1,j−1
= −c∆t
2∆x
1
n
n
= unj + (k1,j
+ k2,j
).
2
(20)
この 2 式は, 整理するとどちらも以下のような式になる.
un+1
= unj − c∆t
j
unj+1 − unj−1 1
unj+2 − 2unj + unj−2
+ (c∆t)2
.
2∆x
2
(2∆x)2
(21)
ここで, j は空間方向の添字である. u を複素数に拡張し, U を振幅, k を波数として
∞
∑
n
Ukn e−ikj∆x を代入すると,
uj =
k=−∞
∞
∑
Ukn+1 e−ikj∆x
=
∞
∑
Ukn e−ikj∆x
k=−∞
k=−∞
1
+
2
(
c∆t
2∆x
∞
)
c∆t ∑ n ( −ik(j+1)∆x
−
Uk e
− e−ik(j−1)∆x
2∆x k=−∞
)2 ∑
∞
(
)
Ukn e−ik(j+2)∆x − 2e−ikj∆x + e−ik(j−2)∆x .
k=−∞
(22)
′
両辺に eik j∆x をかけ, j について総和を取ると,
2015 0402-kawahara.tex
2015/04/02(川原 健史)
GFD ワークノート
J
∞
∑
∑
ルンゲクッタ型公式の安定性
7
′
Ukn+1 e−ikj∆x eik j∆x
j=1 k=−∞
=
J
∞
∑
∑
′
Ukn e−ikj∆x eik j∆x
j=1 k=−∞
J
∞
)
c∆t ∑ ∑ n ( −ik(j+1)∆x ik′ j∆x
−ik(j−1)∆x ik′ j∆x
−
U e
e
−e
e
2∆x j=1 k=−∞ k
(
)2 ∑
J
∞
(
)
∑
1 c∆t
′
′
′
Ukn e−ik(j+2)∆x eik j∆x − 2e−ikj∆x eik j∆x + e−ik(j−2)∆x eik j∆x ,
+
2 2∆x
j=1 k=−∞
∞
∑
Ukn+1
J
∑
j=1
k=−∞
=
′
ei(k −k)j∆x
∞
∑
Ukn
J
∑
′
ei(k −k)j∆x
j=1
k=−∞
J (
)
∑
c∆t
′
′
n
ei(k j−k(j+1))∆x − ei(k j−k(j−1))∆x
−
Uk
2∆x k=−∞
j=1
(
)2 ∑
J (
∞
)
∑
1 c∆t
i(k′ j−k(j+2))∆x
i(k′ j−k)j∆x
i(k′ j−k(j−2))∆x
n
+
e
− 2e
+e
U
2 2∆x k=−∞ k j=1
∞
∑
=
∞
∑
Ukn
J
∑
′
ei(k −k)j∆x
j=1
k=−∞
J
( −ik∆x
)∑
c∆t
′
n
ik∆x
−
ei(k −k)j∆x
Uk e
−e
2∆x k=−∞
j=1
(
)2 ∑
∞
J
(
)∑
1 c∆t
′
+
Ukn e−2ik∆x − 2 + e2ik∆x
ei(k −k)j∆x .
2 2∆x k=−∞
j=1
∞
∑
(23)
三角関数の直交性より,
Ukn+1
′
) 1 ( c∆t )2
)
(
c∆t n ( −ik′ ∆x
′
′
ik′ ∆x
+
=
−
−e
Uk′ e
Ukn′ e−2ik ∆x − 2 + e2ik ∆x
2∆x
2 2∆x
(
)
(
)2
c∆t
1
c∆t
= Ukn′ 1 + i
sin k′ ∆x −
sin2 k′ ∆x
(24)
∆x
2 ∆x
Ukn′
2015 0402-kawahara.tex
2015/04/02(川原 健史)
GFD ワークノート
p=
c∆t
∆x
ルンゲクッタ型公式の安定性
sin k′ ∆x とし, 増幅係数 λ =
Ukn+1
′
Ukn′
8
を導入すると,
p2
λ = 1 + ip −
2
√(
)2
p2
|λ| =
+ p2
1−
2
√
p4
= 1+
≥1
4
(25)
よって, |λ| ≥ 1 となるため, 2 段 2 位のルンゲクッタ型公式は一次元移流方程式に
対して, k∆x = nπ(n は整数) となる波数においては中立, それ以外の波数において
は不安定である.
3 段 3 位のルンゲクッタ型公式の安定性
3 段 3 位のルンゲクッタ型公式のうち, ここではホインの 3 次公式について議論す
る. 時間差分におけるホインの 3 次公式は以下の式で表される.
k1n = ∆tf (un , tn ),
1
k2n = ∆tf (un + k1n , tn +
3
2
k3n = ∆tf (un + k2n , tn +
3
1
un+1 = un + (k1n + 3k3n ).
4
1
∆t),
3
2
∆t),
3
(26)
振動方程式
以下の様な振動方程式を考える.
du(t)
= iωu(t).
dt
(27)
式 (27) に式 (26) を適用すると,
2015 0402-kawahara.tex
2015/04/02(川原 健史)
GFD ワークノート
ルンゲクッタ型公式の安定性
1
1
un+1 = un + iω∆tun − (ω∆t)2 un − i(ω∆t)3 un .
2
6
増幅係数 λ =
un+1
un
9
(28)
を導入すると,,
1
1
λ = 1 + iω∆t − (ω∆t)2 − i(ω∆t)3
2
6
1 2 1 3
= 1 + ip − p − ip .
2
6
(29)
ここで, p = ω∆t とした. ここで, |λ| の大きさを求める.
√
1
1
(1 − p2 )2 + (p − p3 )2
2
6
√
1
1
= 1 − p4 + p6
12
36
|λ| =
(30)
√
√
√
√
− 3 ≤ p ≤ 3 のとき, |λ| ≤ 1 となる. よって, − 3 ≤ p ≤ 3 のとき, ホインの 3
次公式は振動方程式に対して安定である.
一次元移流方程式
以下のような一次元移流方程式を考える.
∂u(x, t)
∂u(x, t)
+c
= 0,
∂t
∂x
c = const.
(31)
また, 領域は x = [0, L] の有限領域であり, 境界条件は周期境界とする. ここで, 空
間差分には 2 次の中心差分, 時間差分にはホインの 3 次光式を用いて差分化する.
このとき, 空間方向に領域を J 分割し, 周期境界より xJ = x0 とする このとき, 一
次元移流方程式は以下のようになる.
2015 0402-kawahara.tex
2015/04/02(川原 健史)
GFD ワークノート
ルンゲクッタ型公式の安定性
n
k1,j
n
k2,j
n
k3,j
un+1
j
unj+1 − unj−1
= −c∆t
,
2∆x
n
n
unj+1 + 13 k1,j+1
− unj−1 − 13 k1,j−1
= −c∆t
2∆x
2 n
n
n
uj+1 + 3 k2,j+1 − unj−1 − 23 k2,j−1
= −c∆t
2∆x
)
(
1
n
n
+ 3k3,j
.
= unj +
k1,j
4
10
(32)
この式を整理すると, 以下のような式になる.
unj+2 − 2unj + unj−2
unj+1 − unj−1 1
+ (c∆t)2
2∆x
2
(2∆x)2
unj+3 − 3unj+1 + 3unj−1 − unj−3
1
.
− (c∆t)3
6
(2∆x)3
un+1
= unj − c∆t
j
(33)
ここで, j は空間方向の添字である. u を複素数に拡張し, U を振幅, k を波数として
∞
∑
n
Ukn e−ikj∆x を代入すると,
uj =
k=−∞
∞
∑
Ukn+1 e−ikj∆x
k=−∞
∞
∑
=
Ukn e−ikj∆x
k=−∞
∞
)
c∆t ∑ n ( −ik(j+1)∆x
−
Uk e
− e−ik(j−1)∆x
2∆x k=−∞
(
)2 ∑
∞
)
(
1 c∆t
+
Ukn e−ik(j+2)∆x − 2e−ikj∆x + e−ik(j−2)∆x
2 2∆x k=−∞
(
)3 ∑
∞
(
)
1 c∆t
−
Ukn e−ik(j+3)∆x − 3e−ik(j+1)∆x + 3e−ik(j−1)∆x − e−ik(j−3)∆x .
6 2∆x k=−∞
(34)
′
両辺に eik j∆x をかけ, j について総和を取ると,
2015 0402-kawahara.tex
2015/04/02(川原 健史)
GFD ワークノート
J
∞
∑
∑
ルンゲクッタ型公式の安定性
11
′
Ukn+1 e−ikj∆x eik j∆x
j=1 k=−∞
=
J
∞
∑
∑
′
Ukn e−ikj∆x eik j∆x
j=1 k=−∞
J
∞
)
c∆t ∑ ∑ n ( −ik(j+1)∆x ik′ j∆x
−ik(j−1)∆x ik′ j∆x
−
U e
e
−e
e
2∆x j=1 k=−∞ k
(
)2 ∑
J
∞
(
)
∑
1 c∆t
′
′
′
Ukn e−ik(j+2)∆x eik j∆x − 2e−ikj∆x eik j∆x + e−ik(j−2)∆x eik j∆x
+
2 2∆x
j=1 k=−∞
(
)3 ∑
J
∞
∑
1 c∆t
−
Ukn (
6 2∆x
j=1 k=−∞
′
′
′
′
e−ik(j+3)∆x eik j∆x − 3e−ik(j+1)∆x eik j∆x + 3e−ik(j−1)∆x eik j∆x − e−ik(j−3)∆x eik j∆x
),
(35)
2015 0402-kawahara.tex
2015/04/02(川原 健史)
GFD ワークノート
∞
∑
Ukn+1
J
∑
12
′
ei(k −k)j∆x
j=1
k=−∞
=
ルンゲクッタ型公式の安定性
∞
∑
Ukn
J
∑
′
ei(k −k)j∆x
j=1
k=−∞
J (
)
∑
c∆t
n
i(k′ j−k(j+1))∆x
i(k′ j−k(j−1))∆x
−
U
e
−e
2∆x k=−∞ k j=1
(
)2 ∑
∞
J (
)
∑
1 c∆t
′
′
′
+
Ukn
ei(k j−k(j+2))∆x − 2ei(k j−k)j∆x + ei(k j−k(j−2))∆x
2 2∆x k=−∞
j=1
(
)3 ∑
∞
J
∑
1 c∆t
n
−
U
(
6 2∆x k=−∞ k j=1
∞
∑
′
′
′
′
ei(k j−k(j+3))∆x − 3ei(k j−k)(j+1)∆x + 3ei(k j−k(j−1))∆x − ei(k j−k(j−3))∆x
)
=
∞
∑
k=−∞
Ukn
J
∑
′
ei(k −k)j∆x
j=1
J
( −ik∆x
)∑
c∆t
′
n
ik∆x
−
ei(k −k)j∆x
Uk e
−e
2∆x k=−∞
j=1
(
)2 ∑
J
∞
( −2ik∆x
)∑
1 c∆t
′
n
2ik∆x
ei(k −k)j∆x
+
Uk e
−2+e
2 2∆x k=−∞
j=1
(
)3 ∑
J
∞
( −3ik∆x
)∑
1 c∆t
′
n
−ik∆x
ik∆x
3ik∆x
−
Uk e
− 3e
+ 3e
−e
ei(k −k)j∆x .
6 2∆x k=−∞
j=1
∞
∑
(36)
三角関数の直交性より,
2015 0402-kawahara.tex
2015/04/02(川原 健史)
GFD ワークノート
ルンゲクッタ型公式の安定性
13
)
c∆t n ( −ik′ ∆x
Ukn+1
= Ukn′ −
Uk′ e
− eik∆x
′
2∆x
(
)2
(
)
1 c∆t
′
′
+
Ukn′ e−2ik ∆x − 2 + e2ik ∆x
2 2∆x
(
)3
(
)
1 c∆t
−3ik′ ∆x
−ik′ ∆x
ik′ ∆x
3ik′ ∆x
n
−
− 3e
+ 3e
−e
Uk′ e
6 2∆x
{
}
(
)2
(
)3
c∆t
c∆t
i
c∆t
1
sin2 k′ ∆x −
sin3 k′ ∆x .
= Ukn′ 1 + i
sin k′ ∆x −
∆x
2 ∆x
6 ∆x
(37)
p=
c∆t
∆x
sin k′ ∆x とし, 増幅係数 λ =
Ukn+1
′
Ukn′
を導入すると,
p2
p3
λ = 1 + ip −
−i
2
6
√(
)
(
)2
2
p2
p3
|λ| =
1−
+ p−
2
6
√
1
1
= 1 − p4 + p6
12
36
(38)
√
√
√
√
− 3 ≤ p ≤ 3 のとき, |λ| ≤ 1 となる. よって, − 3 ≤ p ≤ 3 のとき, ホインの 3
次公式は一次元移流方程式に対して安定である.
4 段 4 位のルンゲクッタ型公式の安定性
4 段 4 位のルンゲクッタ型公式のうち, ここでは古典的ルンゲクッタ公式について
議論する. 時間差分における古典的ルンゲクッタ公式は以下の式で表される.
2015 0402-kawahara.tex
2015/04/02(川原 健史)
GFD ワークノート
ルンゲクッタ型公式の安定性
k1n = ∆tf (un , tn ),
1
1
k2n = ∆tf (un + k1n , tn + ∆t),
2
2
1
1
k3n = ∆tf (un + k2n , tn + ∆t),
2
2
n
n n
n
k4 = ∆tf (u + k3 , t + ∆t),
1
un+1 = un + (k1n + 2k2n + 2k3n + k4 ).
6
14
(39)
振動方程式
以下の様な振動方程式を考える.
du(t)
= iωu(t).
dt
(40)
1
1
1
un+1 = un + iω∆tun − (ω∆t)2 un − i(ω∆t)3 un + (ω∆t)4 un .
2
6
24
(41)
式 (40) に式 (39) を適用すると,
増幅係数 λ =
un+1
un
を導入すると,,
1
1
1
λ = 1 + iω∆t − (ω∆t)2 − i(ω∆t)3 + (ω∆t)4
2
6
24
1 4
1 2 1 3
= 1 + ip − p − ip + p .
2
6
24
(42)
ここで, p = ω∆t とした. ここで, |λ| の大きさを求める.
√
1
1
1
(1 − p2 + p4 )2 + (p − p3 )2
2
24
6
√
1
1 8
= 1 − p6 +
p
72
576
|λ| =
2015 0402-kawahara.tex
(43)
2015/04/02(川原 健史)
GFD ワークノート
ルンゲクッタ型公式の安定性
15
√
√
√
√
−2 2 ≤ p ≤ 2 2 のとき, |λ| ≤ 1 となる. よって, −2 2 ≤ p ≤ 2 2 のとき, 古典
的ルンゲクッタ公式は振動方程式に対して安定である.
一次元移流方程式
以下のような一次元移流方程式を考える.
∂u(x, t)
∂u(x, t)
+c
= 0,
∂t
∂x
c = const.
(44)
また, 領域は有限であり, 境界条件は周期境界とする. ここで, 空間差分には 2 次
の中心差分, 時間差分にはルンゲクッタ型公式を用いて差分化する. 古典的ルンゲ
クッタ公式においては, 以下のようになる.
unj+1 − unj−1
,
2∆x
n
n
unj+1 + 12 k1,j+1
− unj−1 − 12 k1,j−1
= −c∆t
2∆x
n
n
unj+1 + 12 k2,j+1
− unj−1 − 12 k2,j−1
= −c∆t
2∆x
n
n
unj+1 + k3,j+1
− unj−1 − k3,j−1
= −c∆t
2∆x
(
)
1
n
n
n
n
= unj +
k1,j
+ 2k2,j
+ 2k3,j
+ k4,j
.
6
n
k1,j
= −c∆t
n
k2,j
n
k3,j
n
k4,j
un+1
j
(45)
この式を整理すると, 以下のような式になる.
unj+1 − unj−1 1
un − 2unj + unj−2
2 j+2
− c∆t
=
+ (c∆t)
2∆x
2
(2∆x)2
n
n
n
n
uj+3 − 3uj+1 + 3uj−1 − uj−3
1
− (c∆t)3
6
(2∆x)3
n
u − 4unj+2 + 6unj − 4unj−2 + unj−4
1
4 j+4
+ (c∆t)
.
24
(2∆x)4
un+1
j
unj
ここで, u を複素数に拡張し, U を振幅, k を波数として unj =
∞
∑
(46)
Ukn e−ikj∆x を代
k=−∞
入すると,
2015 0402-kawahara.tex
2015/04/02(川原 健史)
GFD ワークノート
∞
∑
ルンゲクッタ型公式の安定性
16
Ukn+1 e−ikj∆x
k=−∞
∞
∑
=
Ukn e−ikj∆x
k=−∞
−
+
−
+
∞
)
c∆t ∑ n ( −ik(j+1)∆x
Uk e
− e−ik(j−1)∆x
2∆x k=−∞
(
)2 ∑
∞
(
)
1 c∆t
Ukn e−ik(j+2)∆x − 2e−ikj∆x + e−ik(j−2)∆x
2 2∆x k=−∞
(
)3 ∑
∞
(
)
1 c∆t
Ukn e−ik(j+3)∆x − 3e−ik(j+1)∆x + 3e−ik(j−1)∆x − e−ik(j−3)∆x
6 2∆x k=−∞
(
)4 ∑
∞
(
)
1 c∆t
Ukn e−ik(j+4)∆x − 4e−ik(j+2)∆x + 6e−ikj∆x − 4e−ik(j−2)∆x + e−ik(j−4)∆x .
24 2∆x k=−∞
(47)
′
両辺に eik j∆x をかけ, j について総和を取ると,
2015 0402-kawahara.tex
2015/04/02(川原 健史)
GFD ワークノート
J
∞
∑
∑
ルンゲクッタ型公式の安定性
17
′
Ukn+1 e−ikj∆x eik j∆x
j=1 k=−∞
=
J
∞
∑
∑
′
Ukn e−ikj∆x eik j∆x
j=1 k=−∞
J
∞
)
c∆t ∑ ∑ n ( −ik(j+1)∆x ik′ j∆x
−ik(j−1)∆x ik′ j∆x
−
U e
e
−e
e
2∆x j=1 k=−∞ k
(
)2 ∑
J
∞
(
)
∑
1 c∆t
′
′
′
Ukn e−ik(j+2)∆x eik j∆x − 2e−ikj∆x eik j∆x + e−ik(j−2)∆x eik j∆x
+
2 2∆x
j=1 k=−∞
(
)3 ∑
J
∞
∑
1 c∆t
−
Ukn (
6 2∆x
j=1 k=−∞
′
′
′
′
e−ik(j+3)∆x eik j∆x − 3e−ik(j+1)∆x eik j∆x + 3e−ik(j−1)∆x eik j∆x − e−ik(j−3)∆x eik j∆x
)
1
+
24
e
(
c∆t
2∆x
)4 ∑
∞
J
∑
Ukn (
j=1 k=−∞
−ik(j+4)∆x ik′ j∆x
e
′
′
− 4e−ik(j+2)∆x eik j∆x + 6e−ikj∆x eik j∆x
′
′
− 4e−ik(j−2)∆x eik j∆x + e−ik(j−4)∆x eik j∆x
),
2015 0402-kawahara.tex
2015/04/02(川原 健史)
GFD ワークノート
∞
∑
Ukn+1
J
∑
18
′
ei(k −k)j∆x
j=1
k=−∞
∞
∑
=
ルンゲクッタ型公式の安定性
Ukn
J
∑
′
ei(k −k)j∆x
j=1
k=−∞
J (
)
∑
c∆t
n
i(k′ j−k(j+1))∆x
i(k′ j−k(j−1))∆x
−
U
e
−e
2∆x k=−∞ k j=1
(
)2 ∑
∞
J (
)
∑
1 c∆t
′
′
′
+
Ukn
ei(k j−k(j+2))∆x − 2ei(k j−k)j∆x + ei(k j−k(j−2))∆x
2 2∆x k=−∞
j=1
(
)3 ∑
∞
J
∑
1 c∆t
n
−
U
(
6 2∆x k=−∞ k j=1
∞
∑
′
′
′
′
′
′
′
ei(k j−k(j+3))∆x − 3ei(k j−k)(j+1)∆x + 3ei(k j−k(j−1))∆x − ei(k j−k(j−3))∆x
(
)4 ∑
J
∞
∑
1 c∆t
n
(
U
+
24 2∆x k=−∞ k j=1
′
)
′
ei(k j−k(j+4))∆x − 4ei(k j−k)(j+2)∆x + 6ei(k j−kj)∆x − 4ei(k j−k(j−2))∆x + ei(k j−k(j−4))∆x
)
∞
∑
=
k=−∞
−
+
−
+
Ukn
J
∑
′
ei(k −k)j∆x
j=1
J
( −ik∆x
)∑
c∆t
′
n
ik∆x
ei(k −k)j∆x
Uk e
−e
2∆x k=−∞
j=1
(
)2 ∑
J
∞
( −2ik∆x
)∑
1 c∆t
′
n
2ik∆x
Uk e
−2+e
ei(k −k)j∆x
2 2∆x k=−∞
j=1
(
)3 ∑
J
∞
( −3ik∆x
)∑
1 c∆t
′
n
−ik∆x
ik∆x
ik∆x
Uk e
− 3e
+ 3e
−e
ei(k −k)j∆x
6 2∆x k=−∞
j=1
(
)4 ∑
∞
J
(
)∑
1 c∆t
′
Ukn e−4ik∆x − 4e−2ik∆x + 6 − 4e2ik∆x + e4ik∆x
ei(k −k)j∆x .
24 2∆x k=−∞
j=1
∞
∑
(48)
三角関数の直交性より,
2015 0402-kawahara.tex
2015/04/02(川原 健史)
GFD ワークノート
ルンゲクッタ型公式の安定性
19
)
c∆t n ( −ik′ ∆x
Ukn+1
= Ukn′ −
Uk′ e
− eik∆x
′
2∆x
(
)2
(
)
1 c∆t
′
′
+
Ukn′ e−2ik ∆x − 2 + e2ik ∆x
2 2∆x
(
)3
(
)
1 c∆t
−3ik′ ∆x
−ik′ ∆x
ik′ ∆x
3ik′ ∆x
n
−
− 3e
+ 3e
−e
Uk′ e
6 2∆x
(
)4
(
)
1 c∆t
′
′
′
′
+
Ukn′ e−4ik ∆x − 4e−2ik ∆x + 6 − 4e2ik ∆x + e4ik ∆x
24 2∆x
= Ukn′ {
(
)2
c∆t
1 c∆t
′
1+i
sin k ∆x −
sin2 k′ ∆x
∆x
2 ∆x
(
)3
(
)4
i c∆t
1 c∆t
3 ′
−
sin k ∆x +
sin4 k′ ∆x
6 ∆x
24 ∆x
}.
p=
c∆t
∆x
(49)
sin k′ ∆x とし, 増幅係数 λ =
Ukn+1
′
Ukn′
を導入すると,
p3 p4
p2
−i +
λ = 1 + ip −
2
6
24
√(
)
(
)2
2
p2 p4
p3
|λ| =
1−
+
+ p−
2
24
6
√
1
1 8
= 1 − p6 +
p
72
576
(50)
√
√
√
√
−2 2 ≤ p ≤ 2 2 のとき, |λ| ≤ 1 となる. よって, −2 2 ≤ p ≤ 2 2 のとき, 古典
的ルンゲクッタ公式は一次元移流方程式に対して安定である.
2015 0402-kawahara.tex
2015/04/02(川原 健史)
© Copyright 2024 ExpyDoc
