Metode Runge-Kutta explicite#

Metodele Runge-Kutta, o întreagă familie de metode, reprezintă generalizări ale metodei lui Euler, fiind utilizate pentru a soluționa ecuațiile diferențiale ordinare cu o acuratețe mult mai bună.

Formularea anatlitică#

Să considerăm o ecuație diferențială ordinară de forma \(\frac{dy}{dt}=f(t,y)\).

Să definim totodată progresia aritmetică \((t_n)_{n\in\mathbb{N}}\) aleasă astfel încât \(t_k=t_0+kh\), \(\forall k\in\overline{1,n}\), unde \(t_0,h\in\mathbb{R}\) sunt constante cunoscute, iar valoarea \(y(t_0)\) se presupune cunoscută.

Vom nota cu \(y_k\) aproximarea pe care o găsim pentru funcția \(y\) în punctul \(t_k\), adică \(y_k\approx y(t_k)\), \(\forall k\in\overline{0,n}\). Metodele Runge-Kutta explicite soluționează această aproximare prin formula:

unde \(\mu\in\mathbb{N}^*\), valorile \(b_1,\dots,b_\mu\in\mathbb{R}\) sunt constante și reprezintă ponderi, iar coeficienții de forma \(\Psi_i\), \(\forall i\in\overline{1,\mu}\), vor fi recalculați pentru fiecare \(k\) astfel:

unde \(c_2,\dots,c_\mu,a_{ij}\in\mathbb{R}\), \(\forall i\in\overline{2,\mu}\), \(j\in\overline{1,i-1}\), sunt mai multe constante.

Explicație detaliată#

Consider că formularea analitică este mult prea complicată pentru a putea fi parcursă așa. Astfel, voi încerca să aduc o explicație mai digerabilă.

Vrem să folosim aceeași idee precum în cazul metodei lui Euler, adică să ne deplasăm dintr-un punct oarecare în direcția potrivită pe distanța \(\Delta t=h\), direcție aleasă în funcție de derivata din punctul în care ne aflăm. Am observat cum această abordare nu aduce rezultate satisfăcătoare numeric.

Ideea principală din spatele metodelor Runge-Kutta (ecuația \(y_{k+1}=y_k+h\sum_{i=1}^\mu b_i\Psi_i\)) este de a folosi mai multe direcții, \(\mu\in\mathbb{N}^*\) la număr, și anume \(\Psi_1,\dots,\Psi_\mu\), fiecare decizând într-o oarecare măsură direcția finală. Impactul fiecărei direcții este dat de ponderea sa, adică de valorile \(b_1,\dots,b_\mu\).

Așadar, ideea este, în mare, aceeași - se calculează o nouă direcție de deplasare pentru fiecare pas. Geniul algoritmului apare însă pentru a se obține direcțiile de forma \(\Psi_i\), \(\forall i\in\overline{1,\mu}\).

Valorile \(\Psi_i\) sunt evaluări ale funcției \(f\) și se folosesc de influența adusă deja de \(\Psi_1,\dots,\Psi_{i-1}\), însă nu în totalitate, ci tot sub forma unor ponderi, anume ponderea \(a_{i1}\) îi revine lui \(\Psi_1\), ponderea \(a_{i2}\) lui \(\Psi_2\) etc. Totodată, evaluările lui \(f\) pentru \(\Psi_i\) nu trebuiesc făcute neapărat din același punct de pornire \(t\), ele pot varia în funcție de unii coeficienți suplimentari, \(c_2,\dots,c_\mu\).

Tot acest cumul de parametri permite metodelor de tip Runge-Kutta să aproximeze mult mai bine o anumită curbă ce respectă o ecuație diferențială dată.

Înainte de a putea vizualiza metoda, vom mai trece prin câteva noțiuni teoretice.

Consistența unei metode#

Considerăm o metodă de tip Runge-Kutta consistentă dacă și numai dacă, folosind notațiile anterioare, are loc relația:

Ne amintim că termenii \(b_1,\dots,b_\mu\) sunt coeficienții utilizați pentru \(y_{k+1}\), \(\forall k\in\overline{0,n-1}\). Atunci, dacă suma lor egalează \(1\), înseamnă că se formează o medie ponderată a valorilor \(\Psi_1,\dots,\Psi_\mu\).

