Populacijska dinamika

Populacijska dinamika

Reševanje sistemov navadnih diferencialnih enačb

Reševanje sistemov se ne razlikuje bistveno od reševanja ene same enačbe. Rešujemo sistem $k$ enačb

\[\begin{array}{lcl} y_1'(t) &=& f_1(y_1(t),y_2(t),\ldots,y_k,t)\\ y_2'(t) &=& f_2(y_1(t),y_2(t),\ldots,y_k,t)\\ \vdots && \vdots\\ y_k'(t) &=& f_k(y_1(t),y_2(t),\ldots,y_k,t)\\ \end{array}\]

z začetnim pogojem $y_i(t_0)=y_{i0}$. Sistem lahko zapišemo kot eno samo vektorsko enačbo. Če označimo z $Y=[y_1,y_2,\ldots,y_k]$ in z $F(Y,t)=[f_1,f_2,\ldots,f_k]$, lahko sistem enačb zapišemo v vektorski obliki

\[Y'(t) = F(Y(t),t).\]

Za reševanje vektorske enačbe lahko uporabimo povsem iste formule kot za eno samo skalarno enačbo.

Numerično reševanje začetnega problema

Iskanju rešitve za diferencialno enačbo

\[y'(t) = f(y(t), t) \]

ki zadošča pogoju $y(t_0)=y_0$ v določenem času $t_0$, pravimo začetni problem za diferencialno enačbo [1]. Numerične metode za reševanje začetnega problema bolj ali manj delujejo po istem principu. Začetni pogoj poda vrednost v točki $t_0$. Najprej izračunamo vrednost $y(t)$ v točki $t_1=t_0+h$, ki je blizu $t_0$. Vrednost $y(t_0+h)$ lahko zapišemo z integralom

\[\begin{equation}\label{eq:integral} y(t_0+h) = y(t_0) + \int_{t_0}^{t_0+h} f(t, y(t))dt. \end{equation}\]

Ker rešitve $y(t)$ ne poznamo, lahko vrednost integrala izračunamo le približno in tako dobimo približek za vrednost $y_1 = y(t_0+h)$. Isti postopek lahko nato ponovimo za začetni problem z začetnimi pogoji $y(t_1)=y_1$ in izračunamo približek za $y(t_1+h_1)$. To ponavljamo, dokler ne dosežemo končne vrednosti $t_n$.

Eulerjeva metoda

Ogledali si bomo preprosto numerično metodo za reševanje začetnega problema imenovano Eulerjeva metoda.

Iščemo rešitev NDE $y'=f(y,t)$ z začetnim pogojem $y(x_0)=y_0$ na intervalu $[t_0,t_k]$. Izberimo si končno mnogo točk

\[x_0< x_1 < \ldots < x_{n+1}=x_k\]

Z Eulerjevo metodo izračunamo zaporedje približkov $y_i$ za $y(x_i)$ z rekurzivno formulo

\[\begin{eqnarray} h_n &=& (x_{n+1}-x_n)\\ y_{n+1} &=& y_n + h_n f(y_n,x_n). \end{eqnarray}\]

Pogosto vzamemo ekvidistančne točke in fiksen korak $h_n=h$.

Opomba: Eulerjeva metoda je 1. reda, kar pomeni, da napaka zelo počasi pada (sorazmerno s $h$), ko manjšamo korak $h$. Zato se v praksi redko uporablja, saj obstajajo metode višjih redov (4,5,6), kjer je hitrost konvergence bistveno večja.

Implicitne metode

Poglejmo si, kaj dobimo, če integral \eqref{eq:integral} izračunamo s trapezno formulo za $t_1 = t_0 +h$

\[\begin{equation}\label{eq:trapez} y(t_1) = y(t_0) + \frac{h}{2}\left(f(y(t_0), t_0) + f(y(t_1), t_1)\right). \end{equation}\]

Vrednost $y_1=y(t_1)$, ki jo želimo izračunati nastopa tudi na desni strani enačbe in je tako ne moremo izračunati eksplicitno. Lahko pa poiščemo rešitev enačbe \eqref{eq:trapez} numerično. Če je $h$ dovolj majhen bo konvergirala že navadna iteracija podana z enačbo \eqref{eq:trapez}.

Implicitne metode

Numerične metode, pri katerih približek ni podan eksplicitno s formulo ali zaporedjem približko, ampak kot rešitev enačbe imenujemo implicitne metode.

[1]

Poleg začetnega problema, lahko za dano diferencialno enačbo rešujemo tudi robni problem, kjer so

pogoji podani v dveh različnih točkah $a$ in $b$. Pri robnih problemih iščemo rešitev znotraj intervala $[a, b]$.

Primer: Lotka-Volterra

Poskusimo rešiti začetni problem za sistem enačb Lotka-Volterra

\[\begin{eqnarray} x' &=& x(\alpha -\beta y)\\ y' = -y(\gamma -\delta x) \end{eqnarray}\]

za parametre $\alpha=\beta=\gamma=1$ in $\delta=2$ in začetnim pogojem $x(0)=y(0)=1$

Red metode

Če rešujemo začetni problem za NDE na fiksnem intervalu $[t_0,t_k]$, napaka približka pada s potenco $h$

\[y_n-y(t_k) = C h^r \]

Eksponent $r$ v potenci $h^r$ imenujemo red metode. Red metode lahko določimo numerično.

Vsak sistem NDE lahko prevedemo na avtonomen sistem

Če za enačbo oziroma sistem enačb $y'(t) = f(y(t), t)$ funkcija $f$ ni eksplicitno odvisno od $t$, se pravi da je $y' = f(y,t) = f(y)$, potem sistem imenujemo avtonomen sistem. V resnici lahko vsak sistem

\[\begin{array}{lcl} y_1'(t) &=& f_1(y_1(t),y_2(t),\ldots,y_k,t)\\ y_2'(t) &=& f_2(y_1(t),y_2(t),\ldots,y_k,t)\\ \vdots && \vdots\\ y_k'(t) &=& f_k(y_1(t),y_2(t),\ldots,y_k,t)\\ \end{array}\]

preprosto razširimo v avtonomen sistem preprosto tako, da dodamo $t$ kot še eno komponento $y_{k+1} = t$ vektorju $y$ in dobimo sistem, kjer $t$ ne nastopa eksplicitno

\[\begin{array}{lcl} y_1'(t) &=& f_1(y_1(t),y_2(t),\ldots,y_k,y_{k+1})\\ y_2'(t) &=& f_2(y_1(t),y_2(t),\ldots,y_k,y_{k+1})\\ \vdots && \vdots\\ y_k'(t) &=& f_k(y_1(t),y_2(t),\ldots,y_k,y_{k+1})\\ y_{k+1}'(t) &=& 1 \end{array}.\]