Kazalo

Uvod¶

Zapis (ene) diferencialne enačbe¶

Predpostavimo, da je mogoče diferencialno enačbo prvega reda zapisati v eksplicitni obliki: $$y'=f(t, y),$$ kjer je $f(t, y)$ podana funkcija in velja $y'=dy/dx$.

Dodatno je podan začetni pogoj: $$y(t_0)=y_0.$$

Cilj reševanja diferencialne enačbe je izračunati funkcijo $y(t)$, ki reši zgoraj definiran začetni problem. Ob določenih pogojih funkcije $f(t, y)$ ima začetni problem enolično rešitev na intevalu, ki vsebuje $t_0$.

Pri numeričnem reševanju vedno računamo tabelo funkcije $y(t_i)$, ki reši dan začetni problem. Pri tem so vozlišča $t_i$ običajno ekvidistantna: $$t_0, t_0+h, t_0+2h,\dots$$ in $h$ imenujemo (časovni) korak (integracije).

Tukaj si bomo pogledali nekatere numerične metode za reševanje diferencialnih enačb pri začetnem pogoju.

Eulerjeva metoda¶

Eksplicitna Eulerjeva metoda temelji na razvoju funkcije $y$ v Taylorjevo vrsto: $$y(t+h)=y(t)+y'(t, y(t))\,h + \mathcal{O}(h^2).$$

Naredimo napako metode $\mathcal{O}(h^2)$, ker zanemarimo odvode drugega in višjih redov; sedaj lahko ob znani vrednosti $y(t)$ in odvodu $y'(t)=f(t,y)$ ocenimo vrednosti pri naslednjem časovnem koraku $t+h$. Ko imamo enkrat znane vrednosti pri $t+h$, ponovimo postopek!

Koraki Eulerjeve metode:

1. Postavimo $i=0$, $t_0$, $y_0=y(t_0)$.
2. Izračun vrednosti funkcije pri $t_{i+1}=t_i+h$: $$y_{i+1}= y_i + f(t_i, y_i)\,h.$$
3. $i=i+1$ in nadaljevanje v koraku 2.

Diferencialno enačbo rešujemo na intervalu $[t_0,t_n]$ in velja $h=(t_n-t_0)/n$. $n$ je število integracijskih korakov (kolikokrat izvedemo korak 2 v zgornjem algoritmu).

Numerična rešitev začetnega problema: $$y_0, y_1, y_2,\dots, y_{n}$$

pri vrednostih neodvisne spremenljivke: $$t_0, t_1, t_2\, \dots \, t_n.$$

Napaka Eulerjeve metode¶

Napaka Eulerjeve metode na vsakem koraku je reda $\mathcal{O}(h^2)$.

Ker na intervalu od $t_0$ do $t_n$ tako napako naredimo $n$-krat, je kumulativna napaka $n\,\mathcal{O}(h^2)=\frac{t_n-t_0}{h}\,\mathcal{O}(h^2)=\mathcal{O}(h)$.

Lokalno je napaka drugega reda, globalno pa je napaka prvega reda in ker je Eulerjeva metoda tako nenatančna jo redko uporabljamo v praksi!

Ocena napake¶

Točna rešitev $y(t_n)$ pri velikosti koraka $h$ je: $$y(t_n)=y_{n,h}+E_h,$$ kjer je $y_{n,h}$ numerični približek in $E_h$ napaka metode. Ker je globalna napaka prvega reda, lahko napako zapišemo kot: $$E_h=k\,h.$$

Podobno lahko za velikost koraka $2h$ zapišemo: $$y(t_n)=y_{n,2h}+E_{2h},$$ kjer je $y_{n,2h}$ numerični približek in $E_{2h}$ napaka metode: $$E_{2h}=k\,2\,h.$$

Ob predpostavki, da je konstanta $k$ pri koraku $h$ in koraku $2h$ enaka, lahko določimo oceno napake pri boljšem približku $E_h$. Očitno velja: $$y_{n,h}+k\,h=y_{n,2h}+2\,k\,h$$ nato določimo oceno napake: $$E_h=k\,h=y_{n,h}-y_{n,2h}.$$

Komentar na implicitno Eulerjevo metodo¶

Pri eksplicitni Eulerjevi metodi računamo rešitev pri $t_{i+1}$ iz izračunane vrednosti pri $t_i$.