Alegerea valorilor \(b_1,\dots,b_\mu\) se face astfel încât să se satisfacă un ordin și un rang, noțiuni pe care nu le vom defini concret întrucât sunt dificil de explicat matematic.

Tabelul Butcher#

Pentru a simplifica notațiile utilizate în cadrul metodelor Runge-Kutta, pentru o anumită metodă, se construiește următorul tabel:

Acest tabel este cunoscut drept un tabel Butcher (termenul în limba română nu este definit, în engleză fiind acceptată forma "Butcher tableau").

Această notație facilitează reținerea diverselor metode de tip Runge-Kutta, motiv pentru care o vom adopta și noi în exemplele ce urmează.

Metoda RK4#

Cea mai cunoscută și utilizată metodă de tip Runge-Kutta este cea de ordin \(4\), notată prescurtat RK4 și definită explicit, \(\forall k\in\overline{0,n-1}\), sub forma:

Alternativ, putem adopta scrierea tabelară:

Aceasă metodă este interesant de vizualizat, așadar o vom transforma într-un studiu de caz.

Intuiția geometrică#

Pentru acest exemplu, vom considera ecuația: \(y=t^2+0.5\Leftrightarrow\dfrac{dy}{dt}=\dfrac{2y-1}{t}\).

Asemănător metodei lui Euler, pornim dintr-un punct \(P_0(t_0,y_0)\) și încercăm să ajungem într-un punct \(P_1(t_0+h,y_1)\) aflat pe graficul unei funcții \(y\), descrisă printr-o ecuație diferențială ordinară.

Pentru a ne muta din punctul \(P_0\) în \(P_1\), vom folosi o sumă ponderată de direcții, decisă de valorile \(\Psi_1,\dots,\Psi_4\). Pentru această însă, avem nevoie să soluționăm întâi valorile de direcție. Vom începe cu \(\Psi_1\), echivalentul direct al metodei lui Euler:

Ecuația \(y=t^2+0.5\Leftrightarrow\dfrac{dy}{dt}=\dfrac{2y-1}{t}\);

Calcularea lui \(\Psi_1\)

Menționăm că am desenat \(\Psi_1\) asemănător unui vector pentru a ilustra influența sa asupra direcției finale ce va fi folosită pentru a calcula \(P_1\). Ponderea a fost și ea respectată în ilustrație.

Cunoaștem faptul că \(\Psi_2\) este calculat din abscisa \(t_0+\frac{h}{2}\), însă vine ca un fel de prelungire a lui \(\Psi_1\) - amintim formula sa analitică:

Observăm deci cum, dacă notăm \(t_0+\frac{h}{2}=x\), parametrii depind de o dreaptă cu panta \(\Psi_1\) - mai exact, primul parametru va fi \(x\), iar al doilea parametru va fi \(x\Psi_1+y_0-t_0\Psi_1\), adică se formează o simplă dreaptă de ecuație \(d:y=x\Psi_1+y_0-t_0\Psi_1\).

Din punct de vedere grafic, descoperirea anterioară implică trasarea unui segment punctat ce pornește din \(P_0\) și se deplasează cu \(\frac{h}{2}\) unități pe \(Ot\), păstrând direcția vectorului \(\overrightarrow{\Psi_1}\). Punctul în care se încheie acest segment va fi locul din care pornește vectorul \(\overrightarrow{\Psi_2}\); abscisa acestui punct va fi folosită pentru a calcula panta lui \(\overrightarrow{\Psi_2}\) asemănător metodei lui Euler, însă panta lui \(\overrightarrow{\Psi_2}\) va fi la rândul său afectată de panta lui \(\overrightarrow{\Psi_1}\):

Calcularea lui \(\Psi_2\)

În aceeași manieră, poate fi explicată geometric și aflarea lui \(\Psi_3\) - se va trasa din \(P_0\) un segment având aceeași pantă precum \(\overrightarrow{\Psi_2}\) și care se desfășoară pe distanța \(\Delta t=\frac{h}{2}\), căci formula analitică a lui \(\Psi_3\) este dată de ecuația:

unde \(x=t_0+\frac{h}{2}\). La capătul segmentului, se va construi vectorul \(\overrightarrow{\Psi_3}\) de pantă determinată prin evaluarea lui \(f\), respectiv prin influența lui \(\overrightarrow{\Psi_2}\):

Calcularea lui \(\Psi_3\)

