Logistična diferencialna enačba DIPLOMSKO DELO

UNIVERZA V LJUBLJANI PEDAGOŠKA FAKULTETA

Sabina Šubic

Logistična diferencialna enačba DIPLOMSKO DELO

Ljubljana, 2015

UNIVERZA V LJUBLJANI PEDAGOŠKA FAKULTETA

Študijski program: Dvopredmetni učitelj

Kandidatka: Sabina Šubic

Mentor: izr. prof. dr. Marko Slapar Somentor: asist. dr. Tadej Starčič

Logistična diferencialna enačba DIPLOMSKO DELO

Ljubljana, 2015

Zahvala

Iskrena hvala mentorju dr. Marku Slaparju in somentorju dr. Tadeju Starčiču, da sta me sprejela pod svoje mentorstvo, mi nudila strokovno pomoč in nasvete pri pisanju diplomskega dela.

Iskrena hvala družini in fantu Blažu za vso podporo in spodbujanje v času študija.

Hvala, da ste me vedno sprejemali, podpirali in razumeli. Vedite, da mi brez vas ne bi uspelo!

Hvala Maji za lektoriranje mojega diplomskega dela in Simoni za pomoč pri prevodu v angleški jezik.

Še enkrat hvala vsem!

Povzetek

V diplomskem delu bomo obravnavali logistično diferencialno enačbo, ki je uporabna pri obravnavi modelov rasti populacij določenih vrst in na številnih drugih področjih, kot bo predstavljeno preko primerov. Začeli bomo z eksponentnim modelom rasti in pokazali pomanjkljivosti le-tega pri rasti populacij zaradi neupoštevanja omeji- tvenih dejavnikov. Tako bomo postopoma prišli do uvedbe logistične diferencialne enačbe, ki upošteva omejitvene dejavnike. V zadnjem delu bomo vpeljali še logi- stično diferencialno enačbo s kritičnim pragom, ki upošteva še minimalno velikost populacije, da se le-ta lahko uspešno širi.

Ključne besede: model rasti, eksponentna rast, logistična diferencialna enačba, logistična rast, stopnja rasti, nosilna zmogljivost okolja, rast populacije, kritični prag, logistična diferencialna enačba s kritičnim pragom.

Abstract

In this diploma thesis we will discuss the logistic differential equation, which is useful when dealing with models of growth of populations of certain species, as well as in many other areas, which will be presented through examples. We will start with the exponential growth model and show its weakness due to disregarding the limiting factors. This is how we will gradually come to the introduction of the logistic differential equation, which takes into account the limiting factors. In the last section we will introduce a logistic differential equation with a critical threshold, which takes into account the minimum population size in order to spread successfully.

Keywords: model of growth, exponential growth, logistic differential equation, logistic growth, rate of growth, environmental carrying capacity, population growth, critical threshold, logistic differential equation with a critical threshold.

Kazalo

Slike

Poglavje 1. Uvod 1

Poglavje 2. Eksponentna rast 3

2.1. Naravna (eksponentna) rast 3

2.2. Primeri eksponentne rasti 4

Poglavje 3. Logistična rast 8

3.1. Logistična diferencialna enačba 8

3.2. Primeri uporabe logistične diferencialne enačbe 13

Poglavje 4. Modeliranje s kritičnim pragom 21

4.1. Kritični prag 21

4.2. Primeri rasti s kritičnim pragom 24

4.3. Logistična diferencialna enačba s kritičnim pragom 24

Literatura 28

Slike

1 Eksponentna rast: y(t) za ^dy_dt =ry (r >0). 4

2 f(y) v odvisnosti od y za ^dy_dt =r1− _K^yy. 9 3 Fazna premica diferencialne enačbe ^dy_dt =r1− _K^yy. 10 4 Graf y(t) logistične diferencialne enačbe ^dy_dt =r1− _K^yy. 10 5 Graf _K^y(t) za naravno rast rib na območju Tihega oceana. 14

6 Graf f(y)za ^dy_dt =−r1− _T^yy. 21

7 Fazna premica rasti s kritičnim pragom. 22

8 Rast s kritičnim pragom: ^dy_dt =−r1− _T^yy. 22

9 Graf f(y)za ^dy_dt =−r1− _T^y1−_K^yy. 25

10 Fazna premica za logistični model rasti s kritičnim pragom. 25 11 Graf y(t) za logistični model rasti s kritičnim pragom: ^dy_dt =

−r1− _T^y1− _K^yy. 26

POGLAVJE 1

Uvod

Področje diferencialnih enačb v matematiki je zelo uporabno na področju eko- logije za modeliranje rasti prebivalstva, v medicini za modeliranje rasti tumorja, v kemiji, fiziki in številnih drugih znanostih. Rast populacije pogosto opišemo z dife- rencialno enačbo z začetnim pogojem. S funkcijot 7−→y(t) povemo, da je v času t velikost populacije enaka y(t). Kako hitro populacija raste v času t, pa pove odvod y⁰(t).

Diplomsko delo bomo začeli z eksponentnim modelom rasti, kjer je hitrost ra- sti populacije v času t kar premo sorazmerna z njeno velikostjo v tem času. Rešili bomo diferencialno enačbo eksponentnega modela rasti ter spoznali pomanjkljivosti tega modela za modeliranje rasti populacije in nujnost upoštevanja omejitvenih de- javnikov, kot so pomanjkanje hrane, prostora, smrtnost, bolezni. Omenjeni model namreč dobro opisuje dejansko rast populacije le za majhne t. Pogledali si bomo še nekatera področja uporabe eksponentnega modela, zlasti tista, kjer eksponentni model precej realno opisuje določeno rast.

V osrednjem delu bo obravnava logistične diferencialne enačbe. Konstanto r premega sorazmerja pri eksponentnem modelu rasti bomo zamenjali s tako funkcijo, da bo dobljena Verhulstova oziroma logistična diferencialna enačba upoštevala prej omenjene omejitvene dejavnike rasti. Populacija, katere rast opisuje Verhulstova diferencialna enačba, se bo namreč približevala vrednosti K, to je nivo nasičenja ali nosilna zmogljivost okolja. Sprva bomo logistično diferencialno enačbo obrav- navali le geometrijsko – že pri slednji bomo prišli do grafa y(t), ki prikazuje rast populacije v odvisnosti od časa, nato pa jo bomo obravnavali še analitično, s čimer bomo potrdili ugotovitve pri geometrijski obravnavi. Pogledali si bomo še področja uporabe logističnega modela rasti, zlasti modeliranje rasti prebivalstva, širjenje nale- zljive bolezni ter uporabo enačbe pri rasti tumorja. Logistična diferencialna enačba namreč zelo dobro opisuje naraščanje populacije ob upoštevanju zgoraj omenjenih omejitvenih dejavnikov.

V zadnjem delu bomo geometrijsko in analitično obravnavali še diferencialno enačbo s kritičnim pragom. To je diferencialna enačba, ki vsebuje kritični prag T; če je začetna vrednost populacije pod kritičnim pragom T, se rast populacije ne pojavi, torej limitira proti 0; če pa je začetna vrednost populacije nad kritičnim pragomT, je rast neomejena. Opisali bomo še primer uporabe diferencialne enačbe s kritičnim pragom. Na koncu bomo modificirali logistično diferencialno enačbo tako, da bomo dobili logistično diferencialno enačbo s kritičnim pragom, ki bo upoštevala tako omejitvene dejavnike kot potrebno velikost začetne populacije, da se lahko pojavi uspešno širjenje le-te.

