Ugrás a tartalomhoz

Impulzív jelenségek modelljei

Karsai János

Typotex

12. fejezet - Impulzusok populációdinamikai modellekben

12. fejezet - Impulzusok populációdinamikai modellekben

12.1. A Malthus modell szabályozása

Tekintsünk egy populációt (állat, növény stb.), jelölje x(t)

x ? ( t ) a populáció egyedeinek számát. Tegyük fel, hogy a populációt egy egyed sem hagyja el, nincs halálozás, és a szaporodásnak a környezet sem szab határt. Az x(t) x ? ( t ) értékei nyilván egész számok, de nagy elemszám esetén folytonosnak és differenciálhatónak tekinthetjük. Tegyük fel továbbá, hogy a szaporodás egy pillanat alatt végbemegy, nincs eltolódás a fogamzás és a szaporodás ideje között. Minden új egyed azonos időnként születik. Ekkor a születések időbeli eloszlása egyenletes. Ezért a népesség számának megváltozása ?t ? t idő alatt ?x=?x?t ? x = ? ? x ? ? t , ahol ? statisztikusan meghatározott arányossági tényező (változási ráta). Ha ?t›0 ?t 0 , az ismert

x'=?x, x ' = ? ? x , (12.1.1)

egyenletet kapjuk.

Ha az adott populációban van születés és halálozás, akkor x(t)

x ? ( t ) változására ?x=?S–?H ? x = ? S - ? H , ahol ?S ? S és ?H ? H a ?t ?t alatti születések és halálozások száma. Innen a ?t idő alatti átlagos változás

(?x/?t)=(?S/?t)–(?H/?t)
? ? x ? ? t = ? ? S ? ? t - ? ? H ? ? t

Ha a populáció elég nagy, és feltesszük, hogy ?S=?x?t

? S = ? ? x ? ? t , ?H=?x?t ? H = ? ? x ? ? t , akkor ?t › 0 ?t 0 esetén kapjuk, hogy

x'=S'–H', S'=?x és H'=µx,
x ' = S ' - H ' , S ' = ? ? x és H ' = µ ? x ,

ahonnan

x'=(?–µ)x.
x ' = ( ? - µ ) ? x .

A fenti érvelések alapján a Malthus modell szerint változó populáció változását ugyanaz a homogén lineáris

x'=kx, x ' = k ? x , (12.1.2)

egyenlet írja le, amely a gyógyszerfelszívódásnál is fontos szerepet játszott. Ha k///>///0

k ///>/// 0 , akkor x x exponenciálisan növekszik, ha k///<///0, k ///</// 0 , akkor csökken. Az egyenlet megoldása x(t)=x_0e^(kt) x ? ( t ) = x 0 ? ? k ? t , ahol x_0?0 x 0 ? 0 az egyedek száma a kezdeti pillanatban. Negatív kezdeti értéknek itt nincs értelme.

Bár jól ismert az egyenlet megoldása, ábrázoljuk az iránymezőt és egy megoldást is az x'=x, x(0)=0.1

x ' = x , x ? ( 0 ) = 0.1 esetben.

Adott populációk esetén a k változási rátát mérési adatok alapján görbeillesztéssel határozzák meg. Lásd a könyv bővített CD-változatának M2.5 fejezetét.

A populáció szabályozása

Ezt a populációt szabályozni kell. A k///>///0

k ///>/// 0 esetben cél x(t) x ? ( t ) adott korlát alatt tartása, azaz az egyedek számát egyenletesen vagy időnként csökkenteni kell. Ezt nevezzük lehalászásnak. A k///<///0 k ///</// 0 esetben pedig x(t) x ? ( t ) -nek adott minimális érték felett kellene maradnia, vagyis a csökkenést pótolni kell. Két kézenfekvő lehetőséggel, az állandó sebességű és impulzív szabályozással foglalkozunk. Előbb áttekintjük a különböző eseteket, majd a további pontokban kísérleteket végzünk.

  • A k///<///0

    k ///</// 0 eset: betelepítés kihaló populáció esetén

