Metode Runge–Kutta

Compararea metodelor Runge-Kutta pentru ecuația diferențială
$\quad y\,'=\sin ^{2}(t)\cdot y$
(cu roșu este soluția exactă)

În analiza numerică metodele Runge–Kutta sunt o familie de metode iterative implicite și explicite^⁠(d), care cuprind metoda Euler^⁠(d), utilizată în discretizarea timpului pentru soluțiile aproximative ale sistemelor de ecuații neliniare^⁠(d). Aceste metode au fost dezvoltate în jurul anului 1900 de matematicienii germani Carl Runge și Martin Wilhelm Kutta.

Metoda Runge–Kutta

Pantele utilizate în metoda Runge–Kutta clasică (RK4)

Cea mai cunoscută metodă dintre metodele Runge–Kutta este cea denumită „RK4”, „metoda clasică Runge–Kutta” sau, simplu, „metoda Runge–Kutta”.

Fie o problemă cu valori inițiale:

${\frac {dy}{dt}}=f(t,y),\quad y(t_{0})=y_{0}$

Aici $y$ este o funcție necunoscută (scalară sau vectorială) de timp, $t,$ care se dorește să fie aproximată. Se știe că viteza de variație ${\frac {dy}{dt}}$ este o funcție de $t$ și de $y.$ La momentul inițial, $t_{0},$ valoarea corespunzătoare a lui $y$ este $y_{0}.$ Sunt date funcția $f$ și condițiile inițiale $t_{0},$ $y_{0}.$

Se alege pasul $h>0$ și se definesc:

${\begin{aligned}y_{n+1}&=y_{n}+{\frac {h}{6}}\left(k_{1}+2k_{2}+2k_{3}+k_{4}\right)\\t_{n+1}&=t_{n}+h\\\end{aligned}}$

pentru n = 0, 1, 2, 3, ..., folosind

${\begin{aligned}k_{1}&=\ f(t_{n},y_{n})\\k_{2}&=\ f\!\left(t_{n}+{\frac {h}{2}},y_{n}+k_{1}{\frac {h}{2}}\right)\\k_{3}&=\ f\!\left(t_{n}+{\frac {h}{2}},y_{n}+k_{2}{\frac {h}{2}}\right)\\k_{4}&=\ f\!\left(t_{n}+h,y_{n}+hk_{3}\right).\end{aligned}}$

(Notă: Diferiți autori folosesc definiții diferite, dar echivalente.)

Aici $y_{n+1}$ este aproximarea RK4 a $y(t_{n+1}),$ a cărei valoare este determinată de valoarea actuală, $y_{n},$ plus media ponderată a patru incremente, unde fiecare increment este produsul dintre mărimea intervalului, h, și o pantă estimată specificată de funcția f în membrul drept al ecuației diferențiale.

• $k_{1}$ este panta la începutul intervalului, folosind $y$ (metoda Euler);
• $k_{2}$ este panta în mijlocul intervalului, folosind $y$ și $k_{1};$
• $k_{3}$ este tot panta la începutul intervalului, dar folosind acum $y$ și $k_{2};$
• $k_{4}$ este panta la sfârșitul intervalului, folosind $y$ și $k_{3}.$

La calcularea mediei celor patru pante se acordă o pondere mai mare pantelor din mijloc. Dacă f este independentă de y, astfel încât ecuația diferențială este echivalentă cu o integrală simplă, atunci RK4 este metoda Simpson^⁠(d).

Metoda RK4 este o metodă de ordinul patru, ceea ce înseamnă că eroarea locală de trunchiere este de ordinul $O(h^{5}),$ în timp ce eroarea totală acumulată este de ordinul $O(h^{4}).$

În multe aplicații practice funcția $f$ este independentă de $t$ (sistem invariant în timp, în special în fizică), iar incrementele lor nu sunt deloc calculate și nu sunt transmise funcției $f,$ care este utilizată doar formula finală pentru $t_{n+1}.$

Metode Runge–Kutta explicite

Familia de metode Runge–Kutta explicite este o generalizare a metodei RK4 menționate mai sus. Este dată de:

$y_{n+1}=y_{n}+h\sum _{i=1}^{s}b_{i}k_{i}$

unde

${\begin{aligned}k_{1}&=f(t_{n},y_{n})\\k_{2}&=f(t_{n}+c_{2}h,y_{n}+(a_{21}k_{1})h)\\k_{3}&=f(t_{n}+c_{3}h,y_{n}+(a_{31}k_{1}+a_{32}k_{2})h)\\&\ \ \vdots \\k_{s}&=f(t_{n}+c_{s}h,y_{n}+(a_{s1}k_{1}+a_{s2}k_{2}+\cdots +a_{s,s-1}k_{s-1})h)\end{aligned}}$

(Notă: Diferiți autori folosesc definiții diferite, dar echivalente.)

Pentru a specifica o anumită metodă trebuie să se furnizeze numărul întreg s (numărul de stadii) și coeficienții a_ij (pentru 1 ≤ j < i ≤ s), b_i (pentru i = 1, 2, ..., s) și c_i (pentru i = 2, 3, ..., s). Matricea [a_ij] se numește matricea Runge–Kutta, în timp ce b_i și c_i sunt cunoscute sub numele de ponderi și noduri. Aceste date sunt de obicei aranjate mnemonic în tabela Butcher (după John C. Butcher):

$0$
$c_{2}$	$a_{21}$
$c_{3}$	$a_{31}$	$a_{32}$
$\vdots$	$\vdots$		$\ddots$
$c_{s}$	$a_{s1}$	$a_{s2}$	$\cdots$	$a_{s,s-1}$
	$b_{1}$	$b_{2}$	$\cdots$	$b_{s-1}$	$b_{s}$

O dezvoltare în serie Taylor arată că metoda Runge–Kutta este consistentă dacă și numai dacă:

$\sum _{i=1}^{s}b_{i}=1$