V diplomskem delu bomo uporabljali diferencialne enačbe 1. reda z ločljivima spremenljivkama, v katerih se neodvisna spremenljivka ne pojavi eksplicitno (avto- nomne enačbe). Zato si poglejmo še nekatere najpomembnejše definicije in izreke v

povezavi z diferencialnimi enačbami, ki jih bomo uporabljali v naslednjih poglavjih.

Definicije in izreki so v večini povzeti po[2], [4] in[10].

Najprej sploh definirajmo diferencialno enačbo. Navadna diferencialna enačba je zveza med neodvisno spremenljivkox, odvisno spremenljivkoy=y(x)in njenimi odvodiy⁰(x), . . ., y⁽ⁿ⁾(x), torej zveza

F(x, y(x), y⁰(x), . . . , y⁽ⁿ⁾(x)) = 0,

kjer jen ∈N, D ⊂ Rⁿ⁺² odprta množica in F :D → R zvezna funkcija. Običajno namestoy(x), y⁰(x), . . . , pišemoy, y⁰, . . .Red diferencialne enačbe je red najvišjega odvoda v enačbi.

Splošno diferencialno enačbo prvega reda podamo v obliki F(x, y, y⁰) = 0.

Njena rešitev je običajno enoparametrična družina krivuljy=f(x, C), kjer jeC∈R parameter. Partikularno rešitev dobimo tako, da predpišemo pogoj, ki mu mora ta rešitev zadoščati. Običajno podamo vrednost v neki točki, torejf(x₀) =y₀.

Diferencialni enačbi skupaj z začetnim pogojem F(x, y, y⁰) = 0, f(x₀) =y₀

pravimo začetni problem. Če je ^∂F_∂y0 6= 0, lahko enačbo F(x, y, y⁰) napišemo v obliki y⁰ =f(x, y). Omenimo še Picard-Lindelöfov izrek, ki pravi:

Naj bo dan naslednji problem z začetno vrednostjo

y⁰(t) = f(t, y(t)), y(t₀) =y₀, t∈[t₀−₀, t₀ +₀],

kjer je f Lipschitzovo zvezna glede na y in zvezna glede na t. Potem za neko vrednost,₀ ≥ >0, obstaja edinstvena rešitevy(t)k problemu začetne vrednosti na intervalu [t₀ −, t₀ +]. Teorem bo za nas pomemben pri obravnavi logistične diferencialne enačbe, saj nam pove, da skozi določeno točko poteka natanko ena rešitev.

V diplomski nalogi bomo uporabljali le diferencialne enačbe z ločljivima spre- menljivkama. EnačbaF(x, y, y⁰) = 0ima ločljivi spremenljivki, če jo lahko zapišemo v oblikiy⁰ =h(x)g(y). Ker jey⁰ = _dx^dy, lahko zapišemo _dx^dy =h(x)g(y) in dobimo

dy

g(y) =h(x)dx.

Splošna rešitev je torej enaka

Z dy g(y) =

Z

h(x)dx+C,

ki jo lahko zapišemo v oblikiG(y) =H(x)+C. Za partikularno rešitev, ki na primer ustreza pogojuy(x₀) = y₀, pa velja G(y₀) = H(x₀) +C, torej je

C =G(y₀)−H(x₀).

Že na začetku smo omenili, da bomo pri konceptu eksponentne in logistične rasti uporabljali diferencialno enačbo prvega reda, v kateri se neodvisna spremenljivka (t) ne pojavi eksplicitno (avtonomne enačbe). Njihova oblika je

dy

dt =f(y). (1)

2

POGLAVJE 2

Eksponentna rast

Thomas Malthus je leta 1798 objavil esej o rasti prebivalstva, v katerem je tr- dil, da prebivalstvo narašča geometrijsko, medtem ko se pridelava hrane povečuje aritmetično. Ta razkorak naj bi vodil v revščino in s tem v propad družbe, zato je podal predlog za omejevanje rasti prebivalstva. Uvod poglavja je povzet po [2] in [6].

Malthusov model rasti, ki ga imenujemo tudi eksponentna ali naravna rast, je model, kjer je hitrost spreminjanja neke količine premo sorazmerna s količino samo.

Če hočemo na primer priti do matematičnega (zveznega) modela za število zajcev na nekem otoku, pri čemer privzamemo, da je zelja in prostora, kolikor ga potrebujejo, s tem pridemo do eksponentne (naravne) rasti. Narava je seveda ne dovoljuje, kot bomo obrazložili pozneje. V tem poglavju se bomo posvetili matematični obravnavi eksponentnega modela rasti ter si pogledali nekaj primerov uporabe tega modela.

2.1. Naravna (eksponentna) rast

Posvetimo se matematičnemu opisu eksponentnega modela rasti. Obravnava je povzeta po[1].

Naj bo y = φ(t) populacija določene vrste ob času t. Najpreprostejša predpo- stavka glede spreminjanja populacije je torej ta, da je stopnja spreminjanja funkcije y premo sorazmerna trenutni vrednosti te funkcije. To pomeni

dy

dt =ry, (2)

kjer konstanto premo sorazmerja r imenujemo stopnja rasti/upadanja, odvisno od tega, ali je konstanta pozitivna ali negativna. Tu privzamemo, da je r > 0, torej populacija raste. Rešimo enačbo (2):

dy

dt = ry;

dy

y = r dt;

ln(y) = rt+c;

y(t) = Ce^rt, C ∈R. Če pri reševanju enačbe (2) upoštevamo začetni pogoj

y(0) =y₀, (3)

dobimo

y=y₀e^rt. (4)

V primeru, ko imamo r >0in y₀ >0, je

t→∞lim y(t) = lim

t→∞y₀e^rt=∞.

Matematični model torej predvideva nenehno rast populacije, kar lahko vidimo tudi na sliki 1, ki prikazuje graf rasti populacije (y) v odvisnosti od časa (t) pri različnih začetnih vrednostihy0.

Slika 1. Eksponentna rast: y(t)za ^dy_dt =ry (r >0).

Pri idealnih pogojih je enačba (4) razmeroma natančna za populacije mnogih vrst, vsaj za določen čas. Seveda je jasno, da takšni idealni pogoji ne vzdržijo za nedolo- čen čas. Sčasoma se pojavijo omejitve glede prostora, hrane oziroma drugi vzroki, ki zmanjšajo stopnjo rasti in ustavijo dotlej neovirano eksponentno rast. Zakon na- ravne rasti zato v daljšem časovnem obdobju ne ustreza nobenemu organizmu ali neki količini v naravi. Omejitvene dejavnike bomo podrobneje opisali v poglavju o logistični rasti. V nadaljevanju so opisani nekateri primeri uporabe eksponentnega modela.

2.2. Primeri eksponentne rasti

V nadaljevanju sta podrobneje opisana dva primera uporabe eksponentne rasti, nato so na kratko opisani še nekateri drugi primeri uporabe modela eksponentne rasti. Primera, ki sta obravnavana podrobneje, sta povzeta po[5].

