システムとモデリング

modelica, Julia, Design Structure Matrix, SysML, 他モデリング全般について。

【Julia入門】Juliaでベクトル場をプロットする。

 

今回も書籍『A Biologist's Guide to Mathematical Modeling in Ecology and Evolution』 の写経をします。

以下の2つのロジスティック方程式について、縦軸に$n_2$、横軸に$n_1$を取ってベクトル場をプロットするとします。

\[ n_{1}(t+1)=n_{1}(t) + r_{1}n_{1}(t)\left(1-\frac{n_{1}(t)+\alpha_{12}n_{2}(t)}{K_{1}}\right)\\ n_{2}(t+1)=n_{2}(t) + r_{2}n_{1}(2)\left(1-\frac{n_{2}(t)+\alpha_{21}n_{1}(t)}{K_{2}}\right) \]

ここで各パラメーターは以下とします。

パラメーター
$r_{1}$ 0.5
$r_{2}$ 0.5
$\alpha_{12}$ 0
$\alpha_{21}$ 0.5
$K_{1}$ 1000
$K_{2}$ 1000

まず、適当な初期値$(n_{1}(0),n_{2}(0))=(50,1200)$を考えます。

この値を上式に代入しt=1のときの値を求めると $(n_{1}(1),n_{2}(1))=(73.75,1065)$になります。

この2点を矢印で結びます。これで1本のベクトルが完成します。

今回は$(n_{1}(0),n_{2}(0))=(50,1200)$としましたが、他の位置を$(n_{1}(0),n_{2}(0))$して上記を 繰り返すことでベクトル場を作成します。

using Plots

r1=0.5
r2=0.5
K1=1000
K2=1000
α12 = 0
α21 =0.5

ParameterN1 = (r1,K1,α12)
ParameterN2 = (r2,K2,α21)

nextgeneration(n1,n2,r,K,α)= n1 + r*n1*(1-(n1+α*n2)/K)

N = zeros(10,10)

global p =plot(xlabel="n1",ylabel="n2",xlims=(0,1300),ylims=(0,1300))

for i in 1:10 ,j in 1:10
global p=plot!([120*i,nextgeneration(120*i,120*j,r1,K1,α12)],[120*j,nextgeneration(120*j,120*i,r2,K2,α21)],arrow = 2, label="",
linecolor = "black")
end
p

無事、ベクトル場をプロットできました。 短いですが今回はここまでにします。