Julia言語で動的モデルのカオス挙動をシミュレーションする
こんばんは!Otepipiです。 今回も前回に引き続きJulia言語を使って数理モデリングとシミュレーションをしていきたいと思います。 例題は下記「数理モデリング入門」のものを使用しています。
- 作者: Mark M. Meerschaert,佐藤一憲,梶原毅,佐々木徹,竹内康博,宮崎倫子,守田智
- 出版社/メーカー: 共立出版
- 発売日: 2015/01/24
- メディア: 単行本
- この商品を含むブログ (1件) を見る
シロナガスクジラの個体数xとナガスクジラの個体数yの関係は以下で表すことができる。
$$ \frac{dx}{dt}=0.05x \left(1-\frac{x}{150,000}\right) -\alpha xy $$
$$ \frac{dy}{dt}=0.08y \left(1-\frac{y}{400,000}\right) -\alpha xy $$
初期値を下記のように仮定する。 $$x(0)=5,000 $$
$$y(0)=70,000$$
$$\alpha = 1.0×10^{-8}$$
未来のxとyを決定せよ。
今回はPlotsパッケージを使用します。
using Plots
前回と同じ連立の常微分方程式ですが、今回はルンゲクッタ法ではなくオイラー法を使用します。下記は差分Δt=1年としたときのものです。 また、ステップ数はN=50回とします。
const α = 1e-8 const N = 50 const Δt = 1 x = zeros(N) y = zeros(N) t = zeros(N) x[1] = 5000 y[1] = 70000 t[1] = 0 for i = 1: N-1 x[i+1]=x[i] + (0.05x[i]*(1-x[i]/150000)-α*x[i]*y[i])*Δt y[i+1]=y[i] + (0.08y[i]*(1-y[i]/400000)-α*x[i]*y[i])*Δt t[i+1] = i*Δt end
この結果をプロットします
plot(t,x,seriestype=:scatter,label="x", title = "Δt=1",xlabel="Year", ylabel="Population"); plot!(t,y,seriestype=:scatter,label="y");
シロナガスクジラの個体数x、ナガスクジラの個体数yともに解が現れています。このタイムスパンでは一定値への収束は確認できません。
次に時間間隔Δtを1年から10年に変更して再計算・結果をプロットしてみます。
解が一定に収束していくのを確認できます。 Δt=24ではどうでしょうか?
ここでも解は収束するものの、序盤に振動が見られます。
では、Δt=27ではどうなるでしょうか?
Δt=27になりますと、個体数は収束せず・発散もせず、周期2の振動を続けます。
Δt=32ではどうでしょう。
ここでも個体数は振動を続けますが、その周期は4に変わっています。
そしてΔt=37を見てみます。
ここではカオスが現れています。
Δtを37より大きくすると、解は無限遠に発散して計算できなくなります。
最後にこれらの画像をgifにしてみます。
まずΔt=1………37に対してオイラー法の計算を繰り返します。
const α = 1e-8 const N = 50 const max = 37 x = zeros(N , max) y = zeros(N , max) t = zeros(N , max) for Δt = 1 : max x[1,Δt] = 5000 y[1,Δt] = 70000 t[1,Δt] = 0 for i = 1: N-1 x[i+1,Δt]=x[i,Δt] + (0.05x[i,Δt]*(1-x[i,Δt]/150000)-α*x[i,Δt]*y[i,Δt])*Δt y[i+1,Δt]=y[i,Δt] + (0.08y[i,Δt]*(1-y[i,Δt]/400000)-α*x[i,Δt]*y[i,Δt])*Δt t[i+1,Δt] = i*Δt end end
次にΔt=1………37までのプロットをアニメーションにします。
anim = @animate for i = 1 : max plot(t[:,i],x[:,i],seriestype=:scatter,label="x", title = string("Δt =", i),xlabel="Year", ylabel="Population") plot!(t[:,i],y[:,i],seriestype=:scatter,label="y") end #gif(anim,"\anoi.gif",fps=30) エラーになってしまいました。
gifで書き出そうと思ったのですがgif()でエラーを吐いてしまいました。 今回は連番で書き出したpingをImageMagicで手動でGIF化することにしました。 完成したGIFがこちらになります。
解が収束した状態からカオスに至るまでが確認できました。