【Julia入門】JuliaでCobweb plotを作成する

今回も書籍『A Biologist's Guide to Mathematical Modeling in Ecology and Evolution』から学習した内容の備忘録になります。

計算モデルを解析する際に使用されるCobweb Plotsを作成してみます。

対象とする計算モデルは以下のロジスティック関数とします。

\[ n(t+1) = n(t) + r_{d}n(t)\left(1-\frac{n(t)}{K}\right) \]

まずn(t)を横軸に、n(t+1)を縦軸にとるグラフを作成します。

using Plots
tmax = 100
r = 0.7
K = 1000
Recursionx  = [i for i = 1:1200]

nextgeneration(r,K,n)= n + r*n*(1-n/K)
Recursiony = nextgeneration.(r,K,Recursionx)

plot(Recursionx,Recursiony,legend=:bottomright,xlabel = "n(t)", ylabel = "n(t+1)",label="Recursion")

次に、このグラフに原点を通り傾き1の直線(Diagonal)を加えます。

plot!([0;maximum(Recursionx)],[0;maximum(Recursionx)],ls=:dash,label = "Diagonal")

ここで、この計算モデルの初期値を$n(0)$としたときに$n(t+1)$がどのように変化していくか考えます。仮に$n(0)=400$とすると、横軸の$n=400$から生じる上向き垂直線と"Recursion"との交点の縦座標が、次の$n(1) = 568$ になります。

n0 = 400
Cobnx = Float64[]
Cobny = Float64[]
push!(Cobnx,n0)
push!(Cobnx,n0)
push!(Cobny,0)
push!(Cobny,nextgeneration(r,K,n0))

plot!(Cobnx,Cobny,label = "Cobweb")

次に、$(400,568)$の点から右に水平線を伸ばし、Diagonalラインとぶつかる点が次の点です。このように、階段状に進めていくことでCobweb Plotは完成します。この階段状のラインは$n=1000$に向けて収束していくように見えます。このことから$n=1000$がこの計算モデルの安定的な平衡値であることがわかります。（このことはK=1000であることからもわかります）

n0 = 400
Cobnx = Float64[n0;n0]
Cobny = Float64[0]


push!(Cobny,nextgeneration(r,K,last(Cobnx)))

for i=1:100
    push!(Cobnx,last(Cobny))
    push!(Cobny,last(Cobny))
    push!(Cobnx,last(Cobnx))
    push!(Cobny,nextgeneration(r,K,last(Cobnx)))
end

plot!(Cobnx,Cobny,label = "Cobweb")

コードをまとめると以下のようになります。

using Plots
tmax = 100
r = 0.7
K = 1000
Recursionx  = [i for i = 1:1200]

nextgeneration(r,K,n)= n + r*n*(1-n/K)
Recursiony = nextgeneration.(r,K,Recursionx)



n0 = 400
Cobnx = Float64[n0;n0]
Cobny = Float64[0]


push!(Cobny,nextgeneration(r,K,last(Cobnx)))

for i=1:100
    push!(Cobnx,last(Cobny))
    push!(Cobny,last(Cobny))
    push!(Cobnx,last(Cobnx))
    push!(Cobny,nextgeneration(r,K,last(Cobnx)))
end

plot(Recursionx,Recursiony,legend=:bottomright,xlabel = "n(t)", ylabel = "n(t+1)",label="Recursion")
plot!([0;maximum(Cobnx)],[0;maximum(Cobnx)],ls=:dash,label = "Diagonal")
plot!(Cobnx,Cobny,label = "Cobweb")

次に、rの値を変えてみます。$r=1.8$で調べてみます。

using Plots
tmax = 100
r = 1.8
K = 1000
Recursionx  = [i for i = 1:1200]

nextgeneration(r,K,n)= n + r*n*(1-n/K)
Recursiony = nextgeneration.(r,K,Recursionx)

n0 = 825
Cobnx = Float64[n0;n0]
Cobny = Float64[0]


push!(Cobny,nextgeneration(r,K,last(Cobnx)))

for i=1:100
    push!(Cobnx,last(Cobny))
    push!(Cobny,last(Cobny))
    push!(Cobnx,last(Cobnx))
    push!(Cobny,nextgeneration(r,K,last(Cobnx)))
end

plot(Recursionx,Recursiony,legend=:bottomright,xlabel = "n(t)", ylabel = "n(t+1)",label="Recursion")
plot!([0;maximum(Cobnx)],[0;maximum(Cobnx)],ls=:dash,label = "Diagonal")
plot!(Cobnx,Cobny,label = "Cobweb")

グラフを形は大きく変わりましたが、安定な平衡点は存在します。

最後は$r=2.5$で調べてみます。

using Plots
tmax = 100
r = 2.5
K = 1000
Recursionx  = [i for i = 1:1200]

nextgeneration(r,K,n)= n + r*n*(1-n/K)
Recursiony = nextgeneration.(r,K,Recursionx)

n0 = 950
Cobnx = Float64[n0;n0]
Cobny = Float64[0]


push!(Cobny,nextgeneration(r,K,last(Cobnx)))

for i=1:100
    push!(Cobnx,last(Cobny))
    push!(Cobny,last(Cobny))
    push!(Cobnx,last(Cobnx))
    push!(Cobny,nextgeneration(r,K,last(Cobnx)))
end

plot(Recursionx,Recursiony,legend=:bottomright,xlabel = "n(t)", ylabel = "n(t+1)",label="Recursion")
plot!([0;maximum(Cobnx)],[0;maximum(Cobnx)],ls=:dash,label = "Diagonal")
plot!(Cobnx,Cobny,label = "Cobweb")