Pythonで学ぶ常微分方程式の数値解

数学の具体的な計算にPythonを使って、数学もPythonも同時に学んでしまいましょう。今回はPythonを使って、常微分方程式の解を数値的に求めてみたいと思います。本記事で記載している図の描画の確認は、JupyterNotebook上で行っています。インターフェイスなどの動作環境の違いによって適宜変更点があるかもしれません。

常微分方程式

常微分方程式

一般に、未知関数とその導関数を含んだ方程式を微分方程式といいます。一変数関数の場合を常微分方程式、多変数関数の場合を偏微分方程式と呼びます。Pythonに実装されているコマンドを使って、力学で現れる具体的な常微分方程式をいくつか解析的に解いているのが次の記事です：

pianofisica.hatenablog.com

本記事では、必ずしも解析的な解が得られない場合も多いことから、初期条件から数値的に常微分方程式の解を求める方法をみていきます。

裳華房の『数値計算』（柳田・中木・三村）では、数値計算の基本的な手法を紹介しています：

<br />

常微分方程式を数値的にいかにして解くかを数学的背景を基礎に解説していて、原理の面から理解したいという人にはうってつけの本だと思います。

Euler法（オイラー法）

１階常微分方程式

$\quad \displaystyle{\frac{dy}{dx}=f(y)}$

を考えましょう。微分の定義から、微小変位 $\Delta x$ に対して

$\quad \displaystyle{\frac{y(x+\Delta x)-y(x)}{\Delta x}\sim f(y(x))+\mathcal{O}(\Delta x^2)}$

が成り立ちます。そこで $h$ を微小な差分間隔として

$\quad \displaystyle{ y(x_{n+1}) = y(x_{n})+h\,f(y(x_{n})) }$

によって各 $x_n=x_0+nh \ (n=0,1,\dots,N)$ での $y$ の値を求めることができます。

具体例として $f(y)=-y$ として

$\quad \displaystyle{\frac{dy}{dx}=-y}$

という微分方程式を考えてみます。この微分方程式の一般解は $y=Ce^{-x}$ で与えられますが、数値的な計算が解析的な解とどれほど違ってくるかを体験してみましょう。初期条件として $C=20$ の場合、差分間隔を $h=0.1$ として計算してみます：

N = 60
h = 0.1
C = 20
X = [0]
Y = [C]
def x(k): return h+X[k-1]
def y(k): return Y[k-1]+h*(-Y[k-1])
for k in range(N):
      X.append(x(k+1))
for k in range(N):
      Y.append(y(k+1))
import matplotlib.pyplot as plt
plt.plot(X,Y)
plt.show()

によって得られるのが次のグラフです：

これと微分方程式の解析的な解 $y=20e^{-x}$ を合わせてグラフにしてみましょう：

import numpy as np
P = 100
xmin = 0
xmax = h*N
p = np.linspace( xmin, xmax, P)
q = 20*np.exp(-p)
plt.plot(X,Y)
plt.plot(p,q)
plt.show()

から作成されるのが次のグラフです：

青線が数値的に求めた解で、オレンジ線が解析的な解です。定性的には解の様子をうまく表していますが、あまり精度が良いようには見えません。

改良Euler法（Runge-Kutta法：ルンゲ=クッタ法）

差分間隔を保ったまま、評価点の個数も変えずに、計算の精度を上げるにはどのようにすれば良いでしょうか？ここでTaylor展開の公式から、微小な差分間隔 $h$ に対して

$\quad \displaystyle{\begin{aligned}y(x_{n+1})&\sim y(x_n)+h\,f(y(x_n))+\frac{h^2}{2}\,y''(x_n)\\&\sim y(x_n)+h\,f(y(x_n))+\frac{h^2}{2}\,f'(y)\,f(y(x_n))\end{aligned}}$

と近似されることを使って $h^2$ の補正項を考慮に入れたのが以下の計算です：

N = 60
h = 0.1
C = 20
X2 = [0]
Y2 = [C]
def x(k): return h+X2[k-1]
def y(k): return Y2[k-1]+h*(-Y2[k-1])+h**2/2*(-1)*(-Y2[k-1])
for k in range(N):
      X2.append(x(k+1))
for k in range(N):
      Y2.append(y(k+1))
import matplotlib.pyplot as plt
plt.plot(X2,Y2)
plt.plot(p,q)
plt.show()

から得られるのが次の図です：

（見づらいですが）解析的な解の曲線とかなりよく一致し、近似が改良されていることがわかります。

２階常微分方程式

運動方程式にみられるように、２階微分方程式

$\quad \displaystyle{\frac{d^2y}{dx^2}=f\left(y,\frac{dy}{dx}\right)}$

は応用上の重要な対象です。これを上のようにして数値的に解析したいわけですが、そのために連立方程式による１階化

$\quad \displaystyle{\frac{dy}{dx}=p\qquad\qquad\frac{dp}{dx}=f\left(y,p\right)}$

を考えます。この連立微分方程式から得られる近似の関係式

$\quad \displaystyle{\begin{aligned} y(x+h) &=y(x)+hp(x)+\frac{h^2}{2}p'(x)+\dots\\ &=y(x)+hp(x)+\frac{h^2}{2}f(y,p)+\dots \\ p(x+h) &=p(x)+hf\left(y,p\right) +\frac{h^2}{2}\left(p\frac{\partial f(y,p)}{\partial y}+f\frac{\partial f(y,p)}{\partial p}\right) +\dots \end{aligned}}$

を用いて上と同様の計算ができます。

もう少し具体的に詳しくみるために単振動の運動方程式

$\quad \displaystyle{\frac{d^2y}{dx^2}=-y}$

を考えましょう。いま連立微分方程式によって線形化して

$\quad \displaystyle{\frac{dy}{dx}=p \qquad\qquad \frac{dp}{dx}=-y}$

です。初期条件として $y(0)=1,\ p(0)=0$ としましょう。Euler法で計算すると

N = 10000
h = 0.01
T = [0]
y_pos = [1]
p_vel = [0]
def t(k): return h+T[k-1]
def y(k): return y_pos[k-1]+h*(p_vel[k-1])
def p(k): return p_vel[k-1]+h*(-y_pos[k-1])
for k in range(N):
      T.append(t(k+1))
      y_pos.append(y(k+1))
      p_vel.append(p(k+1))
import matplotlib.pyplot as plt
Q = np.cos(T)
plt.plot(T,y_pos)
plt.plot(T,Q)
plt.show()

から得られるのが次の図です：

青線は数値計算によって得られる解、オレンジの線が解析解 $y(x)=\cos(x)$ です。数値解は、周期は一致しますが、振幅が徐々に大きくなってしまっていることがわかります。

２次の補正まで入れて（Runge-Kutta法で）計算したのが以下です：

N = 10000
h = 0.01
T = [0]
y2_pos = [1]
p2_vel = [0]
def t(k): return h+T[k-1]
def y2(k): return y2_pos[k-1]+h*(p2_vel[k-1])+h**2/2*(-y2_pos[k-1])
def p2(k): return p2_vel[k-1]+h*(-y2_pos[k-1])+h**2/2*(-p2_vel[k-1])
for k in range(N):
      T.append(t(k+1))
      y2_pos.append(y2(k+1))
      p2_vel.append(p2(k+1))
import matplotlib.pyplot as plt
Q = np.cos(T)
plt.plot(T,y2_pos)
plt.plot(T,Q)
plt.show()