Quit[];

Két fajból álló kompetitív és kooperatív Lotka-Volterra-rendszerek: 2. Kooperatív eset

Farkas Miklós: Dynamical models in biology (Academic Press 2001), 35. old.

A modell

eq={x'[t]== x[t]((Subscript[r, 1])-(Subscript[a, 11])x[t]-(Subscript[a, 12])y[t]),y'[t]== y[t]((Subscript[r, 2])-(Subscript[a, 21])x[t]-(Subscript[a, 22])y[t])};ini={x[0]==x0,y[0]==y0};

f[x_,y_]:=x((Subscript[r, 1])-(Subscript[a, 11])x-(Subscript[a, 12])y);g[x_,y_]:=y((Subscript[r, 2])-(Subscript[a, 21])x-(Subscript[a, 22])y);

eq=={x'[t]==f[ x[t],y[t]],y'[t]==g[x[t],y[t]]}

Jelölések

Ahol:

x(t) az egyik faj létszáma a t időpontban

y(t) a másik faj létszáma a t időpontban

1.  (Subscript[r, 1]),(Subscript[r, 2])>0  fakultatív kooperáció esetén

2.  (Subscript[r, 1]),(Subscript[r, 2])<0  kötelező kooperáció esetén

(Subscript[a, 11]),(Subscript[a, 22])≥0 van (lehet) intraspecifikus verseny

(Subscript[a, 12]),(Subscript[a, 21])<0 van kooperáció

Kiegészítő magyarázat

1. azt jelenti, hogy az egyik faj a másik nélkül fenn tud maradni, ezért lesz a kooperáció fakultatív

2. azt jelenti, hogy az egyik faj a másik nélkül nem tud fennmaradni, ezért lesz a kooperáció kötelező

feltehető (1 valószínűséggel), hogy Det[A]≠0, ahol A:=((Subscript[a, ij])) az interakciós mátrix

Egyensúlyi helyzetek (vagy stacionárius pontok)

stac=Solve[{0==f[x,y],0==g[x,y]},{x,y}]

{Ν3,P3}={x,y}/.stac[[2]]

(E^0):=(0,0)

(E^1):=(((Subscript[r, 1])/(Subscript[a, 11])),0)

(E^2):=(0,((Subscript[r, 2])/(Subscript[a, 22])))

(E^3):=((((Subscript[a, 22]) (Subscript[r, 1])-(Subscript[a, 12]) (Subscript[r, 2]))/(Det[A])),(((Subscript[a, 11]) (Subscript[r, 2])-(Subscript[a, 21]) (Subscript[r, 1]))/(Det[A])))∈((((SuperPlus[ℝ])))^2) az alábbi két esetben:

1. A)   (Subscript[r, 1]),(Subscript[r, 2])>0 és Det[A]>0

2. A)   (Subscript[r, 1]),(Subscript[r, 2])<0 és Det[A]<0

(de nem kizárólag ekkor!)

(E^3):=((((Subscript[a, 22]) (Subscript[r, 1])-(Subscript[a, 12]) (Subscript[r, 2]))/(Det[A])),(((Subscript[a, 11]) (Subscript[r, 2])-(Subscript[a, 21]) (Subscript[r, 1]))/(Det[A])))∉((((SuperPlus[ℝ])))^2) az alábbi két esetben:

1. B)  (Subscript[r, 1]),(Subscript[r, 2])>0 és Det[A]<0

2. B)  (Subscript[r, 1]),(Subscript[r, 2])<0 és Det[A]>0

(de nem kizárólag ekkor!)

Előkészületek

Definiáljuk a különböző eseteket.

kozos:=And[(Subscript[a, 11])≥0, (Subscript[a, 22])≥0,(Subscript[a, 12])<0, (Subscript[a, 21])<0];
a1eset:=And[(Subscript[r, 1])>0,(Subscript[r, 2])>0,(Subscript[a, 11]) (Subscript[a, 22])-(Subscript[a, 12]) (Subscript[a, 21])>0,kozos];
a2eset:=And[(Subscript[r, 1])<0,(Subscript[r, 2])<0,(Subscript[a, 11]) (Subscript[a, 22])-(Subscript[a, 12]) (Subscript[a, 21])<0,kozos];
b1eset:=And[(Subscript[r, 1])>0,(Subscript[r, 2])>0,(Subscript[a, 11]) (Subscript[a, 22])-(Subscript[a, 12]) (Subscript[a, 21])<0,kozos];b2eset:=And[(Subscript[r, 1])<0,(Subscript[r, 2])<0,(Subscript[a, 11]) (Subscript[a, 22])-(Subscript[a, 12]) (Subscript[a, 21])>0,kozos];

