Quit[];

Kolmogorov-rendszerek 1.

Általános, több zsákmányfajból és több ragadozófajból álló modell

Farkas M. 25-30, 38-46. oldalán taláható modellek összefoglalásait adják a Kolmogorov 1.-3. állományok.

A modell

ClearAll[eq,i,j,m,n,ini];eqpray[i_,m_,n_]:=(Subscript[x, i])'[t]== (Subscript[x, i])[t](Subscript[F, i])@@Join[Array[(Subscript[x, #])[t]&,m],Array[(Subscript[y, #])[t]&,n]];
eqpred[j_,m_,n_]:=(Subscript[y, j])'[t]== (Subscript[y, j])[t](Subscript[G, j])@@Join[Array[(Subscript[x, #])[t]&,m],Array[(Subscript[y, #])[t]&,n]];
inipray[i_]:=(Subscript[x, i])[0]==(Subscript[x0, i]);inipred[j_]:=(Subscript[y, j])[0]==(Subscript[y0, j]);

Például 2 ragadozó és 3 zsákmányfaj esetén az egyenletek:

Join[Table[eqpray[i,2,3],{i,2}],Table[eqpred[j,2,3],{j,3}]]//TableForm

és a kezdeti feltételek:

Join[Table[inipray[i],{i,2}],Table[inipred[j],{j,3}]]//TableForm