V kolikor bi nastopala neznana vrednost rešitve pri $t_{i+1}$, to je $y_{i+1}$, tudi na desni strani, bi govorili o implicitni Eulerjevi metodi (ali povratni Eulerjevi metodi):

Ker se iskana vrednost $y_{i+1}$ nahaja na obeh straneh enačbe, moramo za določitev $y_{i+1}$ rešiti (nelinearno) enačbo. Prednost implicitne Eulerjeve metode je, da je bolj stabilna (npr. v primeru togih sistemov, ki jih bomo spoznali pozneje) kakor eksplicitna oblika, vendar pa je numerično bolj zahtevna (zaradi računanja rešitve enačbe).

Numerična implementacija¶

Najprej uvozimo potrebne knjižnice:

Nato definirajmo Eulerjevo metodo:

Pripravimo funkcijo za oceno napake (v numeričnem smislu bi bilo bolje oceno napake vključiti v funkcijo euler, vendar jo zaradi jasnosti predstavimo ločeno):

Numerični zgled¶

Kot primer rešimo diferencialno enačbo, ki opisuje padanje telesa, ki je izpostavljeno sili teže in zračnemu uporu:

Glede na II. Newtonov zakon, lahko zapišemo diferencialno enačbo:

kjer je $m$ masa, $g$ gravitacijski pospešek, $c$ koeficient zračnega upora in $v$ hitrost. Diferencialno enačbo bi hoteli rešiti glede na začetni pogoj:

Funkcija desne strani / prvega odvoda $f(t,y)$ je:

in začetni pogoj: $$y_0=0.$$

Definirajmo funkcijo desnih strani:

Definirajmo začetni pogoj in časovni vektor, kjer nas zanima rezultat:

Kličemo funkcijo euler za izračun vrednosti $y$ (hitrost $v$):

Prikažemo rezultat:

Preverimo sedaj vpliv časovnega koraka:

Opazimo, da se numerična napaka pri spremembi koraka iz 1 na 0,1 bistveno zmanjša!

Ocenimo še napako pri 100 in 1000 odsekih:

Ko smo korak zmanjšali na desetino, se je proporcionalno zmanjšala tudi napaka (prvi red napake).

Poglejmo še primer, ko je zračni upor c argument funkcije euler in je prek **kwargs posredovan v funkcijo f_zračni_upor():

Metoda Runge-Kutta drugega reda¶

Eulerjeva metoda je prvega reda (prvega reda je namreč globalna napaka $\mathcal{O}(h)$). Če bi želeli izpeljati metodo drugega reda napake, bi si morali pomagati z razvojem $y(t+h)$ v Taylorjevo vrsto, kjer bomo zanemarili tretji in višje odvode: $$y(t+h)=y(t)+y'(t)\,h + \frac{1}{2}y''(t)\,h^2+\mathcal{O}(h^3).$$ Lokalna napaka metode bo tako tretjega reda, globalna pa drugega reda.

Uporabimo zamenjavi $y'(t)=f(t,y)$ in $y''(t)=f'(t,y)$: $$y(t+h)=y(t)+f(t,y)\,h + \frac{1}{2}f'(t,y)\,h^2+\mathcal{O}(h^3).$$

Ker je desna stran $f(t,y)$ odvisna od neodvisne $t$ in odvisne spremenljivke $y$, moramo uporabiti implicitno odvajanje: $$f'(t,y)=\frac{\partial f(t,y)}{\partial t}+\frac{\partial f(t,y)}{\partial y}\,\underbrace{\frac{\textrm{d} y}{\textrm{d} t}}_{y'=f(t,y)} =\frac{\partial f(t,y)}{\partial t}+\frac{\partial f(t,y)}{\partial y}\,{f(t, y)}.$$

Vstavimo v izraz za Taylorjevo vrsto: $$y(t+h)_{\textrm{Taylor}}=y(t+h)=y(t)+f(t,y)\,h + \frac{1}{2}\,{\LARGE(} \frac{\partial f(t,y)}{\partial t}+\frac{\partial f(t,y)}{\partial y}\,{f(t, y)} {\LARGE)}\,h^2.$$

Kot je razvidno iz zgornjega izraza, potrebujemo dodatne odvode. To predstavlja določeno težavo, ki se ji lahko izognemo na različne načine; v nadaljevanju si bomo pogledali pristop Runge-Kutta. Ker bomo zgornji izraz pozneje še potrebovali, smo ga tukaj poimenovali $y(t+h)_{\textrm{Taylor}}$.