Există, de asemenea, cerințe însoțitoare dacă se impune ca metoda să aibă un anumit ordin p, ceea ce înseamnă că eroarea locală de trunchiere este O(h^p+1). Acestea pot fi obținute din definiția erorii de trunchiere. De exemplu, o metodă în două stadii are ordinul 2 dacă b₁ + b₂ = 1, b₂c₂ = 1/2 și b₂a₂₁ = 1/2. O condiție populară pentru determinarea coeficienților este:

$\sum _{j=1}^{i-1}a_{ij}=c_{i}\quad {\text{ pentru }}i=2,\ldots ,s$

Totuși, această condiție în sine nu este nici necesară, nici suficientă pentru consistență. Ca un contraexemplu, fie orice schemă Runge–Kutta explicită în 2 stadii cu $b_{1}=b_{2}=1/2$ și $c_{1}$ și $a_{21}$ alese aleatoriu. Această metodă este consistentă și (în general) convergentă de ordinul întâi. Pe de altă parte, metoda într-un singur stadiu cu $b_{1}=1/2$ este inconsistentă și nu reușește să conveargă, chiar dacă $\sum _{j=1}^{i-1}a_{ij}=c_{i}{\text{ pentru }}i=2,\ldots ,s$ este valabilă trivial.

În general, dacă o metodă Runge–Kutta explicită cu s stadii este de ordinul p, atunci se poate demonstra că numărul de stadii trebuie să satisfacă $s\geq p$ și dacă $p\geq 5,$ atunci $s\geq p+1\,.$ Totuși, nu se știe dacă aceste limite sunt exacte în toate cazurile. În unele cazuri se demonstrează că limita nu poate fi atinsă. De exemplu, Butcher a demonstrat că pentru $p>6$ nu există o metodă explicită cu $s=p+1$ stadii. Butcher a demonstrat, de asemenea, că pentru $p>7$ nu există o metodă explicită Runge–Kutta cu $p+2$ stadii. În general numărul minim exact de stadii $s$ pentru ca o metodă explicită Runge–Kutta să fie de ordinul $p$ rămâne o problemă deschisă. Câteva valori cunoscute sunt:

${\begin{array}{c|cccccccc}p&1&2&3&4&5&6&7&8\\\hline \min s&1&2&3&4&6&7&9&11\end{array}}$

Limita demonstrabilă de mai sus implică faptul că nu se pot găsi metode de ordinele $p=1,\,2,\,\ldots ,\,6$ care necesită mai puține stadii decât metodele cunoscute pentru aceste ordine. Lucrarea lui Butcher demonstrează, de asemenea, că metodele de ordinul 7 și 8 au un minim de 9, respectiv 11 stadii. Un exemplu de metodă explicită de ordinul 6 cu 7 stadii poate fi găsit în lucrarea lui Butcher din 1964. Se cunosc și metode explicite de ordinul 7 cu 9 stadii și de ordinul 8 cu 11 stadii. Se găsesc liste cu aceste metode.

Exemple

Metoda RK4 se încadrează în acest cadru. Tabala sa este:

0
1/2	1/2
1/2	0	1/2
1	0	0	1
	1/6	1/3	1/3	1/6

O mică variație din 1901 a metodei Runge–Kutta se datorează tot lui Kutta și se numește regula 3/8. Principalul avantaj al acestei metode este că aproape toți coeficienții de eroare sunt mai mici decât în metoda populară, dar necesită ceva mai multe operații în virgulă mobilă pe un pas de timp. Tabela sa Butcher este:

0
1/3	1/3
2/3	−1/3	1
1	1	−1	1
	1/8	3/8	3/8	1/8

Totuși, cea mai simplă metodă Runge–Kutta este metoda Euler, dată de formula $y_{n+1}=y_{n}+hf(t_{n},y_{n}).$ Aceasta este singura metodă Runge–Kutta explicită consistentă cu un singur stadiu. Tabala corespunzătoare este:

0
	1

Metode de ordinul al doilea cu două stadii

Un exemplu de metodă de ordinul al doilea cu două stadii corespunde metodei explicite a punctului de mijloc^⁠(d):

0
1/2	1/2
	0	1

Metoda punctului de mijloc nu este singura metodă Runge–Kutta de ordinul al doilea cu două stadii; există o familie de astfel de metode, parametrizate de α și date de formula:

$y_{n+1}=y_{n}+h{\bigl (}(1-{\frac {1}{2\alpha }})f(t_{n},y_{n})+{\frac {1}{2\alpha }}f(t_{n}+\alpha h,y_{n}+\alpha hf(t_{n},y_{n})){\bigr )}$

Tabela sa Butcher este:

0
$\alpha$	$\alpha$
	$(1-{\tfrac {1}{2\alpha }})$	${\tfrac {1}{2\alpha }}$

În această familie $\alpha ={\tfrac {1}{2}}$ dă metoda punctului de mijloc, $\alpha =1$ este metoda lui Heun^⁠(d), iar $\alpha ={\tfrac {2}{3}}$ este metoda lui Ralston.

Utilizare

Ca exemplu, fie metoda Runge-Kutta de ordinul al doilea în două stadii cu α = 2/3, cunoscută și sub numele de metoda Ralston. Aceasta este dată de tabela:

	0
	2/3	2/3
		1/4	3/4

cu ecuațiile corespunzătoare:

${\begin{aligned}k_{1}&=f(t_{n},\ y_{n})\\k_{2}&=f(t_{n}+{\frac {2}{3}}h,\ y_{n}+{\frac {2}{3}}hk_{1})\\y_{n+1}&=y_{n}+h\left({\frac {1}{4}}k_{1}+{\frac {3}{4}}k_{2}\right)\end{aligned}}$

Această metodă este utilizată pentru a rezolva problema valorilor inițiale:

${\frac {dy}{dt}}=\tan(y)+1,\quad y_{0}=1,\ t\in [1;\ 1,1]$

cu pasul h = 0,025, astfel metoda necesită patru pași.

