Volterra モデル: 連立常微分方程式の例

Volterra の被食者-捕食者モデル方程式を計算してみよう

これまで単独のスカラー値常微分方程式を対象としてきたが、ここで、連立常微分方程式を対象として数値計算に取り組んでみよう．

tips: その取扱いをみてみると、実は、スカラー値常微分方程式 → 連立常微分方程式 → 偏微分方程式 という流れでつながることが「後で」分かる．今はまあそういうもうのか、と思ってこの流れについていこう．

さて、以前勉強したように、Volterra の被食者-捕食者モデル方程式は、時間変数を $t$, 被食者数を関数 $u(t)$ で、捕食者数を関数 $v(t)$ で表すとして以下のような非線形の常微分微分方程式で表されるものだった．

$$ \left\{\begin{array}{rcl} \displaystyle \frac{du}{dt} & = & - C_1 \, uv + D_1 \, u, \cr \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 \cr \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 & \left(1 + \Delta t D_1 - \Delta t C_1 v_n \right)\, u_n \cr % v_{n+1} & = % & v_n + \Delta t C_2 u_n v_n - \Delta t D_2 v_n & \left(1 - \Delta t D_2 + \Delta t C_2 u_n \right)\, v_n \end{array}\right. $$

というわけで Euler 法による近似式はこれでよさそうだ．

次に、用いる定数だが、とりあえず $C_1 = $ $C_2 = $ $D_1 = $ $D_2 = 1.0$, $u(0) = v(0) = 0.7$, としておこう． あとは、近似のために発生したパラメータ $\Delta t$ だが、まあ最初は $\Delta t = 0.05$ ぐらいでやってみよう．

というわけで、まず定数を設定して、

1
2
3

C1 = C2 = D1 = D2 = 1.0
u0 = v0 = 0.7
Δt = 0.05

次に、Euler スキームによる時間 1ステップ分の発展を関数として定義する．

1
2
3
4
5
6

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

euler (generic function with 1 method)

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

1

euler(u0,v0)

(0.7105, 0.6895)

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

1
2
3
4
5
6
7

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

結果がきちんと出ているか、単純に数字をまずちょっと眺めてみよう．

1

uv_sq

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

うん、数字としては出ているな．では、グラフにしてみよう

1
2

using Plots
gr()

Plots.GRBackend()

1
2
3
4

# plot 命令に行列を与えると、縦ベクトルを別々のグラフとして描くことを利用している．
#  ' をつけて行列を転置していることにも注意．

plot(uv_sq') 

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

横軸を、配列の要素番号ではなく、本来の問題中の時間変数 $t$ の表記にきちんとした方が良いかな？ では例えば次のようにしよう．

1
2

t_sq = [0:500] ＊ Δt
plot(t_sq, uv_sq', xaxis = "time t", yaxis = "u and v", leg = false)

さて、ふ〜む、だんだん振幅が大きくなっているけどこの挙動は信用して良いものだろうか？ ちょっと $\Delta t$ を小さくしてみて、変化したりしないかチェックしてみよう．

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

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13

lastT = 25   # 計算を終える時間 $t\_{last}$ を設定して、

Δt = 0.01    # Δt を小さめに設定し直して、
N = Int(lastT/Δt)   # 最終時間ステップは整数なのでこんな感じで

# あとは、上でやった計算そっくりに. ただし, 計算ステップ数は N に直して．
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

として計算して、グラフを早速描いてみる．

1
2

t_sq = [0:N] * Δt
plot(t_sq, uv_sq', xaxis = "time t", yaxis = "u and v", leg = false)

んん？ 随分変わってしまったな．振幅が余り変わらなくなった． どうやらこの問題でも以前同様、Euler スキームは $\Delta t$ が大きい時は少し信用ならないようだ． あとでもう少し信用できそうな方法で計算し直そう(レポート問題. このページの最後尾を見よう).

相図をみてみよう

今確認したように Euler 法の結果はあまり信用ならないけれども、$\Delta t$ が小さいならばまあまったく外れというわけでもなさそうなので、まずはこの Euler 法の結果を使って相図を作り、様子を見てみよう．なにか見えてくるものがあるかもしれない．

1
2
3
4
5
6
7
8

u_sq = uv_sq[1,:] # とってあるデータの、1行目. これが u の近似値の列のはず．
v_sq = uv_sq[2,:] # とってあるデータの、2行目. これが v の近似値の列のはず．

plot(u_sq, v_sq, aspect_ratio = 1, xaxis = "u", yaxis = "v")
# aspect_ratio は、縦横の軸目盛りの比率を指定する．1 にすると、縦横の縮尺が一致する．

plot!((u0,v0), marker = :circle, lab = "initial value")
# lab に指定した文字列は、legend に現れる

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

蛇足: 相図上の軌道の形

ちなみに、この「ぐるぐる周っている軌道」の形、つまり、軌道を数式で表す方法などは理論的にわかるものなのだろうか？ という疑問が生ずる． 実は、一般に、相図上の軌道を式で表すことは大変難しい問題で、「初等関数では表せないことが証明されている」ようなケースも有るぐらいだ．

そしてさらに、この Volterra モデルの場合は「大変ラッキーなことに」この軌道を数式できちんと書けるのだ！ 詳しいことは次回以降に楽しみに．

漁業の影響を見てみよう

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

まず、モデルの数式は、漁業の影響係数 $E > 0$ を用いて、漁業「有り」だと次のようになると考えるとよさそうだ、と Volterra さんは考えた．

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

そこで、この数式をやはり近似して、計算してみよう．

なお、この係数 $E$ はとりあえず $E = 0.5$ としておこう．

Euler 法で…

先の 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, \cr \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, \cr % 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. $$

実際にその 1ステップ分を関数としてプログラムすると以下のような感じになるだろう．

1
2
3
4
5
6
7
8

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

euler_fishing (generic function with 1 method)

実際に動かしてみよう．

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12

lastT = 25

Δt = 0.01
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

無事に動いたら、これを早速グラフで見てみよう．

1
2

t_sq = [0:N] * Δt
plot(t_sq, uv_f_sq', xaxis = "time t", yaxis = "u and v", leg = false)

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

1
2
3
4
5

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

plot(u_f_sq, v_f_sq, aspect_ratio = 1, xaxis = "u", yaxis ="v")
plot!((u0,v0), marker = :circle, lab = "initial value")

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

1
2
3
4
5
6
7
8

# まず、漁業効果無しのケースをプロットして、
plot(u_sq, v_sq, aspect_ratio = 1, xaxis = "u", yaxis = "v", label = "without fishing")

# 次に plot! 命令で、漁業効果有りのケースを描き足して、
plot!(u_f_sq, v_f_sq, label = "with fishing")

# ついでに初期値も描き足しておこうか．
plot!((u0,v0), marker = :circle, lab = "initial value")

漁業が行われると、相図の上での軌跡の中心が右下にずれるようだ(ついでに、軌跡も大きくなる)、ということがこれから読み取れる． これは、モデルの説明のときにもちょっと示したように、

「漁業を行うと魚の量がかえって増える現象」をうまく説明できている

ようにみえ、Volterra モデルの妥当性を示していると言える．

さて、じつは、もとの Volterra モデルの常微分方程式をじっと眺めると、この「軌跡の中心」とでもいうべき点は幸いにも手計算で求めることが出来ることがわかる．
チャレンジ: 実際にその「中心点」相当をモデルの常微分方程式から導出してみよう．

この中心点相当は、それぞれのケースで、

漁業なし… 中心 $\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)$

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

1
2
3
4
5
6
7

plot(u_sq, v_sq, aspect_ratio = 1, xaxis = "u", yaxis = "v", label = "without fishing")
plot!((D2/C2,D1/C1), marker = :circle, lab = "center of ...")

plot!(u_f_sq, v_f_sq, label = "with fishing")
plot!(((D2 + E)/C2,(D1 - E)/C1), marker = :circle, lab = "center of ...")

plot!((u0,v0), marker = :star, lab = "initial value")

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

リアルタイムに漁業の影響の大きさをいろいろ変えて、相図の様子の変化をみてみよう

グラフの描き方のところでも少し紹介してあるように、Julia の Jupyter 環境などでは Interact package を使うとパラメータをマウスでコントロールしながらリアルタイムに処理を行わせることができる． これを使って、漁業の影響を実際にみてみよう．

まずは Package の使用宣言だ．ちなみに，JuliaBox はこのパッケージはデフォルトでインストールされているので，自分でインストールしなくても良い(しようとすると失敗する)．

1

using Interact

あと、漁業効果の係数 E の影響を明示的にプログラムに反映させておこう．

1
2
3
4
5
6

function euler_fishing(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

euler_fishing (generic function with 2 methods)

あとはこいつを使って、グラフの例と同等に Interact package を使えば良い．

Interact の使い方だが、@manipulate マクロに続けて for 文を書けば良い． 通常だと for 文はその中身がループするが，@manipulate がついているときはループせずに，for で指定されている数量などが画面で制御できるようになる．

まあ，例で見てみたほうがわかりやすいだろうから，実際の以下の使用例を見てみよう．

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

@manipulate for E in 0.0:0.05:0.8

    # ここから end まではこれまでと本質的に同じ．
    u,v = u0, v0
    uv_f_sq = [ u0, v0 ]
    
    for i in 1:N
        u,v = euler_fishing(u,v,E)   # 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, xaxis = "u", yaxis = "v", label = "without fishing")
    plot!((D2/C2,D1/C1), marker = :circle, lab = "center of ...")

    plot!(u_f_sq, v_f_sq, label = "with fishing")
    plot!(((D2 + E)/C2,(D1 - E)/C1), marker = :circle, lab = "center of ...")

    plot!((u0,v0), marker = :star, lab = "initial value")

end

漁業効果の係数 E をマウスで変えられるスライダーが出てくるはずなので、それを掴んで動かして、いろいろ見てみよう． これによって、漁業の影響が多少ならば、魚も増えて良いことばかりのように見えるが、漁業の影響が強くなりすぎるといろいろな点で危ないことなども見えてくるだろう．

さて、@manipulate のところを @gif と変えるだけで、今度はアニメーションを作ることが出来る．やってみよう．

ただし、1,2分待たされるのでじっと待とう．

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20

@gif for E in 0.0:0.05:0.8 # 上の例と、@gif という部分が違うだけ！
    u,v = u0, v0
    uv_f_sq = [ u0, v0 ]
    
    for i in 1:N
        u,v = euler_fishing(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, xaxis = "u", yaxis = "v", label = "without fishing")
    plot!((D2/C2,D1/C1), marker = :circle, lab = "center of ...")

    plot!(u_f_sq, v_f_sq, label = "with fishing")
    plot!(((D2 + E)/C2,(D1 - E)/C1), marker = :circle, lab = "center of ...")

    plot!((u0,v0), marker = :star, lab = "initial value")
end

シチュエーションによっては、こうしたアニメーションのほうが様子を上手く掴めたりするので、上手に使っていくとよいだろう．

Report No.05: 漁業効果あり、漁業効果なし、いずれの場合も、時間方向に対称なスキームを自分で作って動かしてみよう

tips: この問題で時間対称なスキームを作ると、時間ステップを 1ステップ進めるのに連立非線形方程式を解かないといけない形になるはずだ． その解法にもしも Newton 法を使おうとすると、スカラー値常微分方程式よりも少しばかり手強いかもしれない． そうした場合は、Newton 法よりも NLsolve package を使ったほうが楽だろうね．

Report No.06: 漁業効果あり、漁業効果なし、いずれの場合も、Rugne-Kutta 法で数値スキームを作って動かしてみよう

やってみるとわかるが、時間方向に対称なスキームを作ることに比べると格段に容易だ．実感しよう．