2012年度・前期・数理解析・計算機数学２・第１１回 1
● 講義資料
▼ 講義予定
• Symplectic解法
• 常微分方程式の境界値問題の数値解法
• 連立一次方程式の数値解法
● 前回の講義のまとめ
★ 数値積分
• 台形公式・中点公式などの誤差評価としては,以下の式を得ることができる.
Z 1 0
f(x)dx= 1 N
N−1
X
k=0
f(ak) +f′(ξ1) 2N Z 1
0
f(x)dx= 1 N
N−1
X
k=0
f(ak) +f(ak+1)
2 +f′′(ξ2) 12N2 Z 1
0
f(x)dx= 1 N
N−1
X
k=0
f(ak+1/2)−f′′(ξ3) 24N2
• ニュートン・コーツの公式を使って,周期関数をその周期で積分したときには,上の計算誤差 よりも精度の高い計算ができている. その説明はオイラー・マクローリンの和公式
n
X
i=0
f(i) = Z n
0
f(x)dx+1
2(f(0) +f(n)) +
k
X
j=2
Bj
j!(f(j−1)(n)−f(j−1)(0)) +Rk, Rk =(−1)k−1
k!
Z n 0
B˜k(x)f(k)(x)dx
から得ることができる. ただし, ˜Bk(x)は Bk(x)を周期1で拡張した関数であり,Bn は x
ex−1 =X
n=0
Bn
n!xn, Bn(x) =
n
X
j=0
n j
Bjxn−j
で定義される. (Bernoulli数) 特に, Bernoulli 数に関しては, B0(x) = B0= 1, B1=−1/2, B2k+1= 0,k≥1,Bk′(x) =kBk−1(x)が成り立つ.
• 特に,B2k+1= 0となることから,区間の分割幅hに対して,台形公式は,h−2,シンプソンの 公式はh−4 の誤差を与えることが証明できる.
Jun. 27, 2012, Version: 1.0 [email protected]
2012年度・前期・数理解析・計算機数学２・第１１回 2
★ 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 を考える.
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 は,系の全エネ ルギーを与える関数に一致する.
Jun. 27, 2012, Version: 1.0 [email protected]
2012年度・前期・数理解析・計算機数学２・第１１回 3
• 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 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
Jun. 27, 2012, Version: 1.0 [email protected]
2012年度・前期・数理解析・計算機数学２・第１１回 4
２段４次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の消去法（掃出し法）で解くプログラムを書きなさい.
4. （★★）整数係数の4×4行列に対して,それが正則か否かを判定し,もし正則であるときに は,逆行列を有理数係数で求め,正則でないときには, ランク,核の基底,像の基底を求めるプ ログラムを書きなさい,
Jun. 27, 2012, Version: 1.0 [email protected]