Metoda funcționează după cum urmează:

$t_{0}=1\colon$
	$y_{0}=1$
$t_{1}=1,025\colon$
	$y_{0}=1$	$k_{1}=2,557407725$	$k_{2}=f(t_{0}+{\tfrac {2}{3}}h,\ y_{0}+{\tfrac {2}{3}}hk_{1})=2,7138981400$
	$y_{1}=y_{0}+h({\tfrac {1}{4}}k_{1}+{\tfrac {3}{4}}k_{2})={\underline {1,066869388}}$
$t_{2}=1,05\colon$
	$y_{1}=1,066869388$	$k_{1}=2,813524695$	$k_{2}=f(t_{1}+{\tfrac {2}{3}}h,\ y_{1}+{\tfrac {2}{3}}hk_{1})$
	$y_{2}=y_{1}+h({\tfrac {1}{4}}k_{1}+{\tfrac {3}{4}}k_{2})={\underline {1,141332181}}$
$t_{3}=1,075\colon$
	$y_{2}=1,141332181$	$k_{1}=3,183536647$	$k_{2}=f(t_{2}+{\tfrac {2}{3}}h,\ y_{2}+{\tfrac {2}{3}}hk_{1})$
	$y_{3}=y_{2}+h({\tfrac {1}{4}}k_{1}+{\tfrac {3}{4}}k_{2})={\underline {1,227417567}}$
$t_{4}=1,1\colon$
	$y_{3}=1,227417567$	$k_{1}=3,796866512$	$k_{2}=f(t_{3}+{\tfrac {2}{3}}h,\ y_{3}+{\tfrac {2}{3}}hk_{1})$
	$y_{4}=y_{3}+h({\tfrac {1}{4}}k_{1}+{\tfrac {3}{4}}k_{2})={\underline {1,335079087}}$

Soluțiile numerice corespund valorilor subliniate.

Metode Runge–Kutta implicite

Metodele Runge–Kutta explicite sunt în general nepotrivite pentru rezolvarea ecuațiilor rigide^⁠(d) deoarece zona lor de stabilitate absolută este mică; în special, este mărginită. Această problemă este deosebit de importantă în rezolvarea ecuațiilor cu derivate parțiale.

Instabilitatea metodelor Runge–Kutta explicite a motivat dezvoltarea metodelor implicite. O metodă Runge–Kutta implicită are forma:

$y_{n+1}=y_{n}+h\sum _{i=1}^{s}b_{i}k_{i}$

unde

$k_{i}=f\left(t_{n}+c_{i}h,\ y_{n}+h\sum _{j=1}^{s}a_{ij}k_{j}\right),\quad i=1,\,\ldots ,\,s$

Diferența față de o metodă explicită este că, într-o metodă explicită, suma pentru j merge doar până la i − 1. Acest lucru apare și în tabela Butcher: matricea coeficienților $a_{ij}$ a unei metode explicite este triunghiulară inferior. Într-o metodă implicită suma pentru j merge până la s, iar matricea coeficienților nu este strict triunghiulară, rezultând o tabelă Butcher de forma:

${\begin{array}{c|cccc}c_{1}&a_{11}&a_{12}&\dots &a_{1s}\\c_{2}&a_{21}&a_{22}&\dots &a_{2s}\\\vdots &\vdots &\vdots &\ddots &\vdots \\c_{s}&a_{s1}&a_{s2}&\dots &a_{ss}\\\hline &b_{1}&b_{2}&\dots &b_{s}\\\end{array}}={\begin{array}{c|c}\mathbf {c} &A\\\hline &\mathbf {b^{T}} \\\end{array}}$

Consecința acestei diferențe este că la fiecare pas trebuie rezolvat un sistem de ecuații algebrice. Aceasta crește considerabil costul de calcul. Dacă se utilizează o metodă cu s stadii pentru a rezolva o ecuație diferențială cu m componente, atunci sistemul de ecuații algebrice are ms componente. Acest lucru poate fi comparată cu metoda liniară cu pași multipli implicită (cealaltă mare familie de metode pentru EDO): o metodă liniară cu pași multipli implicită cu s pași trebuie să rezolve un sistem de ecuații algebrice cu doar m componente, astfel dimensiunea sistemului nu crește pe măsură ce numărul de pași crește.

Exemple

Cel mai simplu exemplu de metodă Runge–Kutta implicită este metoda Euler înapoi^⁠(d):

$y_{n+1}=y_{n}+hf(t_{n}+h,\ y_{n+1}).\,$

Tabela Butcher pentru asta este simplă:

${\begin{array}{c|c}1&1\\\hline &1\\\end{array}}$

Această tabelă Butcher corespunde formulei:

$k_{1}=f(t_{n}+h,\ y_{n}+hk_{1})\quad {\text{ și }}\quad y_{n+1}=y_{n}+hk_{1}$

care poate fi rearanjată pentru a obține formula pentru metoda Euler înapoi menționată mai sus.

Un alt exemplu de metodă Runge–Kutta implicită este regula trapezului^⁠(d). Tabela sa Butcher este:

${\begin{array}{c|cc}0&0&0\\1&{\frac {1}{2}}&{\frac {1}{2}}\\\hline &{\frac {1}{2}}&{\frac {1}{2}}\\&1&0\\\end{array}}$

Metoda Gauss–Legendre^⁠(d) cu s stadii are ordinul 2s (astfel se pot construi metode de ordin arbitrar de mare). Metoda cu două stadii (deci de ordinul 4) are tabela Butcher:

${\begin{array}{c|cc}{\frac {1}{2}}-{\frac {1}{6}}{\sqrt {3}}&{\frac {1}{4}}&{\frac {1}{4}}-{\frac {1}{6}}{\sqrt {3}}\\{\frac {1}{2}}+{\frac {1}{6}}{\sqrt {3}}&{\frac {1}{4}}+{\frac {1}{6}}{\sqrt {3}}&{\frac {1}{4}}\\\hline &{\frac {1}{2}}&{\frac {1}{2}}\\&{\frac {1}{2}}+{\frac {1}{2}}{\sqrt {3}}&{\frac {1}{2}}-{\frac {1}{2}}{\sqrt {3}}\end{array}}$

