1.4. SciPyを用いたPredator-Preyモデルのシミュレーション

本章の最後に,SciPyを用いた微分方程式の解法例として,predator-preyモデルのシミュレーションについて紹介します.

帯状流と乱流の相互作用は,捕食者—被食者(Predator-Prey)モデルで記述されることが知られており [PP] ,このモデルは1階の連立微分方程式の形をしています. SciPyパッケージのodeintモジュールを使うと1階の常微分方程式の数値解を簡単に得ることが出来ます [1] . odeintはLSODA(Livermore Solver for Ordinary Differential equations with Automatic switching for stiff and non-stiff problems)法を利用した汎用的な積分器ですが, 詳しくはODEPACK Fortran library [ODE] を参照して下さい.

まずは,ソースコードを見てみましょう.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
#!/usr/bin/env python
"""
    Sample Code

    Status
    ------
    Version 1.0

    Authour
    -------
    Shigeru Inagaki                                       
    Research Institute for Applied Mechanics 
    inagaki@riam.kyushu-u.ac.jp  
    
    Revision History
    ----------------
    [11-April-2017] Creation

    Copyright
    ---------
    2017 Shigeru Inagaki (inagaki@riam.kyushu-u.ac.jp)
    Released under the MIT, BSD, and GPL Licenses.

"""
import numpy as np
import scipy.integrate as desol
import matplotlib.pyplot as plt

def predator_prey(f, t, a, b, c, d):
    """ return left hand sides of ordinary differential equations

        model equation:
            dx/dt = ax - bxy
            dy/dt = cxy - dy

    f[0]    - x: Population of prey
    f[1]    - y: Population of predator
    t       - Time
    a,b,c,d - Control parameters
    """
    return [a*f[0]-b*f[0]*f[1], c*f[0]*f[1]-d*f[1]]


#model
#eq1 = r"$\frac{dx}{dt} = ax - bxy$"
#eq2 = r"$\frac{dy}{dt} = cxy - dy$"
eq1 = r"$dx/dt = ax - bxy$"
eq2 = r"$dy/dt = cxy - dy$"

#input parameters
a = 1.0
b = 1.0
c = 1.0
d = 1.0
header = r"$a={0:.1f}, b={1:.1f}, c={2:.1f}, d={3:.1f}$".format(a, b, c, d)

#initial condition
f0 = [1.0, 0.1]

#independent variable
nt = 1000
tmax = 30.0
dt = tmax / nt
t = dt * np.arange(nt)

f = desol.odeint(predator_prey, f0, t, args=(a,b,c,d))

#plot
prey = f[:,0]
predator = f[:,1]

fig = plt.figure()
ax = fig.add_axes([0.15, 0.1, 0.8, 0.8])
ax.plot(t, prey, color='r', label=r"$x$: prey")
ax.plot(t, predator, color='b', label=r"$y$: predator")
handles, labels = ax.get_legend_handles_labels()
ax.legend(handles, labels, loc='best')
ax.text(21, 2.7, eq1, fontsize=16)
ax.text(21, 2.5, eq2, fontsize=16)
ax.set_xlabel("Time")
ax.set_ylabel("Population")
ax.set_title(header)
plt.show()
../_images/predator_prey.png

プログラムの内容は以下のようになっています.

  1. 解析する関数(この場合predator_prey)を定義する
    • 第1引数fが微分方程式中の未知関数
    • 第2引数tが関数のパラメータ(時間に対応)
    • 第3 - 6引数a, b, c, dが定数
    • 戻り値がパラメータtにおけるdx/dt, dy/dtを与える
  2. 微分方程式の定数a, b, c, dを与える
  3. 微分方程式の初期値f0を与える
  4. 未知関数の解析範囲(時間)を与えるパラメータ列tを用意する
  5. 関数 SciPy.integrate.odeint に1 - 4を引数にして呼び出す
  6. 戻り値がパラメータtに対応する未知関数fの各値となる

帯状流とプラズマ乱流の相互作用を当てはめて考えてみると,乱流を餌として発生・成長する帯状流は捕食者の役割を,またプラズマ圧力勾配により発生する線形不安定性を源として成長する乱流は被食者の役割を果たします.

このようにPythonを用いることで,簡単にモデルの計算と可視化をすることができます. コーディングの時間を短縮し,試行錯誤に多くの時間を割けるのがPythonの利点でもありますので,みなさんもまずは簡単なプログラムを作成し,動作を確認してみて下さい.

[PP]小林すみれ 他:プラズマ・核融合学会誌 92, 211 (2016).
[ODE]http://people.sc.fsu.edu/~jburkardt/f77_src/odepack/odepack.html
[1]なお,高階の微分方程式でも,1階の微分方程式に変換することでodeintを用いて計算することができます.