{eqpray[#,2,3],eqpred[#,2,3],inipray[#]}&/@Range[2]

Tömören:

kep[m_,n_]:=Flatten[Reverse[
Partition[Sort[Flatten[Join[{eqpray[#,m,n],inipray[#]}&/@Range[m],{eqpred[#,m,n],inipred[#]}&/@Range[n]]]],m+n]]]

kep[2,3]

Még tömörebben:

ivp[m_,n_]:=Join[Array[eqpray[#,m,n]&,m],Array[eqpred[#,m,n]&,n],
Array[inipray[#]&,m],Array[inipred[#]&,n]]

ivp[2,3]

ivp[1,1]

A Kolmogorov-alak tehát azt jelenti, hogy az a változó kiemelhető a jobb oldalon, amelynek deriváltja szerepel a bal oldalon.

Jelölések

ahol:

(Subscript[x, i])(t) az i-edik zsákmányfaj létszáma a t időpontban (i=1,2,...,m)

(Subscript[y, j])(t) a j-edik ragadozófaj létszáma a t időpontban (j=1,2,...,n)

(Subscript[F, i]) az i-edik zsákmányfaj létszámának változási rátája ((Subscript[F, i])∈(𝒞^1)((ℝ^(m+n)),ℝ);   i=1,2,...(Subscript[,, ii])m)

(Subscript[G, j]) a j-edik  ragadozófaj létszámának változási rátája ((Subscript[G, j])∈(𝒞^1)((ℝ^(m+n))(Subscript[,, i])ℝ);   j=1,2,...,n)

Az (Subscript[F, i]) és (Subscript[G, j]) függvényekre vonatkozó természetes feltételek  (i=1,2,...,m;j=1,2,...,n)

(Subscript[F, i])(0,0)>0, azaz az  i-edik zsákmányfaj önállóan fennmarad

∃ (Subscript[K, i])>0∀(Subscript[x, i])>0&(Subscript[x, i])≠(Subscript[K, i]): ((Subscript[x, i])-(Subscript[K, i]))(Subscript[F, i])(0,...,0,(Subscript[x, i]),0,...,0,0)<0, azaz az  i-edik zsákmányfaj nem növekedhet korlátlanul

(Subscript[D, m+j])(Subscript[F, i])(x,y)<0, (ahol (Subscript[D, m+j]) a j-edik  ragadozófaj létszáma ((Subscript[y, j])) szerinti deriválást jelöli), azaz a j-edik  ragadozófaj létszámának növekedése csökkenti az i-edik zsákmányfaj változási rátáját

(Subscript[G, j])(0,y)<0, azaz  a j-edik  ragadozófaj létszáma táplálék hiányában csökken

(Subscript[D, i])(Subscript[G, j])(x,y)>0, (ahol (Subscript[D, i]) az i-edik zsákmányfaj létszáma ((Subscript[x, i])) szerinti deriválást jelöli), azaz az i-edik  zsákmányfaj létszámának növekedése növeli a j-edik  ragadozófaj változási rátáját

(Subscript[D, m+j])(Subscript[G, j])(x,y)≤0, azaz a j-edik  ragadozófajnál lehet belső (intraspecifikus) verseny

Egy zsákmányfajból és egy ragadozófajból álló modell (m=1,n=1)

A leggyakrabban használt paraméterek

pard={a->4,r->1,μ->16,d->1};
park={a->4,r->1,μ->16,K->5};
par=Intersection[pard,park];

Legyen a továbbiakban

Clear[ϕ,γ];
ϕ[x_,y_,K_,a_,μ_,r_]:=r(1-(x/K))-((μ y)/(a+x));
(Subscript[F, 1])[x_,y_]:=ϕ[x,y,K,a,μ,r]
γ[x_,y_,a_,μ_,d_]:=((μ x)/(a+x))-d;
(Subscript[G, 1])[x_,y_]:=γ[x,y,a,μ,d]

Vegyük észre, hogy γ és (Subscript[G, 1]) valójában nem függ y-tól. Ha

{(Subscript[F, 1])[x,y],(Subscript[G, 1])[x,y]}

ekkor az egyenletek és a kezdeti feltételek:

ivp[1,1]

Jelölések

ahol:

r(1-((Subscript[x, 1])/K)) a zsákmányfaj logisztikus növekedési rátája, és

r>0 a belső növekedési ráta (azaz születési és a halálozási ráta különbsége)

K>0 a környezet fenntartóképessége a zsákmányfaj szempontjából

((μ (Subscript[x, 1])(Subscript[y, 1]))/(a+(Subscript[x, 1]))) az úgynevezett Holling-féle funkcionális válasz, és

μ>0 a ragadozófaj maximális növekedési rátája

a>0 az úgynevezett félszaturációs állandó, ugyanis x=a esetén lesz a ragadozó növekedési rátája a maximális érték fele

a ragadozó létszámával arányosan csökkenti a zsákmányfaj növekedési rátáját

d>0 a ragadozófaj mortalitása

Az Allee-hatás zónája

A zsákmány nullavonala (nullklínája)

ClearAll[preynullcline];
preynullcline:=Function[{x,K,a,μ,r},FullSimplify[Solve[ϕ[x,y,K,a,μ,r]==0,y][[1,1,2]]]]

preynullcline[x,K,a,μ,r]

A fenti függvény a zsákmány nullavonala gyakran használt paraméterértékeknél.

preynullcline[x,K,a,μ,r]/.park

A fenti függvény a zsákmány nullavonala a tipikus paraméterértékek mellett.

ClearAll[maxpy];
maxpy:=Function[{K,a},FullSimplify[Solve[D[preynullcline[x,K,a,μ,r],x]==0,x][[1,1,2]]]]

maxpy[K,a]

A fenti a zsákmány nullavonalának maximumhelye.

maxvalpy:=Function[{K,a,μ,r},preynullcline[maxpy[K,a],K,a,μ,r]]

maxvalpy[K,a,μ,r]

A fenti a zsákmány nullavonalának maximuma.

ClearAll[ply];ply[K_,opts___]:=Plot[Evaluate[preynullcline[x,K,4,16,1]],{x,0,5.1},AxesLabel->{(Subscript[x, 1]),(Subscript[y, 1])},AspectRatio->1,PlotRange->{0,0.4},opts]

p1=ply[3,PlotLabel->Style["K<a",Medium]];
p2=ply[4,PlotLabel->Style["K=a",Medium]];
p3=ply[5,PlotLabel->Style["K>a",Medium],Epilog->{Dashed,Line[{{(1/2),0},{(1/2),((7 6)/(16 9))}}],Text[Style["Allee zóna",8],{0.25,0.1},{0,0},{0,1}]}];
Show[GraphicsRow[{p1,p2,p3}],ImageSize->500,PlotLabel->"A zsákmány tipikus nullklínája"]

Az Allee-zóna különösen az emlős fajokra jellemző. Ha a létszám túl kicsi, akkor a szaporodáshoz nehéz párt találni. Egy küszöbig a létszám növekedése előnyös a növekedési rátára nézve.

A ragadozó nullavonala

ClearAll[prednullcline];
prednullcline:=Function[{a,μ,d},FullSimplify[Solve[γ[x,y,a,μ,d]==0,x][[1,1,2]]]]

prednullcline[a,μ,d]

A fenti függőleges egyenes a ragadozó nullavonala.

prednullcline[a,μ,d]/.pard

A fenti függőleges egyenes a ragadozó nullklínája a tipikus paraméterértékeknél.

prednullcline[4,16,d]

A fenti függőleges egyenes a ragadozó nullklínája d függvényeként, miközben a többi paraméter értékét tipikusnak vesszük.

Ezek után azt vizsgáljuk, hogy milyen d érték mellett esik a ragadozó nullklínája éppen az áldozat nullklínájának maximumhelyére?

dkrit=Solve[prednullcline[a,μ,d]==maxpy[K,a],d][[1,1,2]]//Simplify

Nevezzük ezt a továbbiakban kritikus d-értéknek.

dkrit/.park

A linearizált egyenlet Jacobi-mátrixa

Határozzuk meg először a stacionárius pontokat.

ClearAll[stac];
stac[K_,a_,μ_,r_,d_]:=Sort[FullSimplify[Solve[{ϕ[x,y,K,a,μ,r]==0,γ[x,y,a,μ,d]==0},{x,y}][[1]]]]

stac[K,a,μ,r,d]

Ha tehát d=dkrit, akkor a stacionárius pont első koordinátája éppen a maximumhelyre esik.

{x/.stac[K,a,μ,r,dkrit],maxpy[K,a]}

pp[d_,opts___,igen_:True]:=ParametricPlot[{{10t,preynullcline[10t,5,4,16,1]},{prednullcline[4,16,d],t}},{t,0,0.5},AspectRatio->1,PlotLabel->If[igen,"d = "<>ToString[TraditionalForm[d]],""],Epilog->{Magenta,PointSize[Large],Point[{maxpy[5,4],maxvalpy[5,4,16,1]}],Text["Maximum",{maxpy[5,4],maxvalpy[5,4,16,1]-0.03}],Red,PointSize[Large],Point[{x,y}/.stac[5,4,16,1,d]],Text["Egyensúlyi pont",{x,y+0.03}/.stac[5,4,16,1,d]]},opts]

GraphicsRow[{pp[(16/9),PlotRange->{{-1,5},Automatic},True],pp[5,True]},ImageSize->400]

Maga a Jacobi-mátrix:

MatrixForm[D[{x ϕ[x,y,K,a,μ,r],y γ[x,y,a,μ,d]},{{x,y}}]]//Simplify

És a stacionárius pontban:

ClearAll[J];
J[K_,a_,μ_,r_,d_]:=D[{x ϕ[x,y,K,a,μ,r],y γ[x,y,a,μ,d]},{{x,y}}]/.stac[K,a,μ,r,d]

MatrixForm[J0=Simplify[J[K,a,μ,r,d]]]

Sajátértékei pontosan akkor tiszta képzetes számok, amikor a nyoma nulla:

Reduce[J0[[1,1]]==0,d]//Simplify

A (nemtriviális) harmadik eset éppen dkrit mellett következik be, azaz arra a d értékre, amelyikre a ragadozó nullklínája épp az áldozat nullklínájának maximumán megy át.

Trajektóriák, iránymező

Az Allee-hatás zónájának környezetét kinagyítva Hopf-bifurkáció figyelhető meg.

Manipulate[Show[pp[d,PlotPoints->25,False],StreamPlot[{x ϕ[x,y,5,4,16,1],y γ[x,y,4,16,d]},{x,0,5.1},{y,0,4},StreamStyle->Directive[Orange,Thick]],PlotRange->{0,0.5},AspectRatio->0.7],{{d,(16/9)},0,9,1/10},SaveDefinitions->True]

dkrit/.park

A d=16/9 értéktől távolodva megszűnik a periodikus megoldás.

A következőkben lerajzoljuk a periodikus megoldásokat.

ndsazigazi[d_]:=NDSolve[{x'[t]==x[t]ϕ[x[t],y[t],5,4,16,1],y'[t]==y[t]γ[x[t],y[t],4,16,d],x[0]==0.5,y[0]==0.3},{x,y},{t,0,100}][[1]];

stac[5,a,μ,r,d]/.pard//N

Plot[Evaluate[{x[t],y[t]}/.ndsazigazi[16/9]],{t,0,100},AxesLabel->{"t","x(t),y(t)"}]

ParametricPlot[Evaluate[{x[t],y[t]}/.ndsazigazi[16/9]],{t,0,100},AxesLabel->{"x(t)","y(t)"}]

Plot[Evaluate[{x[t],y[t]}/.ndsazigazi[3]],{t,0,100},PlotLabel->"Másik d értékre nem oszcillál,\nstabilizálódik",AxesLabel->{"t","x(t),y(t)"}]

Plot[Evaluate[{x[t],y[t]}/.ndsazigazi[1]],{t,0,100},PlotLabel->"Másik d értékre nem oszcillál,\ndestabilizálódik",AxesLabel->{"t","x(t),y(t)"}]

Ha a stacionárius ponttól távolról indítjuk, akkor d=(16/9) esetén sem oszcillál.

Manipulate[Module[{nds,x,y},nds=NDSolve[{x'[t]==x[t]ϕ[x[t],y[t],5,4,16,1],y'[t]==y[t]γ[x[t],y[t],4,16,d],x[0]==.5,y[0]==.3},{x,y},{t,0,100}][[1]];
ParametricPlot3D[Evaluate[{{0,x[t],y[t]},{t,x[t],0},{t,4,y[t]}}/.nds],{t,0,100},BoxRatios->{1,1,1},PlotPoints->100,PlotRange->{{0,100},{0,4},{0,1}},AxesLabel->{Style["t",Italic],Row[{Style["x",Italic],"(",Style["t",Italic],")"}],Row[{Style["y",Italic],"(",Style["t",Italic],")"}]}, PlotStyle->Directive[Thick,RGBColor[.49,0,0]],ImageSize->{400,400},PlotLabel->"Trajectory and coordinate functions of the solutions\n\n"]],
{{d,1.77778},0,5,0.1,Appearance->"Labeled"},SaveDefinitions->True]

Ehhez vegyük például az (5,3) kezdeti feltételt.

Manipulate[Module[{nds,x,y},nds=NDSolve[{x'[t]==x[t]ϕ[x[t],y[t],5,4,16,1],y'[t]==y[t]γ[x[t],y[t],4,16,d],x[0]==5,y[0]==3},{x,y},{t,0,100}][[1]];
ParametricPlot3D[Evaluate[{{0,x[t],y[t]},{t,x[t],0},{t,4,y[t]}}/.nds],{t,0,100},BoxRatios->{1,1,1},PlotPoints->100,PlotRange->{{0,100},{0,4},{0,1}},AxesLabel->{Style["t",Italic],Row[{Style["x",Italic],"(",Style["t",Italic],")"}],Row[{Style["y",Italic],"(",Style["t",Italic],")"}]}, PlotStyle->Directive[Thick,RGBColor[.49,0,0]],ImageSize->{400,400},PlotLabel->"Trajectory and coordinate functions of the solutions\n\n"]],
{{d,1.77778},0,5,0.1,Appearance->"Labeled"},SaveDefinitions->True]

Hivatkozás

Farkas, M.: Dynamical models in biology, Academic Press, New York, 2001.