Stabilitate

Avantajul metodelor Runge–Kutta implicite față de cele explicite constă în stabilitatea lor mai mare, în special atunci când sunt aplicate la ecuații rigide. Fie ecuația de test liniară $y\,'=\lambda y.$ O metodă Runge–Kutta aplicată acestei ecuații se reduce la iterația $y_{n+1}=r(h\lambda )\,y_{n}$, cu r dat de:

$r(z)=1+zb^{T}(I-zA)^{-1}e={\frac {\det(I-zA+zeb^{T})}{\det(I-zA)}}$

unde e este vectorul unitate. Funcția r este funcția de stabilitate. Din formulă rezultă că r este câtul a două polinoame de gradul s dacă metoda are s stadii. Metodele explicite au o matrice triunghiulară strict inferioară A, ceea ce implică faptul că det(I − zA) = 1 și că funcția de stabilitate este un polinom.

Soluția numerică a ecuației liniare de test descrește la zero dacă $|r(z)|<1\ {\text{ cu }}z=h\lambda .$ Mulțimea unor astfel de z este domeniul stabilității absolute. Metoda este absolut stabilă dacă toate z cu Re(z) < 0 se află în domeniul stabilității absolute. Funcția de stabilitate a unei metode Runge–Kutta explicite este un polinom, deci metodele Runge–Kutta explicite nu pot fi niciodată A-stabile.

Dacă metoda este de ordinul p, atunci funcția de stabilitate satisface $r(z)={\textrm {e}}^{z}+O(z^{p+1})$ când $z\to 0\,.$ Prin urmare, este interesant să se studieze coeficienții polinoamelor de grad dat care aproximează cel mai bine funcția exponențială. Acestea sunt cunoscute sub numele de aproximare Padé^⁠(d). O aproximare Padé cu numărător de gradul m și numitor de gradul n este A-stabilă dacă și numai dacă m ≤ n ≤ m + 2.

Metoda Gauss–Legendre cu s stadii este de ordinul al doilea, deci funcția sa de stabilitate este aproximarea Padé cu m = n = s. Rezultă că metoda este A-stabilă. Aceasta arată că metoda Runge–Kutta A-stabilă poate avea un ordin oricât de mare. În schimb, ordinul metodelor A-stabile liniare cu mai mulți pași nu poate fi mai mare ca doi.

Metode Runge–Kutta adaptive

Metodele adaptive sunt concepute pentru a produce o estimare a erorii locale de trunchiere a unui singur pas Runge–Kutta. Acest lucru se realizează prin utilizarea a două metode, una de ordinul $p$ și una de ordinul $p-1\,.$ Aceste metode au pași intermediari comuni. Datorită acestui fapt, estimarea erorii are un cost computațional mic sau neglijabil în comparație cu un pas din metoda de ordin superior.

În timpul integrării dimensiunea pasului este adaptată astfel încât eroarea estimată să rămână sub un prag definit de utilizator: Dacă eroarea este prea mare, un pas este repetat cu un pas mai mic. Dacă eroarea este mult mai mică, pasul poate fi crescut pentru a economisi timp. Acest lucru are ca rezultat o dimensiune a pasului (aproape) optimă, ceea ce economisește timp de calcul. Mai mult, utilizatorul nu trebuie să consume timpul său pentru a găsi o dimensiune a pasului adecvată.

Pasul de ordin inferior este dat de:

$y_{n+1}^{*}=y_{n}+h\sum _{i=1}^{s}b_{i}^{*}k_{i}$

unde $k_{i}$ sunt aceliași ca la metoda de ordin superior. Atunci eroarea este:

$e_{n+1}=y_{n+1}-y_{n+1}^{*}=h\sum _{i=1}^{s}(b_{i}-b_{i}^{*})k_{i}$

care este $O(h^{p}).$ Tabela Butcher pentru acest tip de metodă este extins pentru a da valorile lui $b_{i}^{*}$: ${\begin{array}{c|cccc}c_{1}&a_{11}&a_{12}&\dots &a_{1s}\\c_{2}&a_{21}&a_{22}&\dots &a_{2s}\\\vdots &\vdots &\vdots &\ddots &\vdots \\c_{s}&a_{s1}&a_{s2}&\dots &a_{ss}\\\hline &b_{1}&b_{2}&\dots &b_{s}\\&b_{1}^{*}&b_{2}^{*}&\dots &b_{s}^{*}\\\end{array}}$

Metoda Runge–Kutta–Fehlberg^⁠(d) are două variante, de ordinul 5 și 4. Tabela sa Butcher extinsă este:

	0
	1/4	1/4
	3/8	3/32	9/32
	12/13	1932/2197	−7200/2197	7296/2197
	1	439/216	−8	3680/513	-845/4104
	1/2	−8/27	2	−3544/2565	1859/4104	−11/40
		16/135	0	6656/12825	28561/56430	−9/50	2/55
		25/216	0	1408/2565	2197/4104	−1/5	0

Totuși, cea mai simplă metodă adaptivă Runge–Kutta combină metoda lui Heun, care este de ordinul al doilea, cu metoda lui Euler, care este de ordinul întâi. Tabela sa Butcher extinsă este:

	0
	1	1
		1/2	1/2
		1	0

Alte metode adaptive Runge–Kutta sunt metoda Bogacki–Shampine^⁠(d) (de ordinele 3 și 2), metoda Cash–Karp^⁠(d) și metoda Dormand–Prince^⁠(d) (ambele de ordinele 5 și 4).

Metode Runge-Kutta neconfluente

Despre o metodă Runge-Kutta se spune că este „neconfluentă” dacă $c_{i},\,i=1,2,\ldots ,s$ sunt distincte.

