• /Users/bob/Desktop/maple_ode/python_ode.ipynb
• origin git@github.com:daddygongon/maple_ode.git (fetch)

Table of Contents

Euler法による落下運動¶

重力場中の運動¶

重力場中のボールの落下を考えて，１軸で考えた運動方程式を立てます．

$$ \begin{aligned} v &= \frac{dx}{dt} \\ a &= \frac{d^2 x}{dt^2} \end{aligned} $$

質量を$m$, 重力加速度を$g$として，働く力が$F=-mg$であるとすると，ニュートンの運動方程式$F=ma$は，

$$ -mg = m \frac{d^2 x}{dt^2} $$

となります．

Euler法¶

1次の微分方程式の一般形は

$$ \frac{dx}{dt}=f(x,t) $$

と書けます．この微分方程式を簡単な近似から求めるオイラー法を示します． $x(t+\delta t)$をテイラー級数展開すると，

$$ x(t+\delta t) \simeq x(t) + \frac{dx}{dt} \delta t $$

となります．これらを代入すると，計算アルゴリズムはつぎのようになります，

$$ x_{i+1} = x_i + f_i\,\delta t $$

ここで，$f_i$は点($x_i, t_i$)における関数の値です．このアルゴリズムを適用して，$t+\delta t$の座標$x_{i+1}$を一つ前の時間の座標$x_i$から導くことができます．これを重力場中の運動方程式に適用します．

重力場中の運動をEuler法で解いたら¶

Euler法は一階の微分方程式に対する定式化をしています．ところが，重力場中の運動は2階の微分方程式です．このようなときには媒介変数を導入して1次連立方程式に置き直します．

媒介変数として速度$v$を使って，2階の運動方程式 $$ \frac{d^2x}{dt^2} = -g $$ が，1階の連立方程式 $$ \begin{aligned} \frac{dv}{dt} & = -g \\ \frac{dx}{dt} & = v \end{aligned} $$ で置き換えられると考えることに相当します．アルゴリズムにすると，

$$ \begin{aligned} v_{i+1} &= v_i - g\, dt \\ x_{i+1} &= x_i + v_i\, dt \end{aligned} $$

なる連立方程式を解くことに置き換わります．これをpythonで関数にして，さらに計算結果を表示させてみます．

Eulerは$x_i$, $v_i$を受け取って，先ほど導いた簡単な計算によって，$v_{i+1}$, $x_{i+1}$を順次計算して返します．結果は，

位置(height:$x$)と速度($v$)の時間($t$)変化を表示させています．

時間とともに位置は放物線状に変化し，速度は一定の傾きで増加していく，等加速度運動を再現しています．Euler法ではこのように非常に簡単なcodeによって微分方程式で表される現象をシミュレートできることがわかるでしょう．

空気抵抗がある水滴の落下¶

ballの落下ではわかりにくいですが，より小さな質量の水滴では，速度に比例する空気抵抗が効いてきます．この様子を見ましょう．微分方程式には，

$$ F_x = - C v_x $$

項が付与されます．そうすると運動方程式は

$$ m \frac{dv_x}{dt} = - C v_x - mg $$

となります．これにともなったv_xの時間変化に対して，今までは単純に重力加速していたのが，v_xに比例する空気抵抗を記述する項が付与されます．この変化をEuler2に入れ込むと少しの修正ですが，結果は劇的に変化します．

