07. 連立常微分方程式にチャレンジ
Photo by Marcelo Cidrack on Unsplash
Volterra の被食者-捕食者モデル方程式の計算にチャレンジ
前回まで単独のスカラー値常微分方程式を対象としてきたが、ここで、連立常微分方程式を対象として数値計算に取り組んでみよう.
ターゲットの取扱いをみてみると、実は、スカラー値常微分方程式 → 連立常微分方程式 → 偏微分方程式 という流れでつながることが「後で」分かる.今はまあそういうもうのか、と思ってこの流れについていこう.
さて、以前勉強したように、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} & = & \left(1 + \Delta t \left( D_1 - C_1 v_n \right) \right)\, u_n , \cr % v_{n+1} & = & \left(1 + \Delta t \left( - D_2 + C_2 u_n \right) \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$ ぐらいでやってみよう.
では、プログラムだ. まず定数を設定して、
|
|
次に、Euler スキームによる時間 1ステップ分の発展を関数として定義する.
|
|
早速、初期値を試しに入れて、動作するかどうか確かめてみよう
|
|
(0.7104999999999999, 0.6895)
とりあえず動きそうなので,思い切って 500ステップほど動かしてみようか.
|
|
結果がきちんと出ているか、単純に数字をまずちょっと眺めてみよう.
|
|
501×2 Array{Float64,2}:
0.7 0.7
0.7105 0.6895
0.721531 0.679519
0.733092 0.670058
…略…
0.355027 1.07067
0.353772 1.03614
0.353133 1.00266
0.353086 0.970234
うん、結果が数字として出ていて,エラーにはなっていないようだ. では、早速グラフで様子を見てみよう.
|
|
青線が $u_n$ (魚に相当)で、赤線が $v_n$ (サメに相当)だね.
横軸を、配列の要素番号ではなく、本来の問題中の時間変数 $t$ の表記にきちんとした方が良いかな? 後学のために,label もきちんとした図を一回作っておこう(面倒なのでこの後は単純なものに戻すけどね).
|
|
さて、計算結果の内容のチェックに戻ろう. 図をみるとだんだん振幅が大きくなっているけどこの挙動は信用して良いものだろうか?
ちょっと $\Delta t$ を小さくして、全体の様子が変化しないかチェックしてみよう.
上のグラフと比べたいので、同じ時間( $t=25.0$ )ぐらいまで計算することにしよう.すると次のような感じかな.
|
|
として計算して、グラフを早速描いてみる.
|
|
んん? 随分変わってしまったな.振幅が余り変わらなくなった. どうやらこの問題でも以前同様、Euler スキームは $\Delta t$ が大きい時は少し信用ならないようだ. あとでもう少し信用できそうな方法で計算し直そう(レポート問題. このページの最後尾を見よう).
相図をみてみよう
今確認したように Euler 法の結果はあまり信用ならないけれども、$\Delta t$ が小さいならばまあまったく外れというわけでもなさそうなので、まずはこの Euler 法の結果を使って相図を作り、様子を見てみよう.なにか見えてくるものがあるかもしれない.
|
|
どうやら、相図の上では、ぐるぐると周っている状態のようだ.円とはすこし違った歪んだ形のようではあるが.
相図上の軌道の形
ちなみに、この「ぐるぐる周っている軌道」の形、つまり、軌道を数式で表す方法などは理論的にわかるものなのだろうか? という疑問が生ずる. 実は、一般に、相図上の軌道を式で表すことは大変難しい問題で、「初等関数では表せないことが証明されている」ようなケースも有るぐらいだ.
そしてさらに、この 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} & = & \left(1 + \Delta t \left( D_1 - E - C_1 v_n \right) \right)\, u_n, \cr % v_{n+1} & = & \left(1 + \Delta t \left( - D_2 - E + C_2 u_n \right) \right)\, v_n . \end{array}\right. \]
実際にその 1ステップ分を関数としてプログラムすると以下のような感じになるだろう.
|
|
実際に動かしてみよう.
|
|
無事に動いたら、これを早速グラフで見てみよう.
|
|
んんん? これは随分変わったな.早速相図を描いてみよう.
|
|
だいぶ変わった…が,比べてみないとよくわからないな.
変化を見たいので、漁業効果の入っていないデータと一緒に描いてみよう.
|
|
以前の授業で解説したように,漁業が行われると、相図の上での軌跡の中心が右下にずれる ことがわかる. ついでに、軌跡も大きくなるようだ.
これは、モデルの説明のときにも示したように、
「漁業を行うと魚の量がかえって増える現象」をうまく説明できている
ようにみえ、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)$
ということがわかっている.そこでこれも相図に描き込んでみよう.
|
|
なんとなく、いろいろと分かってきたようなきがするのではないだろうか.
リアルタイムに漁業の影響の大きさをいろいろ変えて、相図の様子の変化をみてみよう
グラフの描き方のところでも紹介したように、Julia の Jupyter 環境などでは Interact package を使うとパラメータをマウスでコントロールしながらリアルタイムに処理を行わせることができる. これを使って、漁業の影響を実際にみてみよう.
まずは Package の使用宣言だ.
|
|
あとはこいつを使って、グラフの例と同等に Interact package を使えば良い.
以前示したように,intearct package はうまく動かない場合がある.
そのときは「@manipulate
命令がうまく働かない場合」に対処法が書いてあるので読もう.
Interact の使い方はいぜん示したように @manipulate マクロに続けて for 文を書けば良い. 通常だと for 文はその中身がループするが,@manipulate がついているときはループせずに,for で指定されている数量などが画面で制御できるようになる.
今回の場合は以下のようにするのが目的だ.
|
|
漁業効果の係数 E をマウスで変えられるスライダーが出てくるはずなので、それを掴んで動かして、いろいろ見てみよう. これによって、漁業の影響が多少ならば、魚も増えて良いことばかりのように見えるが、漁業の影響が強くなりすぎるといろいろな点で危ないことなども見えてくるだろう.
さて、@manipulate のところを @gif と変えるだけで、今度はアニメーションを作ることが出来る.やってみよう.
|
|
┌ Info: Saved animation to
│ fn = c:\home\julia-programs\v1.5\tmp.gif
└ @ Plots c:\julia\PKG\v1.5\packages\Plots\vsE7b\src\animation.jl:104
tmp.gif という名前のファイルで保存したよ! とメッセージが出ているので、jupyter で web browser を開いた最初のタブを見てみよう.そうしたファイルがあるはずだ.それをダウンロードしてみよう.すると次のような動画になっているはずだ.
シチュエーションによっては、こうしたアニメーションのほうが様子を上手く掴めたりするので、上手に使っていくとよいだろう.
レポート
下記要領でレポートを出してみよう.
- e-mail にて,
- 題名を 2020-numerical-analysis-report-07 として,
- 教官宛(アドレスは web の "TOP" を見よう)に,
- 自分の学籍番号と名前を必ず書き込んで,
- 内容はテキストでも良いし,pdf などの電子ファイルを添付しても良いので,
下記の内容を実行して,結果や解析,感想等をレポートとして提出しよう.
- 漁業効果無しの Volterraモデルに対して、時間方向に対称なスキームを自分で作って計算してみよう.
この問題で時間対称なスキームを作ると、時間ステップを 1ステップ進めるのに連立非線形方程式を解かないといけない形になるはずだ. その解法にもしも Newton 法を使おうとすると、スカラー値常微分方程式よりも少しばかり手強いかもしれない. そうした場合は、Newton 法よりも NLsolve package を使ったほうが楽だろうね. - 上同様だが漁業効果有りにして、時間方向に対称なスキームを自分で作って計算してみよう.
- 漁業効果無しで、Runge-Kutta スキームを自分で作って計算してみよう.
やってみるとわかるが、時間方向に対称なスキームを作ることに比べると Runge-Kutta 法を用いることは格段に容易だ.実感しよう. - 漁業効果有りで、Runge-Kutta スキームを自分で作って計算してみよう.