■2 階常微分方程式  ′′  y (x) = −y(x) , x ∈ (0, L] y(0) = 1  ′ y (0) = 0 (1) z(x) = y ′ (x) とおくと、 z ′ (x) = y ′′ (x) = −y(x) y(x) = t(y(x), z(x)) とすると、(1) は、 y ′ (x) = ( ( = ( y ′ (x) z ′ (x) z(x) −y(x) ) ) )( ) 0 1 y(x) −1 0 z(x) ( ( )) 0 1 = Ay(x) A = −1 0 = と書き換えることができる。よって、未知関数 1 つの 2 階常微分方程式は、未知関数を成分に持つ 2 次ベクトルの 1 階常微分方程式に帰着された。     y ′ (x) =    y(0) = ( Ay(x) , x ∈ (0, L] ( ) 1 0 ( A= 0 1 −1 0 )) (2) (2) の数値計算は、これまでやった未知関数 1 つの 1 階常微分方程式と同様にして出来る。 ■オイラー法 [0, L] を均等な幅 h で N 分割する。すなわち h = L/N 。よって、離散点 xi (i = 0, · · · , N ) は、xi = ih と表せる。y(xi )(= (y(xi ), z(xi )) の数値解を yi (= (yi , zi )) とする。(2) をオイラー法により離散化する。 yi+1 − yi = Ayi h yi+1 = yi + hAyi = (E + hA)yi これを各成分について書き出すと、 { yi+1 zi+1 = yi + hzi = zi − hyi 初期条件は、y0 = 1, z0 = 0 である。よって、これを解けばよい。 1 ■リープ・フロッグ法 [0, L] を均等な幅 h で 2N 分割する。すなわち h = L/2N 。よって、離散点 xi (i = 0, · · · , N ) は、xi = ih と表せる。y(xi ) の数値解を yi とする。中心差分法 y ′ (xi ) を次のように近似することを中心差分法という。 y ′ (xi ) ≃ y(xi+1 ) − y(xi−1 ) 2h 中心差分法を用いると (2) の離散式は次のようになる。 yi+1 − yi−1 = Ayi 2h yi+1 = yi−1 + 2hAyi (3) これを各成分について書き出すと、 { yi+1 zi+1 = yi−1 + 2hzi = zi−1 − 2hyi (4) (3) において i = 0 とすると、 y1 = y−1 + 2hAy0 である。しかし、y−1 が分からないので (3) により y1 を与えるこはできない。そこで、オイラー法により、 y1 を与える。 y1 = (E + hA)y0 これを各成分について書き出すと、 { y1 z1 = y0 + hz0 = z0 − hy0 (4) を次のようにして解く方法をリープ・フロッグ法という。 y0, z1 を初期値にして解く場合 { y2k z2k+1 = y2k−2 + 2hz2k−1 = z2k−1 − 2hy2k (k = 1, · · · , N ) (k = 1, · · · , N − 1) (5) y1, z0 を初期値にして解く場合 { z2k y2k+1 = z2k−2 − 2hy2k−1 = y2k−1 + 2hz2k (k = 1, · · · , N ) (k = 1, · · · , N − 1) (6) (5) では、y は偶数番目、z は奇数番目の離散点での値を使って、y は次の偶数番目、z は次の奇数番目の離散点での値を求めている。 (6) では、y は奇数番目、z は偶数番目の離散点での値を使って、y は次の奇数番目、z は次の偶数番目の離散点での値を求めている。 2 ■数値計算例 L = 10, N = 50 として、オイラー法とリープ・フロッグ法により求めた数値解を比較する。厳密解との比較それぞれの方法で求めた数値解を厳密解と比較する。 2 1.5 1 0.5 0 -0.5 -1 -1.5 -2 -2.5 -3 0 2 4 6 8 10 cos(x) eular method leap-frog method 3 位相図数値解 (yi , zi ) を図示する。厳密解は y(x) = cos (x) である。つまり、(y(x), z(x)) = (cos (x), − sin (x)) であるから、これは原点を中心にもつ半径 1 の円を表す。それぞれの方法で求めた数値解 (yi , zi ) をプロットして、解の軌道が半径 1 の円周上にあるか調べる。 3 2 1 0 -1 -2 -3 -3 -2 eular method -1 0 1 2 leap-frog method 4 3 エネルギーの保存 (1) を積分すると次のようになる。 ∫ y ′′ (x) + y(x) = 0 ∫ ) d ( ′ 2 2 (y (x)) + y(x) = 0 dx 1 1 ′ 2 (y (x)) + y(x)2 = C (C は定数) 2 2 (7) つまり、(7) の左辺の値は任意の x で一定である。これは、(7) 左辺の第 1 項は運動エネルギー、第 2 項は弾性エネルギーであり、力学的エネルギーの保存側により、これらのエネルギーの和は場所によらず一定であることを意味している。よって、それぞれの方法で求めた数値解より、各場所（離散点）でのエネルギーの和を計算して、それが保存しているか比較する。初期条件より、x = 0 でのエネルギーの和は、1/2 であることが分かるので、すべての離散点 xi でこれが保たれているか調べる。 4 3.5 3 2.5 2 1.5 1 0.5 0 0 2 eular method 4 6 8 10 leap-frog method 5