その前に，ここで示した修正がどれほど普通の微分方程式の解法としてややこしいかを少し示しておきます． $$ m \frac{dv_x}{dt} = - C v_x - mg $$ を整理すると， $$ \frac{dv_x}{dt} = - \frac{C}{m} v_x - g $$ です．右辺を$v_x$の係数でくくって，変形していくと， $$ \begin{aligned} \frac{dv_x}{dt} & = - \frac{C}{m}\left(v_x - \frac{m}{C}g\right) \\ \frac{1}{\left(v_x - \frac{m}{C}g\right)}dv_x & = - \frac{C}{m}dt \end{aligned} $$ これで変数分離されているので積分をとって， 公式$\int \frac{1}{x+a} dx = \log |x+a| + C_1$を使うと $$ \log \left|v_x - \frac{m}{C}g \right| = - \frac{C}{m} t + C_1 $$ なんですが，両辺の指数をとって「適当」に変形すると $$ \begin{aligned} \left| v_x - \frac{m}{C}g\right| & = \exp(- \frac{C}{m}t + C_1 ) \\ v_x &= \frac{gm}{C} + \exp(- \frac{C}{m}t + C_1 ) \end{aligned} $$ となります．最後の絶対値の変形は怪しいんであまり信頼しないでね．

これをpythonが用意している数式処理libraryのsympyを使うと次の通りとなります．

ちょっと気をつけたいのは，$-Cv$の向きです．vは例えば自由落下の時には下向きのベクトルなんで，抵抗は反対向き，すなわち上むき，数値的には正，になります．

高精度計算¶

バネの運動¶

今度はバネの運動です．空気抵抗との違いはほんの少しで， $$ F_x = -k x $$ と今度は，位置xに力が比例することです．そうすると運動方程式は，

$$ m \frac{dv_x}{dt} = - k x $$

となります．

連立方程式は $$ \begin{aligned} \frac{dv}{dt} & = - \frac{k}{m}x \\ \frac{dx}{dt} & = v \end{aligned} $$ となるんで，アルゴリズムに置き換えると，

$$ \begin{aligned} v_{i+1} & = v_i - \frac{k}{m} x_i\, dt \\ x_{i+1} & = x_i + v_i\, dt \end{aligned} $$

なる連立方程式を解くことに置き換わります．これをpythonで関数にして，さらに計算結果を表示させてみます．

kはあらかじめ$m$で割られて正規化されているとします．これをEuler法で計算すると上のような結果が得られます．

徐々に発散していく様子がわかると思います．本来，摩擦のないバネは定常的に振動します．この発散の原因は，Euler法の計算誤差が大きいせいです．そこで，より精度の高いRunge-Kutta法を導入します．

2次のRunge-Kuttaの導出¶

一般にRunge-Kutta法と呼ばれる手法は，4次の古典的Runge-Kutta法を指します．導出は意外と面倒なので，2次の場合のさわりを紹介して，そこからの類推としましょう．$x,v,t$などのパラメータ名を変更しますのでご注意あれ．

(「ANSI Cによる数値計算法」堀之内聰一，酒井幸吉，榎園茂，森北2002, p.133) テイラー展開により，$h^2$の精度まで展開する． $$ y(x_0+h) = y(x_0) +y_0'h+\frac{1}{2!}y_0''h^2 + O(h^3) $$

この式において，$y_0' = f(x_0,y_0)$は既知とする．一方，$y_0''$は$f(x_0,y_0)$から直接的には求められない．したがって，この式の右辺と$h^2$の項まで一致する近似値を，$f(x_0,y_0)$だけを既知として算出する方法を考えよう．

平均値の定理より， $$ \Delta y = y(x_0+h)-y(x_0) = hy'(x_0 + \theta h), \, 0<\theta<1 $$

