Python3で数理生態

常微分方程式モデル

ここでは、logistic growth model, prey-predator population dynamics model等のお馴染みの古典個体群動態モデルを使って、以下の数値計算のアルゴリズムの解説とコーディングを行う．

[1] 陽的オイラー法（一次元）

[2] 4次のルンゲクッタ法（一次元）とオイラー法との精度の違い

[3] 4次のルンゲクッタ法（多次元）：汎用的な関数として定義する（関数を引数とする関数）

[4] C, Rとの速度比較

[5] Pythonの高速化

[6] クラス、メソッドとしてオイラー法とルンゲクッタ法を実装する

[1] 陽的オイラー法（一次元）

一番単純な場合として，一次元の常微分方程式（1D ODE）を考える．下の一つ目の式はN(t)に関する方程式、二つ目の式は初期条件（t = 0における個体群サイズ）．f(N(t), t)は、微分方程式の右辺がN(t)とtの関数であることを表す．右辺が陽に時刻tに依存しないとき、たとえば、通常のlogistic growth model: dN(t)/dt = r*N(t)*(1 - N(t)/K) 、二つ目のパラメータtは必要なく、f(N(t))とすればよい(autonomous ODE)．

陽的オイラー法は，以下のような方法（式１）である．数学的な厳密解とオイラー法による近似解を区別するために，大文字N(t)は厳密解、小文字n(t)は近似解とする．式１は，時刻t = t_1のときの近似解n(t_1)から，時間hだけ経過後の近似解n(t_1 + h)を求めるアルゴリズムである．式１を変形すると式２が得られる．この式2は「差分モデル」のページで扱った数式と同じ形式なので，for loopを使って式２を時間（インターバル）hのステップごとに逐次的に解いていけば上記の常微分方程式の近似解が必要な期間だけ(0 <= t <= T) 得られる.

陽的オイラー法は精度の点から全く実用的ではないが，常微分方程式を「差分化」して数値的に説く方法のエッセンスが詰まっているのでここで解説する．