1. primer: biologija – naravna rast populacije poljskih miši

Populacija poljskih miši naseljuje določeno podeželsko območje. Če ne upošte- vamo plenilcev, se vsak mesec populacija miši poveča za 50 %. Vendar na tem območju živi mnogo sov, ki ubijejo 15 miši na dan. Poiščimo diferencialno enačbo, ki opisuje velikost populacije v odvisnosti od časa.

Najprej identificirajmo spremenljivke: y naj bo velikost populacije miši int čas v mesecih. Eksponentno rast nam opisuje diferencialna enačba

dy

dt = 0,5y.

Ob upoštevanju plenilcev moramo odšteti še mesečno izgubo miši. Ker je vsak dan ubitih 15 miši, to pomeni 15·30 = 450 ubitih miši v enem mesecu. Dobimo torej

4

modificirano diferencialno enačbo eksponentne rasti dy

dt = 0,5y−450. (5)

Rešimo enačbo (5). Najprej ločimo spremenljivki dy

0,5y−450 =dt;

integriramo obe strani enačbe

2 ln|0,5y−450|=t+c in odtod dobimo

y=Ce^2t + 900. (6)

Če enačbo (6) analiziramo, ugotovimo naslednje:

• če je začetna vrednost y(0)>900, se y stalno povečuje;

• če je začetna vrednost y(0) = 900, bo y konstantno 900;

• če je začetna vrednost y(0)<900, bo y upadal proti 0 z naraščajočim t.

y= 900 je torej nestabilna ravnotežna rešitev.

Opomba: izračun velikosti populacije po eksponentnem modelu rasti je realen za majhne vrednosti t. Za večje vrednosti ta model ni primeren, ker ne upošteva šte- vilnih omejitvenih dejavnikov, kot bomo videli v poglavju 3.

2. primer: fizika – gravitacija

Po principu eksponentne rasti s časom narašča tudi hitrost predmeta, ki pada proti tlom. Poglejmo si konkreten primer, kjer moramo poleg teorije diferencialnih enačb upoštevati tudi Newtonov zakon. Denimo, da predmet z neke višine pada proti morski gladini. Predpostavimo, da je zračni upor premo sorazmeren hitrosti s koeficientom zračnega upora 2 ^kg_s . Po kolikšnem času bo predmet dosegel tla in kakšna je njegova hitrost ob tem času, če je masa predmeta 10 kg in ga spustimo z višine 300 m?

Najprej identificirajmo spremenljivke: naj bo v hitrost in t čas. Priti hočemo do diferencialne enačbe, ki opisuje padanje predmeta. Fizikalni zakon, ki opisuje gibanje, je drugi Newtonov zakon F = ma. Vemo, da je pospešek a definiran kot a = ^dv_dt, torej je F =m^dv_dt. Na predmet delujeta dve sili: sila gravitacije Fg = mg, kjer je g = 9,81_s^m₂ , in sila upora, ki je po predpostavki enaka 2v. Ker ti dve sili delujeta v nasprotnih smereh, je rezultanta delujočih sil enaka razliki teh dveh sil.

Dobimo torej enačbo

mdv

dt =mg−2v.

Če upoštevamo še dane podatke, torej maso predmeta in silo gravitacije, dobimo dv

dt = 9,81−v

5. (7)

Ker je postopek reševanja diferencialne enačbe (7) analogen kot pri diferencialni enačbi (5), zapišimo le rešitev:

v = 49,05−Ce^−t5. (8)

Ob upoštevanju začetnega pogojav(0) = 0, saj predmet spustimo z določene višine, dobimo

v = 49,05−49,05e^−5t. (9) Če analiziramo dobljeno rešitev (9), ugotovimo, da se dolgoročno hitrost približuje vrednosti49,05^m_s. Če pa analiziramo splošno rešitev (8), ugotovimo naslednje:

• če je začetna hitrost med 0^m_s in49,05^m_s, bo hitrost s časom naraščala proti vrednosti 49,05 ^m_s;

• če je začetna hitrost49,05^m_s, bo hitrost imela konstantno vrednost49,05^m_s;

• če je začetna hitrost večja od 49,05 ^m_s, bo hitrost s časom padala proti vrednosti 49,05 ^m_s.

Torej je v = 49,05stabilna ravnotežna rešitev, ker se v primeru, ko t → ∞, hitrost približuje tej vrednosti. Da bomo lahko poiskali čas padca in hitrost ob tem času, moramo poiskati formulo, ki opisuje razdaljo x v odvisnosti od časa t. Vemo, da je v = ^dx_dt, zato iz enačbe (9) dobimo

dx= (49,05−49,05e^−5t)dt;

x= 49,05t+ 245,25e^−5t +c.

Ker koordinatni sistem izberemo tako, da gravitacija deluje v pozitivni smeri, torej navzdol, jexrazdalja predmeta od njegovega začetnega položaja in zato veljax(0) = 0. Ob upoštevanju začetnega pogoja dobimo iskano formulo

x= 49,05t+ 245,25e^−5t −245,25, (10) iz katere lahko s pomočjo katerega izmed matematičnih računalniških programov dobimo čas, pri katerem je prišlo do padca predmeta. Spomnimo se, da predmet spustimo z višine 300 m.

300 = 49,05t+ 245,25e^−5t −245,25;

t ∼= 10,51.

Hitrost ob tem času je

v(10,51)∼= 43,06 m s .

Ostali primeri eksponentne rasti

V nadaljevanju so na kratko opisani še nekateri primeri eksponentne rasti z raz- ličnih področij. Opisi so povzeti po [2] in[7].

6

• BIOLOGIJA: po principu naravne rasti se obnaša količina lesa v drevesu, dokler ne naleti na omejitve, ki te rasti ne dovolijo več. Rast lesa je nepre- kinjena (do omejitev), saj ne moremo reči, da na pol leta ali vsak mesec ali celo vsak dan v hipu zraste nekaj novega lesa. Vendar model dobro opisuje rast lesa le za manjše t. Vemo namreč, da so številni omejitveni dejavniki, ki take rasti ne dovoljujejo. Poleg tega se rast lesa enkrat ustavi, česar eksponentna rast ne upošteva.

• BIOLOGIJA: število mikroorganizmov se bo eksponentno povečevalo pri optimalnih pogojih (dovolj velika količina hranila, prostora . . . ). Prvi or- ganizem se namreč razdeli v dva hčerinska organizma, ki se zopet delita naprej – v naslednji fazi torej dobimo 4 mikroorganizme in tako dalje.

• BIOLOGIJA: virus (na primer črne koze) se eksponentno širi, če ni na voljo umetnega cepiva. Vsak človek lahko okuži več drugih ljudi.

• BIOLOGIJA: številni odzivi živih bitij na dražljaje, vključno s človeškim zaznavanjem, so logaritemski odzivi, ki so inverzni eksponentnim odzivom.

Glasnost in frekvenco zvoka zaznavajo logaritmično (celo komaj zaznavne dražljaje) v mejah zaznavanja. To je tudi razlog, da eksponentno povečanje svetlosti vizualnih dražljajev ljudje sprejemajo bolj kot linearno in ne kot eksponentno povečanje, kar je zelo pomembno za preživetje.

