2015年度・前期・数理解析・計算機数学３・第１１回 1
● 講義資料
▼ 講義予定
• Symplectic解法
• 連立一次方程式の数値解法
● 前回の講義のまとめ
★ 数値積分
• Runge-Kutta型公式をx′(t) =f(t),x(0) = 0に適用して,x(1)の近似値を求めることは,積 分I=
Z 1 0
f(t)dtの近似値を求めることに他ならない.
• Runge-Kutta型公式をx′(t) =f(t),x(0) = 0に適用すると, xn+1=xn+h
s
X
i=1
biki, ki=f(t+hci), ci=X
j
aij
と書ける. 特に,前進オイラー法,ホインの２次公式, Runge-Kutta法,陰的中点法を使うと, ti=ih, 1 =N h とおいて,それぞれ,
Z 1 0
f(t)dt∼h
N−1
X
i=0
f(ti), Z 1
0
f(t)dt∼h 2
N−1
X
i=0
(f(ti) +f(ti+1)) （複合）台形公式 Z 1
0
f(t)dt∼h 6
N−1
X
i=0
(f(ti) + 4f(ti+1/2) +f(ti+1)) （複合）シンプソンの公式 Z 1
0
f(t)dt∼h
N−1
X
i=0
f(ti+1/2) （複合）中点則 を得る.
• 台形公式,シンプソンの公式,中点則の誤差は,以下のように評価できる Z 1
0
f(t)dt∼h 2
N−1
X
i=0
(f(ti) +f(ti+1))∼h2(f′(1)−f′(0)) Z 1
0
f(t)dt−h 6
N−1
X
i=0
(f(ti) + 4f(ti+1/2) +f(ti+1))∼h4(f(3)(1)−f(3)(0)) Z 1
0
f(t)dt−h
N−1
X
i=0
f(ti+1/2)∼h2(f′(1)−f′(0))
となる. すなわち, 台形公式と中点則は, １次以下の多項式に対して正しい積分値を与え, シンプソンの公式は, ３次以下の多項式に対して正しい積分値を与える. この誤差評価は,
Euler-MacLaurinの総和公式を用いて証明することができる.
Jun. 29, 2015, Version: 1.0 [email protected]
2015年度・前期・数理解析・計算機数学３・第１１回 2
• 一般に, [0,1]上の相異なる N+ 1点{ti}Ni=0 を与えた時, (ti, f(ti))を通るN 次多項式fN
がただ一つ存在する. このfN は,f を, サンプル点{ti} を用いてLegrandre補間した多項 式である. このとき,
Z 1 0
f(t)dt= Z 1
0
fN(t)dt=
N
X
i=0
WifN(ti), Wi= Z 1
0
Q
j6=i(t−tj) Q
j6=i(ti−tj)dt
が成り立つ. よって,一般にはN+ 1点のサンプル点を使った数値積分公式は,N 次以下の 多項式に対して正しい積分値を与える.
• N+ 1次 Legendre多項式PN+1 の零点を{τi}Ni=0 とすると,
N
X
i=0
WifN(τi), Wi= Z 1
0
Q
j6=i(t−τj) Q
j6=i(τi−τj)dt は, 2N+ 1次以下の多項式f に対して正しい積分値
Z 1 0
f(t)dtを与える.（ガウスの積分法）
● 今回の講義資料
★ Symplectic法
• Runge-Kutta法は,広範囲の常微分方程式の初期値問題に適用可能な汎用的な方法であるが,
力学に由来するエネルギー保存則が成り立つ系では, Runge-Kutta法による数値解はエネル ギーを保存しない.
たとえば,x′′=−xの解に沿って,エネルギー E(t) =1
2((x(t))2+ (x′(t))2)
は一定値を持つ. この証明はE′(t) = 0を示す方法もあるが,方程式の両辺にx′ をかけると 0 =x′′(t)x′(t) +x(t)x′(t) =1
2 d
dt((x′(t))2+ (x(t))2) となることからも証明できる.
しかし,前進オイラー法ではEn = (1/2)(x2n+vn2)は単調増加となり,古典的Runge-Kutta 法ではEn は単調減少となってしまい,そのような数値解は「よい」数値解と言うことは難 しい. そこで,エネルギーまたはそれに近いものが保存するような数値解を構成することが 望まれる.
• ここでは, Hamilton系と呼ばれる常微分方程式を考える.
未知関数を(q1, . . . , qN, p1, . . . , pn), qi: [0, T]−→R, pi: [0, T]−→R として, Hamiltonian H:R2N −→Rが与えられているととき,常微分方程式
dq dt =∂H
∂p, dp
dt =−∂H
∂q を考える.
Jun. 29, 2015, Version: 1.0 [email protected]
2015年度・前期・数理解析・計算機数学３・第１１回 3
x′′=−xに対応する問題としては,q=x,p=x′ と考え,H(q, p) = (1/2)(q2+p2)とすれば, dq
dt = ∂H
∂p =p, dp
dt =−∂H
∂q =−q
より d
dt dq dt = d
dtp=−q
となるので,x′′=−xを考えることと,H(q, p) = (1/2)(q2+p2)に対するHamilton系を考 えることは同じである.
また, HamlitonianH は解に沿って一定値をとる. すなわち, dH
dt = 0 が成り立つ. （実際に計算してみればよい.）
なお,系の運動エネルギーが(1/2)P
q2i で書けているとき, HamiltonianH は,系の全エネ ルギーを与える関数に一致する.
• Hamilton系の持つ重要な性質は, Hamiltonianが保存することのほかに,次のLiouville の 定理が成り立つことである.
φt:R2N −→R2N を初期条件(q(0), p(0))に対する解(q(t), p(t))によるt 秒後の位置とす る. （この写像を相流,または1-parameter 変換群と呼ぶ.）この時,φt は R2N の任意の 領域の面積を保存する. すなわち,φ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=pn−hdU
dq(qn+1),
をSymplectic Euler 法とよび, 1次のSymplectic解法となっていることが証明できる.
• Symplectic Euler法を拡張した方法のパラメータ
St¨olmer-Verlet (p= 2)
1 2
bi 0 1
ci 1/2 1/2
Jun. 29, 2015, Version: 1.0 [email protected]
2015年度・前期・数理解析・計算機数学３・第１１回 4
Ruth (p= 3)
1 2 3
bi 7/24 3/4 −1/24 ci 2/3 −2/3 −1
Sanz-Serna (p= 4)
1 2 3 4 5 6
bi 7/48 3/8 −1/48 −1/48 3/8 7/48
ci 1/3 −1/3 1 −1/3 1/3 0
• 陰的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. （★★★）６回目資料の「単振動」と「単振り子」の常微分方程式の初期値問題の数値解を, Symplectic Euler法で計算しなさい.
2. （★★★）６回目資料の「単振動」と「単振り子」の常微分方程式の初期値問題の数値解を, 陰的中点法で計算しなさい.
3. （★★★）4×4 行列の積を計算するプログラムを書きなさい.
4. （★★★）4×4 行列の転置行列を計算するプログラムを書きなさい.
5. （★★★）連立一次方程式
3x+ 12y+ 9z= 3 2x+ 5 y+ 4z= 4
−x+ 3 y+ 2z= 5
をGauss-Jordanの消去法（掃出し法）で解くプログラムを書きなさい. また,この方程式を
Gaussの消去法で解くプログラムを書きなさい.
Jun. 29, 2015, Version: 1.0 [email protected]
