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 , 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 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}.$