Numeričke metode rješavanja diferencijalnih jednadžbi

Eulerova metoda

Eulerova metoda je vrlo popularna metoda za rješavanje diferencijalnih jednadžbi, ali ima nedostatak zbog velikih pogrešaka prilikom zaokruživanja rezultata. Ova metoda je zasnovana na korištenju Taylorovog teorema reda  na intervalu , gdje je , te se za primjer diferencijalne jednadžbe prvog reda:

može pisati:

za neki . U ovom slučaju je  pravo rješenje diferencijalne jednadžbe u trenutku . Ako iz jednadžbe (13) izbacimo ostatak, dobivamo aproksimaciju točnog rješenja:

Jednadžba (14) predstavlja Eulerovu metodu.

Jednadžba za Eulerovu metodu je nesimetrična, što znači da nalazi rješenje za vremenski interval , ali pri tome koristi samo informacije o derivaciji s početka intervala (slika 3).

Slika 3.  Eulerova metoda za izračun slijedeće vrijednosti funkcije koristi derivaciju na početku intervala

Prilikom analize kvalitete numeričke metode bitno je odrediti najveću moguću vrijednost pogreške, tj. vrijednosti . Činjenica je da se pogreška metode povećava s korakom .

Pošto parcijalna derivacija po  funkcije  u nekom četverokutu koji sadrži početnu točku  mora biti kontinuirana, mora vrijediti i nejednakost:

za neku konstantu  i za svaki  i  unutar četverokuta domene. Uz pretpostavku da postoji  za koji vrijedi  za sve  za koje se traži rješenje diferencijalne jednadžbe, može se odrediti gornja vrijednost izraza za pogrešku.

Ako se definira , tada nejednakost (16) prelazi u oblik , za . Početna greška je . Greška pri koraku  je:

Greška pri koraku je:

Greška pri koraku  je:

Općenito za vrijedi:

Koristeći činjenicu da je , za  i , te supstitucijom , može se gornja suma (20) dalje raspisati:

gdje je krajnja granica  konstanta koja sadrži  i , te koristi činjenicu da je eksponencijalna funkcija omeđena na zadanom intervalu , tj . Izvod (21) je pokazao da je gornja granična vrijednost pogreške ovisna o konstanti , te veličini koraka . Može se zaključiti da će manji korak povećati preciznost daljnjeg proračuna, ali će sukladno s time, zbog većeg broja potrebnih iteracija da se dođe do krajnjeg , i vrijeme računanja porasti, te se pri upotrebi algoritma mora naći kompromis između preciznosti i vremena trajanja izračuna.

Dosadašnja analiza se zasnivala na činjenici da greška u izračunu nastaje samo uslijed zanemarivanja članova višeg stupnja u aproksimaciji zadane funkcije Taylorovim redom, međutim grešku u proračun unosi i aritmetika brojeva s pomičnim zarezom implementirana u današnjim računalima. Zbog nedostataka u implementaciji aritmetike realnih brojeva u računalima, prilikom izračuna vrijednosti  u rezultantnu vrijednost se unosi sitna greška  nastala prilikom zaokruživanja brojeva s pomičnim zarezom. Rezultat dobiven Eulerovom metodom priliko korištenja aritmetike brojeva s pomičnim zarezom može se prikazati kao:

I u ovom slučaju izraz za pogrešku se prikaže na sličan način kao prije, te se dobiva:

gdje je gornja granica svih  vrijednosti.

Grešku u iteraciji  definiramo kao:, ali se u ovom slučaju početna vrijednost pogreške ne može uzeti kao , jer nije sigurno da je početna vrijednost  u obliku broja s pomičnim zarezom potpuno jednaka vrijednosti . Za iteraciju  vrijedi:

Za iteraciju  vrijedi:

Za iteraciju  vrijedi:

Uz  i neke konstante , , , općeniti oblik izraza za pogrešku je

Bitna stvar za uočiti u gornjoj nejednakosti je postojanje izraza  po kojem se ova nejednakost razlikuje od izraza za gornju graničnu grešku u kojem nije u obzir uzeta pogreška prilikom rada s brojevima s pomičnim zarezom (21). Dok je prijašnji rezultat (21) ukazivao na činjenicu da se gornja granična greška smanjuje smanjenjem intervala , izvod (27) ukazuje na mogućnost povećanja gornje granične pogreške do beskonačnosti u slučaju smanjenja intervala  na nulu. Da bi se postigla najmanja vrijednost pogreške potrebno je odabrati vrijednost  takvu da je funkcija  u svom globalnom minimumu, što je zadovoljeno za .

Runge-Kutta metode

Runge-Kutta metode su zamišljene da imaju aproksimacijske greške višeg stupnja, ali i da za njih nije potrebno računati parcijalne derivacije funkcije , što je u fizikalnim simulacijama neizvedivo, zbog toga što se unaprijed ne zna buduće stanje sustava.

Za izvod Runge-Kuta metoda potrebno je proširiti Taylorov teorem na funkciju s dvije varijable.