Nu în ultimul rând, se va calcula \(\Psi_4\), trasându-se un segment ce pornește din \(P_0\), are aceeași pantă precum \(\overrightarrow{\Psi_3}\) și se desfășoară pe distanța \(\Delta t = h\), datorită formulei analitice:

unde \(x=t_0+h\). În completarea segmentului, răsare vectorul \(\overrightarrow{\Psi_4}\), ce are panta influențată de evaluarea lui \(f\) în punctul respectiv și de vectorul \(\overrightarrow{\Psi_3}\):

Calcularea lui \(\Psi_4\)

Motivul pentru care ne-am însușit o scriere vectorială a noțiunilor este pentru a facilita soluționarea sistemului - se adună vectorii, iar vârful găsit va fi chiar punctul \(P_1\):

Însumarea vectorilor \(\overrightarrow{\Psi_1},\dots,\overrightarrow{\Psi_4}\) pentru a se obține poziția lui \(P_1\)

Acest rezultat poate fi verificat și din punct de vedere matematic; nu vom elabora prea mult, însă calculele converg către valorile:

Vizualizarea algoritmului#

Să ne amintim exemplul ecuației diferențiale:

Considerând \(h=0.5\), vom calcula curba folosind algoritmul RK4 și îl vom compara cu metoda lui Euler când \(h=0.5\), respectiv când \(h=0.25\) - vom trece direct la grafic:

Curba calculată prin RK4 (\(h=0.5\)) și comparată cu metoda lui Euler, atât pentru \(h=0.5\), cât și pentru \(h=0.25\)

Este inutilă ilustrarea erorii, întrucât, în cel mai rău caz, valoarea sa absolută nu depășește \(0.007\), ceea ce face algoritmul foarte puternic. Prin comparație, metoda lui Euler are o eroare supraunitară, chiar și pentru \(h=0.25\).

Cod ilustrativ#

Algoritmul metodei Runge-Kutta explicită de ordin 4 poate fi scris în Matlab sub forma funcției RK4. Așadar, dacă:

• a - capătul de la care se pornește, echivalentul lui \(t_0\);

• b - capătul la care se dorește să se ajungă, echivalentul lui \(t_n\);

• n - numărul de pași pe care îi va face algoritmul (tehnic, impune pasul \(h=\frac{b-a}{n}\));

• y0 - evaluarea funcției pe care o căutăm în punctul \(t_0\), echivalentul lui \(y_0=y(t_0)\);

• f - partea dreaptă a ecuației diferențiale \(\frac{dy}{dt}=f(t,y)\).

Atunci funcția va returna y, soluția ecuației diferențiale.

function y = RK4(a, b, n, y0, f)
    h = (b - a) / n
    y = zeros(n+1, 1);
    y(1) = y0;
    for k = [1 : n]
        t = a + (k - 1) * h;
        psi1 = f(t,       y(k))
        psi2 = f(t + h/2, y(k) + h*psi1/2)
        psi3 = f(t + h/2, y(k) + h*psi2/2)
        psi4 = f(t + h,   y(k) + h*psi3)
        y(k+1) = y(k) + (h/6) * (psi1 + 2*psi2 + 2*psi3 + psi4);
    end
end

Alte metode RK#

Vom prezenta foarte pe scurt alte câteva metode Runge-Kutta pe care le considerăm interesante.

Metoda lui Euler#

După cum spuneam la începutul capitolului, metodele Runge-Kutta sunt generalizări ale metodei Euler de rezolvare a unei ecuații diferențiale ordinare. De aceea, metoda lui Euler poate fi scrisă sub următoarea formă tabelară:

Metode RK2#

La modul general, o metodă RK2 (ordin \(2\), rang \(2\)) se va putea scrie folosind următorul tabel, unde \(\alpha\in\mathbb{R}^*\) este o constantă oarecare:

Metoda lui Heun#

Este metoda lui Euler impresionantă sub forma unei metode RK? Nu foarte, însă demonstrează generalitatea acestor metode. Interesantă este însă îmbunătățirea sa, cunoscută sub denumirea de metoda lui Heun:

Metoda punctului de mijloc#

O ultimă metodă RK2 pe care o prezentăm este metoda punctului de mijloc:

Metoda RK3#

Încheiem prezentând tabelul pentru metoda RK3 propusă de însuși Kutta:

Licență#

The book "Metode Numerice", written by Valentin-Ioan Vintilă, is licensed under CC BY-NC-SA 4.0