• EKONOMIJA: kot eksponentni model rasti se obnaša rast glavnice pri obre- stnem obrestovanju. Denimo, da vložimo na banko G₀ denarja, banka pa nam ponujap%letne obresti. Pon-letih imamo naslednjo količino denarja:

Gn=G0

1 + p 100

n

.

• FIZIKA: eksperimenti so pokazali, da je količina radioaktivne snovi, ki razpade na časovno enoto, premo sorazmerna s količino snovi, ki je še na razpolago v danem trenutku. Konstanta premega sorazmerja k je običajno določena z meritvami.

POGLAVJE 3

Logistična rast

Uvod poglavja je povzet po [3] in[7].

Eksponentno (naravno) rast populacije poljubne vrste sčasoma začne omejevati ko- ličina hrane, vode in drugih naravnih virov, pomanjkanje prostora za rast, plenilci, bolezni ali kakšen drug ekološki dejavnik. Če je rast omejena z viri, se eksponentna rast začne upočasnjevati, ko konkurenca za prej omenjene vire narašča. Sčasoma se rast populacije upočasni skoraj na 0, ko populacija doseže nosilnost okolja (K) za dano vrsto.

Leta 1838 je Pierre François Verhulst predstavil logistično rast prebivalstva, ki upošteva omejitve, ki onemogočajo eksponentno naraščanje populacije. Verhulstov model rasti pravi, da je rast populacije sorazmerna številu osebkov v populaciji in razpoložljivim virom, ki omejujejo rast.

V nadaljevanju bomo postopoma prišli do logističnega modela rasti, ki upo- števa zgoraj opisane omejitve. Do rešitev diferencialne enačbe bomo sprva prišli s pomočjo geometrijske obravnave diferencialne enačbe in tako pridobili kvalitativno pomembne informacije. Z njimi bomo nato prišli do iskanega grafa y(t), ki nam pove, kako se število populacije spreminja s časom. Diferencialno enačbo bomo na koncu še rešili in jo analitično obravnavali.

3.1. Logistična diferencialna enačba

Obravnava in reševanje logistične diferencialne enačbe je povzeto po [1].

Če upoštevamo dejstvo, da je stopnja rasti odvisna od populacije, zamenjamo kon- stanto r v enačbi (2) s funkcijo h(y) in s tem dobimo spremenjeno enačbo

dy

dt =h(y)y. (11)

Sedaj hočemo izbrati tako funkcijo h(y), da h(y) ∼= r > 0 za majhne vrednosti y, h(y) upada pri povečevanju y in h(y) < 0, ko je y dovolj velik. Najpreprostejša funkcija, ki vsebuje opisane lastnosti, je h(y) = r −ay, kjer je a neka pozitivna konstanta. Če uporabimo to funkcijo v enačbi (11), dobimo

dy

dt = (r−ay)y. (12)

Enačba (12) je znana kot Verhulstova enačba ali logistična enačba. Običajno je bolj priročno, da logistično enačbo zapišemo v ekvivalentni obliki

dy dt =r

1− y K

y, (13)

kjerK = ^r_a. Konstanto r imenujemo intrinzična (bistvena) stopnja rasti – to je sto- pnja rasti, če ni omejitvenih dejavnikov. Konstanto K bomo interpretirali kasneje, prav tako bomo kasneje raziskovali natančne rešitve enačbe (13).

8

3.1.1. Geometrijsko reševanje logistične diferencialne enačbe.

Poglejmo si, kako lahko enostavno narišemo kvalitativno pravilno skico rešitve dane diferencialne enačbe. Enake metode se uporabljajo tudi za bolj splošne enačbe, kot je na primer enačba (1).

Poiščimo najprej najbolj preproste rešitve enačbe (13), to so konstantne funkcije.

Za slednji tip funkcij je ^dy_dt = 0 za vsak t, torej mora konstantna rešitev enačbe (13) zadostiti naslednji enačbi:

r

1− y K

y= 0.

Vidimo, da sta konstantni rešitvi y = φ₁(t) = 0 in y = φ₂(t) = K. Ta tip rešitev imenujemo ravnotežne rešitve enačbe (13), saj vrednostiyne variirajo, kotnarašča.

Analogno lahko poiščemo ravnotežne rešitve bolj splošne enačbe ^dy_dt =f(y)tako, da nastavimo enačbof(y) = 0. Ničle funkcijef(y)imenujemo tudi kritične točke.

Da bi lažje prišli še do drugih rešitev enačbe (13) in skicirali grafe teh rešitev, začnemo z risanjem grafa f(y) = r1− _K^yy. Dobimo neko parabolo, kot lahko vidimo na spodnji sliki.

Slika 2. f(y) v odvisnosti od y za ^dy_dt =r1− _K^yy.

Presečišči z vodoravno osjo (0,0) in (K,0) ustrezata kritičnima točkama enačbe (13), vrh parabole je^K₂,^rK₄ .

Na grafu na sliki 2 opazimo, da je ^dy_dt > 0 za 0 < y < K; to pomeni, da y(t) narašča, ko je y v tem intervalu, kar lahko označimo tudi z desno usmerjenimi puščicami na omenjenem grafu. Če y > K, potem ^dy_dt < 0, zato je y(t) padajoča funkcija za y > K, kar lahko označimo z levo usmerjeno puščico na istem grafu. V omenjenem kontekstu zato y-os pogosto imenujemo fazna premica in jo ponavadi predstavimo v navpičnem položaju, kot kaže slika 3. Piki pri y = 0 in y = K predstavljata kritični točki oziroma ravnotežni rešitvi. Puščice označujejo (kot je že zgoraj omenjeno), da y(t) narašča, ko 0< y < K ter pada, ko y > K.

Slika 3. Fazna premica diferencialne enačbe ^dy_dt =r1− _K^yy.

Opazimo še (slika 2), da je naklonf(y)blizu0, ko jeyblizu0aliK. Naklon postaja večji, ko se vrednost y oddaljuje od okolice 0ali K.

Sedaj lahko skiciramo rešitve enačbe (13). Začnemo z ravnotežnima rešitvama y= 0 in y=K; nato narišemo ostale krivulje, ki naraščajo v intervalu

0 < y < K, padajo, ko y > K, ter se približujejo konstantni vrednosti, ko se y približuje vrednosti 0 ali K. Torej imajo grafi rešitev enačbe (13) splošno obliko, kot jo lahko vidimo na sliki 4, ne glede na vrednost r aliK.

Slika 4. Graf y(t) logistične diferencialne enačbe ^dy_dt =r1−_K^yy.

Na grafu se zdi, kot da bi druge rešitve lahko sekale ravnotežno rešitevy=K, vendar to ni mogoče zaradi Picard-Lindelöfovega izreka, ki smo ga omenili v uvodu. Čeprav so druge rešitve asimptotično približujejo ravnotežni rešitvi y = K, ko t → ∞, je ne morejo sekati za noben končen t. Posledično rešitev, ki se začne v intervalu 0< y < K, ostane v tem intervalu, podobno velja za rešitev, ki se začne v intervalu K < y <∞.

