Lotka Volterra 捕食モデルの数値解

Lotka Volterra捕食モデル

モデル式の定義

deq = {h'[t]==rH h[t] - a h[t]p[t],
p'[t]==-rP p[t]+b h[t]p[t],
h[0]==h0,p[0]==p0}

パラメーター値の設定

para = {rH->1,a->1, rP->1.2, b->1}

初期値の設定 1

deqN1 =deq/.para/.{h0->0.2, p0->0.3}

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

tEnd=50;
sol1 = NDSolve[deqN1,{h[t],p[t]},{t,0,tEnd}]

Plot[ Evaluate[{h[t],p[t]}/.sol1],{t,0,tEnd},PlotRange->{0,2},AxesLabel->{"t","h[t], p[t]"}]

g1 = ParametricPlot[ Evaluate[{h[t],p[t]}/.sol1],{t,0,tEnd},AspectRatio->1,PlotRange->{{0,5},{0,5}},AxesLabel->{"h[t]", "p[t]"},PlotStyle->RGBColor[0,0,0]]

初期値の設定 2

deqN2 =deq/.para/.{h0->0.5, p0->0.5}

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

tEnd=50;
sol2 = NDSolve[deqN2,{h[t],p[t]},{t,0,tEnd}]

Plot[ Evaluate[{h[t],p[t]}/.sol2],{t,0,tEnd},PlotRange->{0,2},AxesLabel->{"t","h[t], p[t]"}]

g2 = ParametricPlot[ Evaluate[{h[t],p[t]}/.sol2],{t,0,tEnd},AspectRatio->1,PlotRange->{{0,5},{0,5}},AxesLabel->{"h[t]", "p[t]"},PlotStyle->RGBColor[0,0,0]]

初期値の設定 3

deqN3 =deq/.para/.{h0->0.9, p0->0.7}

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

tEnd=50;
sol3 = NDSolve[deqN3,{h[t],p[t]},{t,0,tEnd}]

Plot[ Evaluate[{h[t],p[t]}/.sol3],{t,0,tEnd},PlotRange->{0,2},AxesLabel->{"t","h[t], p[t]"}]

g3 = ParametricPlot[ Evaluate[{h[t],p[t]}/.sol3],{t,0,tEnd},AspectRatio->1,PlotRange->{{0,5},{0,5}},AxesLabel->{"h[t]", "p[t]"},PlotStyle->RGBColor[0,0,0]]

Show[g1,g2,g3]

isoH=rH h - a h p
isoP = -rP p + b h p

isoHN = isoH/.para
isoPN = isoP/.para

gisoH=ContourPlot[isoHN==0,{h,-0.1,5},{p,-0.1,5},ContourStyle->{Thick,RGBColor[0,0,1]},PlotPoints->100]

gisoP=ContourPlot[isoPN==0,{h,-0.1,5},{p,-0.1,5},ContourStyle->{Thick,RGBColor[1,0,0]},PlotPoints->100]

Show[gisoH,gisoP]

Show[gisoH,gisoP,g1,g2,g3]

保存量

fn = b H - rP Log[H]+a P -rH Log[P]

para={rH->1,rP->1,a->1, b->1.2}

fnN=fn/.para

Plot3D[fnN,{H,0,6},{P,0,6}]

ContourPlot[fnN,{H,0,6},{P,0,6}]

Lotka Volterra捕食モデル + Logistic growth of the prey

deq = {h'[t]==rH( 1-h[t]/k)h[t] - a h[t]p[t],
p'[t]==-rP p[t]+b h[t]p[t],
h[0]==h0,p[0]==p0}

para = {rH->2,a->1, rP->0.8, b->0.4,k->8}

deqN1 =deq/.para/.{h0->1, p0->0.5}

tEnd=50;
sol1 = NDSolve[deqN1,{h[t],p[t]},{t,0,tEnd}]

Plot[ Evaluate[{h[t],p[t]}/.sol1],{t,0,tEnd},PlotRange->{0,10},AxesLabel->{"t","h[t], p[t]"},PlotStyle->Thick]

g1 = ParametricPlot[ Evaluate[{h[t],p[t]}/.sol1],{t,0,tEnd},AspectRatio->1,PlotRange->{{0,8},{0,8}},AxesLabel->{"h[t]", "p[t]"},PlotStyle->{Thick,RGBColor[0,0,0]}]

isoH=rH( 1-h/k)h - a h p
isoP = -rP p + b h p

isoHN=isoH/.para
isoPN=isoP/.para

gisoH=ContourPlot[isoHN==0,{h,-0.1,10},{p,-0.1,10},ContourStyle->{Thick,RGBColor[0,0,1]},PlotPoints->100]

gisoP=ContourPlot[isoPN==0,{h,-0.1,10},{p,-0.1,10},ContourStyle->{Thick,RGBColor[1,0,0]},PlotPoints->100]

Show[gisoH,gisoP]

Show[gisoH,gisoP,g1]

Lotka Volterra捕食モデル+ Logistic growth of the 
prey + Satiation of the predator

Plot[{a h/(1+c h),a h}/.{a->1,c->0.2},{h,0,40},PlotRange->{0,6},PlotStyle->Thick,AxesLabel->{"H","a H/(1+c H)"}]

deq = {h'[t]==rH( 1-h[t]/k)h[t] - a h[t]p[t]/(1+c h[t]),
p'[t]==-rP p[t]+b h[t]p[t]/(1+c h[t]),
h[0]==h0,p[0]==p0}

para = {rH->1,a->1, rP->1, b->8,k->1,c->2}

Hのヌルクライン （放物線） の軸

(1-1/c)/2/.para

Pのヌルクラインの軸

1/(b-c)/.para

deqN1 =deq/.para/.{h0->0.2, p0->0.1}

tEnd=100;
sol1 = NDSolve[deqN1,{h[t],p[t]},{t,0,tEnd}]

Plot[ Evaluate[{h[t],p[t]}/.sol1],{t,0,tEnd},PlotRange->{0,5},AxesLabel->{"t","h[t], p[t]"}]

gDyna = ParametricPlot[ Evaluate[{h[t],p[t]}/.sol1],{t,0,tEnd},AspectRatio->1,PlotRange->{{0,1},{0,4}},AxesLabel->{"h[t]", "p[t]"},PlotStyle->{Thick,RGBColor[0,0,0]}]

isoH=rH h(1-h/k)-a h p/(1+c h);
isoP=-rP p+b h p/(1+c h);

isoHN=isoH/.para
isoPN=isoP/.para

gisoH=ContourPlot[isoHN==0,{h,-0.1,1.5},{p,-0.1,5},ContourStyle->{Thick,RGBColor[0,0,1]},PlotPoints->100]

gisoP=ContourPlot[isoPN==0,{h,-0.1,1.5},{p,-0.1,5},ContourStyle->{Thick,RGBColor[1,0,0]},PlotPoints->100]

Show[gisoH,gisoP]

Show[gisoH,gisoP,gDyna]