ClearAll[nds];nds[x0_,y0_,ize_]:=NDSolve[{x'[t]==x[t]((Subscript[r, 1])-(Subscript[a, 11])x[t]-(Subscript[a, 12])y[t]),y'[t]==y[t]((Subscript[r, 2])-(Subscript[a, 21])x[t]-(Subscript[a, 22])y[t]),x[0]==x0,y[0]==y0}/.ize,{x[t],y[t]},{t,0,1000}][[1]]

xnull=Solve[f[x,y]==0,y][[1,1,2]]

ynull=Reduce[{g[x,y]==0,y>0,(Subscript[a, 21])>0,(Subscript[a, 22])>0,(Subscript[r, 2])>0},y,Reals][[5,2]]

Részletes vizsgálata az 1. A) esetben

ClearAll[par1a];
par1a={(Subscript[r, 1])->1,(Subscript[r, 2])->1,(Subscript[a, 11])->2,(Subscript[a, 12])->-1,(Subscript[a, 21])->-2,(Subscript[a, 22])->2};

a1eset/.par1a

Valóban idetartozó esetet írnak le a paraméterek.

Következzenek a nullavonalak a stacionárius pontokkal együtt..

nullclines=Plot[{xnull,ynull}/.par1a,{x,0,2.5},PlotRange->{-0.1,3.5},AspectRatio->(1/GoldenRatio),Epilog->{Red,PointSize[0.02],Point/@({x,y}/.stac/.par1a)},PlotLabel->"1. A) eset"]

Magyarázó feliratokat is elhelyezünk.

nullclinesfull=Plot[{xnull,ynull}/.par1a,{x,0,2.5},PlotRange->{-0.1,3.5},AspectRatio->(1/GoldenRatio),Epilog->{Red,PointSize[0.02],Point/@({x,y}/.stac/.par1a),{Black,Text[((Subscript[r, 1])/(Subscript[a, 11])),{0.53,0.5}]},{Black,Text[((Subscript[r, 2])/(Subscript[a, 22])),{0.15,0.5}]}},PlotLabel->"1. A) eset"]

vp=VectorPlot[{f[x,y],g[x,y]}/.par1a,{x,0,2.5},{y,0,4},AxesOrigin->{0,0},VectorScale->0.1,AspectRatio->(1/GoldenRatio),PlotLabel->"1. A) eset"];

Igazából a StreamPlot eredménye használhatób a legtöbb esetben.

sp=StreamPlot[{f[x,y],g[x,y]}/.par1a,{x,0,2.5},{y,0,3.5},AxesOrigin->{0,0},AspectRatio->(1/GoldenRatio),VectorScale->0.1,PlotLabel->"1. A) eset"]

Show[nullclines,sp,PlotLabel->"1. A) eset"]

pp[x0_,y0_,ize_]:=ParametricPlot[Evaluate[{x[t],y[t]}/.
nds[x0,y0,ize]],{t,0,1000},PlotRange->All,PlotPoints->150,AspectRatio->(1/GoldenRatio),PlotStyle->Directive[Thick,Red]];

p1=pp[0.2,0.4,par1a];p2=pp[0.2,3,par1a];p3=pp[1.2,0.3,par1a];p4=pp[2,0.3,par1a];p5=pp[2,2,par1a];p6=pp[2,4,par1a];

Show[nullclines,sp,p1,p2,p3,p4,p5,p6,PlotLabel->"1. A) eset"]

Próbáljuk meg a Locator felhasználásával megmutatni a kezdeti értékektől való függést.

Megjegyzések

Látszik, hogy (E^3) aszimptotikusan stabilis (globálisan is), (E^0),(E^1),(E^2) pedig instabilis

Szemléletesen nyilvánvaló, hogy Det[A]>0, azaz -(Subscript[a, 12]) (Subscript[a, 21])+(Subscript[a, 11]) (Subscript[a, 22])>0 azt jelenti, hogy az intraspecifikus verseny erősebb, mint az interspecifikus kooperáció.

Összegezve és a későbbiek kedvéért:

A fenti számolásokat automatikusan végeztethetjük el az alábbi programmal, amely már csak az eset nevét kéri.

abra[adat_]:=Module[{ize=ToExpression["par"<>ToString[adat]]},nullclines=Plot[{xnull,ynull}/.ize,{x,-0.1,3},PlotRange->{-0.1,3.5},Epilog->{Red,PointSize[0.02],Point/@({x,y}/.stac/.ize)}];
sp=StreamPlot[{f[x,y],g[x,y]}/.ize,{x,-0.1,3},{y,-0.1,3.5},AxesOrigin->{0,0},VectorScale->0.1];
Show[nullclines,sp,PlotLabel->ToString[adat]<>") eset"]]

