python入門: 18. Euler法で微分方程式を解こう

微分方程式は自然科学・社会科学の広い分野で使われる, 強力な数学的ツールである。微分方程式を解析的に解くには数学の訓練が必要で大変だが、これまで学んだ方法を使うと, 微分方程式を数値的に解くことができる。その最も簡単な方法であるEuler法(Eulerはオイラーと読む)を学ぼう。

考え方はシンプルである。方程式に基づいて関数の増加分の値を逐次計算して、それを数値積分すればよい。それを最も単純に行うのがEuler法である。

例として, ロジスティック方程式

dN/dt = a N - b N²

を解いてみよう(資源に限りがある場での生物の繁殖モデル)。初期条件をN(0)=10とし, a=1, b=0.05としよう。

まず, いつものお約束で, 数学計算ライブラリのnumpyと, グラフを描くためのモジュールmatplotlib.pyplotを用意する:

import numpy as np                                                      
import matplotlib.pyplot as plt 

次にt (時刻)を配列（numpy.ndarrayのインスタンス）として用意する。ここでは0以上9未満とし, 刻みdtは0.2にしておこう:

t=np.arange(0, 9, 0.2)

次に, Nの配列を用意する。といってもまだNの値は求まっていないので, とりあえず全部0で埋めておこう:

N=np.zeros(t.size)

np.zerosというのは, 項がぜんぶ0であるような配列を作る関数である。Nはいずれtと対応して値が定まるので, tと同じ数の項を用意しておくために, Nの項数としてt.size（=tの項数）を指定している。

次に初期条件を設定する:

N[0]=10

ところでこの微分方程式は,

dN = ( a N - b N² )dt

すなわち,

N(t+dt) = N(t)+( a N - b N² )dt

と書ける（N(t)をNと省略されていたのである）。ここでdtはさきほど0.2と決めた。すると, ある時刻tでのN(t)がわかっていれば, この右辺を計算してN(t+dt)が求まるのだ。それを繰り返していけば, tの全ての値について, N(t)が求まる。pythonで表現すれば, たとえば

N[1]=N[0]+(a*N[0]-b*N[0]*N[0])*0.1

となる, ということだ。ここでN[0]はN(t₀)を意味し, N[1]はN(t₁)を意味することに注意せよ（つまり数学の関数N()と, pythonのコードのN[]は, 括弧の中の意味・扱いが違う。前者はtの値で, 後者は0から始まる項番号である。）

これをfor文で繰り返せばよい:

a, b = 1, 0.05
for k in range(t.size-1):
    N[k+1] = N[k] + (a*N[k]-b*N[k]*N[k])*0.2

これで終了である。ではグラフに描いてみよう:

plt.plot(t, N)
plt.show()

よく見慣れた, 個体群動態のグラフ（Logistic曲線）が現れるだろう。

課題18-1: 以上の操作を再現せよ。

さて, これでもよいのだが, ちょっと工夫すると, 良いプログラムになる:

import numpy as np                                                      
import matplotlib.pyplot as plt 
t0, t1, dt = 0, 9, 0.2  # 時刻の開始, 終了, 刻みをここで設定する。
a, b = 1, 0.05   # ロジスティック方程式のパラメータをここで設定する。
y0 = 0.05        # 初期条件をここで設定する。
t=np.arange(t0, t1, dt)
y=np.zeros(t.size)
y[0]=y0
for k in range(t.size-1):
    y[k+1]=y[k]+(a*y[k]-b*y[k]*y[k])*dt
ya=1/((b/a) + (1/y0 - b/a)*np.exp(-a*t)) # 比較のために解析解も作る。
plt.plot(t, y)
plt.plot(t, ya)
plt.show()

↑このプログラムのポイントは, 時刻の範囲や刻み, パラメータ, 初期条件など, 条件設定に関する数値をプログラムの冒頭部分にまとめたことである。こうすることで, 条件を変えたらどうなる? というときに変更や調整がやりやすくなる。

課題18-2: 以下の微分方程式をEuler法で解け:

dy/dt = -y t

初期条件はt =0のときy =1とする。解析的に解ける人は、解析解と比較せよ。

課題18-3: Lotka-Volterra方程式

ある島に羊と狼がそれぞれ何匹かいる状況を考えよう。羊は草を食べて生活・繁殖するが、狼は羊を食べて生活・繁殖する。

羊の頭数の変化＝( 羊の数×自然繁殖率 - 狼の数×羊の数×捕食による死亡率)×時間間隔 狼の頭数の変化＝(-狼の数×自然死亡率＋狼の数×羊の数×捕食による繁殖率)×時間間隔

注：　狼と羊がばったり出くわす頻度（つまり狼が羊をつかまえて食べることのできる頻度）は、狼の数と羊の数の積に比例する。

時刻をt, 羊の頭数をS, 狼の頭数をW, 羊の自然繁殖率をp, 捕食による羊の死亡率をq, 狼の自然死亡率をr, 捕食による狼の繁殖率をsとすると、

dS/dt = pS - qWS

dW/dt = -rW + sWS

が成り立つ（と考える）。あるいは同じことだが,

dS = （ pS - qWS）dt

dW = （-rW + sWS）dt

である。この連立微分方程式を、Lotka-Volterra(ロトカ・ヴォルテラ方)程式と呼ぶ。p, q, r, sは定数（パラメータ）である。

Lotka-Volterra方程式を、p=1.1, q=0.6, r=1.3, s=0.6のもとで解け。初期条件(t=0のときの値)はS=2.0, W=1.0とする。刻みdtは0.01とし, t=0からt=20程度まで解け。結果をグラフにしてみよ。（以下にプログラムの例を示すが, これを見ないで自分で組んでみよ。うまくいったら, プログラムを比較してみよ）

# Lotka-Volterra方程式のプログラムの例
import numpy as np
import matplotlib.pyplot as plt
t0, t1, dt = 0, 20, 0.001    # 時刻範囲と刻み
p, q, r, s = 1.1, 0.6, 1.3, 0.6  # パラメータ設定
S0, W0=2.0, 1.0             # 初期条件
t=np.arange(0, 20, dt)
S=np.zeros(t.size)
W=np.zeros(t.size)
S[0], W[0] = S0, W0
for k in range(t.size-1):
    S[k+1]=S[k]+( p*S[k] - q*S[k]*W[k])*dt
    W[k+1]=W[k]+(-r*W[k] + s*S[k]*W[k])*dt
plt.plot(t, S)
plt.plot(t, W)
plt.show()

参考になる教材： 「大学生のための応用数学入門」P37~39。

課題18-4:　Lotka-Volterra方程式で、初期条件をいろいろ変えると解はどのように変わるか観察せよ。

課題18-5:　パラメータをいろいろ変えると解はどのように変わるか観察せよ。

（ブラウザの「戻る」ボタンで戻ってください。）