● 講義資料

2013年度・後期・数理解析・計算機数学３・第１１回 1

● 講義資料

▼ 講義予定

• 連立一次方程式の数値解法

● 前回の講義のまとめ

★ Symplectic法

• Runge-Kutta法は,広範囲の常微分方程式の初期値問題に適用可能な汎用的な方法であるが,

力学に由来するエネルギー保存則が成り立つ系では, Runge-Kutta法による数値解はエネル ギーを保存しない.

たとえば,x^′′=−xの解に沿って,エネルギー E(t) =1

2((x(t))²+ (x^′(t))²)

は一定値を持つ. この証明はE^′(t) = 0を示す方法もあるが,方程式の両辺にx^′ をかけると 0 =x^′′(t)x^′(t) +x(t)x^′(t) =1

2 d

dt((x^′(t))²+ (x(t))²) となることからも証明できる.

しかし,前進オイラー法ではEⁿ = (1/2)(x²n+vn²)は単調増加となり,古典的Runge-Kutta 法ではEⁿ は単調減少となってしまい,そのような数値解は「よい」数値解と言うことは難 しい. そこで,エネルギーまたはそれに近いものが保存するような数値解を構成することが 望まれる.

• ここでは, Hamilton系と呼ばれる常微分方程式を考える.

未知関数を(q1, . . . , qN, p1, . . . , pn), qi: [0, T]−→R, pi: [0, T]−→R として, Hamiltonian H:R^2N −→Rが与えられているととき,常微分方程式

dq dt =∂H

∂p, dp

dt =−∂H

∂q を考える.

x^′′=−xに対応する問題としては,q=x,p=x^′ と考え,H(q, p) = (1/2)(q²+p²)とすれば, dq

dt = ∂H

∂p =p, dp

dt =−∂H

∂q =−q

より d

dt dq dt = d

dtp=−q

となるので,x^′′=−xを考えることと,H(q, p) = (1/2)(q²+p²)に対するHamilton系を考 えることは同じである.

Dec. 18, 2013, Version: 1.0 [email protected]

2013年度・後期・数理解析・計算機数学３・第１１回 2

また, HamlitonianH は解に沿って一定値をとる. すなわち, dH

dt = 0 が成り立つ. （実際に計算してみればよい.）

なお,系の運動エネルギーが(1/2)P

q²i で書けているとき, HamiltonianH は,系の全エネ ルギーを与える関数に一致する.

• Hamilton系の持つ重要な性質は, Hamiltonianが保存することのほかに,次のLiouville の 定理が成り立つことである.

φ^t:R^2N −→R^2N を初期条件(q(0), p(0))に対する解(q(t), p(t))によるt 秒後の位置とす る. （この写像を相流,または1-parameter 変換群と呼ぶ.）この時,φ^t は R^2N の任意の 領域の面積を保存する. すなわち,φ^tの微分Dφ^tに対して,

|det(Dφ^t)|= 1 が成り立つ. （φtは面積保存写像となる.）

• このことから, Hamilton 系に対する数値解法として, (qn, pn)7→(qn+1, pn+1)が面積保存写 像となるような解法が望ましいことがわかる. そのような解法をSymplectic 解法と呼ぶ.

• HamiltonianH が

H(q, p) =T(p) +U(q) と分解できている場合に,以下の計算スキーム

qn+1 =qn+hdT dp(pn), pn+1=pⁿ−hdU

dq(qn+1),

をSymplectic Euler 法とよび, 1次のSymplectic解法となっていることが証明できる.

• Symplectic Euler法を拡張した方法のパラメータ St¨olmer-Verlet (p= 2)

1 2

bi 0 1

cⁱ 1/2 1/2 Ruth (p= 3)

1 2 3

bⁱ 7/24 3/4 −1/24 cⁱ 2/3 −2/3 −1

Sanz-Serna (p= 4)

1 2 3 4 5 6

bⁱ 7/48 3/8 −1/48 −1/48 3/8 7/48

cⁱ 1/3 −1/3 1 −1/3 1/3 0

Dec. 18, 2013, Version: 1.0 [email protected]

2013年度・後期・数理解析・計算機数学３・第１１回 3

• 陰的Runge-Kutta法で, Symplectic解法となるものの例

陰的中点法（１段２次Gauss-Legrandre法）

1/2 1/2 1

２段４次Gauss-Legrandre法 3−√

3 6

1 4

3−2√ 3 3 +√ 12

3 6

3 + 2√ 3 12

1 4

1/2 1/2

３段６次Gauss-Legrandre法 5−√

15 10

5 36

80−24√ 15 360

50−12√ 15 1 360

2

50 + 15√ 15 360

2 9

50−15√ 15 5 +√ 360

15 10

50 + 12√ 15 360

80 + 24√ 15 360

5 36

5/18 4/9 5/18

● 実習内容

1. （★★★）4×4 行列の積を計算するプログラムを書きなさい.

2. （★★★）4×4 行列の転置行列を計算するプログラムを書きなさい.

3. （★★★）連立一次方程式











3x+ 12y+ 9z= 3 2x+ 5 y+ 4z= 4

−x+ 3 y+ 2z= 5

をGauss-Jordanの消去法（掃出し法）で解くプログラムを書きなさい. また,この方程式を

Gaussの消去法で解くプログラムを書きなさい.

4. （★★）整数係数の4×4行列に対して,それが正則か否かを判定し,もし正則であるときに は,逆行列を有理数係数で求め,正則でないときには, ランク,核の基底,像の基底を求めるプ ログラムを書きなさい,

Dec. 18, 2013, Version: 1.0 [email protected]