Sedaj lahko raziščemo še konkavnost oziroma konveksnost krivulj ter prevoje, in sicer z iskanjem drugega odvoda ^d_dt^2y2. Iz diferencialne enačbe (1) dobimo (s pomočjo verižnega pravila):

d²y dt² = d

dt dy dt = d

dtf(y) = f⁰(y)dy

dt =f⁰(y)f(y). (14) Graf funkcije y v odvisnosti od časa t je konveksen, ko y⁰⁰ > 0. To pomeni, da imata f in f⁰ enak predznak. Podobno je graf konkaven, ko y⁰⁰ <0, torej imata f

10

in f⁰ nasproten predznak. Predznake funkcij f in f⁰ lahko enostavno določimo iz grafa f(y) v odvisnosti od y. Do prevojev lahko pride, ko f⁰(y) = 0. V primeru enačbe (13) so krivulje rešitev konveksne za 0 < y < ^K₂, ker je f v tem intervalu pozitivna in narašča (slika 2), zato sta f in f⁰ pozitivna. Krivulje rešitev so prav tako konveksne za y > K, kjer je f negativna in padajoča funkcija, torej sta f in f⁰ obe negativni. V intervalu ^K₂ < y < K pa so krivulje rešitev konkavne, ker je f pozitivna in padajoča funkcija, torej je f pozitivna in f⁰ negativna. Do prevoja pride, ko graf y v odvisnosti od t prečka črto y = ^K₂, saj na grafu f(y) opazimo, da pri y = ^K₂ funkcija preide iz naraščajoče v padajočo – to pomeni, da na grafu y(t)funkcija preide iz konveksne krivulje v konkavno, torej je to točka prevoja. Vse opisane lastnosti grafaf(t)lahko opazimo na sliki 4.

Na koncu opazimo še, da je vrednost K zgornja meja, ki se ji vrednosti y(t) približujejo, a je ne dosežejo. V primeru, da je začetna vrednost y₀ < K, vrednosti y(t) naraščajo in se približujejo vrednosti K. V primeru, ko y₀ > K, pa vrednosti y(t)padajo protiK. ZatoK imenujemo nivo nasičenja ali nosilna zmogljivost okolja za dano vrsto.

Če primerjamo eksponentno (slika 1) in logistično rast (slika 4), ugotovimo, da so rešitve nelinearne enačbe (13) bistveno drugačne od rešitev linearne enačbe (1), še posebno za večje vrednostit. Ne glede na vrednostK (torej ne glede na to, kako majhen je nelinearen člen v enačbi (13)) rešitve te enačbe dosežejo končno vrednost, ko t → ∞, medtem ko rešitve enačbe (1) eksponentno rastejo brez zgornje meje.

Torej ima še tako majhen nelinearen člen v diferencialni enačbi (13) odločilen vpliv na rešitve za velike vrednosti t.

3.1.2. Analitično reševanje logistične diferencialne enačbe – logistična funkcija.

Velikokrat zadostuje, da imamo kvalitativno informacijo o rešitviy=φ(t)enačbe (13), kot na sliki 2. Te informacije smo pridobili zgolj iz grafaf(y), ne da bi pri tem rešili diferencialno enačbo (13). Če hočemo povsem natančen opis logistične rasti, na primer, če hočemo vedeti vrednost populacije ob točno določenem času, moramo rešiti enačbo (13) glede na začetni pogoj (3). Pri sklepu, da y6= 0 in y6=K, lahko zapišemo enačbo (13) v obliki

dy

1−_K^yy =r dt.

Najprej uporabimo razcep na parcialne ulomke, da bomo lahko integrirali ulomek

1

(¹⁻K^y)^y:

1

1− _K^yy = A

1− _K^y +B y; 1

1− _K^yy = Ay+B1−_K^y

1− _K^yy ; 1

1− _K^yy = B+A−_K^By

1− _K^yy .

Torej jeB = 1 inA = _K¹.

Sedaj lahko rešimo preoblikovano diferencialno enačbo (13):

1 y +

1 K

1− _K^y

!

dy=r dt;

z integriranjem obeh strani dobimo ln|y| −ln

1− y K

=rt+c. (15)

Odstranimo absolutne vrednosti, saj vemo, da če 0 < y₀ < K, je y ves čas v temu intervalu. Dobimo

y

1− _K^y =Ce^rt, (16)

kjer C = e^c in c poljubna konstanta integriranja, ki jo določimo iz začetnega po- goja y(0) = y₀. Sedaj moramo ugotoviti še vrednost konstante C, da zadostimo začetnemu pogojuy(0) =y₀:

y₀

1− ^y_K⁰ =Ce⁰; C = y₀

1− ^y_K⁰. (17)

Če upoštevamo (17) in izpostavimoy, dobimo

y= y₀K

y₀+ (K −y₀)e^−rt. (18)

Rešitev (18) smo dobili ob predpostavki, da 0 < y₀ < K, zato rešimo enačbo (15) še ob predpostavki, day₀ > K:

ln|y| −ln

1− y

K

=rt+c;

odstranimo absolutne vrednosti in pri tem upoštevajmo, day₀ > K:

ln(y)−ln

y K −1

=rt+c, kar je ekvivalentno

y

y

K −1 =Ce^rt. Upoštevamo še začetni pogojy(0) =y₀ in dobimo

C = y₀

y0

K −1. (19)

Ko izpostavimoy in vstavimo konstanto C (19), pridemo do enačbe (18). Torej je enačba (18) veljavna za oba intervala 0 < y < K ter y > K. Enačba (18) vsebuje tudi ravnotežni rešitvi y = φ₁(t) = 0 in y = φ₂(t) = K, ki zaporedoma ustrezata začetnim pogojemy₀ = 0iny₀ =K. Torej lahko vse kvalitativne informacije, ki smo jih na začetku pridobili s pomočjo geometrijske utemeljitve, potrdimo z analitično proučitvijo rešitve (18).

Če y₀ = 0, potem iz enačbe (18) dobimo y(t) = 0za vsak t. Če y₀ >0int → ∞ v enačbi (18), potem dobimo

t→∞lim y(t) = y₀K y₀ =K.

12

Za vsaky₀ >0se rešitev torej asimptotično približuje ravnotežni rešitviy =φ₂(t) = K, kot → ∞. Sedaj res vemo, da jeK meja, ki se ji populacija približuje ne glede na njeno začetno vrednost. Pravimo, da jeK asimptotično stabilna rešitev enačbe (13), oziroma da je točka y = K asimptotično stabilno ravnotežje ali kritična točka. Po daljšem obdobju je ne glede na začetno velikost populacije (dokler je pozitivna) le-ta blizu nivoja nasičenja K. Ostale rešitve z rastočim r hitreje dosegajo ravnotežno rešitev.

Situacija za ravnotežno rešitev y = φ1(t) = 0 pa je precej drugačna. Rešitve, ki se začnejo blizu 0, z rastočim t naraščajo in dosegajo K, ko t → ∞. Pravimo, da je φ₁(t) = 0 nestabilna ravnotežna rešitev ali da je y = 0 nestabilno ravnotežje ali kritična točka. Torej je edini način, da rešitev ostane blizu 0, ta, da je začetna vrednost natanko 0.

3.2. Primeri uporabe logistične diferencialne enačbe