Egyszerű szabályozás a folyamatos B0(t,x)

B0 ( t , x ) sebességű betelepítés. Ha ez függ az x-től, akkor a szabályozás figyeli az egyedek számát, egyébként pedig rögzített szabály szerint történik. Ekkor a változás egyenlete:

x'=–|k| x+B0(t,x). x ' = - | k | x + B0 ( t , x ) . (12.1.3)

Impulzív szabályozás esetén valamely t_i

t i időpillanatokban (i=1,2,...) ( i = 1 , 2 , ... ) növeljük meg az egyedek számát valamely B_i(x(t_i–0))///>///0 B i ( x ? ( t i - 0 ) ) ///>/// 0 értékkel. Ekkor az

x^'=–|k|x, t?t_i, x(t_i+0)=x(t_i–0)+B_i x ' = - | k | x , t ? t i , x ? ( t i + 0 ) = x ? ( t i - 0 ) + B i (12.1.4)

impulzív rendszert kapjuk. Ilyen eset például az ivadékok rendszeres kihelyezése halastavakba.

Végezhetünk állapotfüggő szabályozást is. Ha az egyedek száma a kritikus x_(min)

x min értéket eléri, akkor történjék betelepítés. Ekkor a rendszer az alábbi alakú:

x^'=–|k|x, (x(t)///>///x_(min)) x(t+0)=x(t–0)+B, ha x(t–0)?x_(min). x ' = - | k | x , ( x ( t ) ///>/// x min ) x ? ( t + 0 ) = x ? ( t - 0 ) + B , ha x ( t - 0 ) ? x min . (12.1.5)

Vegyük észre, hogy pontosan ugyanezeket a rendszereket kaptuk az egyszerű egyrekeszes gyógyszeradagolási modelleknél is.

  • A k///>///0

    k ///>/// 0 eset: lehalászás exponenciálisan növekvő populáció esetén

A k///<///0

k ///</// 0 esethez hasonló szabályozási stratégiákat követünk. A túlszaporodás elkerülésére csökkentsük az egyedek számát folyamatosan B0(t,x) B0 ( t , x ) sebességgel. A tavakban való folyamatos horgászat ilyen szabályozás. Ekkor a változás egyenlete:

x'=kx–B0(t,x) x ' = k ? x - B0 ( t , x ) (12.1.6)

Impulzív szabályozás esetén valamely t_n

t n időpillanatokban (n=1,2,...) ( n = 1 , 2 , ... ) csökkentsük az egyedek számát valamely B_n(x(t_n–0))///>///0 B n ( x ? ( t n - 0 ) ) ///>/// 0 értékkel. A halastavak időszakonkénti lehalászása felel meg ennek a modellnek. Ha B_n(x(t_n–0))?x(t_n–0) B n ( x ? ( t n - 0 ) ) ? x ? ( t n - 0 ) , akkor a populáció kihal. Ekkor az

x^'=kx, t?t_i, x(t_i+0)=(x(t_i–0)–B_i(x(t_i–0)))_+ x ' = k ? x , t ? t i , x ? ( t i + 0 ) = ( x ? ( t i - 0 ) - B i ( x ? ( t i - 0 ) ) ) + (12.1.7)

impulzív rendszerhez jutunk, ahol (.)_+

( . ) + valamely szám pozitív része. Emiatt a rendszer nemlineárissá válik.

Az állapot figyelése esetén akkor történik szabályozás, ha az egyedek száma elér egy kritikus x_(max)

x max értéket. Ekkor a rendszer az alábbi alakú:

x^'=kx, x(t)///<///x_(max), x(t+0)=(x(t–0)–B(x(t–0)))_+, ha x(t–0)?x_(max). x ' = k ? x , x ? ( t ) ///</// x max , x ? ( t + 0 ) = ( x ? ( t - 0 ) - B ? ( x ? ( t - 0 ) ) ) + , ha x ( t - 0 ) ? x max . (12.1.8)

Egyenletes lehalászás

Most tekinstük az egyenletes sebességű lehalászást, amelyet az

x^'=kx–B0
x ' = k ? x - B0

inhomogén lineáris egyenlet ír le, ahol k///>///0

k ///>/// 0 és B0///>///0. B0 ///>/// 0. Az x(0)=x0 x ? ( 0 ) = x0 kezdőértékhez tartozó megoldás:

xdot={kx–B0};

xdot = { k ? x - B0 } ;

thsol[t_]=x[t]/.DSolve[ODEICGen[xdot,{x0},{0,t},{x}],x[t],t][[1]]

thsol [ t_ ] = x [ t ] /. ? DSolve [ ODEICGen [ xdot , { x0 } , { 0 , t } , { x } ] , x [ t ] , t ] [ 1 ]

(B0–B0e^(kt)+e^(kt)kx0/k)

B0 - B0 ? ? k ? t + ? k ? t ? k ? x0 k

Az egyenlet egyensúlyi helyzete pedig:

Solve[xdot==0,x]

Solve [ xdot == 0 , x ]

{x›(B0/k)}

{ x B0 k }

Ez az egyensúlyi helyzet instabilis (használjuk az (x–(B0/k))^2

( x - B0 k ) 2 Ljapunov függvényt). Az x0///>///(B0/k) x0 ///>/// B0 k tulajdonságú megoldások exponenciálisan növekednek. Ha x0///<///(B0/k) x0 ///</// B0 k , akkor a megoldások exponenciálisan csökkennek, és a

{t›–(1/k)log(1–(kx0/B0))}
{ t - 1 k ? log ? ( 1 - k ? x0 B0 ) }

pillanatban a populáció kihal. Ez a szabályozási mód nem működik elvárásainknak megfelelően. Azt sem várhatjuk, hogy az ezzel analóg impulzív rendszer rögzített pillanatokban ható impulzusokkal is hatékony legyen. Ábrázoljuk a fenti program segítségével az iránymezőt és néhány megoldást az alábbi paraméterekkel:

parm={k›1,B0›1}; x0rul={{x0›0.9},{x0›1.1}};

parm = ? { k 1 , B0 1 } ; x0rul = { { x0 0.9 } , { x0 1.1 } } ;

plv0=PlotVectorField[{1,xdot[[1]]}/.parm,{t,0,2},{x,0,2}];

plv0 = PlotVectorField [ { 1 , xdot [ 1 ] } /. ? parm , { t , 0 , 2 } , { x , 0 , 2 } ] ;

pls0=Plot[Evaluate[thsol[t]/.parm/.x0rul],{t,0,2}];

pls0 = Plot [ Evaluate [ thsol [ t ] /. ? parm /. ? x0rul ] , { t , 0 , 2 } ] ;

Show[plv0,pls0,Axes›True,AspectRatio›0.5];

Show [ plv0 , pls0 , Axes True , AspectRatio 0.5 ] ;

Kvadratikus lehalászás: logisztikus növekedés

Láttuk, hogy a konstans sebességű lehalászás nem alkalmas a Malthus modell szerint változó populáció szabályozására. Ismerve a lineáris inhomogén egyenletek megoldásainak általános alakját, könnyű látni, hogy hasonló problémával találkoznánk, ha a szabályozás az idő függvénye lenne (ellenőrizzük kísérletekkel is). Ha a szabályozó függvény B0 x(t)

B0 ? x ? ( t ) alakú, akkor az egyenlet alakja x^'=(k–B0)x x ' = ( k - B0 ) ? x , ami szintén nem képes szabályozni. Most tekinstünk kvadratikus sebességű (vagy magasabb hatvány szerinti) lehalászást, amelyet az

x^'=kx–B0x^2
x ' = k ? x - B0 ? x 2

egyenlet ír le (k///>///0 és B0///>///0)

( k ///>/// 0 és B0 ///>/// 0 ) . Ez az egyenlet

x^'=kx(1–(B0/k)x) x ' = k ? x ? ( 1 - B0 k ? x ) (12.1.9)

alakban is írható. Ezzel olyan szabályozást definiálunk, amely matematikailag azonos a logisztikus szaporodási folyamatot leíró rendszerrel. Az ilyen rendszerekkel a következő pontban foglalkozunk. Írjuk fel az általános megoldást, és majd a következő pontban megkeressük az ezzel analóg impulzív szabályozást. Az előző megoldóprogram utasításait használjuk.

xdot={kx–B0x^2};

xdot = { k ? x - B0 ? x 2 } ;

(e^(kt)kx0/k+B0e^(kt)x0–B0 x0)

? k ? t ? k ? x0 k + B0 ? ? k ? t ? x0 - B0 ? x0

Az egyenlet egyensúlyi helyzetei pedig:

{{x›0},{x›(k/B0)}}.
{ { x 0 } , { x k B0 } } .

Ábrázoljuk az iránymezőt és néhány megoldást az alábbi paraméterekkel:

parm={k›1,B0›1}; x0rul={{x0›0.2},{x0›1.2}};

parm = ? { k 1 , B0 1 } ; x0rul = { { x0 0.2 } , { x0 1.2 } } ;

Láthatjuk, és a stabilitási kritériumokkal igazolhatjuk, hogy az x=0

x = 0 egyensúlyi helyzet instabilis, az x=(k/B0) x = k B0 pedig aszimptotikusan stabilis, és minden pozitív megoldás ehhez tart, ha t›?. t ? . Éppen ezt szerettük volna elérni a szabályozással.

Impulzív lehalászás: korlátos növekedés

Keressünk olyan impulzív rendszert, mely az előző pontban leírt szabályozó rendszerhez hasonló eredményt ad. A kx–B0x^2

k ? x - B0 ? x 2 sebesség két egymást követő impulzus között eredményezzen ugyanakkora változást a közönséges differenciálegyenlet megoldásainál, mint az impulzív rendszer. Ha t_i=iT t i = i ? T és T elegendően kicsi, akkor az

?_(t_i)^(t_(i+1))(kx(t)–B0x^2(t))dt??_(t_i)^(t_(i+1))kx(t)dt–B_i(x(t_i))
? t i t i + 1 ( k ? x ? ( t ) - B0 ? x 2 ( t ) ) ? ? t ? ? t i t i + 1 k ? x ? ( t ) ? ? t - B i ( x ? ( t i ) )

közelítő egyenlőségnek kellene teljesülni a B_i

B i impulzusra. Ha az impulzusok közt eltelt idő elegendően kicsi és állandó, akkor némi számolás és becslés után az alábbi rendszerhez jutunk:

x^'=kx, t?iT, x(iT+0)=x(iT–0)–B1x^2(iT–0),
x ' = k ? x , t ? i ? T , x ? ( i ? T + 0 ) = x ? ( i ? T - 0 ) - B1 ? x 2 ( i ? T - 0 ) ,

valamilyen B1(?B0)

B1 ? ( ? B0 ) konstanssal. Az impulzus alakjából azonnal látható, hogy nagy x értékek esetén negatív megoldást eredményezne, ami ebben a modellben a kihalást jelentené. Ez a példa is mutatja, hogy nemlineáris rendszerek esetén adott közönséges modellhez nem mindig egyszerű analóg impulzív rendszert találni.

Próbáljunk olyan impulzív szabályozást bevezetni, ami a logisztikus növekedéshez hasonló megoldásokat ad. Vegyük észre, hogy a közönséges (12.1.9) egyenletben kis x-ek esetén a kvadratikus szabályozás hatását nagy x-ek esetén fejti ki úgy, hogy minden pozitív megoldás korlátos marad és egy aszimptotikusan stabilis egyensúlyi helyzethez tart.

Tegyük fel, hogy egyenlő időközökben ugyanaz az impulzusfüggvény hat, vagyis B_i(.)=B(.)

B i ( . ) = B ? ( . ) . Biztosan korlátos megoldásokat kapunk, ha a megoldások nagy értékeire teljesül az

x(iT+0)=B(x(iT–0))?x((i–1)T+0)
x ? ( i ? T + 0 ) = B ( x ? ( i ? T - 0 ) ) ? x ? ( ( i - 1 ) ? T + 0 )

feltétel, amit egy adott impulzus függvény esetén ellenőrizni tudunk, mivel

x(iT–0)=x((i–1)T+0)e^(kT).
x ? ( i ? T - 0 ) = x ? ( ( i - 1 ) ? T + 0 ) ? ? k ? T .

Az alábbiakban két speciális alakú impulzust tekintünk. A rendszer tulajdonságainak formális ellenőrzését az olvasóra bízzuk.

  • Hatványimpulzus: x(iT+0)=B0 x^ß(iT–0)

    Hatvány ? impulzus : x ( i ? T + 0 ) = B0 ? x ß ( i ? T - 0 )

Alkalmazhatunk

x(iT+0)=B0x^ß(iT–0)
x ? ( i ? T + 0 ) = B0 ? x ß ( i ? T - 0 )

alakú impulzust, ahol 0///<///ß///<///1

0 ///</// ß ///</// 1 . Az olvasó ellenőrizheti és kísérleteink is megerősítik, hogy ez az impulzus teljesíti a megoldások korlátosságát biztosító feltételt. Ekkor kis egyedszám esetén is beleavatkozunk a folyamatba. Legyen ß=0.5 ß = 0.5 . Az egyenlet és a megoldás paraméterei:

eqnparm={k›1.};Imax=50;Iparm={T›N[1.],B0›0.5};

eqnparm = { k 1. } ; Imax = 50 ; Iparm = { T N [ 1. ] , B0 0.5 } ;

xvar={x};xdot={kx};

xvar = { x } ; xdot = { k ? x } ;

tn=Table[iT,{i,0,Imax}];

tn = Table [ i ? T , { i , 0 , Imax } ] ;

Bn=Table[{B0},{i,0,Imax}]/.Iparm;

Bn = Table [ { B0 } , { i , 0 , Imax } ] /. ? Iparm ;

Imp[n_,tn_List,x_List]:=Bn[[n]]?x;

Imp [ n_ , tn_List , x_List ] := Bn [ n ] ? x ;

Oldjuk meg a rendszert és ábrázoljuk a megoldásokat:

x0tab={{0.01},{0.5},{2.5}};

x0tab = { { 0.01 } , { 0.5 } , { 2.5 } } ;

t0=0; t1=30; x1=0;x2=3;

t0 = 0 ; t1 = 30 ; x1 = 0 ; x2 = 3 ;

Sol=IDESolve[xdot/.eqnparm,xvar,tn/.Iparm,Imp, x0tab,{t,t0,t1}];

Sol = IDESolve [ xdot /. eqnparm , xvar , tn /. Iparm , Imp , x0tab , { t , t0 , t1 } ] ;

PlotSol=Plot[Evaluate[Sol],{t,t0,t1}];

PlotSol = Plot [ Evaluate [ Sol ] , { t , t0 , t1 } ] ;

Az ábra alapján sejthető, hogy a megoldások egy periodikus megoldáshoz tartanak, ha t›?

t ? (lásd a fejezet végén a feladatokat).

  • Logaritmikusimpulzus: x(iT+0)=B0Log(x(iT–0)+1)

    Logaritmikus ? impulzus : x ? ( i ? T + 0 ) = B0 ? Log ( x ? ( i ? T - 0 ) + 1 )

Alkalmazhatunk olyan impulzust is, amelynek hatása kis x értékekre elhanyagolható a szaporodás saját törvényei mellett, nagy x értékekre pedig erősen konkáv. Ilyen például a logaritmikus

x(iT+0)=B0 log_?(x(iT–0)+1)
x ? ( i ? T + 0 ) = B0 ? log ? ( x ? ( i ? T - 0 ) + 1 )

alakú impulzus, ahol 1///<///?

1 ///</// ? . Ez az impulzus is teljesíti a megoldások korlátosságát biztosító feltételt (bizonyítandó). Oldjuk meg a rendszert az előző eset kezdeti értékeivel. A paraméterek, a megoldó és ábrázoló utasítás azonosak az előzővel, csak az impulzusfüggvény alakja más:

Imp[n_,tn_List,x_List]:=Bn[[n]]Log[2,x+1];

Imp [ n_ , tn_List , x_List ] := Bn [ n ] ? Log [ 2 , x + 1 ] ;

Ábrázoljuk a megoldásokat:

Folyamatos figyelés, diszkrét szabályozás

Próbálkozzunk a legegyszerűbb technikával, a folyamatos állapotfigyeléssel. Ha az egyedek száma elér egy maximális megengedett értéket, akkor csökkentjük valamekkora 0///<///B///<///x_(max)

0 ///</// B ///</// x max mennyiséggel. A rendszer alakja:

x^'=kx, x(t)///<///x_(max), x(t+0)=(x(t–0)–B(x(t–0)))_+, ha x(t–0)?x_(max)
x ' = k ? x , x ( t ) ///</// x max , x ? ( t + 0 ) = ( x ? ( t - 0 ) - B ? ( x ? ( t - 0 ) ) ) + , ha x ( t - 0 ) ? x max

Ez az eljárás műszaki berendezéseknél könnyen automatizálható, az ökológiában azonban szervezeti hátteret kell kiépíteni. Sokkal egyszerűbben megvalósíthatók a rögzített algoritmusú, folyamatos figyelést nem igénylő eljárások. Írjuk le a rendszert számítógéppel.

var={x}; xdot={kx}; parm={k›1.,xmax›2., B›1.};

var = { x } ; xdot = { k ? x } ; parm = { k 1. , xmax 2. , B 1. } ;

Impulse={{x–xmax,{x–B},0}};

Impulse = { { x - xmax , { x - B } , 0 } } ;

t0=0;t1=5;dt=0.01; x0list={{0.5}};

t0 = 0 ; t1 = 5 ; dt = 0.01 ; x0list = { { 0.5 } } ;

sol=IDERKSolve[xdot/.parm,Impulse/.parm,var,x0list, {t,t0,t1,dt}];

sol = IDERKSolve [ xdot /. ? parm , Impulse /. ? parm , var , x0list , { t , t0 , t1 , dt } ] ;

plv1=PlotVectorField[Flatten[{1,xdot}]/.parm,{t,t0,t1},{x,0,2.2}];

plv1 = PlotVectorField [ Flatten [ { 1 , xdot } ] /. ? parm , { t , t0 , t1 } , { x , 0 , 2.2 } ] ;

pls1=ListPlot1[sol,PlotJoined›True];

pls1 = ListPlot1 [ sol , PlotJoined True ] ;

Show[plv1,pls1,Axes›True];

Show [ plv1 , pls1 , Axes True ] ;

Nyilvánvaló, hogy a változás periodikus az első szabályozó pillanattól kezdve. Innen kezdve az egyedek átlagos száma

atlag=Simplify[(k/Log[xmax]–Log[xmax–B])?_((Log[xmax–B]/k))^((Log[xmax]/k))e^(kt)dt]

atlag = Simplify [ k Log [ xmax ] - Log [ xmax - B ] ? ? Log [ xmax - B ] k Log [ xmax ] k ? k ? t ? ? t ]

(B/Log[xmax]–Log[–B+xmax])

B Log [ xmax ] - Log [ - B + xmax ]

Az adott paraméterek mellett:

atlag/.parm

atlag /. parm

1.4427

1.4427

Feladatok

12.1.1. Feladat.

Bizonyítsuk be, hogy a rögzített pillanatokban végzett impulzív szabályozások esetén minden megoldás korlátos és aszimptotikusan T-periodikus. A periódussal való eltolástól eltekintve létezik pontosan egy aszimptotikusan stabilis pozitív periodikus megoldás. Határozzuk meg a periodikus megoldás átlagos értékét (lásd a 11.2. fejezetben az ismételt gyógyszeradagolásnál a periodikus megoldás létezésére vonatkozó gondolatmenetet).

12.1.2. Feladat.

Keressünk más hasonló hatású impulzusokat a Malthus modell szabályozására.