Volterra の、被食者-捕食者関係のモデル方程式を扱ってみよう

これは、被食者数を $u(t)$ で、捕食者数を $v(t)$ で表すとして、 以下のような非線形の常微分微分方程式で表される．

$\left\{\begin{array}{rcl} \displaystyle \frac{du}{dt} & = & - C_1 uv + D_1 u \\ \displaystyle \frac{dv}{dt} & = & C_2 uv - D_2 v \end{array}\right.$

ただし、係数 $C_1$, $C_2$, $D_1$, $D_2$ は全て定数で正としておこう．

相図等についてはまた話すとして、まずはこいつに従って計算してみよう．

まあまずは簡単な Euler 法で

離散時間の近似解を $u_n \cong u(n\Delta t)$, $v_n \cong v(n\Delta t)$, と書くとして、Euler スキームは次のようになるだろう．

$\left\{\begin{array}{rcl} \displaystyle \frac{u_{n+1} - u_n}{\Delta t} & = & - C_1 u_n v_n + D_1 u_n \\ \displaystyle \frac{v_{n+1} - v_n}{\Delta t} & = & C_2 u_n v_n - D_2 v_n \end{array}\right.$

これを分かりやすく書き直すと次のようになる．

$\left\{\begin{array}{rcl} u_{n+1} & = & u_n - \Delta t C_1 u_n v_n + \Delta t D_1 u_n = (1 + \Delta t D_1 - \Delta t C_1 v_n ) u_n \\ % v_{n+1} & = & v_n + \Delta t C_2 u_n v_n - \Delta t D_2 v_n = (1 - \Delta t D_2 + \Delta t C_2 u_n ) v_n \end{array}\right.$

定数などを、 $C_1 = $ $C_2 = $ $D_1 = $ $D_2 = 1.0$, $u(0) = v(0) = 0.7$, $\Delta t = 0.05$ としてやってみよう．

まず定数を設定して、

In [2]:
C1 = C2 = D1 = D2 = 1.0
u0 = v0 = 0.7
Δt = 0.05
Out[2]:
0.05

Euler スキームを関数として定義する．

In [3]:
function euler(u,v)
    u_new = (1 + Δt*D1 - Δt*C1*v) * u
    v_new = (1 - Δt*D2 + Δt*C2*u) * v
    u,v = u_new, v_new
    return u, v
end
Out[3]:
euler (generic function with 1 method)

初期値を試しに入れて、動作するかどうか確かめてみよう

In [4]:
euler(u0,v0)
Out[4]:
(0.7105, 0.6895)

とりあえず動きそうなので,思い切って 500ステップほど動かしてみようか．

In [5]:
u,v = u0, v0
uv_sq = [ u0, v0 ]

for i in 1:500
    u,v = euler(u,v)
    uv_sq = hcat(uv_sq, [u,v])
end

結果をまずちょっと眺めてみて、

In [6]:
uv_sq
Out[6]:
2×501 Array{Float64,2}:
 0.7  0.7105  0.721531  0.733092  0.745186  …  0.353772  0.353133  0.353086
 0.7  0.6895  0.679519  0.670058  0.661116     1.03614   1.00266   0.970234

グラフにしてみよう

In [7]:
using Plots
gr()
Out[7]:
Plots.GRBackend()
In [8]:
plot(uv_sq') 
# 行列を与えると、縦ベクトルを別々のグラフとして描いてくれる．ちなみに,' は行列の転置記号．
Out[8]:
<?xml version="1.0" encoding="utf-8"?> 100 200 300 400 500 0.5 1.0 1.5 2.0 y1 y2

青線が $u_n$ で、赤線が $v_n$ だね．

横軸を、配列の要素番号ではなくてもとの問題の時間 $t$ にあわせたければ、次のようにすれば良い．

In [9]:
t_sq = [0:Δt:500*Δt]
plot(t_sq, uv_sq')
Out[9]:
<?xml version="1.0" encoding="utf-8"?> 0 5 10 15 20 25 0.5 1.0 1.5 2.0 y1 y2

ふ〜む、だんだん振幅が大きくなっているけど本当かな？ ちょっと $\Delta t$ を小さくしてみて、変化したりしないかチェックしてみよう．

上のグラフと比べたいので、同じ時間($t=25.0$)まで計算することにしよう．すると次のような感じかな．

In [10]:
Δt = 0.01
lastT = 25
N = Int(lastT/Δt)

u,v = u0, v0
uv_sq = [ u0, v0 ]

for i in 1:N
    u,v = euler(u,v)
    uv_sq = hcat(uv_sq, [u,v])
end
In [11]:
t_sq = 0:Δt:lastT
plot(t_sq, uv_sq')
Out[11]:
<?xml version="1.0" encoding="utf-8"?> 0 5 10 15 20 25 0.75 1.00 1.25 1.50 y1 y2

んん？ 随分変わったぞ？ 振幅が余り変わらなくなった．どうやら Euler スキームはこの問題だと少し信用ならないようだな．あとでもう少し信用できそうな方法で計算し直そう．

相図をみてみよう

あまり信用ならないけれども、まあまったく外れというわけでもないだろうから、まずは Euler 法の結果を相図にして、様子を見てみよう．なにか見えてくるものがあるかもしれない．

In [12]:
u_sq = uv_sq[1,:]
v_sq = uv_sq[2,:]

plot(u_sq, v_sq, aspect_ratio = 1)
# aspect_ratio は、縦横の軸目盛りの比率を指定する．1 にすると、縦横の縮尺が一致する．
plot!((u0,v0), marker = :circle, lab = "initial value")
# lab に指定した文字列は、legend に現れる
Out[12]:
<?xml version="1.0" encoding="utf-8"?> 0.75 1.00 1.25 1.50 0.75 1.00 1.25 1.50 y1 initial value

どうやら、相図の上では、ぐるぐると周っている状態のようだ．円とはすこし違った歪んだ形のようではあるが．

時間対称なスキームを作ってみよう(レポート問題)

自力で作ってみよう！ そして、プログラムして、動かしてみよう！

注意として、この問題は非線形なので、時間対称なスキームを作ると当然ながらそれは陰的なスキームになる．つまり、非線形方程式を説かないと、時間ステップを進めることが出来ないのだ．

というわけで、スキームを作ったら、Newton 法を自力で用いるか、NLsolve パッケージを使うかしないといけないな．そういう風にプログラムをして、自分で動かしてみよう．

漁業の影響を見てみよう

漁業の影響がどう出るか、Volterra モデルで見てみよう． これがこのモデルの「キモ」であるので、注意深く考えながらやろう．

まず、モデルの数式は、漁業の影響係数 $E > 0$ を用いて、次のようになる．

$\left\{\begin{array}{rcl} \displaystyle \frac{du}{dt} & = & - C_1 uv + D_1 u - E u\\ \displaystyle \frac{dv}{dt} & = & C_2 uv - D_2 v - E v \end{array}\right.$

そこで、この係数 $E$ を $E = 0.5$ として以下、計算してみよう．

Euler 法で…

スキームを一部変えるだけでいいはずだ．次のようになるだろう．

$\left\{\begin{array}{rcl} \displaystyle \frac{u_{n+1} - u_n}{\Delta t} & = & - C_1 u_n v_n + D_1 u_n - E u_n\\ \displaystyle \frac{v_{n+1} - v_n}{\Delta t} & = & C_2 u_n v_n - D_2 v_n - E v_n \end{array}\right.$

これを分かりやすく書き直すと次のようになる．

$\left\{\begin{array}{rcl} u_{n+1} & = & u_n - \Delta t C_1 u_n v_n + \Delta t D_1 u_n - \Delta E u_n = (1 + \Delta t D_1 - \Delta t E - \Delta t C_1 v_n ) u_n \\ % v_{n+1} & = & v_n + \Delta t C_2 u_n v_n - \Delta t D_2 v_n - \Delta t E v_n = (1 - \Delta t D_2 - \Delta t E + \Delta t C_2 u_n ) v_n \end{array}\right.$

In [13]:
E = 0.5

function euler_fishing(u,v)
    u_new = (1 + Δt*D1 - Δt*E - Δt*C1*v) * u
    v_new = (1 - Δt*D2 - Δt*E + Δt*C2*u) * v
    u,v = u_new, v_new
    return u, v
end
Out[13]:
euler_fishing (generic function with 1 method)
In [14]:
Δt = 0.01
lastT = 25
N = Int(lastT/Δt)

u,v = u0, v0
uv_f_sq = [ u0, v0 ]

for i in 1:N
    u,v = euler_fishing(u,v)
    uv_f_sq = hcat(uv_f_sq, [u,v])
end
In [15]:
t_sq = 0:Δt:lastT
plot(t_sq, uv_f_sq')
Out[15]:
<?xml version="1.0" encoding="utf-8"?> 0 5 10 15 20 25 0.5 1.0 1.5 2.0 2.5 y1 y2

んんん？ これは随分変わったな．早速相図を描いてみよう．

In [16]:
u_f_sq = uv_f_sq[1,:]
v_f_sq = uv_f_sq[2,:]

plot(u_f_sq, v_f_sq, aspect_ratio = 1)
plot!((u0,v0), marker = :circle, lab = "initial value")
Out[16]:
<?xml version="1.0" encoding="utf-8"?> 1.0 1.5 2.0 2.5 3.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 y1 initial value

変化を見たいので、漁業効果の入っていないデータと一緒に描いてみよう．

In [17]:
plot(u_sq, v_sq, aspect_ratio = 1)
plot!(u_f_sq, v_f_sq)
plot!((u0,v0), marker = :circle, lab = "initial value")
Out[17]:
<?xml version="1.0" encoding="utf-8"?> 0.5 1.0 1.5 2.0 2.5 3.0 0.3 0.6 0.9 1.2 1.5 y1 y2 initial value

どうやら、漁業が行われると、相図の上での軌跡の中心が右下にずれるようだ．

じつはこの「軌跡の中心」は幸いにも手計算で求めることが出来て、

漁業なし… 中心 $\displaystyle = \left(\frac{D_2}{C_2}, \frac{D_1}{C_1}\right)$,

漁業あり… 中心 $\displaystyle = \left(\frac{D_2 + E}{C_2}, \frac{D_1 - E}{C_1}\right)$

ということがわかっている．そこでこれも相図に描き込んでみよう．

In [18]:
plot(u_sq, v_sq, aspect_ratio = 1)
plot!(u_f_sq, v_f_sq)
plot!((u0,v0), marker = :circle, lab = "initial value")
plot!((D2/C2,D1/C1), marker = :circle, lab = "center without fishing")
plot!(((D2 + E)/C2,(D1 - E)/C1), marker = :circle, lab = "center with fishing")
Out[18]:
<?xml version="1.0" encoding="utf-8"?> 0.5 1.0 1.5 2.0 2.5 3.0 0.3 0.6 0.9 1.2 1.5 y1 y2 initial value center without fishing center with fishing

なんとなく、いろいろと分かってきたようなきがするのではないだろうか．

漁業の影響の大きさをいろいろ変えてみてみよう

Interact パッケージを使えば良いわけだ．

In [19]:
using Interact

係数 E の影響を計算がきちんと受けるように、関数の引数に明示しておこう．

In [20]:
function euler_fishing_mod(u,v,E)
    u_new = (1 + Δt*D1 - Δt*E - Δt*C1*v) * u
    v_new = (1 - Δt*D2 - Δt*E + Δt*C2*u) * v
    u,v = u_new, v_new
    return u, v
end
Out[20]:
euler_fishing_mod (generic function with 1 method)

あとはこいつを使って、これまでと同等のことをやれば良い．

In [21]:
@manipulate for E in 0.0:0.01:0.99
    u,v = u0, v0
    uv_f_sq = [ u0, v0 ]
    
    for i in 1:N
        u,v = euler_fishing_mod(u,v,E)
        uv_f_sq = hcat(uv_f_sq, [u,v])
    end
    
    u_f_sq = uv_f_sq[1,:]
    v_f_sq = uv_f_sq[2,:]
    
    plot(u_sq, v_sq, aspect_ratio = 1)
    plot!(u_f_sq, v_f_sq)
    plot!((u0,v0), marker = :circle, lab = "initial value")
    plot!((D2/C2,D1/C1), marker = :circle, lab = "center without fishing")
    plot!(((D2 + E)/C2,(D1 - E)/C1), marker = :circle, lab = "center with fishing")
end
Out[21]:
<?xml version="1.0" encoding="utf-8"?> 0.5 1.0 1.5 2.0 2.5 0.3 0.6 0.9 1.2 1.5 y1 y2 initial value center without fishing center with fishing

漁業の影響が多少ならば、魚も増えて良いことばかりのように見えるが、漁業の影響が強くなりすぎるといろいろな点で危ないことが見えてくるだろう．

漁業効果ありの場合も、時間対称なスキームを自分で作って動かしてみよう(レポート問題)

もう難しくないだろう．チャレンジしてみよう．