V nadaljevanju je obravnavanih nekaj najbolj značilnih primerov uporabe logi- stične diferencialne enačbe na različnih področjih.

1. primer: modeliranje rasti skupne mase rib na določenem območju

Logistični model rasti je bil uporabljen za rast rib na določenih območjih Tihega oceana. Obravnava primera je povzeta po[1].

Naj bo y (v kilogramih) skupna masa rib ob časut. Parametri logistične enačbe so ocenjeni na naslednje vrednosti: r = 0,71/leto in K = 80,5×10⁶ kg. Naj bo začetna skupna masa rib y₀ = 0,25K in pokažimo, kolikšna je skupna masa rib 2 leti kasneje. Poleg tega ugotovimo, po kolikšnem časut je y(t) = 0,75K.

Za računanje bo lažje, če enačbo (18) zapišemo v naslednji obliki y

K =

y0

K

_y

0

K

+^h1− ^y_K⁰ⁱe^−rt. Uporabimo podatke, ki so nam na voljo in dobimo

y(2)

K = 0,25

0,25 + 0,75e^−1,42,

torej jey(2) = 46,7×10⁶ kg. Izračunati moramo še, po kolikšnem času t masa rib naraste na y(t) = 0,75K. Najprej iz enačbe (18) izpostavimo t:

e^−rt =

y0

K

1− _K^y

y K

1− ^y_K⁰; t = −1

r ln

y0

K

1− _K^y

y K

1− ^y_K⁰. Vstavimo vrednosti za r, ^y_K⁰ ter upoštevamo, da je _K^y = 0,75:

t=− 1

0,71ln0,25(1−0,75)

0,75(1−0,25) =− 1 0,71ln1

9 = 1

0,71ln 9 ∼= 3,095 let.

Na sliki 4 lahko vidimo graf _K^y(t)za dane vrednosti parametrov ter različne začetne pogoje. Črna krivulja predstavlja primer, ki smo ga ravnokar rešili(y₀ = 0,25K).

Slika 5. Graf _K^y(t) za naravno rast rib na območju Tihega oceana.

2. primer: modeliranje rasti prebivalstva

Leta 1800 je bilo v ZDA 5,3 milijona prebivalcev, leta 1850 pa 23,1 milijona.

Napovejmo, kakšno je bilo število prebivalcev leta 1950 ter kakšen je nivo nasičenja, t.j. maksimalna vrednost populacije, če je konstanta r = 0,031476. Obravnava primera je povzeta po [8].

Upoštevamo, da populacija narašča logistično, zato jo lahko opišemo s funkcijo (18)

y= y0K

y₀+ (K −y₀)e^−rt, iz katere izrazimoK:

K = yy₀e^−rt−yy₀ ye^−rt−y₀ , ki ga lahko izračunamo iz danih podatkov:

K = (23,1)10⁶(5,3)10⁶e(−0,031476)(50)−(23,1)10⁶(5,3)10⁶ (23,1)10⁶e(−0,031476)(50)−(5,3)10⁶ ;

K ∼= 189⁰419⁰781.

Sedaj lahko izračunamo še število prebivalcev leta 1950:

y(150) = (5,3)10⁶(189⁰419⁰781)

(5,3)10⁶+ (189⁰419⁰781−(5,3)10⁶)e(−0,031476)(150); y(150)∼= 144⁰675⁰180.

Leta 1950 naj bi bilo torej v ZDA glede na naš model okrog 144,7 milijona prebival- cev. V resnici jih je bilo 150,7 milijona, torej malo več od pričakovanega števila, in sicer zaradi tako imenovanega "baby boom-a".

14

3. primer: širjenje nalezljive bolezni

Opis primera uporabe je povzet po [1].

Uporaba matematičnih metod pri proučevanju širjenja nalezljivih bolezni sega vsaj že v leto 1760, ko je Daniel Bernoulli metode uporabljal pri proučevanju širjenja črnih koz. V zadnjih letih je bilo proučevanih in predlaganih več različnih matematičnih modelov za različne bolezni. Poglejmo si enega od predlaganih modelov.

Domnevamo, da dano populacijo lahko razdelimo na dva dela: tiste, ki imajo dano bolezen in lahko okužijo druge, ter tiste, ki bolezni nimajo, ampak so vseeno dovzetni za okužbo. Naj bo x delež posameznikov, dovzetnih za okužbo, iny delež okuženih posameznikov; potem jex+y= 1. Privzemimo, da se bolezen širi ob kon- taktu med zdravim in okuženim posameznikom populacije in da je stopnja širjenja

dy

dt premo sorazmerna številu takih kontaktov. Privzemimo še, da se posamezniki obeh skupin populacije prosto gibljejo, torej je število kontaktov premo sorazmerno produktu xin y. Ker je x= 1−y, dobimo naslednji problem z začetno vrednostjo:

dy

dt =αy(1−y), y(0) =y₀, (20)

kjer jeα pozitiven faktor premo sorazmerja iny₀ začetni delež okuženih posamezni- kov.

Ravnotežni rešitvi diferencialne enačbe (20) sta y =φ₁(t) = 0 in y=φ₂(t) = 1.

ZaK = 1v enačbi (13) dobimo ravno diferencialno enačbo (20). Ravnotežna rešitev y = φ1(t) = 0 je torej nestabilna, y = φ2(t) = 1 pa stabilna. Rešitev diferencialne enačbe (20) ob upoštevanju začetne vrednostiy(0) =y0 je

y= y₀

y0+ (1−y0)e^−αt (21)

in vsebuje tudi ravnotežni rešitvi y = φ₁(t) = 0 in y = φ₂(t) = 1, ki zaporedoma ustrezata začetnima pogojema y₀ = 0 in y₀ = 1. Iz enačbe (21) je razvidno, da delež okuženih posameznikov gre proti 1, ko t → ∞, torej se sčasoma okuži celotna populacija.

4. primer: modeliranje rasti tumorja

Opis primera uporabe je povzet po [3] in [9].

Ena od uporab logistične krivulje v medicini je uporaba logistične diferencialne enačbe pri modeliranju rasti tumorja. Ta uporaba se lahko obravnava kot razširitev zgoraj omenjene uporabe pri ekologiji. Če zy(t) označimo velikost tumorja ob času t, njegovo dinamiko opisuje diferencialna enačba

dy dt =r

1− y K

y, ki je enačba tipa

y⁰ =f(y)y, f⁰(y)≤0,

kjer jef(y)stopnja širjenja tumorja. Če se prične zdravljenje s kemoterapijo (log-kill effect), je enačbo morda treba spremeniti v naslednjo obliko

dy dt =r

1− y K

y−c(t)y,

kjer je c(t) stopnja smrtnosti (celic), ki jo povzroča terapija. V idealiziranem pri- meru zelo dolgega zdravljenja c(t) lahko modeliramo bodisi kot periodično funkcijo (obdobja T), v primeru neprekinjenega zdravljenja z infuzijo pa kot konstantno funkcijo, in če je

1 T

Z T 0

c(t)dt > r, je mogoče pokazati, da

t→+∞lim y(t) = 0.