ここで注意しなければいけないのは，厳密解N(t)と近似解n(t)の差（＝局所的誤差）である．もっと具体的には，「時刻t = t_1時点で近似解が厳密解と同一だとしたとき(n(t_1) = N(t_1))，一つのステップでどれだけ誤差が生じるか」について考える必要がある．厳密解に関して, 時刻t = t_1のまわりでテイラー展開を用いれば時刻t = t_1 + hにおける厳密解を式３のように書き下すことができる．ここで式３の右辺の最後の項は，hに関する３次以降の高次項をまとめたものである（高次の無限小）．ここでテイラー展開のhの一次の係数(t = t_1での偏微分係数）は元の微分方程式の右辺をt=t_1で評価した値となるため，実は式２の右辺の第二項の係数f(n(t_1), t_1)と一致する．したがって厳密解と近似解の差は，式３から式２を引き算すればよく，以下のように計算できる．

エラーの大きさはステップ幅hの二乗に比例する．そこで「精度(accuracy)」の点からはステップ幅hの一乗に比例するので、first order accuracyを持つことになる．ちなみに式２，３、誤差の関係は以下のように図示できる．

では実際に，以下のモデルを具体例に数値計算をしてみる．以下のモデルはロジスティック成長モデル(dN/dt= rN(1-N/K)の環境収容力Kに周期的な時間変動項(a*sin(2*pi*t))を入れたモデルとなっている．

このモデルを陽的オイラー法で解いて結果を図示するには以下のようなスクリプトを書けばよい．

まずはいくつか必要なモジュール、関数をインポートする．

import numpy as np   #for numerical calculations

import pandas as pd  #read with pandas as dataframe

import math    #for mathematical functions

import matplotlib.pyplot as plt   #for graphics

from numpy import exp

from numpy import sin

from numpy import pi

それからモデルに必要なパラメータ(r, K,a)と数値計算に必要なパラメータ(end_time, h)を用意する．

#Logistic growth model with seasonally-fluctuating carrying capacity

r = 1.0  #growth parameter

K = 1.0 #constant part of carrying capacity

a = 0.5 #amplitude of fluctuation

end_time = 10.0  #end time of simulation

h = 0.1 #time interval for integration

注目する微分方程式に特有の関数を定義する．

#definition of differential equation

def f1(x, t):

    temp = r*x*(1.0 - x/(K + K*a*sin(2.0*pi*t)))

    return temp

初期条件（個体群サイズnと時刻time），データ書き込み間隔を調整するパラメータを用意する．

n = 0.1 #initial condition (initial population size)

time = 0.0 #initial condition (initial time, 0)

write_index = 1 #need for counting the numerical step

write_interval = 0.1 #with this interval, the result will be saved

初期条件（数値）を文字列に変換後，ファイルに書き込む．

#write initial condition

text=str(time)+','+str(n)+'\n'

with open('result_ode.csv', 'w') as f:

    f.write(text)

forループを用いて微分方程式の数値解を求め、write_intervalごとにファイルに書き込む．すべての間隔hごとに結果をベクトルに保存したりファイルに保存したりするのはメモリのムダ使いになるため注意．

for i in range(0, int(end_time/h)):

    n = n + h*f1(n, time)

    time += h

    if(write_index < int(write_interval/h)):

        write_index += 1

    else:

        text = str(time) + ',' + str(n) + '\n'

        with open('result_ode.csv', 'a') as f:

            f.write(text)

        write_index = 1

結果を再度ファイルから読み出し後，plot関数で線グラフを描く．

#load result as a dataframe

table=pd.read_csv('result_ode.csv', header=None)

#line plot

plt.plot(table.iloc[:,0], table.iloc[:,1], color='blue', linestyle='solid')

下のグラフはh=0.1のときのもの．

さっきの精度の話題の例としてh=0.01とセットしなおして，それ以降の部分について保存先を'result_ode2.csv'と変えた上で実行し直し，グラフを重ねて描いて数値計算の結果を視覚的に比較してみる．赤い線がh=0.01のときの結果となる．

#comparison

table2=pd.read_csv('result_ode2.csv', header=None)

plt.xlim([7,10])

plt.ylim([0.65, 0.95])

plt.plot(table.iloc[:,0], table.iloc[:,1], color='blue', linestyle='solid')

plt.plot(table2.iloc[:,0], table2.iloc[:,1], color='red', linestyle='solid')

plt.legend(['h=0.1', 'h=0.01'])

これは些細な違いとは言えない．一般にhを小さくすれば精度はよくなっていきますが，hを10倍細かくすれば当然計算時間は10倍かかることになる．オイラー法以来，多くの手法が開発されてきたがその中でも古典的なものの一つとして４次の(陽的)Runge-Kutta(ルンゲクッタ)法について次に解説，オイラー法との比較を行う．

[2] 4次のルンゲクッタ法（一次元）とオイラー法との精度の違い

ここでは，（固定幅）（陽的）４次のルンゲクッタ法（4th order explicit Runge-Kutta method with fixed interval）について解説する．アイデアとしては，t=t_1で一度限り微分係数を計算して次ステップの解を近似するという陽的オイラー法のアルゴリズム（２）を変更し，何度も微分係数を計算しそれらを重み付けして「平均値」を計算して次ステップの解を近似することにする．具体的には、以下のk1からk4までを計算して重み付けして平均値(1/6)*(k1 + 2*k2 + 2*k3 + k4)を求め，それを使うことにする．

もしも「平均値」の計算にk1しか使わなかった場合は陽的オイラー法と同一となることに注意．オイラー法のときと同様にテイラー展開を使えば（計算は省略）以下のように誤差が刻み幅hの５次のオーダーで生じることを示すことができる．「４次」と名前に付くゆえんは、精度がhの４次のオーダーとなるからである．

このアルゴリズムを特に深く考えることもなく関数を使わずに実装するのは非常に簡単で以下のようにやればよい．先ほどのオイラー法の結果を重ね合わせて結果もグラフにしてみる．

#simple way of numerical calculation for 4th-order Runge Kutta for 1D ODE

n = 0.1 #initial condition (initial population size)

time = 0.0 #initial condition (initial time, 0)

write_index = 1 #need for counting the numerical step

write_interval = 0.1 #with this interval, the result will be saved

h = 0.1

#write initial condition

text=str(time)+','+str(n)+'\n'

with open('result_ode3.csv', 'w') as f:

    f.write(text)

for i in range(0, int(end_time/h)):

    k1 = f1(n, time)

    k2 = f1(n + (h/2.0)*k1, time + h/2.0)

    k3 = f1(n + (h/2.0)*k2, time + h/2.0)

    k4 = f1(n + h*k3, time + h)

    n = n + (h/6.0)*(k1 + 2.0*k2 + 2.0*k3 + k4)

    time += h   #update time

    ###recored the result every "write_interval" only

    if(write_index < int(write_interval/h)):

        write_index += 1

    else:

        text = str(time) + ',' + str(n) + '\n'

        with open('result_ode3.csv', 'a') as f:

            f.write(text)

        write_index = 1

#load result as a dataframe

table3=pd.read_csv('result_ode3.csv', header=None)

#comparison

plt.xlim([7,10])

plt.ylim([0.65, 0.95])

plt.plot(table.iloc[:,0], table.iloc[:,1], color='blue', linestyle='solid')

plt.plot(table2.iloc[:,0], table2.iloc[:,1], color='red', linestyle='solid')

plt.plot(table3.iloc[:,0], table3.iloc[:,1], color='green', linestyle='solid')

plt.legend(['Euler: h=0.1', 'Euler: h=0.01', 'rk4:h=0.1'])

グラフは以下の感じになる．ルンゲクッタ法(rk4)を使えば，刻み幅をh=0.1と大きくしてもオイラー法のh=0.01の時の結果と大差ないことが見て取れるだろう．

では以下のコードでシステマティックに誤差の大きさを比較してみる．刻み幅hの初期値を0.2にセット後，ステップごとに3分の１に小さくしていき、オイラー法とルンゲクッタでどれだけ誤差に差が出るかを計算する．一番小さいhの値でルンゲクッタ方を用いたときの時刻t = 20の個体群サイズを「真の」値と仮定し，それぞれの方法、それぞれのhの値を用いたときの数値解と「真の値」との絶対値を対数スケールでプロットした．

#Compare the performance of Euler and rk4:

vec_comp = np.zeros([20,3])  #vector to record the result

h = 0.2

end_time = 20.0

for i in range(0, 13):

    n_ee = 0.1

    n_rk4 = 0.1

    time = 0.0

    for j in range(0, int(end_time/h)):

        n_ee = n_ee + h*f1(n_ee, time)  #Explicit Euler one step

        #for rk4

        k1 = f1(n_rk4, time)

        k2 = f1(n_rk4 + (h / 2.0) * k1, time + h / 2.0)

        k3 = f1(n_rk4 + (h / 2.0) * k2, time + h / 2.0)

        k4 = f1(n_rk4 + h * k3, time + h)

        n_rk4 = n_rk4 + (h / 6.0) * (k1 + 2.0 * k2 + 2.0 * k3 + k4)

        time += h   #update time

    vec_comp[i,0] = h

    vec_comp[i,1] = n_ee #record the value n from Euler method  at t = end_time

    vec_comp[i,2] = n_rk4 #record the value n from rk4 method at t = end_time

    h /= 3.0

plt.plot(vec_comp[0:i,0], abs(vec_comp[i,2] - vec_comp[0:i,1]), color='red', linestyle='solid', marker="o")

plt.plot(vec_comp[0:i,0], abs(vec_comp[i,2] - vec_comp[0:i,2]), color='green', linestyle='solid', marker="o")

plt.xscale('log')

plt.yscale('log')

plt.ylim([1.0e-9,0.1])

plt.grid(which='major', color='black', linestyle='-')

plt.legend(['Euler', '4th-RK'])

4次のルンゲクッタ法ではh=0.01程度で真の値との差が限界まで小さくなっている(1.0e-8以下)のにオイラー法は仮にh=0.000001としても同じ誤差まで下げることはできていない．オイラー法は話にならないことが分かるだろう．４次のルンゲクッタ法だって多くの実用的・工学的場面によっては役に立たないのだが（たとえ基礎科学であってもカオス的挙動の解析などでは信頼できない），ここではこれ以上は深入りしない．「刻み幅自動調節」、「精度保証付き」などのキーワードでググれば，微分方程式の数値解法の深淵を覗くことができるだろう．

[3] 4次のルンゲクッタ法（多次元）：汎用的な関数として定義する（関数を引数とする関数）

配列に関して要素ごとにfor loopで繰り返す場合と、配列全体での計算式を使う場合では大きく計算時間が異なることに注意．当然ループを使うほうが遅い．参照元：http://japanichaos.appspot.com/NumpyArray.html

しかし注意しないといけないのは、関数の引数として配列を用いたときの挙動である．まずは「プログラミングの基礎」の”関数，値渡し，参照渡し ”のセクションを復習しておく．

さて，ルンゲクッタ法を多次元化するのは実は簡単で，上のk1, k2, k3, k4を使ったアルゴリズムにおいてすべてをベクトル化（配列化）すればよい．微分方程式も多次元化するので，以下のような配列を引数とする関数で微分方程式の右辺を定義する必要がある．ただしまずは上で扱った一次元のロジスティック成長モデルに適用するので実質的に一次元のままであることに注意する．

#definition of differential coefficients, should be finally a part of method

def diff(in_vec, out_vec, t):

    out_vec[0] = f1(in_vec[0], t)

そして多次元ルンゲクッタは、直感的には以下のような関数を用意すればよいのだが・・・．実は引数としての配列の更新がうまくいかない．

def rk4_test(in_vec, out_vec, time, h_interval, dim, diff_vec):

   #in_vec: the vector with the value at "time"

   #out_vec: the vector with which the updated vector value after one time step is saved

   #time: time value

   #h_interval: time step for discritizing ODE

   #dim: dimension of the ODE

   #diff_vec: function that determines the r.h.s. of the ODE

    k1 = np.zeros([dim])

    k2 = np.zeros([dim])

    k3 = np.zeros([dim])

    k4 = np.zeros([dim])

    temp_vec = np.zeros([dim])

    h_half = h_interval / 2.0

    t_h = time + h_half

    diff_vec(in_vec, k1, time)  # calculate k1

    temp_vec = in_vec + h_half * k1

    diff_vec(temp_vec, k2, t_h)  # calcualte k2

    temp_vec = in_vec + h_half * k2

    diff_vec(temp_vec, k3, t_h)  # calcualte k3

    temp_vec = in_vec + h_interval * k3

    diff_vec(temp_vec, k4, time + h_interval)  # calcualte k4

    out_vec = in_vec + (h_interval / 6.0) * (k1 + 2.0 * k2 + 2.0 * k3 + k4)

これは任意の微分方程式を解くための情報をパラメータにした関数である．目新しいところは関数もパラメータにしていることであり，微分方程式の右辺を定義した関数を引数としてdiff_vecに渡せば任意の微分方程式が解ける．さらにforループを使わず計算速度も速いはずである．しかし，これを使って計算してみると以下のスクリプトでは一次元ベクトルnvの値が全く更新されないことが分かる．

nv = np.zeros([1])

nv[0] = 0.1  # initial condition (initial population size)

time = 0.0  # initial condition (initial time, 0)

write_index = 1  # need for counting the numerical step

write_interval = 0.1  # with this interval, the result will be saved

h = 0.01

# write initial condition

text = str(time) + ',' + str(nv[0]) + '\n'

with open('result_ode3.csv', 'w') as f:

    f.write(text)

for i in range(0, int(end_time / h)):

    rk4_test(nv, nv, time, h, 1, diff)  #update vector value

    #rk4(nv, nv, time, h, 1, diff) #update vector value

    time += h  # update time

    ###recored the result every "write_interval" only

    if (write_index < int(write_interval / h)):

        write_index += 1

    else:

        text = str(time) + ',' + str(nv[0]) + '\n'

        with open('result_ode3.csv', 'a') as f:

            f.write(text)

        write_index = 1

# load result as a dataframe

table = pd.read_csv('result_ode3.csv', header=None)

print(table)

実は関数rk_4の最後の行だけが問題で（この行では、関数内部で新しい変数out_vecが引数のout_vecとは異なるIDで生成されてしまい、引数の中身が更新されない），ここをforループで書けばちゃんと引数の中身の値が変更させる．

for i in range(0, dim):

    out_vec[i] = in_vec[i] + (h_interval / 6.0) * (k1[i] + 2.0 * k2[i] + 2.0 * k3[i] + k4[i])

したがって、以下のような関数を作って上のスクリプトのrk4_testの呼び出しをrk4で置き換えればうまくいく（コメントアウトしているrk4の行を復活させればよい）．

def rk4(in_vec, out_vec, time, h_interval, dim, diff_vec):

    #in_vec: the vector with the value at "time"

    #out_vec: the vector with which the updated vector value after one time step is saved

    #time: time value

    #h_interval: time step for discritizing ODE

    #dim: dimension of the ODE

    #diff_vec: function that determines the r.h.s. of the ODE

    k1 = np.zeros([dim])

    k2 = np.zeros([dim])

    k3 = np.zeros([dim])

    k4 = np.zeros([dim])

    temp_vec = np.zeros([dim])

    h_half = h_interval / 2.0

    t_h = time + h_half

    diff_vec(in_vec, k1, time)  # calculate k1

    temp_vec = in_vec + h_half * k1

    diff_vec(temp_vec, k2, t_h)  # calcualte k2

    temp_vec = in_vec + h_half * k2

    diff_vec(temp_vec, k3, t_h)  # calcualte k3

    temp_vec = in_vec + h_interval * k3

    diff_vec(temp_vec, k4, time + h_interval)  # calcualte k4

    for i in range(0, dim):

        out_vec[i] = in_vec[i] + (h_interval / 6.0) * (k1[i] + 2.0 * k2[i] + 2.0 * k3[i] + k4[i])

さて、一次元でうまくいったのでこの関数を以下のprey-predatorモデル（二次元）に使ってみよう．x， yはそれぞれprey，predatorの個体群密度を表す．

まずはモデルのパラメータを用意する．

#define parameters for ODE for prey predator model

r = 1.0

a = 1.0

b = 0.5

H_t = 1.0

m = 0.1

K  = 1.0

次に，x, yというモデル変数とそれを保存するベクトルの各要素(0, 1)とを対応付けるように整数インデックスを用意する．このようなトリックは二次元程度のモデルでは必要性を感じないかもしれないが、植物、動物、微生物等何十何百もの種をモデリングする際などモデル式の変数とプログラムコードの中のベクトル要素を一対一に体操させるときに便利である．

#index to link the model variables (x, y) and the elements of the array(vecotor)

j_x = 0

j_y = 1

上の微分方程式dx/dt, dy/dtの右辺に対応する関数をそれぞれ用意する．

#definition of differential equation

def dxdt(in_vec, t):

    growth = r*in_vec[j_x]*(1.0 - in_vec[j_x]/K)

    mortality = a*in_vec[j_x]*in_vec[j_y]/(1.0 + H_t*a*in_vec[j_x])

    return growth-mortality

def dydt(in_vec, t):

    growth = b*in_vec[j_x]*in_vec[j_y]/(1.0 + H_t*a*in_vec[j_x])

    mortality = m*in_vec[j_y]

    return growth-mortality

さらに微分方程式の右辺の値をベクトル値in_vec, 時刻tですべて計算してベクトルout_vecに格納する関数を用意する．

#function to calculate all coefficients, two dimensional

def diff2(in_vec, out_vec, t):

    out_vec[j_x] = dxdt(in_vec, t)

    out_vec[j_y] = dydt(in_vec, t)

最後に，モデル変数x,yを格納するための二次元配列nvを用意し，初期値を入れ，あとはロジスティック成長モデルのときと同様にfor_loopによって逐次的に微分方程式の数値解を求める．

nv = np.zeros([2])

nv[j_x] = 0.1  # initial condition (initial population size)

nv[j_y] = 0.1  # initial condition (initial population size)

time = 0.0  # initial condition (initial time, 0)

write_index = 1  # need for counting the numerical step

write_interval = 0.1  # with this interval, the result will be saved

h = 0.01

end_time = 100.0

# write initial condition

text = str(time) + ',' + str(nv[j_x]) + ',' + str(nv[j_y]) + '\n'

with open('result_ode4.csv', 'w') as f:

    f.write(text)

for i in range(0, int(end_time / h)):

    rk4(nv, nv, time, h, 2, diff2) #update vector value

    time += h  # update time

    ###recored the result every "write_interval" only

    if (write_index < int(write_interval / h)):

        write_index += 1

    else:

        text = str(time) + ',' + str(nv[j_x]) + ',' + str(nv[j_y]) + '\n'

        with open('result_ode4.csv', 'a') as f:

            f.write(text)

        write_index = 1

グラフは横軸を時間にするなら以下のようなスクリプトで描ける．

# load result as a dataframe

table = pd.read_csv('result_ode4.csv', header=None)

plt.plot(table.iloc[:, 0], table.iloc[:, j_x+1], color='blue', linestyle='solid')

plt.plot(table.iloc[:, 0], table.iloc[:, j_y+1], color='red', linestyle='solid')

plt.xlabel('time')

plt.ylabel('x, y')

plt.legend(['prey x', 'predator y'])

横軸にx，縦軸にyをプロットするなら以下のような感じである．ちなみにK=1.0, K=1.2, K=2.0の三つの場合を重ね合わせた．これが有名な"Paradox of enrichment"である．

plt.close()

table = pd.read_csv('result_ode4.csv', header=None)

table2 = pd.read_csv('result_ode5.csv', header=None)

table3 = pd.read_csv('result_ode6.csv', header=None)

plt.plot(table.iloc[:, j_x+1], table.iloc[:, j_y+1], color='blue', linestyle='solid')

plt.plot(table2.iloc[:, j_x+1], table2.iloc[:, j_y+1], color='green', linestyle='solid')

plt.plot(table3.iloc[:, j_x+1], table3.iloc[:, j_y+1], color='red', linestyle='solid')

plt.xlabel('prey x')

plt.ylabel('predator y')

plt.legend(['K = 1.0','K = 1.2','K = 2.0'])

[4] C, Rとの速度比較

ここではC言語及びR言語との速度比較をしてみる．まず，ターミナルから実行するための独立したpythonファイルtest_ode2.pyを用意する. 結果を出力する関数（print）の速度で律速されてしまうと良い比較にならないので，t=0での初期値とt=1000での最終値だけ出力する．

import numpy as np   #for numerical calculations

import pandas as pd  #read with pandas as dataframe

import math    #for mathematical functions

import time

#import matplotlib.pyplot as plt   #for graphics

#from numpy import exp

#from numpy import sin

#from numpy import pi

start = time.time()

def rk4(in_vec, out_vec, time, h_interval, dim, diff_vec):

    #in_vec: the vector with the value at "time"

    #out_vec: the vector with which the updated vector value after one time step is saved

    #time: time value

    #h_interval: time step for discritizing ODE

    #dim: dimension of the ODE

    #diff_vec: function that determines the r.h.s. of the ODE

    k1 = np.zeros([dim])

    k2 = np.zeros([dim])

    k3 = np.zeros([dim])

    k4 = np.zeros([dim])

    temp_vec = np.zeros([dim])

    h_half = h_interval / 2.0

    t_h = time + h_half

    diff_vec(in_vec, k1, time)  # calculate k1

    temp_vec = in_vec + h_half * k1

    diff_vec(temp_vec, k2, t_h)  # calcualte k2

    temp_vec = in_vec + h_half * k2

    diff_vec(temp_vec, k3, t_h)  # calcualte k3

    temp_vec = in_vec + h_interval * k3

    diff_vec(temp_vec, k4, time + h_interval)  # calcualte k4

    for i in range(0, dim):

        out_vec[i] = in_vec[i] + (h_interval / 6.0) * (k1[i] + 2.0 * k2[i] + 2.0 * k3[i] + k4[i])

#define parameters for ODE for prey predator model

r = 1.0

a = 1.0

b = 0.5

H_t = 1.0

m = 0.1

K  = 1.0

#index to link the model variables (x, y) and the elements of the array(vecotor)

j_x = 0

j_y = 1

#definition of differential equation

def dxdt(in_vec, t):

    growth = r*in_vec[j_x]*(1.0 - in_vec[j_x]/K)

    mortality = a*in_vec[j_x]*in_vec[j_y]/(1.0 + H_t*a*in_vec[j_x])

    return growth-mortality

def dydt(in_vec, t):

    growth = b*in_vec[j_x]*in_vec[j_y]/(1.0 + H_t*a*in_vec[j_x])

    mortality = m*in_vec[j_y]

    return growth-mortality

#function to calculate all coefficients, two dimensional

def diff2(in_vec, out_vec, t):

    out_vec[j_x] = dxdt(in_vec, t)

    out_vec[j_y] = dydt(in_vec, t)

nv = np.zeros([2])

nv[j_x] = 0.1  # initial condition (initial population size)

nv[j_y] = 0.1  # initial condition (initial population size)

t = 0.0  # initial condition (initial time, 0)

write_index = 1  # need for counting the numerical step

write_interval = 0.1  # with this interval, the result will be saved

h = 0.01

end_time = 1000.0

K = 1.0

# write initial condition

print(t, nv[j_x], nv[j_y])

for i in range(0, int(end_time / h)):

    rk4(nv, nv, t, h, 2, diff2) #update vector value

    t += h  # update time

    ###recored the result every "write_interval" only

    if (write_index < int(write_interval / h)):

      write_index += 1

    else:

      #print(time, nv[j_x], nv[j_y])

      write_index = 1

print(t, nv[j_x], nv[j_y])

end = time.time()

print("elapsed_time:{0}".format(end - start) + "[sec]")

つづいて，R言語のスクリプトwith_r_ode.Rを用意する．

start<-proc.time()

rk4<-function(in_vec, time, h_interval, dim, diff_vec)

  #in_vec: the vector with the value at "time"

  #out_vec: the vector with which the updated vector value after one time step is saved

  #time: time value

  #h_interval: time step for discritizing ODE

  #dim: dimension of the ODE

  #diff_vec: function that determines the r.h.s. of the ODE

  k1 <- numeric(dim)

  k2 <- numeric(dim)

  k3 <- numeric(dim)

  k4 <- numeric(dim)

  temp_vec <- numeric(dim)

  h_half <- h_interval / 2.0

  t_h <- time + h_half

  k1 <- diff_vec(in_vec, time, dim)  # calculate k1

  temp_vec <- in_vec + h_half * k1

  k2 <- diff_vec(temp_vec, t_h, dim)  # calcualte k2

  temp_vec <- in_vec + h_half * k2

  k3 <- diff_vec(temp_vec, t_h, dim)  # calcualte k3

  temp_vec <- in_vec + h_interval * k3

  k4 <- diff_vec(temp_vec, time + h_interval, dim)  # calcualte k4

  temp_vec <- in_vec + (h_interval / 6.0) * (k1 + 2.0 * k2 + 2.0 * k3 + k4)

  temp_vec

#define parameters for ODE for prey predator model

r <- 1.0

a <- 1.0

b <- 0.5

H_t <- 1.0

m <- 0.1

K  <- 1.0

#index to link the model variables (x, y) and the elements of the array(vecotor)

j_x <- 1

j_y <- 2

#definition of differential equation

dxdt <-function(in_vec, t)

  growth <- r*in_vec[j_x]*(1.0 - in_vec[j_x]/K)

  mortality <- a*in_vec[j_x]*in_vec[j_y]/(1.0 + H_t*a*in_vec[j_x])

  return(growth-mortality)

dydt <-function(in_vec, t)

  growth <- b*in_vec[j_x]*in_vec[j_y]/(1.0 + H_t*a*in_vec[j_x])

  mortality <- m*in_vec[j_y]

  return(growth-mortality)

#function to calculate all coefficients, two dimensional

diff2 <-function(in_vec, t, dim){

  temp_vec <- numeric(dim)

  temp_vec[j_x] <- dxdt(in_vec, t)

  temp_vec[j_y] <- dydt(in_vec, t)

  temp_vec

nv <- numeric(2)

nv[j_x] <- 0.1  # initial condition (initial population size)

nv[j_y] <- 0.1  # initial condition (initial population size)

time <- 0.0  # initial condition (initial time, 0)

write_index <- 1  # need for counting the numerical step

write_interval <- 0.1  # with this interval, the result will be saved

h <- 0.01

end_time <- 1000.0

K <- 1.0

# write initial condition

cat(time, nv, "\n")

for(i in 1:as.integer(end_time / h)){

  nv <- rk4(nv, time, h, 2, diff2)

  time <- time + h # update time

  ###recored the result every "write_interval" only

  if(write_index < as.integer(write_interval/h)){

    write_index <- write_index + 1

  else {

    #cat(time, nv,"\n")

    write_index <- 1

cat(time, nv, "\n")

proc.time()-start

最後にC言語のソースコードtest_ode.cを用意する．

//this is the model that considered the temperature-dependent movement only.

//from 20170418

#include <stdio.h>

#include <math.h>

#include <stdlib.h>

#include <time.h>

void rk4(double yin[], int n, double t, double h, double yout[], void (*derivs) (double, double [], double []));

void differential(double time, double in[], double out[]);

double dxdt(double time, double vr[]);

double dydt(double time, double vr[]);

int j_x, j_y;

double r = 1.0;

double a = 1.0;

double b = 0.5;

double H_t = 1.0;

double m = 0.1;

double K = 1.0;

int main(int arg, char *argv[])

  clock_t start,end;

  start = clock();

  double t, tzero, end_time, deltat, write_interval;

  double v[2];

  long n_variable = 2;

  long i, j, k;

  long write_index;

  tzero = 0.0;

  end_time = 1000.0;

  deltat = 0.01;

  write_interval = 0.1;

  j_x = 0;

  j_y = 1;

  write_index = 1;

  //initial conditions

  v[j_x] = 0.1;

  v[j_y] = 0.1;

  printf("%lf\t%.10lf\t%.10lf\n", t, v[j_x], v[j_y]);

  for(t = tzero; t <= end_time;) {

    rk4(v, n_variable, t, deltat, v, differential);

    t += deltat;

    if(write_index < write_interval/deltat){

      write_index += 1;

    else {

      //printf("%lf\t%.10lf\t%.10lf\n", t, v[j_x], v[j_y]);

      write_index = 1;

  }//end of for t

  printf("%lf\t%.10lf\t%.10lf\n", t, v[j_x], v[j_y]);

  end = clock();

  printf("The calculation time was:%.6f [sec]\n",(double)(end-start)/CLOCKS_PER_SEC);

  return 0;

void differential(double time, double in[], double out[])

  out[j_x] = dxdt(time, in);

  out[j_y] = dydt(time, in);

double dxdt(double time, double vr[])

  double temp;

  temp = r*vr[j_x]*(1.0 - vr[j_x]/K);

  temp -= a*vr[j_x]*vr[j_y]/(1.0 + vr[j_x]);

  return temp;

double dydt(double time, double vr[])

  double temp;

  temp =  b*vr[j_x]*vr[j_y]/(1.0 + H_t*a*vr[j_x]);

  temp -= m*vr[j_y];

  return temp;

void rk4(double yin[], int n, double t, double h, double yout[], void (*derivs) (double, double [], double []))

  int i;

  double th, hh, h6;

  double k1[n], k2[n], k3[n], k4[n];

  double temp_y[n];

  hh = h*0.5;

  h6 = h/6.0;

  th = t + hh;

  (*derivs)(t, yin, k1);

  for (i=0;i<n;i++) temp_y[i] = yin[i] + hh*k1[i];

  (*derivs)(th, temp_y, k2);

  for (i=0;i<n;i++) temp_y[i] = yin[i] + hh*k2[i];

  (*derivs)(th, temp_y, k3);

  for (i=0;i<n;i++) temp_y[i] = yin[i] + h*k3[i];

  (*derivs)(t+h, temp_y, k4);

  for (i=0;i<n;i++) yout[i] = yin[i] + h6*(k1[i] + 2.0*k2[i] + 2.0*k3[i] + k4[i]);

これらのコードを使って速度比較してみる.

icc -fastでコンパイルした場合のC：0.0067 sec

gcc -O2でコンパイルした場合のC：0.0088 sec (iccの1.3倍)

Rscript with_r_ode.R: 4.3 sec (iccの642倍, gccの489倍)

python3 test_ode2.py: 3.0 sec (iccの448倍, gccの341倍)

RほどではないにせよPython遅すぎ．

[5] Pythonの高速化

高速化参照リンク：

https://qiita.com/Kuma_T/items/f9572682fd4fac700c0b

いろいろ科学計算に便利なライブラリを一括したanacondaディストリビューションをインストールする.

https://www.anaconda.com/download/

Win, MacではインストーラをクリックすればOK. Linuxではインストーラをダウンロード後，以下のようにすればよい．

sudo bash ./Anaconda3-5.1.0-Linux-x86_64.sh

ここでは@numba.jitという仕組みを使う．pythonではforループが遅いので高速化したいループを含む計算を関数として定義し直して@numba.jitをその前に追加すればよい(test_ode4.py)．単純に１から順に自然数の和を求めるようなスクリプトは以下のように書ける．

import numpy as np   #for numerical calculations

import pandas as pd  #read with pandas as dataframe

import math    #for mathematical functions

import time

import numba

#from numba.decorators import jit

#import matplotlib.pyplot as plt   #for graphics

#from numpy import exp

#from numpy import sin

#from numpy import pi

#@numba.jit #for fast calclations

#start = time.time()

@numba.jit

def test_loop(n):

  sum_num = 0;

  for i in range(n):

    sum_num += i

  return sum_num

start = time.time()

k = 100000000

test_loop(k)

end = time.time()

print("elapsed_time:{0}".format(end - start) + "[sec]")

def test_loop(n)の前の@numba.jitをコメントアウトした場合の実行時間は6.2 sec，コメントアウトしない場合は0.13 secと50倍近い高速化が実現する．このような単純なループでは異なるiの間の計算に相互依存性がないからこそ高速化できる．しかし問題は，微分方程式の逐次解法におけるforループのように一つ前のiに依存して回るようなforループは高速化ができないか非常に難しい．

ただし非常に多くの種がいる場合とか空間差分化した反応拡散方程式のような場合には、常微分方程式の次元が上がるだけなのでこの部分の繰り返しは高速化できるはず. しかしいろいろやってもうまくいかないのでとりあえず保留．いろいろやってみたサンプルファイルはこれ（test_ode2.1.py, test_ode3.py）．

さらに重要なのは型指定。http://yutori-datascience.hatenablog.com/entry/2014/12/10/123157

たとえば、上のコードで、@numba.jitの代わりに@numba.jit('i8(i8)')と爆速になる．

[6] クラス、メソッドとしてオイラー法とルンゲクッタ法を実装する

オイラー法、ルンゲクッタ法の部分をクラス，メソッド化してみる．numeric_odeという名前のクラスを定義しその中でパラメータのいくつかを自動的に初期化するためのinitメソッド(コンストラクタ）を定義し，オイラー法およびルンゲクッタ法の一ステップ分の計算をしてくれる関数をそれぞれee，rk4という名前のメソッドで定義する．普通の関数との違いは一つ目の引数は必ずselfとすること．あと、アンダーバー二つ(__)で始まるものはprivate変数．

#definition of class for the numerical integration, including Explicit Euler and 4th-order RK

class numeric_ode:

    def __init__(self, h_interval, dim):

        self.h = h_interval

        self.dim = dim

    # definition of explicit Euler

    def ee(self, in_vec, time, diff_vec):

        self.__k1 = np.zeros([self.dim])

        self.__out_vec = np.zeros([self.dim])

        # calculate the next step vector

        diff_vec(in_vec, self.__k1, time)  # calculate k1

        self.__out_vec = in_vec + self.h * self.__k1  #update vector

        return self.__out_vec

    def rk4(self, in_vec, time, diff_vec):

        self.__k1 = np.zeros([self.dim])

        self.__k2 = np.zeros([self.dim])

        self.__k3 = np.zeros([self.dim])

        self.__k4 = np.zeros([self.dim])

        self.__temp_vec = np.zeros([self.dim])

        self.__h_half = self.h/2.0

        self.__t_h = time + self.__h_half

        diff_vec(in_vec, self.__k1, time)  # calculate k1

        self.__temp_vec = in_vec + self.__h_half * self.__k1

        diff_vec(self.__temp_vec, self.__k2, self.__t_h)  # calcualte k2

        self.__temp_vec = in_vec + self.__h_half * self.__k2

        diff_vec(self.__temp_vec, self.__k3, self.__t_h)  # calcualte k3

        self.__temp_vec = in_vec + self.h * self.__k3

        diff_vec(self.__temp_vec, self.__k4, time + self.h)  # calcualte k4

        self.__out_vec = in_vec + (self.h / 6.0) * (self.__k1 + 2.0 * self.__k2 + 2.0 * self.__k3 + self.__k4)

        return self.__out_vec

上の例で関数(メソッド）は関数を引数にすることが可能なことに注意．ここではeeおよびrk4の最後の引数diff_vecは関数である．

rk4の方を例に，使い方は以下の感じ．

まずは，これまでの例のようにprey-predator modelの方程式を以下のように定義する．

#define parameters for ODE for prey predator model

r = 1.0

a = 1.0

b = 0.5

H_t = 1.0

m = 0.1

K  = 1.0

#index to link the model variables (x, y) and the elements of the array(vecotor)

j_x = 0

j_y = 1

#definition of differential equation

def dxdt(in_vec, t):

    growth = r*in_vec[j_x]*(1.0 - in_vec[j_x]/K)

    mortality = a*in_vec[j_x]*in_vec[j_y]/(1.0 + H_t*a*in_vec[j_x])

    return growth-mortality

def dydt(in_vec, t):

    growth = b*in_vec[j_x]*in_vec[j_y]/(1.0 + H_t*a*in_vec[j_x])

    mortality = m*in_vec[j_y]

    return growth-mortality

#function to calculate all coefficients, two dimensional

def diff2(in_vec, out_vec, t):

    out_vec[j_x] = dxdt(in_vec, t)

    out_vec[j_y] = dydt(in_vec, t)

そして初期条件の設定，ファイルへの書き込みをやっておく．

nv = np.zeros([2])

nv[j_x] = 0.1  # initial condition (initial population size)

nv[j_y] = 0.1  # initial condition (initial population size)

t = 0.0  # initial condition (initial time, 0)

write_index = 1  # need for counting the numerical step

write_interval = 0.1  # with this interval, the result will be saved

h = 0.01

end_time = 100.0

K = 1.0

#write initial condition

text=str(t)+','+str(nv[j_x])+','+str(nv[j_y])+'\n'

with open('result_ode5.csv', 'w') as f:

    f.write(text)

ここで，クラスを実体化させたもの（＝インスタンス）をパラメータを与えて「生成」する．ここではnum_odeという名前のインスタンスを生成する．

#generate an instance

num_ode = numeric_ode(dim=2, h_interval=h)

次にfor loopの中でrk4メソッドを呼び出して使う。

for i in range(0, int(end_time / h)):

  nv = num_ode.rk4(nv, t, diff2) #call the rk4 method from the instance

  t += h  # update time

  ###recored the result every "write_interval" only

  if (write_index < int(write_interval / h)):

    write_index += 1

  else:

      text = str(t) + ',' + str(nv[j_x]) + ',' + str(nv[j_y]) + '\n'

      with open('result_ode5.csv', 'a') as f:

          f.write(text)

      write_index = 1

結果が書き込まれているのでグラフ化するには以下の通り（グラフは省略）．

#load result as a dataframe

table=pd.read_csv('result_ode5.csv', header=None)

plt.plot(table.iloc[:, 0], table.iloc[:, j_x+1], color='blue', linestyle='solid')

plt.plot(table.iloc[:, 0], table.iloc[:, j_y+1], color='red', linestyle='solid')

plt.xlabel('time')

plt.ylabel('x, y')

plt.legend(['prey x', 'predator y'])

以上、常微分方程式モデルについての解説でした．

Google Sites

Report abuse