$y'(x_0 +\theta h)$の近似値として，$\theta=0, \theta=1$の場合を考えると， $$ \begin{aligned} \Delta y &\simeq hy'(x_0) \, &{\textrm where}\, &\theta=0\\ \Delta y &\simeq hy'(x_0+h) \, &{\textrm where}\, &\theta=1 \end{aligned} $$ これらの値は単独では$\Delta y$に対して$h^2$の精度をもつ近似値にならないが，これらの一次結合$\alpha h y'(x_0)+\beta hy'(x_0+h)$を$\alpha, \beta$をうまく定めることによって，$\Delta y$の$h^2$の精度をもつ近似値にすることができる．実際， $$ \begin{aligned} \alpha h y'(x_0)+\beta hy'(x_0+h) & =\alpha h y'(x_0)+\beta h \{y'(x+0)+y''(x_0)h+O(h^2)\} \\ & =(\alpha+\beta)hy'_0 + \beta h^2 y_0'' + O(h^3) \end{aligned} $$

したがって，テイラー展開式と係数を比較して， $$ \begin{aligned} \alpha + \beta = 1, \, \beta = \frac{1}{2}\\ \alpha = \frac{1}{2}, \, \beta =\frac{1}{2} \end{aligned}$$ となり，

$$ \Delta y = \frac{1}{2}hy'(x_0) + \frac{1}{2}hy'(x_0+h)+O(h^3) $$

いま， $$ k_1 =hy'(x_0) =hf(x_0,y_0) $$

とおこう．上式に代入して， $$ \begin{aligned} hy'(x_0+h) &= hf(x_0+h,y(x_0+h)) \\ & =hf(x_0+h, y(x_0)+y'(x_o)h+O(h^2)) \\ & =hf(x_0+h, y_0+k_1+O(h^2)) \\ & =hf(x_0+h, y_0+k_1)+O(h^2) \end{aligned}$$

したがって， $$ k_2 = hf(x_0+h, y_0+k_1) $$ とおけば， $$ \Delta y = \frac{1}{2}k_1 + \frac{1}{2}k_2 + O(h^3) $$ となる．これより， $$ k = \frac{1}{2}(k_1+k_2), y_1 = y_0 +k $$

とおくと，$y(x_0+h) = y_1+O(h^3)$となり，$y_1$は$h^2$の精度の近似値となる．

Runge-Kutta2次公式¶

こうして得られたRunge-Kuttaの２次の公式を定義すると次の通りです．

微分方程式 $$ \frac{dy}{dx} = f(x,y), \, where \, y(x_0)=y_0 $$

の数値解は，刻み幅を$h$，$x_n=x_0+nh$として，次の漸化式 $$ y_{n+1} = y_n +k (n=0,1,2,\cdots) $$ で与えられる．ここに，$k$は次で定める． $$ \begin{aligned} k_1 & = hf(x_n,y_n), \\ k_2 & = hf(x_n+h, y_n+k_1), \\ k & = \frac{1}{2}(k_1+k_2) \end{aligned} $$

Runge-Kutta4次公式¶

2次と同様の考え方で，$h^4$の精度を持つRunge-Kutta４次公式を作ることができる．

微分方程式 $$ \frac{dy}{dx} = f(x,y), \, {\textrm where}y(x_0)=y_0 $$

の数値解は，刻み幅を$h$，$x_n=x_0+nh$として，次の漸化式 $$ y_{n+1} = y_n +k (n=0,1,2,\cdots) $$ で与えられる．ここに，$k$は次で定める． $$ \begin{aligned} k_1 & = hf(x_n,y_n), \\ k_2 & = hf(x_n+\frac{h}{2}, y_n+\frac{k_1}{2}), \\ k_3 & = hf(x_n+\frac{h}{2}, y_n+\frac{k_2}{2}), \\ k_4 & = hf(x_n+h, y_n+k_3), \\ k & = \frac{1}{6}(k_1+2k_2+2k_3+k_4) \end{aligned} $$

連立方程式にRunge-Kutta4次公式を¶

連立微分方程式 $$ \begin{aligned} \frac{dy}{dx} &= f(x,y,z), &\, \, where \, y(x_0)=y_0 \\ \frac{dz}{dx} &= g(x,y,z), &\, \, where \, z(x_0)=z_0 \\ \end{aligned} $$