abra["1a"]

Ugyanazt adja, ezt fogjuk használni.

Részletes vizsgálata az 2. A) esetben

ClearAll[par2a];
par2a={(Subscript[r, 1])->-1,(Subscript[r, 2])->-1,(Subscript[a, 11])->1,(Subscript[a, 12])->-1,(Subscript[a, 21])->-2,(Subscript[a, 22])->1};

a2eset/.par2a

Valóban idetartozó esetet írnak le a paraméterek.

nullclines=Plot[{xnull,ynull}/.par2a,{x,0,3},PlotRange->{-0.1,3.5},AspectRatio->(1/GoldenRatio),Epilog->{Red,PointSize[0.02],Point/@({x,y}/.stac/.par2a)},PlotLabel->"2. A) eset"];

nullclinesfull=Plot[{xnull,ynull}/.par2a,{x,0,3},PlotRange->{-0.1,3.5},AspectRatio->(1/GoldenRatio),Epilog->{Red,PointSize[0.02],Point/@({x,y}/.stac/.par2a),{Black,Text[((Subscript[r, 2])/(Subscript[a, 21])),{0.52,0.45}]},{Black,Text[((Subscript[r, 1])/(Subscript[a, 12])),{0.12,0.94}]}},PlotLabel->"2. A) eset"]

vp=VectorPlot[{f[x,y],g[x,y]}/.par2a,{x,0,2},{y,0,4},AxesOrigin->{0,0},VectorScale->0.1,AspectRatio->(1/GoldenRatio),PlotLabel->"2. A) eset"];

sp=StreamPlot[{f[x,y],g[x,y]}/.par2a,{x,0,3},{y,0,3.5},AxesOrigin->{0,0},AspectRatio->(1/GoldenRatio),VectorScale->0.1,PlotLabel->"2. A) eset"];

Ahhoz, hogy az abra program helyesen működjön, újra ki kellett számolni az aktuális paraméterekkel az alkatrészeket.

abra["2a"]

ClearAll[nds];nds[x0_,y0_,ig_,ize_]:=NDSolve[{x'[t]==x[t]((Subscript[r, 1])-(Subscript[a, 11])x[t]-(Subscript[a, 12])y[t]),y'[t]==y[t]((Subscript[r, 2])-(Subscript[a, 21])x[t]-(Subscript[a, 22])y[t]),x[0]==x0,y[0]==y0}/.ize,{x[t],y[t]},{t,0,ig}][[1]]

pp[x0_,y0_,ig_,ize_]:=ParametricPlot[Evaluate[{x[t],y[t]}/.
nds[x0,y0,ig,ize]],{t,0,ig},PlotRange->All,PlotPoints->150,AspectRatio->(1/GoldenRatio),PlotStyle->Directive[Thick,Red]];

p1=pp[0.7,3,1000,par2a];p2=pp[0.1,3,1000,par2a];p3=pp[1.2,0.3,1000,par2a];p4=pp[2,0.3,1000,par2a];p5=pp[2,2,1000,par2a];p6=pp[2,3,1000,par2a];

Show[nullclines,sp,p1,p2,p3,p4,p5,p6,PlotLabel->"2. A) eset"]

Megjegyzések

Látszik, hogy (E^0) aszimptotikusan stabilis, globálisan is.

Szemléletesen nyilvánvaló, hogy Det[A]>0, azaz (Subscript[a, 11]) (Subscript[a, 22])-(Subscript[a, 12]) (Subscript[a, 21])>0 azt jelenti, hogy túl erős az intraspecifikus verseny az interspecifikus kooperációhoz képest, ezért a fajok kihalnak.

Az itt adott magyarázatokat az alábbi három esetben nem ismételjük meg.

Részletes vizsgálata az 1. B) esetben

ClearAll[par1b];
par1b={(Subscript[r, 1])->1,(Subscript[r, 2])->1,(Subscript[a, 11])->1,(Subscript[a, 12])->-2,(Subscript[a, 21])->-2,(Subscript[a, 22])->2};

b1eset/.par1b

Valóban idetartozó esetet írnak le a paraméterek.

nullclinesfull=Plot[{xnull,ynull}/.par1b,{x,0,4},PlotRange->{-0.1,3.5},AspectRatio->(1/GoldenRatio),Epilog->{Red,PointSize[0.02],Point/@({x,y}/.stac/.par1b),{Black,Text[((Subscript[r, 1])/(Subscript[a, 11])),{1.03,0.45}]},{Black,Text[((Subscript[r, 2])/(Subscript[a, 22])),{0.25,0.45}]}},PlotLabel->"1. B) eset"]