Ideja pristopa Runge-Kutta¶

Zgornjo dilemo metoda Runge-Kutta (razvita leta 1901) rešuje z idejo, ki smo jo sicer že srečali pri Gaussovi integraciji: točnejšo rešitev poskuša najti z uteženo dodatno vrednostjo funkcije $f$:

kjer so $c_0$, $c_1$, $p$ in $q$ neznane konstante (načeloma od 0 do vključno 1). Če bi v zgornjem izrazu uporabili $c_1=0$, bi izpeljali metodo prvega reda; z dodatno funkcijsko vrednostjo ($A$) pa se bo izkazalo, da bomo izpeljali metodo drugega reda.

Iskanje neznanih konstant $c_0$, $c_1$, $p$, $q$ nadaljujemo z zapisom $A$ v obliki Taylorjeve vrste prvega reda: $$f{\large(}t+p\,h,y+q\,h\,f(t,y){\large)}= \underbrace{ f{\large(}t,y{\large)}+ \frac{\partial f(t,y)}{\partial t}\,\left(p\,h\right)+ \frac{\partial f(t,y)}{\partial y}\,\left(q\,h\,f(t, y)\right) }_{B}.$$

Vstavimo sedaj izpeljani $B$ nazaj v izraz za $y(t+h)_{\textrm{Runge-Kutta}}$: $$y(t+h)_{\textrm{Runge-Kutta}}=y(t)+c_0\,f\left(t,y\right)\,h +c_1 {\LARGE(} f{\large(}t,y{\large)}+ \frac{\partial f(t,y)}{\partial t}\,\left(p\,h\right)+ \frac{\partial f(t,y)}{\partial y}\,\left(q\,h\,f(t, y)\right) {\LARGE)} \,h.$$

Nadaljujemo z izpeljevanjem in enačbo preoblikujemo, da bo podobna zgoraj izpeljani s Taylorjevo vrsto $y(t+h)_{\textrm{Taylor}}$:

Primerjajmo sedaj z zgoraj izpeljanim izrazom: $$y(t+h)_{\textrm{Taylor}}=y(t)+f(t,y)\,h + \frac{1}{2}\,{\LARGE(} \frac{\partial f(t,y)}{\partial t}+\frac{\partial f(t,y)}{\partial y}\,{f(t, y)} {\LARGE)}\,h^2.$$

Ugotovimo, da za enakost mora veljati:

Imamo torej tri enačbe in štiri neznanke. Eno od konstant si tako lahko poljubno izberemo, ostale tri pa izračunamo. Če na primer izberemo $c_0=0$, bi to imenovali spremenjena Eulerjeva metoda in bi ostali parametri bili: $c_1=1$, $p=q=1/2$. Izbira parametrov nima bistvenega vpliva na rešitev. Sicer pa velja omeniti, da tudi metodo Runge-Kutta drugega reda redko uporabljamo, saj obstajajo boljše metode.

Parametre $c_0$, $c_1$, $p$ in $q$ vstavimo v prvo enačbo tega poglavja. Ko je definiran začetni čas $t_0$ in začetni pogoj $y_0$, uporabimo metodo Runge-Kutta drugega reda:

Metoda Runge-Kutta četrtega reda¶

Podobno kot smo izpeljali metodo Runge-Kutta drugea reda, se izpelje metodo Runge Kutta četrtega reda. Tudi pri metodi četrtega reda obstaja več različic in kot metoda Runge-Kutta četrtega reda razumemo naslednjo metodo:

kjer so: $$\begin{align} k_0&=h\,f(t_i,y_i)\\ k_1&=h\,f\left(t_i+\frac{h}{2},y_i+\frac{k_0}{2}\right)\\ k_2&=h\,f\left(t_i+\frac{h}{2},y_i+\frac{k_1}{2}\right)\\ k_3&=h\,f\left(t_i+h,y_i+k_2\right). \end{align} $$

Koraki metode Runge-Kutta četrtega reda so:

1. Določitev $i=0$ in $t_0$, $y_0=y(t_0)$,
2. Izračun koeficintov: $k_0$, $k_1$, $k_2$, $k_3$,
3. Izračun vrednosti rešitve diferencialne enačbe pri $t_{i+1}=t_i+h$: $$y_{i+1}=y_i+\frac{1}{6}(k_0+2\,k_1+2\,k_2+k_3),$$
4. $i=i+1$ in nadaljevanje v koraku 2.

