Lotka Volterra 競争モデルの数値解

モデル式の定義

deq = {n1'[t]==r1(1-(n1[t] + a12 n2[t])/k1)n1[t],
n2'[t]==r2(1-(a21 n1[t] + n2[t])/k2)n2[t],
n1[0]==n10,n2[0]==n20}

パラメーター値の設定

(* para = {r1->1, k1->1, a12->0.8, r2->1, k2->1, a21->0.6} <--2種共存 *)
(* para = {r1->1, k1->1, a12->0.8, r2->1, k2->1, a21->1.2} <--1種のみ *)
(* para = {r1->1, k1->1, a12->1.2, r2->1, k2->1, a21->0.6} <--1種のみ *)
(* para = {r1->1, k1->1, a12->1.2, r2->1, k2->1, a21->1.1} <--1種のみ *)

para = {r1->1, k1->1, a12->0.8, r2->1, k2->1, a21->0.6}

初期値1の設定

deqN1 =deq/.para/.{n10->0.1, n20->0.1}

NDSolveで微分方程式を数値的に解く

tEnd=50;
sol1 = NDSolve[deqN1,{n1[t],n2[t]},{t,0,tEnd}]

gDyna1=Plot[ Evaluate[{n1[t],n2[t]}/.sol1],{t,0,tEnd},AxesLabel->{"t","n1, n2"},PlotRange->{0,1}]

g1 = ParametricPlot[ Evaluate[{n1[t],n2[t]}/.sol1],{t,0,50},PlotRange->All,AxesLabel->{"n1","n2"},AxesOrigin->{0,0},PlotStyle->RGBColor[1,0,0]]

異なる初期値2を設定

deqN2 =deq/.para/.{n10->1, n20->0.1}

sol2 = NDSolve[deqN2,{n1[t],n2[t]},{t,0,tEnd}]

gDyna2=Plot[ Evaluate[{n1[t],n2[t]}/.sol2],{t,0,tEnd},AxesLabel->{"t","n1, n2"}]

g2 = ParametricPlot[ Evaluate[{n1[t],n2[t]}/.sol2],{t,0,50},PlotRange->All,AxesLabel->{"n1","n2"},AxesOrigin->{0,0},PlotStyle->RGBColor[1,0,0]]

異なる初期値3を設定

deqN3 =deq/.para/.{n10->0.1, n20->1}

sol3 = NDSolve[deqN3,{n1[t],n2[t]},{t,0,tEnd}]

gDyna3=Plot[ Evaluate[{n1[t],n2[t]}/.sol3],{t,0,tEnd},AxesLabel->{"t","n1, n2"}]

g3 = ParametricPlot[ Evaluate[{n1[t],n2[t]}/.sol3],{t,0,50},PlotRange->All,AxesLabel->{"n1","n2"},AxesOrigin->{0,0},PlotStyle->RGBColor[1,0,0]]

異なる初期値3を設定

deqN4 =deq/.para/.{n10->1, n20->1}

sol4 = NDSolve[deqN4,{n1[t],n2[t]},{t,0,tEnd}]

gDyna4=Plot[ Evaluate[{n1[t],n2[t]}/.sol4],{t,0,tEnd},AxesLabel->{"t","n1, n2"}]

g4 = ParametricPlot[ Evaluate[{n1[t],n2[t]}/.sol4],{t,0,50},PlotRange->All,AxesLabel->{"n1","n2"},AxesOrigin->{0,0},PlotStyle->RGBColor[1,0,0]]

Show[g1,g2,g3,g4,PlotRange->{{0,1.5},{0,1.5}}]

アイソクラインの描画

iso1 = r1(1-(n1 + a12 n2)/  k1)n1 ;
iso2 = r2(1-(n2 + a21 n1 )/ k2)n2 ;

Solve[{iso1==0,iso2==0},{n1,n2}]

iso1N = iso1/.para
iso2N = iso2/.para

giso1=ContourPlot[iso1N==0,{n1,-0.1,2},{n2,-0.1,2},ContourStyle->RGBColor[0,0,1],PlotPoints->100]

giso2=ContourPlot[iso2N==0,{n1,-0.1,2},{n2,-0.1,2},ContourStyle->RGBColor[0,1,0],PlotPoints->100]

Show[giso1,giso2,g1,g2,g3,g4]

gVector=VectorPlot[{iso1N,iso2N},{n1,0,2},{n2,0,2}]

gVector=VectorPlot[{iso1N,iso2N},{n1,0,2},{n2,0,2},VectorScale->{Automatic,Automatic,None}]

Show[giso1,giso2,g1,g2,g3,g4,gVector,PlotRange->{{-0.1,1.5},{-0.1,1.5}}]

Show[GraphicsGrid[{{gDyna1,gDyna2},{gDyna3,gDyna4}}]]