vp=VectorPlot[{f[x,y],g[x,y]}/.par1b,{x,0,4},{y,0,4},AxesOrigin->{0,0},VectorScale->0.1,AspectRatio->(1/GoldenRatio),PlotLabel->"1. B) eset"];

StreamPlot[{f[x,y],g[x,y]}/.par1b,{x,0,3},{y,0,3},AxesOrigin->{0,0},VectorScale->0.1,AspectRatio->(1/GoldenRatio),PlotLabel->"1. B) eset"];

abra["1b"]

p1=pp[0.4,.8,1.1,par1b];
p2=pp[0.2,3,0.9,par1b];
p3=pp[1.2,0.3,0.9,par1b];
p4=pp[2,0.3,0.6,par1b];
p5=pp[0.8,2.5,0.45,par1b];
p6=pp[2,3,0.15,par1b];

Show[nullclines,sp,p1,p2,p3,p4,p5,p6,PlotLabel->"1. B) eset"]

Megjegyzések

Látszik, hogy (E^0),(E^1),(E^2) instabilis

Szemléletesen nyilvánvaló, hogy Det[A]<0, azaz -(Subscript[a, 12]) (Subscript[a, 21])+(Subscript[a, 11]) (Subscript[a, 22])<0 azt jelenti, hogy az intraspecifikus verseny kisebb, mint az interspecifikus kooperáció, azaz a fajok túlságosan segítik egymást, ezáltal túlszaporodnak.

Részletes vizsgálata az 2. B) esetben

ClearAll[par2b];
par2b={(Subscript[r, 1])->-1,(Subscript[r, 2])->-1,(Subscript[a, 11])->4,(Subscript[a, 12])->-2,(Subscript[a, 21])->-2,(Subscript[a, 22])->2};

b2eset

b2eset/.par2b

Valóban idetartozó esetet írnak le a paraméterek.

stac/.par2b

nullclines=Plot[{xnull,ynull}/.par2b,{x,0,4},PlotRange->{-0.1,3.5},AspectRatio->(1/GoldenRatio),Epilog->{Red,PointSize[0.02],Point/@({x,y}/.stac/.par2b)},PlotLabel->"2. B) eset"];

nullclinesfull=Plot[{xnull,ynull}/.par2b,{x,0,4},PlotRange->{0,4},AspectRatio->(1/GoldenRatio),Epilog->{Red,PointSize[0.02],Point/@({x,y}/.stac/.par2b),{Black,Text[((Subscript[r, 2])/(Subscript[a, 21])),{0.7,0.4}]},{Black,Text[((Subscript[r, 1])/(Subscript[a, 12])),{0.18,0.6}]}},PlotLabel->"2. B) eset"]

vp=VectorPlot[{f[x,y],g[x,y]}/.par2b,{x,0,4},{y,0,4},AxesOrigin->{0,0},VectorScale->0.1,AspectRatio->(1/GoldenRatio),PlotLabel->"2. B) eset"];

sp=StreamPlot[{f[x,y],g[x,y]}/.par2b,{x,0,4},{y,0,4},AxesOrigin->{0,0},VectorScale->0.1,AspectRatio->(1/GoldenRatio),PlotLabel->"2. B) eset"];

abra["2b"]

p2=pp[0.4,0.2,0.2,par2b];
p3=pp[1.2,0.3,4,par2b];
p4=pp[2,0.3,4,par2b];
p5=pp[2,2,2.1,par2b];
p6=pp[.3,1.,1.5,par2b];
p7=pp[0.4,3,2.2,par2b];
p8=pp[2.4,1.2,1.8,par2b];
p9=pp[2.4,1.5,2.3,par2b];
p10=pp[3.4,1.4,1.4,par2b];

Show[nullclines,sp,p1,p2,p3,p4,p5,p6,p7,p8,p9,p10,PlotLabel->"2. B) eset"]

stac/.par2b

Megjegyzések

Látszik, hogy (E^0) aszimptotikusan stabilis, de nem globálisan, a többi egyensúlyi helyzet pedig kívül esik az első negyeden.

Szemléletesen nyilvánvaló, hogy Det[A]<0, azaz (Subscript[a, 11]) (Subscript[a, 22])-(Subscript[a, 12]) (Subscript[a, 21])<0 azt jelenti, hogy akkor halnak ki a fajok, ha az egyik kooperáló létszáma eredendően kevés ahhoz, hogy a másik faj számára segítséget nyújtson, illetve túlszaporodnak, ha eredendően mindkét faj létszáma elég nagy volt, mert ekkor túlságosan segítik egymást.