To pomeni, da v primeru, ko je povprečna stopnja smrtnosti (celic) zaradi terapije večja od stopnje širjenja tumorja, sledi izkoreninjenje bolezni. Seveda smo prikazali poenostavljen model rasti tumorja in terapije. Rast tumorja bolj realno opisuje Gompertzova diferencialna enačba

dy

dt =rln K y

!

y,

ki je limitni primer splošne logistične diferencialne enačbe.

Leta 1960 je A. K. Laird prvič uspešno uporabil Gompertzovo krivuljo kot apro- ksimacijo podatkov o rasti tumorja. Tumorji so populacije celic, ki rastejo v ome- jenih prostorih, kjer je razpoložljivost hranil omejena. Pri tumorju število rakastih celic narašča najprej eksponentno, nato se rast umiri in sčasoma ustavi. S poskusi so ugotovili, da krivulja naraščanja ni logistična, temveč Gompertzova krivulja. Ena od vidnih razlik je, da pri njej prevoj nastopi precej prej kot pri logistični. Število rakastih celic narašča sorazmerno z velikostjo tumorja, vendar se sorazmernostni faktor spreminja s časom.

Rešimo Gompertzovo diferencialno enačbo dy

dt =rln K y

!

y,

kjer sta r in K pozitivni konstanti ter y pozitivna funkcija. Diferencialna enačba je avtonomna z ločljivima spremenljivkama. Enačbo delimo s K in jo zapišemo v obliki

d dt

y K

=−r

y K

ln

y K

. Sedaj naredimo substitucijo z = _K^y >0 in dobimo

dz

dt =−rzln(z).

Ob predpostavki, da z6= 1, integriramo

Z dz

zln z =−

Z

r dt;

ln|ln z|=−rt+c, kar je ekvivalentno

lnz =±e^ce^−rt.

Če združimo to s primeromz = 1 (ravnotežna rešitev), dobimo lnz =Ce^−rt,

kar je ekvivalentno

z =e^Ce−rt,

16

kjer je C poljubna realna konstanta. Torej je y=Kz =Ke^Ce

−rt

.

Če upoštevamo začetni pogojy(0) =y₀, dobimo y₀ =y(0) =Ke^Ce

0

=Ke^C, torej jee^C = ^y_K⁰ in dobimo

y=K

y0

K

e^−rt

. Ponekod lahko vidimo tudi ekvivalentno obliko

y=Ke^ln(^yK⁰)^e−rt.

Vendar je stopnja širjenja celic omejena s časom delitve celic, zato Gompertzova diferencialna enačba ni dobra za modeliranje rasti majhnih tumorjev. Gompertzova enačba poleg tega tudi ne upošteva možnosti imunskega nadzora bolezni (tumorja).

Zaradi teh pomislekov so predlagali nov model rasti tumorja, ki se imenuje model Gomp-Ex, ki nekoliko spreminja Gompertzov zakon. Model Gomp-Ex predposta- vlja, da se na začetku, ko še ni konkurence za vire, populacija celic širi eksponentno.

Vendar pa vsebuje kritično mejno vrednost T, tako da za y > T rast sledi Gomper- tzovemu zakonu.

5. primer: rast tkiva

Opis primera je povzet po [5].

Naj bo y(t) območje tkiva (v cm²) v času t (v dnevih). Ko je rast popolna, je končno območje tkiva veliko 8 cm². Večina delitev celic se pojavi na obodu tkiva, kjer je število celic premo sorazmerno√

y. Smiseln model rasti tkiva torej dobimo ob predpostavki, da je hitrost rasti tkiva premo sorazmerna produktu√

yin8−y. Če je konstanta premega sorazmerja ¹₃ in je začetna površina tkiva 1,5 cm², izračunajmo velikost površine tkiva po desetih dneh.

Ker je hitrost rasti tkiva premo sorazmerna √

y in 8−y, je diferencialna enačba, ki opisuje rast tkiva, naslednja

dy

dt =k√

y(8−y). (22)

Ravnotežni rešitvi sta y = φ₁(t) = 0 ter y =φ₂(t) = 8. Če bi narisali graf f(y) = k√

y(8−y), ki ima definicijsko območje le za pozitivne y, bi dobili graf podoben paraboli na sliki 2, ki seka vodoravno os pri y= 0 in y= 8. Torej je y =φ₁(t) = 0 nestabilna ravnotežna rešitev in y = φ₂(t) = 8 stabilna, zato je v tem primeru nosilnost okoljaK = 8. Zaradi tega pričakujemo, da bodo rešitve v pasu0< y < 8 naraščale proti vrednosti y = 8, ko y > 8, pa bodo padale proti y = 8. Lahko bi analizirali še konveksnost/konkavnost ter prevoje funkcije y(t), vendar je postopek podoben kot v podpoglavju 3.1. pri geometrijskem reševanju logistične diferencialne enačbe, kjer smo iz grafa f(y) razbrali kvalitativno pomembne informacije za graf y(t).

Da bi dobili natančno vrednost velikosti tkiva po desetih dneh, rešimo diferen- cialno enačbo še analitično. Ob upoštevanju vrednosti konstante k = ¹₃ in začetne vrednostiy(0) = ³₂ dobimo diferencialno enačbo

dy dt = 1

3

√y(8−y). (23) Ko ločimo spremenljivki, dobimo

√ dy

y(8−y) = 1 3dt.

Uvedemo novo spremenljivko

√y = z;

√dy

y = 2dz in dobimo

Z 2dz 8−z² =

Z 1

3dt. (24)

Uporabimo razcep na parcialne ulomke:

−2 (z−√

8)(z+√

8) = A

z−√

8+ B z+√

8;

−2 (z−√

8)(z+√

8) = (A+B)z+ (A−B)√ 8 (z−√

8)(z+√ 8) . Torej jeA =−^√1

8 inB = ^√1

8 in dobimo

− 1

√8

Z 1 z−√

8− 1 z+√

8

!

dz =

Z 1 3 dt.

Integriramo in dobimo

− 1

√8

lnz−√

8−lnz+√

8= 1 3t+c, kar je ekvivalentno

ln

z+√ 8 z−√

8

=

√8 3 t+c₁. Upoštevamo, da √

y=z in y <8, torej √ y <√

8, ter dobimo ln

√y+√

√ 8 8−√

y

!

=

√8 3 t+c₁. Končna rešitev diferencialne enačbe (23) je

y= 8

Ce

√ 8 3 t−1

2

Ce

√ 8 3 t+ 1

2 ,

kjer je vrednost konstante C, ki jo dobimo iz začetne vrednosti y(0) = ³₂, enaka C =

√1,5 +√

√ 8 8−√

1,5.

18

Če sedaj izračunamo vrednosty pri t= 10, dobimoy∼= 7,99898, kar pomeni, da se površina tkiva zelo hitro približa končni vrednosti 8 cm².

6. primer: hitrost kemijske reakcije drugega reda

Opis primera je povzet po [1].

Kemijska reakcija drugega reda je reakcija, katere hitrost je odvisna od spremembe dveh koncentracij. Če sta to koncentraciji reaktantov P in Q, iz katerih dobimo novo snov X, to zapišemo kot P +Q → X. Predpostavimo, da sta p in q, kjer p6= q, zaporedoma začetni koncentraciji snovi P in Q in naj bo x(t) koncentracija snoviX ob časut. Potem stap−x(t)inq−x(t)koncentraciji snoviP inQob času t in potem je hitrost, pri kateri se kemijska reakcija odvija, podana z diferencialno enačbo