Napaka metode Runge-Kutta četrtega reda¶

Metodo Runge-Kutta četrtega reda imenujemo tako zato, ker ima lokalno napako petega reda $\mathcal{O}(h^5)$, vendar pa to napako naredimo $n$-krat, zato je globalna napaka četrtega reda $\mathcal{O}(h^4)$.

Ocena napake¶

Točen rezultat $y(t_n)$ pri velikosti koraka $h$ je:

kjer je $y_{n,h}$ numerični približek rešitve in $E_h$ napaka metode. Ker je globalna napaka četrtega reda, lahko napako zapišemo tako:

Podobno lahko za velikost koraka $2h$ zapišemo: $$y(t_n)=y_{n,2h}+E_{2h},$$

kjer je $y_{n,2h}$ numerični približek rešitve in $E_{2h}$ napaka metode: $$E_{2h}=k\,(2\,h)^4=16\,k\,h^4.$$

Ob predpostavki, da je konstanta $k$ pri koraku $h$ in koraku $2h$ enaka, lahko izračunamo oceno napake pri boljšem približku $E_h$.

Najprej je res:

sledi:

in nato določimo oceno napake natančnejše rešitve:

Numerična implementacija¶

Funkcija za oceno napake:

Numerični zgled¶

Poglejmo sedaj primer izračuna hitrosti padajoče mase:

Podajmo začetni pogoj in časovni vektor, kjer nas zanima rezultat:

Za primerjavo izračunajmo rešitev s funkcijo euler ter runge_kutta_4:

Prikažemo rezultat:

Poglejmo še numerično napako:

Pri zmanjšanju koraka na desetino, se je napaka zmanjšala za približno $10^4$-krat (kar ustreza pričakovanjem za metodo četrtega reda).

Uporaba scipy za reševanje navadnih diferencialnih enačb¶

Paket scipy ima implementiranih veliko numeričnih metod za reševanje začetnih problemov navadnih diferencialnih enačb. Tukaj si bomo ogledali funkcijo scipy.integrate.solve_ivp, ki je primerna za večino začetnih problemov.

scipy.integrate.solve_ivp¶

Sintaksa za uporabo (IVP - angl. Initial Value Problem):

scipy.integrate.solve_ivp(fun, t_span, y0, method='RK45', 
                          t_eval=None, dense_output=False, 
                          events=None, vectorized=False, **options)

Pojasnilo vseh argumentov je v dokumentaciji, tukaj bomo izpostavili nekatere:

• fun je desna stran (func(t, y...)),
• t_span terka (t0, tf), ki definira začetni t0 in končni čas tf,
• y0 seznam začetne(ih) vrednosti,
• method definira numerično metodo (privzeta je eksplicitna 'RK45', za toge sisteme pa so boljše: 'Radau', 'BDF', 'LSODA'),
• t_eval uporabimo, če želimo rešitve ob določenih vrednosti neodvisne spremenljivke,
• dense_output ali se pripravi tudi zvezna rešitev, privzeto False,
• events za sledenje dogodkov, ki ustavijo reševanje začetnega problema.

Rezultat klicanja solve_ivp je objekt z atributi (izbrani):

• t vrednosti neodvisne spremenljivke pri katerih je izračunan rezultat,
• y rezultat,
• sol funkcija zvezne rešitve (samo v primeru dense_output=True)
• success je True, če je bila rešitev uspešno izračunana ali prekinjena zaradi dogodka (events).

Numerični zgled¶

Definirajmo začetni pogoj in časovni vektor, kjer nas zanima rezultat (prikažemo samo prvih deset elementov):

Argumente (npr.: c=0.4) v funkcijo fun vpeljemo prek lambda izraza (dokumentacija):

Zgoraj uporabljen lambda izraz lambda t, y: f_zračni_upor(t, y, c=0.4) je ekvivalenten:

def ime_funkcije(t, y):
    return f_zračni_upor(t, y, c=0.4)

Poglejmo rezultat:

Ker je solve_ivp že pripravljen za sistem navadnih diferencialnih enačb prvega reda, je rezultat podan v dveh dimenzijah. Iskan rezultat je torej: y_ivp.y[0].

Za primerjavo še rezultat z lastno implementacijo Runge-Kutta 4. reda:

