3. 常微分方程式(Ordinary Differential Equation, ODE)

　常微分方程式(Ordinary Differential Equation, ODE)は物理のあらゆる場面で現れる。この微分方程式を数値計算で解いてみよう。

3.1. odeint

　微分方程式(ODE)のアルゴリズムには、オイラー法、Runge-Kutta法などがある。一度くらいはこれらの解法の仕組みを知っておく必要があるが、それは後回しにして、SciPyに用意されているODEを解くfunctionであるodeintを使って、実際に計算してみることから始めよう。なお、ODEを解くfunctionのことをODEソルバーと言う。

(例題) 微分方程式、$\frac{dx}{dt}=-x$を初期条件$x(0)=1$で解き、$0\leq t \leq 10$の範囲で図示せよ。

　まずは、必要なモジュールをインポートする。

In [1]:
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

　次に$\frac{dx}{dt}$をfunctionとして定義する。

In [2]:
def dxdt(x, t):
    dxdt = -x
    return dxdt

　$t$の範囲を指定し、初期値$x(0)$はx0とする。

In [3]:
t = np.arange(0, 10., 0.1)
x0 = 1

　最後にodeintを使って$x$を求めてプロットする。

In [20]:
x = odeint(dxdt, x0, t)
plt.plot(t, x)
plt.xlabel('$t$')
plt.ylabel('$x$');

　実行すると指数関数的に減少するグラフが書ける。縦軸をlogスケールにしたければ、plt.yscale('log')を追加するか、plt.semilogyでプロットするとよい。

In [21]:
plt.plot(t, x);
plt.xlabel('$t$')
plt.ylabel('$x$')
plt.yscale('log');
In [11]:
plt.semilogy(t, x)
plt.xlabel('$t$')
plt.ylabel('$x$');

　このように直線で減少していくことが分かる。

　以上のように、解きたいODEを関数で定義すれば、ODEを解くアルゴリズムを知らなくても簡単に解ける。コードを書くのに慣れるために以下の練習問題を試してみることをお勧めする。

練習問題

(3.1.1.) 微分方程式、$at^2 \frac{dx}{dt}+b+x^2=0$を$a=2$、$b=5$、初期条件$x(1)=2$で解き、$1\le t \le 10$の範囲で図示せよ。

(3.1.2.) 素粒子が崩壊して別の素粒子に変わるとき、粒子数が減る速さは残存する原子数$N$に比例する。つまり、比例定数を$\lambda$として、

\begin{equation*}
\frac{dN}{dt} = -\lambda N
\end{equation*}

と書くことができる。年代測定に使われる炭素C14の場合、$\lambda=1.2096 \times 10^{-4}$年である。初期値$N(0)=10000$の場合に$N$がどう変化するか、$2$万年まで図示するコードを書け。また、粒子数が半分になる時間(半減期)を求めよ。

解答例

(3.1.1.) この式のままだとfunctionに直せないので、式変形して、$\frac{dx}{dt}=-\frac{b+x^2}{at^2}$にする。

In [22]:
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
In [23]:
def dxdt(x, t, a, b):
    dxdt = -(b + x**2)/(a*t**2)
    return dxdt

　今の例では、a, bは定数なのでこの時点で具体的に数字を入れてもいいのだが、汎用性を持たせるためにdxdtの中に具体的な値は示さずに、odeintの引数に設定する。具体的には、odeintの引数の最後にargs=(2,5)と書けばよい。argsは引数(argument)の略。a, bの順番に()の中に並べて書く。詳細はSciPyの公式参照。

In [26]:
x0 = [2, 1]
t = np.arange(1, 10., 0.1)
x = odeint(dxdt, x0, t, args=(2,5))
plt.plot(t, x[:,0])
plt.xlabel('$t$')
plt.ylabel('$x$');

　一階の微分方程式なので、最初の例と同じように指数関数的に減少していく。

(3.1.2.) 物理っぽい問題にしているが、数学的には特に難しいことはないだろう。

In [27]:
def dNdt(N, t):
    dNdt = -lamb * N
    return dNdt
In [28]:
lamb = 1.2096e-4 
N_init = 10000
N_0 = [N_init, 0]
t = np.arange(0, 20000, 1)
N = odeint(dNdt, N_0, t)
plt.plot(t, N[:,0])
plt.xlabel('t')
plt.ylabel('N');

　半減期は存在する粒子数が半分になるまでの時間のことなので、そのまま素直に数値計算の結果から求めるなら、argminとabsを使って粒子数が半分N_init/2になるN配列番号を取得する。

In [29]:
ind, n = np.argmin(abs(N-N_init/2), axis=0)
print(t[ind])

5730

　と求められる。1行目のabs(N-N_init/2)とは、N_init/2からNまでの距離(絶対値abs)を要素に持つ配列である。この配列の要素のうち最小となる要素が、粒子数が半分に一番近いことになるので、その要素番号をnp.argminを使ってindという変数に格納する、というやり方で求めている。

　実は数値計算の計算結果と関係なく$\lambda$から$\tau = \frac{1}{\lambda}\log 2$として半減期は求められるので、直接求めると、

In [30]:
1/lamb * np.log(2)
Out[30]:
5730.383437168861

となる。

当サイトのテキスト・画像の無断転載・複製を固く禁じます。

Unauthorized copying and replication of the contents of this site, text and images are strictly prohibited.
© 2019 Go Yusa