Neka su funkcija i njene parcijalne derivacije do uključujući stupnja  kontinuirane na domeni . Neka je . Za svaki  postoje  i  takvi da vrijedi , gdje je  Taylorov polinom stupnja :

a  ostatak razvoja u Taylorov red:

Metode drugog stupnja

Primjena Taylorovog teorema na rezultat  dovodi do jednadžbe  , gdje je  ostatak stupnja . Koristeći jednakosti  i  dobiva se:

Za računanje gornjeg izraza je potrebno znati parcijalne derivacije prvoga reda funkcije , što se želi izbjeći primjenom ove metode. Novom aproksimacijom se mogu zamijeniti derivacije funkcije  s njenom trenutnom vrijednošću, tj. želi se dobiti izraz u obliku

za neke ,,, te ostatke  stupnja . Korištenjem Taylorovog teorema za funkcije dvije varijable dobiva se:

gdje je  stupnja . Desna strana jednadžbe (36) uz uvrštenje gornje jednakosti ima isti oblik kao jednadžba (34), pa se mogu odrediti koeficijenti , i , a izraz za grešku je . Uvrštenjem , i  u jednadžbu (36) dobiva se:

Iz gornje jednadžbe se dobiva numerička metoda nazvana metoda srednje točke (eng. midpoint method):

Metoda srednje točke je nastala kao rezultat korištenja Taylorovog polinoma stupnja 1 da se aproksimira , te numerička metoda ima grešku stupnja . Kod ove metode se derivacija na početku intervala koristi da bi se našla derivacija na sredini intervala. Nakon toga se ta vrijednost derivacije koristi kao derivacija za izračun funkcije na kraju intervala.

Slika 4.  Metoda srednje točke

Moglo bi se pomisliti da bi se istim načinom, ali korištenjem Taylorovog polinoma stupnja 2 mogla postići greška stupnja , međutim ovaj način konstruiranja numeričke metode može dovesti jedino do metoda s greškama stupnja . Razvoj Taylorovog reda stupnja 2 daje:

gdje je  stupnja . Da se zamjene izrazi derivacija funkcije s vrijednostima funkcije koristi se izraz:

uz neke koeficijente , ,  i , te ostatak  stupnja . Koristeći Taylorov teorem prošire se izrazi koji uključuju funkciju :

gdje je  stupnja . Izrazi (40) i (41), uz uvrštenu gornju jednakost (42), skoro odgovaraju, osim elementa . Ako se ti izrazi pokušaju izjednačiti i izračunati odgovarajući koeficijenti, može se odrediti ,  i  uz izraze , i . Svi ostali izrazi iz jednadžbe (40) moraju postati dio ostatka , čiji stupanj se time smanjuje na . Iz tri dobivene jednadžbe postoji više različitih rješenja za koeficijente. Jedno od rješenja može biti: ,  i . Time se dobiva numerička metoda poznata pod nazivom modificirana Eulerova metoda:

Druga metoda je poznata pod nazivom Heunova metoda. Dobije se za rješenja jednadžbi , ,  i  i glasi:

Metode trećeg stupnja

Ako se prilikom zamjene izraza derivacije funkcije s vrijednostima funkcije umjesto jednostruko ugniježđene funkcije  koristi dvostruko ugniježđena funkcija

tada se mogu podudarati i izrazi višeg stupnja. Da bi se dobila Runge-Kutta numerička metoda trećeg stupnja potrebno je proširiti gornji izraz korištenjem Taylorovog teorema. To se postiže tako da se prvo proširi izraz . Rezultat se supstituira u ukupni izraz, te se takav proširi upotrebom Taylorovog teorema.

Jedna od metoda koja se dobije takvim postupkom je:

Dok je druga:

Metode četvrtog stupnja

Izvod Runge-Kutta numeričke metode četvrtog stupnja kreće od aproksimacije funkcije  po Taylorovom teoremu:

Zamjena izraza derivacija funkcije s vrijednostima same funkcije se provodi kao u prethodnim slučajevima, ali će ovdje postojati po jedan izraz bez ugnježđivanja, te po jedan izraz s jednostrukim, dvostrukim, te trostrukim ugnježđivanjem. Krajnji stupanj izraza za pogrešku će biti .

Jedna od rezultirajućih metoda je:

Među češće korištenim metodama je:

i Gillova metoda:

Na slici 5 je prikazan način izračuna vrijednosti funkcije korištenjem Runge-Kutta metode četvrtog stupnja. Pri izračunu svakog koraka gleda se derivacija četiri puta: jedan put u početnoj točci, dva puta u ispitnoj srednjoj točci, te jedan put u ispitnoj krajnjoj točci. Za svaku od tih vrijednosti se računa krajnja vrijednost funkcije označena praznim krugovima, te se korištenjem tih rezultata dobije konačno rješenje označeno punim krugom.

Slika 5. Runge-Kutta metoda četvrtog stupnja