Julia言語で動的モデルのカオス挙動をシミュレーションする

こんばんは！Otepipiです。 今回も前回に引き続きJulia言語を使って数理モデリングとシミュレーションをしていきたいと思います。 例題は下記「数理モデリング入門」のものを使用しています。

数理モデリング入門 ―ファイブ・ステップ法― 原著第4版

シロナガスクジラの個体数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になりますと、個体数は収束せず・発散もせず、周期２の振動を続けます。

Δ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がこちらになります。

解が収束した状態からカオスに至るまでが確認できました。