Prikažemo rezultat:

Sistem navadnih diferencialnih enačb¶

Zgoraj smo si pogledali reševanje začetnega problema ene navadne diferencialne enačbe; sedaj bomo reševanje posplošili na začetni problem za sistem $m$ navadnih diferencialnih enačb prvega reda.

Takšen sistem zapišemo: $$\mathbf{y}'=\mathbf{f}(t, \mathbf{y}),$$ kjer so podani dodatni (začetni) pogoji: $$\mathbf{y}(t_0)=\mathbf{y}_0.$$ S $t$ smo označili neodvisno spremenljivko (ni nujno, da je to vedno čas) in $\mathbf{f}$ vektor desnih strani.

Računamo rešitev sistema $m$ diferencialnih enačb $\mathbf{y}'=\mathbf{f}(t, \mathbf{y})$ pri $a=t_0, t_1,\dots,t_{n-1}=b$, to je $m$ funkcijskih vrednosti na vsakem koraku $y_k(t_i)$, $k=0,1,\dots,m-1$ in $i=0,1,\dots,n-1$. Lahko se dokaže, da lahko uporabimo vsako metodo, ki smo jo izpeljali za diferencialne enačbe prvega reda tudi za sistem navadnih diferencialnih enačb prvega reda, če zamenjamo skalarne veličine z ustreznimi vektorskimi.

Numerična implementacija¶

Metodi Euler in Runge-Kutta četrtega reda potrebujeta zgolj malenkostne popravke (y0 je numerično polje f vrne seznam vrednosti odvodov):

Numerični zgled¶

Padanje mase nadgradimo v ravninsko gibanje; velikost sile upora zraka naj bo definirana kot: $$|\mathbf{F}(\mathbf{v})|=c\,|\mathbf{v}|^2.$$ Sila upora zraka deluje v nasprotno stran kot kaže vektor hitrosti.

Ob pomoči slike, definiramo silo v $x$ smeri: $$F_x=-c\,\left(v_x^2+v_y^2\right)\,\cos(\alpha)=-c\,\left(v_x^2+v_y^2\right)\,\frac{v_x}{\sqrt{v_x^2+v_y^2}}=-c\,v_x\,\sqrt{v_x^2+v_y^2}.$$

Podobno je sila v $y$ smeri: $$F_y=-c\,v_y\,\sqrt{v_x^2+v_y^2}.$$

Glede na drugi Newtonov zakon zapišemo sistem dveh (vezanih) diferencialnih enačb prvega reda:

$m$ je masa, $g$ gravitacijski pospešek, $c$ koeficient zračnega upora ter $v_x$ in $v_y$ hitrost v $x$ oz $y$ smeri. Diferencialno enačbo bi želeli rešiti glede na začetni pogoj:

Definirajmo seznam desnih strani:

Definirajmo začetni pogoj in časovni vektor, kjer nas zanima rezultat:

Za primerjavo izračunajmo rešitev s funkcijo runge_kutta_4 ter solve_ivp:

Poglejmo rezultat:

Prikažemo rezultat:

Preoblikovanje diferencialne enačbe višjega reda v sistem diferencialnih enačb prvega reda¶

Pogledali si bomo, kako navadno diferencialno enačbo poljubnega reda $n$:

pri začetnih pogojih: $$y(t_0)=k_0,\quad y'(t_0)=k_1,\quad\dots\quad y^{(n-1)}(t_0)=k_{n-1}$$

preoblikujemo v sistem diferencialnih enačb prvega reda.

Najprej namesto odvodov vpeljemo nove spremenljivke:

Diferencialna enačba $n$-tega reda zapisana z novimi spremenljivkami je:

Nove funkcije odvajamo po neodvisni spremenljivki:

in $$y_{n-1}'=f(t, y_0, y_1, y_2,\dots,y_{n-1}).$$

Dobili smo sistem navadnih diferencialnih enačb prvega reda: $$ \begin{array}{rcl} y_0'&=&y_1\\ y_1'&=&y_2\\ &\dots\\ y_{n-1}'&=&f(t, y_0, y_1, y_2,\dots,y_{n-1})\\ \end{array} $$

pri začetenih pogojih: $$y_0(t_0)=k_0,\quad y_1(t_0)=k_1,\quad\dots\quad y_{(n-1)}(t_0)=k_{n-1}.$$

Numerični zgled¶