の数値解は，刻み幅を$h$，$x_n=x_0+nh$として，次の漸化式 $$ \begin{aligned} y_{n+1} & = y_n +k &\, (n=0,1,2,\cdots) \\ z_{n+1} & = z_n +l &\, (n=0,1,2,\cdots) \\ \end{aligned} $$ で与えられる．ここに，$k,l$は次で定める． $$ \begin{aligned} k_1 &= hf(x_n,y_n,z_n), \, &l_1 &= hg(x_n,y_n,z_n), \\ k_2 &= hf(x_n+\frac{h}{2}, y_n+\frac{k_1}{2}, z_n+\frac{l_1}{2}), \, &l_2 &= hg(x_n+\frac{h}{2}, y_n+\frac{k_1}{2}, z_n+\frac{l_1}{2}), \\ k_3 &= hf(x_n+\frac{h}{2}, y_n+\frac{k_2}{2}, z_n+\frac{l_2}{2}), \, &l_3 &= hg(x_n+\frac{h}{2}, y_n+\frac{k_2}{2}, z_n+\frac{l_2}{2}), \\ k_4 &= hf(x_n+h, y_n+k_3, z_n+l_3), \, &l_4 &= hg(x_n+h, y_n+k_3, z_n+l_3), \\ k &= \frac{1}{6}(k_1+2k_2+2k_3+k_4), \, &l &= \frac{1}{6}(l_1+2l_2+2l_3+l_4) \end{aligned} $$

Runge-Kuttaの４次公式をそのままcodingすると上のようになります．これを先ほどのバネ運動の問題に当てはめてみましょう．

先ほど導出した運動方程式の漸化式 $$ \begin{aligned} v_{i+1} & = v_i - \frac{k}{m} x_i\, dt \\ x_{i+1} & = x_i + v_i\, dt \end{aligned} $$

とRunge-Kutta4次公式を示した連立微分方程式 $$ \begin{aligned} \frac{dy}{dx} &= f(x,y,z), \, where \, y(x_0)=y_0 \\ \frac{dz}{dx} &= g(x,y,z), \, where \, z(x_0)=z_0 \\ \end{aligned} $$ とを比べて，変数の表記の違いと関数$f,g$を具体的に書き下します．

4次の公式	運動方程式
x	t
y	x
z	v
f(x,y,z)	f(t,x,v) = v
g(x,y,z)	g(t,x,v) = -k x

この変数の書き換えを吸収する中間関数としてEuler3を書き換えます．RungeKutta4の仮引数を上の表に従って置き換えて，数値を渡しています．また，関数$f,g$も定義しておきます．

これを前と同様に走らせると 発散も収束もすることなく，定常的に振動を繰り返していることが見て取れます．

太陽系¶

分子動力学=Verlet(ベルレ)法¶

惑星の軌道は摩擦のない宇宙空間での２体の重力だけで考えることができる単純な動きなので，比較的容易に再現することができます． この運動を記述する方程式がニュートンの運動方程式 $$ F = ma $$ です．ここで$F$はForce（力）, $m$はmass（質量）, $a$はaccelleration（加速度）です．

$a$ を微分で明示すると， $$ F = m \frac{d^2}{dt^2}r(t) $$ となります． $r(t)$がある時間での惑星の位置をさしているとします． これを微小時間$h$で級数展開して $$ r(t + h) = r(t) + \frac{dr}{dt} h + \frac{1}{2} \frac{d^2r}{dt^2} h^2 + O(h^3) $$ と $$ r(t - h) = r(t) - \frac{dr}{dt} h + \frac{1}{2} \frac{d^2r}{dt^2} h^2 + O(h^3) $$ での和をとって高次項を無視して，$r(t+h)$について解いて，$d^2r/dt^2=F/m$を代入すると， $$ r(t + h) = 2r(t) - r(t-h) + \frac{F}{m} h^2 $$ となり，$r(t+h)$の位置を$r(t - h)$と$r(t)$での情報から予測することができます．

このように微分方程式を差分方程式に変換し，数値的に解いていく方法がVerlet法で，そのほかの種々の分子動力学法の手法の基本となる考え方です．