dx

dt =α(p−x)(q−x), (25)

kjer je α pozitivna konstanta.

Če x(0) = 0, najprej določimo mejno vrednost za x(t), ko t → ∞, ne da bi rešili diferencialno enačbo. Nato rešimo problem z začetno vrednostjo in poiščimox(t)za poljubent.

Ravnotežni rešitvi diferencialne enačbe (25) sta y =φ₁(t) =p in y=φ₂(t) =q.

Za graff(y) =α(p−x)(q−x)dobimo parabolo, ki je glede na parabolo na sliki 2 pre- zrcaljena čez vodoravno os, in tu imamo namesto ravnotežnih rešitev y=φ₁(t) = 0 iny=φ₂(t) = K prej omenjeni ravnotežni rešitvi. Fazna premica ima v primerjavi s fazno premico na sliki 3 namesto vrednosti0inK vrednostipinq ter ravno obratno obrnjene puščice. Ker p6=q, moramo ločiti 2 primera:

• p < q: v tem primeru je rešitev y = φ₁(t) = p asimptotično stabilna in y=φ2(t) =q nestabilna rešitev. Torej x(t)→p, ko t→ ∞.

• p > q: v tem primeru je rešitev y = φ₁(t) = p asimptotično nestabilna in y=φ₂(t) =q stabilna rešitev, zato x(t)→q, ko t → ∞.

Sedaj z reševanjem diferencialne enačbe (25) poiščimo x(t) za poljuben t:

dx

(p−x)(q−x) =α dt.

Ulomek na levi strani enačbe razdelimo na parcialne ulomke, integriramo obe strani enačbe ter dobimo

Z 1

(q−p)(p−x)+ 1 (p−q)(q−x)

!

dx = αt+c;

1 p−q

Z

− 1

p−x + 1 q−x

!

dx = αt+c;

1

p−q(ln|p−x| −ln|q−x|) = αt+c.

Absolutne vrednosti lahko spustimo, saj smo že na začetku premislili, da v primeru, kop < q, x(t)→p, ko t→ ∞. V primeru, ko p > q, pa x(t)→q, ko t → ∞. Torej

sta vrednosti (p−x) in(q−x) vedno pozitivni. Dobimo ln p−x

q−x

!

=α(p−q)t+c₁, kar je ekvivalentno

p−x

q−x =Ce^α(p−q)t; izrazimo šex in dobimo

x= pqe^α(p−q)t−1 pe^α(p−q)t−q .

Z analizo dobljene enačbe lahko potrdimo trditve, ki smo jih izpeljali zgolj iz dife- rencialne enačbe (25):

• p < q:

t→∞lim

pqe^α(p−q)t−1

pe^α(p−q)t−q = pq(0−1)

0−q = −pq

−q =p;

• p > q:

t→∞lim

pqe^α(p−q)t−1

pe^α(p−q)t−q = pq1−e^{−α(p−q)t}

p−qe^{−α(p−q)t} = pq−pqe^{−α(p−q)t}

p−qe^{−α(p−q)t} = pq−pq(0) p−q(0) =q.

Sklepi, ki smo jih izpeljali že iz diferencialne enačbe (25) na začetku, so torej pravilni.

20

POGLAVJE 4

Modeliranje s kritičnim pragom

Logistični model, ki upošteva omejitvene dejavnike rasti, je pomanjkljiv v smislu, ker ne upošteva dejstva, da se rast populacije ne pojavi, če je le-ta premajhna. V nadaljevanju bomo zato najprej obravnavali diferencialno enačbo s kritičnim pragom ter primere njene uporabe, nato pa bomo modificirali logistično diferencialno enačbo tako, da bo upoštevan tudi kritični prag, pod katerim se rast populacije ne pojavi.

4.1. Kritični prag Imamo diferencialno enačbo

dy dt =−r

1− y T

y, (26)

kjer sta r in T dani pozitivni konstanti. Opazimo, da se od enačbe (13) razlikuje le v predznaku na desni strani enačbe in v zamenjavi parametra K s T. Pozneje pa bomo videli, da se rešitve enačbe (26) precej razlikujejo od rešitev enačbe (13).

Obravnava in reševanje diferencialne enačbe s kritičnim pragom je povzeto po [1].

4.1.1. Geometrijsko reševanje diferencialne enačbe s kritičnim pra- gom.

Podobno kot pri logistični diferencialni enačbi tudi tu diferencialno enačbo (26) rešimo najprej geometrijsko. Za enačbo (26) dobimo za graf f(y) = −r1−_T^yy parabolo, ki je prikazana na sliki 6.

Slika 6. Graff(y) za ^dy_dt =−r1−_T^yy.

Presečišči z y-osjo sta kritični točki y = 0 in y = T, ki ustrezata ravnotežnima rešitvamay =φ₁(t) = 0 iny =φ₂(t) =T. Če 0< y < T, potem ^dy_dt <0in y pada, ko t narašča. Za y > T pa je ^dy_dt > 0 in y narašča z rastočim t. Torej je φ₁(t) = 0 asimptotično stabilna ravnotežna rešitev in φ₂(t) =T nestabilna.

Poglejmo si še predznake drugega odvoda^d_dt²2^y

, ki smo ga analizirali v obravnavi (14): f⁰(y)je negativen za0< y < ^T₂ in pozitiven za intervala ^T₂ < y < T ter y > T. Graf y(t) je torej konveksen na intervalih 0 < y < ^T₂ ter y > T in konkaven na intervalu ^T₂ < y < T. Torej je v točki (^T₂,−r^T₄) prevoj. Na sliki 7 lahko vidimo fazno premico (os y) enačbe (26). Piki pri y = 0 in y = T predstavljata kritični točki ali ravnotežni rešitvi, puščice pa označujejo, kje so rešitve naraščajoče oziroma padajoče.

Slika 7. Fazna premica rasti s kritičnim pragom.

Sedaj lahko narišemo približne krivulje rešitev enačbe (26): najprej narišemo rav- notežni rešitviy= 0 iny=T. Nato narišemo krivulje na intervalu 0< y < T, ki so padajoče z rastočimtin spremenijo konkavnost, ko sekajo črtoy = ^T₂. Krivulje nad črtoy=T naraščajo bolj in bolj strmo, koyintnaraščata. Ko seypribližuje bodisi 0ali T, krivulje postajajo vse bolj položne. Dobimo torej kvalitativno precej točno skico rešitev enačbe (26) za poljuben r in T. Opazimo, da se z rastočim t y pribli- žuje bodisi0ali pa neomejeno raste, odvisno od tega, ali je bila začetna vrednosty0

pod ali nad vrednostjo T. T je torej mejna vrednost, pod katero ne pride do rasti.

Imenujemo jo tudi mejni prag ali kritični prag. Z geometrijskim sklepanjem smo torej prišli do pomembnih informacij in s tem do grafa y(t) za enačbo (26), ki ga vidimo na sliki 8.

Slika 8. Rast s kritičnim pragom: ^dy_dt =−r1− _T^yy.

22