Vrnemo se k padajoči masi:

Vendar tokrat drugi Newtonov zakon zapišimo glede na pomik (diferencialna enačba drugega reda):

$m$ je masa, $g$ gravitacijski pospešek, $c$ koeficient zračnega upora ter $x''$ in $y''$ pospešek v izbranem koordinatnem sistemu. Začetni pogoji:

Imamo sistem dveh diferencialnih enačb drugega reda. Z uvedbo novih spremenljivk $y_i$:

Pripravimo sistem diferencialnih enačb prvega reda: $$ \begin{array}{rcl} y_0'&=&y_1\\ y_1'&=&F_x/m\\ y_2'&=&y_3\\ y_3'&=&g+F_y/m \end{array} $$

Definirajmo Pythonovo funkcijo desnih strani / prvih odvodov:

Definirajmo začetni pogoj in časovni vektor, kjer nas zanima rezultat:

Izračunamo rešitev:

Poglejmo rezultat ($[x, x', y, y']$):

Prikažemo rezultate; najprej hitrost, nato lego ($y$ koordinata je pozitivna navzdol)!

Stabilnost reševanja diferencialnih enačb*¶

Numerično reševane diferencialnih enačb je izpostavljeno napaki metode in zaokrožitveni napaki. Te napake so pri različnih numeričnih metodah različne.

Reševanje diferencialne enačbe je stabilno, če majhna sprememba začetnega pogoja vodi v majhno spremembo izračunane rešitve; sicer govorimo o nestabilnosti reševanja.

Stabilnost je odvisna od diferencialne enačbe, od uporabljene numerične metode in od koraka integracije $h$.

Primer preprostega nihala¶

Poglejmo si najprej primer reševanja diferencialne enačbe preprostega nihala.

Slika (vir: Slavič, Dinamika, mehanska nihanja in mehanika tekočin, 2017) prikazuje dinamski sistem (masa $m$, togost $k$), katerega diferencialna enačba je

Tako diferencialno enačbo preoblikujemo v standardno obliko lastnega nihanja:

kjer je lastna krožna frekvenca:

in pričakujemo odziv oblike:

Če so začetni pogoji:

je rešitev začetnega problema:

Numerični zgled¶