Metode Runge–Kutta–Nyström

Metodele Runge–Kutta–Nyström (RKN) sunt o familie de metode bazate pe aceleași principii ca și metodele Runge–Kutta, dar pentru probleme cu valori inițiale de ordinul al doilea, adică pentru probleme de forma:

${\frac {d^{2}y}{dt^{2}}}=f(t,{\frac {dy}{dt}},y),\quad y(t_{0})=y_{0},\quad {\frac {dy}{dt}}(t_{0})=y\,'_{0}$

Există două variante de aproximare. O metodă Runge–Kutta–Nyström utilizează două matrici Runge–Kutta, $a_{ij},a'_{ij},$ și două seturi de ponderare $b_{i},b'_{i},$ dar are nevoie de doar un singur set de noduri, $c_{i}.$ Aceasta produce o tabelă Butcher de forma:

${\begin{array}{c|cccc}c_{1}&a_{11}&a_{12}&\dots &a_{1s}\\c_{2}&a_{21}&a_{22}&\dots &a_{2s}\\\vdots &\vdots &\vdots &\ddots &\vdots \\c_{s}&a_{s1}&a_{s2}&\dots &a_{ss}\\\hline &a'_{11}&a'_{12}&\dots &a'_{1s}\\&a'_{21}&a'_{22}&\dots &a'_{2s}\\&\vdots &\vdots &\ddots &\vdots \\&a'_{s1}&a'_{s2}&\dots &a'_{ss}\\\hline &b_{1}&b_{2}&\dots &b_{s}\\&b'_{1}&b'_{2}&\dots &b'_{s}\\\end{array}}={\begin{array}{c|c}\mathbf {c} &\mathbf {A} \\\hline &\mathbf {A'} \\\hline &\mathbf {b} ^{\top }\\&\mathbf {b'} ^{\top }\end{array}}$

Se presupune că aproximările au fost efectuate până la $t_{n},$ cu $y_{n}$ aproximarea lui $y(t_{n})$ și $y\,'_{n}$ aproximarea lui ${\frac {dy}{dt}}(t_{n}).$ Aproximările $y_{n+1},y\,'_{n+1}$ la $t_{n+1}=t_{n}+h$ sunt soluțiile următorului sistem:

${\begin{cases}g_{i}=y_{n}+c_{i}hy\,'_{n}+h^{2}\sum _{j=1}^{s}a_{ij}f(t_{n}+c_{j}h,g\,'_{j},g_{j}),&i=1,2,\ldots ,s\\g\,'_{i}=y\,'_{n}+h\sum _{j=1}^{s}a\,'_{ij}f(t_{n}+c_{j}h,g\,'_{j},g_{j}),&i=1,2,\ldots ,s\\\\y_{n+1}=y_{n}+hy\,'_{n}+h^{2}\sum _{j=1}^{s}b_{j}f(t_{n}+c_{j}h,g\,'_{j},g_{j})\\y\,'_{n+1}=y\,'_{n}+h\sum _{j=1}^{s}b\,'_{j}f(t_{n}+c_{j}h,g\,'_{j},g_{j})\end{cases}}$

unde $g_{i},g\,'_{i}$ sunt aproximările intermediare ale lui $y$ și ${\frac {dy}{dt}}.$ Este strict echivalent să se lucreze cu valorile $k_{j}=f(t_{n}+c_{j}h,g\,'_{j},g_{j})$ unde $g_{j},g\,'_{j}$ au fost înlocuite cu formula lor, în loc să se lucreze cu $g_{i},g\,'_{i},$ similar cu metodele Runge–Kutta anterioare, dar sistemul este mai ușor de scris în acest fel.

Metoda Runge–Kutta–Nyström este explicită dacă ambele $A,A'$ sunt strict triunghiulare inferior, și în acest caz, sumele ${\textstyle \sum _{j=1}^{s}}$ din expresiile lui $g_{i},g\,'_{i}$, pot fi înlocuite cu ${\textstyle \sum _{j=1}^{i-1}}$. În plus, o metodă Runge–Kutta–Nyström este de ordinul $p$ dacă eroarea locală de trunchiere a ambelor $y_{n+1},y\,'_{n+1}$ este $O(h^{p+1}).$

Dacă funcția $f$ a problemei cu valori inițiale considerate este independentă de ${\frac {dy}{dt}},$ nu este nevoie să se aproximeze valorile intermediare $g\,'_{i}$ pentru a calcula aproximările, prin urmare, ponderile $a\,'_{ij}$ sunt inutile și, în schimb, se poate scrie o metodă creată doar pentru acest caz particular folosind o tabelă de forma:

${\begin{array}{c|cccc}c_{1}&a_{11}&a_{12}&\dots &a_{1s}\\c_{2}&a_{21}&a_{22}&\dots &a_{2s}\\\vdots &\vdots &\vdots &\ddots &\vdots \\c_{s}&a_{s1}&a_{s2}&\dots &a_{ss}\\\hline &b_{1}&b_{2}&\dots &b_{s}\\&b'_{1}&b'_{2}&\dots &b'_{s}\\\end{array}}={\begin{array}{c|c}\mathbf {c} &\mathbf {A} \\\hline &\mathbf {b} ^{\top }\\&\mathbf {b'} ^{\top }\end{array}}$

Acest caz particular este deosebit de interesant, deoarece permite un ordin mai mare decât ceea ce poate realiza o metodă Runge–Kutta–Nyström în general. De exemplu, două metode RKN explicite de ordinul patru sunt date de următoarea tabelă Butcher:

${\begin{array}{c|ccc}c_{i}&&a_{ij}&\\{\frac {3+{\sqrt {3}}}{6}}&0&0&0\\{\frac {3-{\sqrt {3}}}{6}}&{\frac {2-{\sqrt {3}}}{12}}&0&0\\{\frac {3+{\sqrt {3}}}{6}}&0&{\frac {\sqrt {3}}{6}}&0\\\hline b_{i}&{\frac {3-2{\sqrt {3}}}{12}}&{\frac {1}{2}}&{\frac {3+2{\sqrt {3}}}{12}}\\\hline b'_{i}&{\frac {5-3{\sqrt {3}}}{24}}&{\frac {3+{\sqrt {3}}}{12}}&{\frac {1+{\sqrt {3}}}{24}}\\\end{array}}$

${\begin{array}{c|ccc}c_{i}&&a_{ij}&\\{\frac {3-{\sqrt {3}}}{6}}&0&0&0\\{\frac {3+{\sqrt {3}}}{6}}&{\frac {2+{\sqrt {3}}}{12}}&0&0\\{\frac {3-{\sqrt {3}}}{6}}&0&-{\frac {\sqrt {3}}{6}}&0\\\hline b_{i}&{\frac {3+2{\sqrt {3}}}{12}}&{\frac {1}{2}}&{\frac {3-2{\sqrt {3}}}{12}}\\\hline b'_{i}&{\frac {5+3{\sqrt {3}}}{24}}&{\frac {3-{\sqrt {3}}}{12}}&{\frac {1-{\sqrt {3}}}{24}}\\\end{array}}$

Aceste două scheme au și proprietăți de conservare simplectică atunci când ecuația originală este obținută dintr-un sistem mecanic clasic conservativ, adică atunci când

$f_{i}(x_{1},\ldots ,x_{n})={\frac {\partial V}{\partial x_{i}}}(x_{1},\,\ldots ,\,x_{n})$

pentru funcția scalară $V.$

B-stabilitate

Conceptul de A-stabilitate la rezolvarea ecuațiilor diferențiale este legat de ecuația liniară autonomă $y\,'=\lambda y.$ Dahlquist a propus cercetarea stabilității schemelor numerice atunci când este aplicată sistemelor neliniare care satisfac o condiție de monotonie. Conceptele corespunzătoare au fost definite ca G-stabilitate pentru metodele cu mai mulți pași (și metodele într-un singur pas aferente) și B-stabilitate pentru metodele Runge–Kutta. O metodă Runge–Kutta aplicată sistemului neliniar $y\,'=f(y),$ care verifică $\langle f(y)-f(z),y-z\rangle \leq 0,$ este B-stabilă dacă această condiție satisface relația $\|y_{n+1}-z_{n+1}\|\leq \|y_{n}-z_{n}\|$ pentru două soluții numerice.

Fie $B,$ $M$ și $Q$ trei matrici $s\times s$ definite astfel:

${\begin{aligned}B&=\operatorname {diag} (b_{1},b_{2},\ldots ,b_{s})\\[4pt]M&=BA+A^{T}B-bb^{T}\\[4pt]Q&=BA^{-1}+A^{-T}B-A^{-T}bb^{T}A^{-1}\end{aligned}}$

O metodă Runge–Kutta este stabilă algebric dacă matricile $B$ și $M$ sunt ambele pozitiv definite. O condiție suficientă pentru B-stabilitate este: $B$ și $Q$ sunt pozitiv definite.

Obținerea metodei Runge-Kutta de ordinul patru

În general, o metodă Runge–Kutta de ordinul s poate fi scrisă astfel:

$y_{t+h}=y_{t}+h\cdot \sum _{i=1}^{s}a_{i}k_{i}+{\mathcal {O}}(h^{s+1})$

unde:

$k_{i}=\sum _{j=1}^{s}\beta _{ij}f(k_{j},\ t_{n}+\alpha _{i}h)$

sunt incremente obținute calculând derivatele lui $y_{t}$ de ordinul $i$.

Se aplică procedura pentru metoda Runge–Kutta de ordinul patru folosind formula generală cu $s=4$ evaluată, așa cum s-a explicat mai sus, în punctul de pornire, la mijloc și în punctul final al oricărui interval $(t,\ t+h).$ Se alege:

${\begin{aligned}&\alpha _{i}&&\beta _{ij}\\\alpha _{1}&=0&\beta _{21}&={\frac {1}{2}}\\\alpha _{2}&={\frac {1}{2}}&\beta _{32}&={\frac {1}{2}}\\\alpha _{3}&={\frac {1}{2}}&\beta _{43}&=1\\\alpha _{4}&=1&&\\\end{aligned}}$

și $\beta _{ij}=0$ altfel. Se definesc următoarele mărimi:

${\begin{aligned}y_{t+h}^{1}&=y_{t}+hf\left(y_{t},\ t\right)\\y_{t+h}^{2}&=y_{t}+hf\left(y_{t+h/2}^{1},\ t+{\frac {h}{2}}\right)\\y_{t+h}^{3}&=y_{t}+hf\left(y_{t+h/2}^{2},\ t+{\frac {h}{2}}\right)\end{aligned}}$

unde $y_{t+h/2}^{1}={\dfrac {y_{t}+y_{t+h}^{1}}{2}}$ și $y_{t+h/2}^{2}={\dfrac {y_{t}+y_{t+h}^{2}}{2}}.$ dacă se definește:

${\begin{aligned}k_{1}&=f(y_{t},\ t)\\k_{2}&=f\left(y_{t+h/2}^{1},\ t+{\frac {h}{2}}\right)=f\left(y_{t}+{\frac {h}{2}}k_{1},\ t+{\frac {h}{2}}\right)\\k_{3}&=f\left(y_{t+h/2}^{2},\ t+{\frac {h}{2}}\right)=f\left(y_{t}+{\frac {h}{2}}k_{2},\ t+{\frac {h}{2}}\right)\\k_{4}&=f\left(y_{t+h}^{3},\ t+h\right)=f\left(y_{t}+hk_{3},\ t+h\right)\end{aligned}}$

iar pentru relațiile anterioare se poate demonstra că următoarele egalități sunt valabile pentru ${\mathcal {O}}(h^{2})$:

${\begin{aligned}k_{2}&=f\left(y_{t+h/2}^{1},\ t+{\frac {h}{2}}\right)=f\left(y_{t}+{\frac {h}{2}}k_{1},\ t+{\frac {h}{2}}\right)\\&=f\left(y_{t},\ t\right)+{\frac {h}{2}}{\frac {d}{dt}}f\left(y_{t},\ t\right)\\k_{3}&=f\left(y_{t+h/2}^{2},\ t+{\frac {h}{2}}\right)=f\left(y_{t}+{\frac {h}{2}}f\left(y_{t}+{\frac {h}{2}}k_{1},\ t+{\frac {h}{2}}\right),\ t+{\frac {h}{2}}\right)\\&=f\left(y_{t},\ t\right)+{\frac {h}{2}}{\frac {d}{dt}}\left[f\left(y_{t},\ t\right)+{\frac {h}{2}}{\frac {d}{dt}}f\left(y_{t},\ t\right)\right]\\k_{4}&=f\left(y_{t+h}^{3},\ t+h\right)=f\left(y_{t}+hf\left(y_{t}+{\frac {h}{2}}k_{2},\ t+{\frac {h}{2}}\right),\ t+h\right)\\&=f\left(y_{t}+hf\left(y_{t}+{\frac {h}{2}}f\left(y_{t}+{\frac {h}{2}}f\left(y_{t},\ t\right),\ t+{\frac {h}{2}}\right),\ t+{\frac {h}{2}}\right),\ t+h\right)\\&=f\left(y_{t},\ t\right)+h{\frac {d}{dt}}\left[f\left(y_{t},\ t\right)+{\frac {h}{2}}{\frac {d}{dt}}\left[f\left(y_{t},\ t\right)+{\frac {h}{2}}{\frac {d}{dt}}f\left(y_{t},\ t\right)\right]\right]\end{aligned}}$

unde:

${\frac {d}{dt}}f(y_{t},\ t)={\frac {\partial }{\partial y}}f(y_{t},\ t){\dot {y}}_{t}+{\frac {\partial }{\partial t}}f(y_{t},\ t)=f_{y}(y_{t},\ t){\dot {y}}_{t}+f_{t}(y_{t},\ t):={\ddot {y}}_{t}$

este derivata totală a lui $f$ în funcție de timp.

Dacă se exprimă formula generală folosind cele de mai sus se obține:

${\begin{aligned}y_{t+h}={}&y_{t}+h\left\lbrace a\cdot f(y_{t},\ t)+b\cdot \left[f(y_{t},\ t)+{\frac {h}{2}}{\frac {d}{dt}}f(y_{t},\ t)\right]\right.+\\&{}+c\cdot \left[f(y_{t},\ t)+{\frac {h}{2}}{\frac {d}{dt}}\left[f\left(y_{t},\ t\right)+{\frac {h}{2}}{\frac {d}{dt}}f(y_{t},\ t)\right]\right]+\\&{}+d\cdot \left[f(y_{t},\ t)+h{\frac {d}{dt}}\left[f(y_{t},\ t)+{\frac {h}{2}}{\frac {d}{dt}}\left[f(y_{t},\ t)+\left.{\frac {h}{2}}{\frac {d}{dt}}f(y_{t},\ t)\right]\right]\right]\right\rbrace +{\mathcal {O}}(h^{5})\\={}&y_{t}+a\cdot hf_{t}+b\cdot hf_{t}+b\cdot {\frac {h^{2}}{2}}{\frac {df_{t}}{dt}}+c\cdot hf_{t}+c\cdot {\frac {h^{2}}{2}}{\frac {df_{t}}{dt}}+\\&{}+c\cdot {\frac {h^{3}}{4}}{\frac {d^{2}f_{t}}{dt^{2}}}+d\cdot hf_{t}+d\cdot h^{2}{\frac {df_{t}}{dt}}+d\cdot {\frac {h^{3}}{2}}{\frac {d^{2}f_{t}}{dt^{2}}}+d\cdot {\frac {h^{4}}{4}}{\frac {d^{3}f_{t}}{dt^{3}}}+{\mathcal {O}}(h^{5})\end{aligned}}$

și comparând cu seria Taylor a lui $y_{t+h}$ în jurul lui $t$:

${\begin{aligned}y_{t+h}&=y_{t}+h{\dot {y}}_{t}+{\frac {h^{2}}{2}}{\ddot {y}}_{t}+{\frac {h^{3}}{6}}y_{t}^{(3)}+{\frac {h^{4}}{2}}4y_{t}^{(4)}+{\mathcal {O}}(h^{5})=\\&=y_{t}+hf(y_{t},\ t)+{\frac {h^{2}}{2}}{\frac {d}{dt}}f(y_{t},\ t)+{\frac {h^{3}}{6}}{\frac {d^{2}}{dt^{2}}}f(y_{t},\ t)+{\frac {h^{4}}{24}}{\frac {d^{3}}{dt^{3}}}f(y_{t},\ t)\end{aligned}}$

se obține sistemul constrângerilor asupra coeficiențior:

${\begin{cases}&a+b+c+d=1\\[6pt]&{\frac {1}{2}}b+{\frac {1}{2}}c+d={\frac {1}{2}}\\[6pt]&{\frac {1}{4}}c+{\frac {1}{2}}d={\frac {1}{6}}\\[6pt]&{\frac {1}{4}}d={\frac {1}{24}}\end{cases}}$

care rezolvat dă: $~a={\frac {1}{6}},~b={\frac {1}{3}},~c={\frac {1}{3}},~d={\frac {1}{6}},~$ cum se spune mai sus.

Bibliografie

• en Press, William H.; Teukolsky, Saul A.; Vetterling, William T.; Flannery, Brian P. (2007), „Section 17.1 Runge-Kutta Method”, Numerical Recipes: The Art of Scientific Computing (ed. 3rd), Cambridge University Press, ISBN 978-0-521-88068-8 . Also, Section 17.2. Adaptive Stepsize Control for Runge-Kutta.
• en Süli, Endre; Mayers, David (2003), An Introduction to Numerical Analysis, Cambridge University Press, ISBN 0-521-00794-1 .
• en Atkinson, Kendall A. (1989), An Introduction to Numerical Analysis (ed. 2nd), New York: John Wiley & Sons, ISBN 978-0-471-50023-0 .
• en Hairer, Ernst; Nørsett, Syvert Paul; Wanner, Gerhard (1993), Solving ordinary differential equations I: Nonstiff problems, Berlin, New York: Springer-Verlag, ISBN 978-3-540-56670-0 .
• en Kaw, Autar; Kalu, Egwu (2008), Numerical Methods with Applications (ed. 1st), autarkaw.com, arhivat din original la 13 aprilie 2009, accesat în 5 august 2026 .
• en Stoer, Josef; Bulirsch, Roland (2002), Introduction to Numerical Analysis (ed. 3rd), Berlin, New York: Springer-Verlag, ISBN 978-0-387-95452-3 .
• en Ascher, Uri M.; Petzold, Linda R. (1998), Computer Methods for Ordinary Differential Equations and Differential-Algebraic Equations, Philadelphia: Society for Industrial and Applied Mathematics, ISBN 978-0-89871-412-8 .
• en Butcher, John C. (2008), Numerical Methods for Ordinary Differential Equations, New York: John Wiley & Sons, ISBN 978-0-470-72335-7 .
• en Iserles, Arieh (1996), A First Course in the Numerical Analysis of Differential Equations, Cambridge University Press, Bibcode:1996fcna.book.....I, ISBN 978-0-521-55655-2 .
• en Butcher, J.C. (1996), „A History of Runge-Kutta Methods”, Applied Numerical Mathematics, 20 (3): 247–260, doi:10.1016/0168-9274(95)00108-5 .
• en Butcher, J.C. (1985), „The non-existence of ten stage eighth order explicit Runge-Kutta methods”, BIT Numerical Mathematics, 25 (3): 521–540, doi:10.1007/BF01935372 .
• en Butcher, J.C. (1965), „On the attainable order of Runge-Kutta methods”, Mathematics of Computation, 19 (91): 408–417, doi:10.1090/S0025-5718-1965-0179943-X .
• en Curtis, A.R. (1970), „An eighth order Runge-Kutta process with eleven function evaluations per step”, Numerische Mathematik, 16 (3): 268–277, doi:10.1007/BF02219778 .
• en Dahlquist, Germund (1963), „A special stability problem for linear multistep methods”, BIT, 3: 27–43, doi:10.1007/BF01963532, hdl:10338.dmlcz/103497 , ISSN 0006-3835 .
• en Butcher, John C. (mai 1964), „On Runge-Kutta processes of high order”, Journal of the Australian Mathematical Society, 4 (2): 179–194, doi:10.1017/S1446788700023387  
• en Hairer, Ernst; Wanner, Gerhard (1996), Solving ordinary differential equations II: Stiff and differential-algebraic problems (ed. 2nd), Berlin, New York: Springer-Verlag, ISBN 978-3-540-60452-5 .
• en Lambert, J.D (1991), Numerical Methods for Ordinary Differential Systems. The Initial Value Problem, John Wiley & Sons, ISBN 0-471-92990-5 
• en Butcher, John C. (1975), „A stability property of implicit Runge-Kutta methods”, BIT, 15 (4): 358–361, doi:10.1007/bf01931672 .

Lectură suplimentară

• de Runge, Carl David Tolmé (1895), „Über die numerische Auflösung von Differentialgleichungen”, Mathematische Annalen, Springer, 46 (2): 167–178, doi:10.1007/BF01446807 .
• de Kutta, Wilhelm (1901), „Beitrag zur näherungsweisen Integration totaler Differentialgleichungen”, Zeitschrift für Mathematik und Physik, 46: 435–453 .
• en Butcher, John C. (mai 1963), „Coefficients for the study of Runge-Kutta integration processes”, Journal of the Australian Mathematical Society, 3 (2): 185–201, doi:10.1017/S1446788700027932  .
• en Butcher, John C. (2000), „Numerical methods for ordinary differential equations in the 20th century”, J. Comput. Appl. Math., 125 (1–2): 1–29, Bibcode:2000JCoAM.125....1B, doi:10.1016/S0377-0427(00)00455-6  .
• en Cellier, F.; Kofman, E. (2006), Continuous System Simulation, Springer Verlag, ISBN 0-387-26102-8 .
• en Forsythe, George E.; Malcolm, Michael A.; Moler, Cleve B. (1977), Computer Methods for Mathematical Computations, Prentice-Hall  (see Chapter 6).
• en Tan, Delin; Chen, Zheng (2012), „On A General Formula of Fourth Order Runge-Kutta Method” (PDF), Journal of Mathematical Science & Mathematics Education, 7 (2): 1–10 .
• en John C. Butcher: "B-Series : Algebraic Analysis of Numerical Methods", Springer(SSCM, volume 55), ISBN: 978-3030709556 (April, 2021).
• en Cooper, G.J.; Verner, J.H. (1972), „Some Explicit Runge–Kutta Methods of High Order”, SIAM Journal on Numerical Analysis, 9 (3): 389–405, Bibcode:1972SJNA....9..389C, doi:10.1137/0709037 .

Legături externe

• Materiale media legate de metoda Runge–Kutta la Wikimedia Commons
• en Hazewinkel, Michiel, ed. (2001), „Runge-Kutta method”, Encyclopaedia of Mathematics, Kluwer Academic Publishers, ISBN 978-1556080104 
• en Runge–Kutta 4th-Order Method Arhivat în 1 septembrie 2006, la Wayback Machine.
• en Tracker Component Library Implementation in Matlab — Implements 32 embedded Runge Kutta algorithms in RungeKStep, 24 embedded Runge-Kutta Nyström algorithms in RungeKNystroemSStep and 4 general Runge-Kutta Nyström algorithms in RungeKNystroemGStep.