Najprej definirajmo vektor začetnih pogojev in funkcijo desnih strani / prvih odvodov (diferencialno enačbo drugega reda pretvorimo v sistem diferencialnih enačb prvega reda $\mathbf{y}'=\mathbf{f}(t, \mathbf{y})$), kjer velja $y_0=x, y_1=x'$:

Definirajmo podatke in analitično rešitev:

Rešitev s pomočjo metod Euler in Runge-Kutta četrtega reda:

Prikažimo rezultate:

Opazimo, da je Eulerjeva metoda nestabilna in če bi povečali korak, bi postala nestabilna tudi metoda Runge-Kutta četrtega reda.

Zakaj je Eulerjeva metoda tako nestabilna?¶

Spomnimo se Eulerjeve metode:

ki nam pove, da pomik $x(t+h)$ določimo glede na lego $x(t)$ in hitrost $x'(t)$. Začetni pogoji izhajajo iz skrajnega odmika $x(t=0)=1$ in takrat je hitrost $x'(t=0)=0$, kar pomeni, da bo $x(t+h)=1$. Že v prvem koraku torej naredimo razmeroma veliko napako. Vendar zakaj potem začne vrednost alternirajoče naraščati?

Spomnimo se, da je analitična rešitev $x(t)=x_0\,\cos(\omega_0\,t)$ in je torej $x'(t)=-\omega_0\,x_0\,\sin(\omega_0\,t)$.

Vstavimo pripravljena izraza v Eulerjevo metodo in uredimo:

Predpostavimo, da gledamo stanje ob takem času $t=\pi/(2\omega_0)$, ko velja $\cos(\omega_0\,t)=0$ in $\sin(\omega_0\,t)=1$:

V kolikor bo absolutna vrednost izraza $A$ večja kot 1, bo pri času $t+h$ vrednost večja kot v predhodnem koraku in v sledečem verjetno spet. Sledi, da lahko pride do nestabilnosti. Da se je izognemo, mora veljati: $$|A|<1\qquad\rightarrow\qquad h<\frac{1}{\omega_0}.$$

Opomba: v nekaterih knjigah boste videli tudi vrednost $h<2/\omega_0$; enoliče meje za vse diferencialne enačbe ni mogoče definirati; v splošnem pa velja, da je korak definiran relativno glede na najkrajšo periodo $T$ v diferencialni enačbi (npr.: $h<2/\omega_0$ je v bistvu enako $h<2/(2\pi/T)$ oziroma $h<T/\pi$). Perioda $T$ je definirana glede na najvišjo lastno frekvenco sistema $T=1/f_{\textrm{max}}$, ki jo izračunamo iz lastne vrednosti sistema.

Nekaj vprašanj za razmislek!¶

1. Na sliki (vir: Numerical Methods in Engineering With Python 3, 3rd Ed, Jaan Kiusalaas) je prikazan izstrelek mase $m$, ki ga izstrelimo s hitrosjo $v_0$ pod kotom $\alpha$. Če je sila upora zraka: $F=c\,v^{3/2}$, potem sta gibalni enačbi: $$\ddot x(t)=-F\,\cos(\alpha)/m\qquad\ddot y(t)=-F\,\sin(\alpha)/m-g.$$ Opomba: $v=\sqrt{\dot x^2+\dot y^2}.$ Ustrezne parametre si izmislite. Sistem dveh diferencialnih enačb drugega reda zapišite v sistem diferencialnih enačb prvega reda.
2. Določite vektor začetnih pogojev, ki smo ga zgoraj označili z $\mathbf{y}$.
3. Določite funkcijo desnih strani, c naj bo parameter.
4. Definirajte začetne pogoje in rešite nalogo s poljubnimi podatki.
5. Prikažite $(x, y)$ lego masne točke, spreminjajte koeficient upora $c$.
6. Prikažite hitrost v odvisnosti od časa. Določite minimum hitrosti in čas, pri katerem nastane.

Primer Van der Polovega nihala¶

Namen tega primera je pokazati, kako lahko izbira integratorja vpliva na hitrost reševanja problema! Van der Polovo nihalo je opisano tukaj.

Definirajmo seznam odvodov:

Rešitev po metodi RK45 (gre za eksplicitno shemo, ki ni primerna za toge sisteme diferencialnih enačb; reševanje je zelo počasno, zato rešitev računamo samo do t1/100):

Implicitna shema BDF (angl. Backward Differentiation Formulas) se tukaj izkaže kot bistveno bolj primerna. Zaradi stabilnosti, so koraki lahko bistveno večji in zato je reševanje bistveno hitrejše:

Simbolno reševanje diferencialne enačbe drugega reda¶

Pogledali si bomo primer, prikazan na sliki, kjer je masa $m$ na klancu naklona $\alpha$. Koeficient trenja je $\mu$, težnostni pospešek pa $g$. Začetna hitrost je $\dot y(0\,\textrm{s})=v_0$, pomik $y(0\,\textrm{s})=0\,\textrm{m}$. Gibalna enačba (samo za smer $y$) je definirana glede na II. Newtonov zakon (glejte diagram sil na prosto telo).

Izpeljava gibalne enačbe¶

Rešitev enačbe je:

Da določimo $C_1$ in $C_2$, vstavimo $t=0\,\textrm{s}$:

Nato odvajamo po času in ponovno vstavimo $t=0\,\textrm{s}$:

Glede na začetne pogoje smo torej določili konstante:

Sledi rešitev:

Pripravimo si funkciji za numerični klic:

Pripravimo prikaz:

Simbolno preoblikovanje diferencialne enačbe v sistem diferencialnih enačb prvega reda¶

Spomnimo se izvorne diferencialne enačbe:

Definirajmo nove spremenljivke in pripravimo funkcijo $f$:

Povežimo sedaj nove spremenljivke.

$d y_0/dt$ naj bo enako $y_1$:

Odvod $d y_1/dt$ (v bistvu je to $y''$) naj bo enak funkciji $f$:

Zgornje izraze zapišemo v vektorski obliki: $$\mathbf{y}'=\mathbf{f}(t, \mathbf{y}).$$

Spomnimo se sedaj f_vec:

Če rešujemo numerično, potem funkcijo $\mathbf{f}(t, \mathbf{y})$ zapišemo:

Preverimo funkcijo pri začetnem času $t=0\,$s in pri začetnih pogojih $[y_0, y_1]=[0, v_0]$:

Uporabimo sedaj Eulerjevo metodo:

Prikažemo in primerjamo z analitično rešitvijo: