METODE NUMERICE

PENTRU
ECUAȚII DIFERENȚIALE ORDINARE
NOTE DE CURS - versiunea 1.0, august 2014

Alexandru I. Nicolin

Departamentul de Fizică Computațională și Tehnologii

Informaționale
Institutul Național de Cercetare-Dezvoltare pentru Fizică și
Inginerie Nucleară „Horia Hulubei”
Măgurele, 2014
Turning to ordinary differential equations, we feel that no one understands
the real nature of a differential equation dy/dx = f (x, y) until he has experi-
mented with computational or graphical procedures for its solution. Such
a simple process as Euler’s method, while not accurate, is both useful and
worthy of study. What normally passes for numerical methods in differential
equations courses bears little resemblance to how such equations are solved
in practice. Why not mix realistic practice right with the theory, where they
can influence and strengthen each other?

The Role of Numerical Analysis in an Undergraduate Program, George E.

Forsythe, The American Mathematical Monthly 66, 651 (1959)

3
Cuprins

Cuvânt liminar 9

1 Introducere 11
1.1 Limbaje de programare și sisteme de calcul . . . . . . . . . . . . . . . . . 12
1.2 Bibliografie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

2 Metode multipas 17
2.1 Metode Adams-Bashforth . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
2.2 Metode Adams-Moulton . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
2.3 Implementări software disponibile . . . . . . . . . . . . . . . . . . . . . . . 25
2.4 Exerciții . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
2.5 Bibliografie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26

3 Metode unipas 27
3.1 Metoda Euler . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
3.2 Metode Runge-Kutta . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
3.2.1 Metode Runge-Kutta explicite . . . . . . . . . . . . . . . . . . . . 30
3.2.2 Metode Runge-Kutta implicite . . . . . . . . . . . . . . . . . . . . 33
3.2.3 Metode Runge-Kutta împerecheate. Controlul erorii . . . . . . . . 38
3.3 Implementări software disponibile . . . . . . . . . . . . . . . . . . . . . . . 43
3.4 Exerciții și probleme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
3.5 Bibliografie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45

5
Listă de tabele

3.1 Metode Runge-Kutta explicite. Tabel Butcher general. . . . . . . . . . . . 30

3.2 Metoda Runge-Kutta explicită (generală) cu s = p = 2 și c2 parametru
liber. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
3.3 Metoda Runge-Kutta explicită (generală) cu s = p = 3. Clasa I de soluții
corespunde metodelor pentru care c2 ̸= 0 ̸= c3 ̸= c2 ̸= 23 . . . . . . . . . . . 32
3.4 Metoda Runge-Kutta explicită (generală) cu s = p = 3. Clasa II de soluții
corespunde metodelor pentru care c2 = c3 = 23 , b3 ̸= 0. . . . . . . . . . . 32
3.5 Metoda Runge-Kutta explicită (generală) cu s = p = 3. Clasa III de
soluții corespunde metodelor pentru care c2 = 23 , c3 = 0, b3 ̸= 0. . . . . . 32
3.6 Metoda Runge-Kutta de ordin 4. . . . . . . . . . . . . . . . . . . . . . . . 34
3.7 Metoda Runge-Kutta de ordin 4, cunoscută drept „regula 3/8”. . . . . . 34
3.8 Metodă Runge-Kutta de ordin 5 determinată de Kutta. . . . . . . . . . . 35
3.9 Metodă Runge-Kutta de ordin 5 determinată de Kutta. Metoda este
incorectă, variantă corectă a coeficienților fiind determintă de Nyström în
1925 (a se vedea tabelul Butcher (3.10)). . . . . . . . . . . . . . . . . . . 35
3.10 Metodă Runge-Kutta de ordin 5 corectată de Nyström. . . . . . . . . . . 36
3.11 Metodă Runge-Kutta bazată pe cuadraturi gaussiene cu s = 1 și p = 2. . 37
3.12 Metodă Runge-Kutta bazată pe cuadraturi gaussiene cu s = 2 și p = 4. . 37
3.13 Metodă Runge-Kutta bazată pe cuadraturi gaussiene cu s = 3 și p = 4. . 37
3.14 Metode Runge-Kutta împerecheate. Prima metodă, cea superioară, este
de ordin p = 2 în timp ce a doua, cea inferioară, este de ordin p = 3. . . 39
3.15 Metode Runge-Kutta-Merson împerecheate. Prima metodă, cea superi-
oară, este de ordin p = 3 în timp ce a doua, cea inferioară, este de ordin
p = 4. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
3.16 Metode Runge-Kutta-Zonneveld împerecheate. Prima metodă, cea supe-
rioară, este de ordin p = 4 în timp ce a doua, cea inferioară, este de ordin
p = 5. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
3.17 Metode Runge-Kutta triplu împerecheate făcute pe baza tabelului But-
cher (3.3). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
3.18 Metode Runge-Kutta triplu împerecheate făcute pe baza tabelului But-
cher (3.4). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
3.19 Metode Runge-Kutta triplu împerecheate făcute pe baza tabelului But-
cher (3.5). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42

7
Cuvânt liminar

Cursurile tradiționale de analiză matematică discută în partea de început aspectele ge-

nerale ce țin de teoria funcțiilor continue și de limita unei funcții într-un punct, pentru
ca mai apoi sa fie introdusă derivarea și integrarea funcțiilor. Până și titlul generic Cal-
cul diferențial și integral induce studentului ideea că derivarea vine, într-un fel, înaintea
integrării, chiar daca dezvoltarea istorica a analizei matematice înregistreaza întâi con-
tribuțiilor remarcabile ale lui Kepler cu privire la calculul volumului butoaielor de vin
cuprinse în Nova stereometria doliorum vinariorum (1615) și abia apoi pe acelea ale lui
Newton din Annotations from Wallis (1665) cu privire la derivate. De fapt, predarea cu-
rentă a analizei matematice urmează firul invers al dezvoltarii ei istorice, subiectele din
finalul cursurilor de profil fiind, de obicei, cele care au fost investigate cel mai devreme.
Așa se face că metodele numerice și, într-un cadru mai larg, analiza numerică ajung sa
fie studiate la finalul cursurilor de analiză matematică, deseori de-o maniera superficială,
chiar dacă o abordare diacronică a domeniului arată că dezvoltarea acestuia ar fi fost
imposibilă fără utilizarea calcului numeric.

Prezentele note de curs prezintă principalele metode folosite în determinarea soluțiilor

numerice ale ecuațiilor diferențiale ordinare, îmbinând aspectele istorice (extrem de uti-
le întelegerii intime a domeniului) cu cele aplicative și, cumpătat, cu cele de natură
teoretică. Cursul nu se adreseaza, așadar, matematicienilor și, în general, acelora care
doresc o tratare riguroasă a domeniului, ci fizicienilor de toate specializările care doresc
să înțeleagă la nivel calitativ proprietățile metodelor numerice folosite în studiul ecuați-
ilor diferențiale și să poată utiliza cu ușurință informațiile dobândite pentru rezolvarea
ecuațiilor ce descriu sisteme fizice reale. Capitolele care urmează prezintă în detaliu cele
mai folosite metode multipas și unipas, acceptul fiind pus pe înțelegerea proprietăților
generale ale acestora (în special pe acelea legate de stabilitate și calcului erorilor) și pe
particularitățile ce privesc implementarea numerică a acestora.

În final autorului mulțumește colegilor din Departamentul de Fizică Computațională

și Tehnologii Informaționale al Institutului Național pentru Fizică și Inginerie Nuclea-
ră „Horia Hulubei” pentru comentariile și sugestiile lor asupra textului și studenților
Facultății de Fizică a Universității din București ale căror observații au sporit conside-
rabil claritatea textului. Mulțumirile autorului se îndreaptă, de asemenea, către Prof.

9
Cuvânt liminar

Univ. Dr. Virgil BĂRAN, șeful departamantului de Fizică teoretică, Matematici, Op-
tică, Spectroscopie, Plasmă, Laseri, pentru oportunitatea de a preda studenților înscriși
în programul masteral Fizică teoretică și computațională prezentele capitole de metode
numerice.

Orice comentariu cu privire la prezentele note de curs pot fi trimise autorului pe adresa
[Link]@[Link].

Alexandru Nicolin,
Măgurele, august, 2014

10
1 Introducere
Ideea de bază a acestui curs este că pentru o ecuație diferențială de tipul

y ′ = f (x, y), y (x0 ) = y0 (1.1)

există două tipuri distincte de metode numerice cu ajutorul cărora obținem aproximația
numerică a soluției ecuației precedente. Pentru o discretizare a lui x de tipul x0 , x1 =
x0 + h, x2 = x1 + h = x0 + 2h, x3 = x2 + h = x1 + 2h = x0 + 3h, etc., ambele tipuri
de metode numerice determină valorile lui y pentru precedentele valori ale lui x, anume
y (x1 ) = y1 , y (x2 ) = y2 , y (x3 ) = y3 , etc., însă diferă calitativ prin tipul polinomului care
interpolează funcția y. Astfel, cunoscută fiind valoarea funcției y de la punctul inițial x0
până la un punct oarecare xn și fiind necesară aflarea funcției în punctul xn+1 , ecuația
(1.1) se rescrie de obicei sub forma

ˆ xn+1
yn+1 = yn + dxf (x, y(x)) , (1.2)
xn

determinarea lui yn+1 depinzând acum de calculul integralei din partea dreaptă a ecua-
ției.
Primele metode folosite pentru calculul numeric al acestei integrale porneau de la
construirea unui polinom de interpolare F (x), ce coincide cu f atât în xn , cât și într-un
număr de puncte anterioare lui xn , anume xn−1 , xn−2 , etc., valoarea aproximativă a lui
yn+1 fiind dată de
ˆ xn+1
yn+1 ≈ yn + dxF (x; xn , xn−1 , xn−2 ...)
xn
= yn + h (βn fn + βn−1 fn−1 + βn−2 fn−2 ....) , (1.3)

unde βn , βn−1 , βn−2 , etc., sunt coeficienți numerici a căror valoare depinde de ordinul
polinomului de interpolare. Deoarece polinomul de interpolare folosit în calcul integralei
este determinat folosind informații cu privire la soluția ecuației în puncte precedente lui
xn aceste metode sunt denumite în mod curent metode multi-pas sau metode cu memo-
rie. Aceste metode apar în a doua jumătate a secolului al XIX-lea, fiind introduse într-o
carte celebră scrisă de Francis Bashforth și John Couch Adams asupra acțiunii capilare
(Bashforth și Adams (1883)) și sunt apoi rafinate de Forest Ray Moulton, un important
astronom american, în a cărui carte New methods in exterior ballistics (Moulton (1926))
sunt cuprinse o serie de noi metode numerice de determinare a traiectoriei proiectile-
lor, dezvoltate de autor în perioada Primului Război Mondial pentru Armata Statelor

11
1 Introducere

Unite ale Americii1 . Particularitatea acestor metode numerice este aceea că polinomul
de interpolare F este determinat folosind și punctul xn+1 , nu doar punctul xn și cele
precedente, valoarea aproximativă a lui yn+1 fiind

yn+1 ≈ yn + h (βn+1 fn+1 + βn fn + βn−1 fn−1 + βn−2 fn−2 ....) . (1.4)

Aceste metode sunt așadar implicite, valoarea numerice a lui yn+1 fiind obținută după
rezolvarea unei ecuații al cărei grad de neliniaritate depinde de dependența funcțională
a lui f de y. Capitolul 2 al acestui volum cuprinde o discuție detaliată asupra acestor
metode, fiind prezentată pe larg o serie largă de formule de calcul (atât implicite cât și
explicite).
Cea de-a doua mare categorie de metode folosite pentru determinarea numerică a
soluțiilor ecuației (1.1) este reprezentată de așa-numitele metode unipas, în cadrul cărora
polinomul de interpolare F care aproximează pe f este determinat fără a folosi puncte
anterioare lui xn . Aceste metode au fost introduse într-un articol devenit acum clasic al
lui Carl David Tolmé Runge în care este propusă următoare schemă pentru determinarea
lui yn+1
( )
h h
yn+1 = yn + hf xn + , yn + f (xn , yn ) . (1.5)
2 2
După cum se observă cu ușurință schema de calcul este una explicită, partea dreaptă
a ecuației precedente fiind calculată exclusiv pe baza valorii funcției y în punctul xn .
Aceste metode explicite au fost dezvoltate apoi de către Karl Heun, Martin Wilhelm
Kutta, Nystrom și Erwin Fehlberg, acesta din urmă propunând o serie de metode simple
de calcul al erorii locale. Un rol aparte în dezvoltarea domeniului l-a avut John C.
Butcher care a introdus conceptul metodelor implicite și a schematizat metodele prin așa-
numitul tabel Butcher. Cea mai simplă asemenea metodă este metoda Euler implicită

yn+1 = yn + hf (xn+1 , yn+1 ) (1.6)

ce necesită rezolvarea unei ecuații algebrice la fiecare pas de integrare. Capitolul 3 cu-
prinde o prezentare detaliata a metodelor unipas, atât cele explicite, cât și cele implicite.

1.1 Limbaje de programare și sisteme de calcul

Dezvoltarea metodelor numerice a avut loc în strânsă legătură cu dezvoltarea limbajelor
de programare. Astfel, primul limbaj de programare care a avut librarii importante
1
Implicarea matematicienilor americani în Al Doilea Război Mondial este magistral surprinsă într-o
recenzie (nu lipsită de critici) la cartea lui Moulton: „Toward the close of the World War, and
even during the period of hostilities, every person who made any pretense of being interested in
contributions to mathematics longed to secure a book on ballistics from which he could compute a
trajectory; not by the old antiquated methods which, he had heard, had been entirely discarded, but
by the methods that had been devised by American mathematicians who had patriotically devoted
their talents to the solution of these important problems” (a se vedea Moulton on exterior ballistics,
J. E. Rowe, Bulletin of the American Mathematical Society 34, 229 (1928)).

12
 1.1 Limbaje de programare și sisteme de calcul

dedicate rezolvării numerice a ecuațiilor diferențiale ordinare a fost FORTRAN, al cărui

nume înseamnă FORmula TRANslation (System). Limbajul a fost creat în 1957 de că-
tre John Warner Backus, iar succesul metodelor numerice implementate în FORTRAN
este cel mai bine ilustrat de faptul că NAG Toolbox for MATLAB, unanim recunoscută
ca cea mai mare și mai vastă librarie de rutine numerice pentru MATLAB, utilieaza
încă (în 2014!) implementari Fortran. Prima varianta a FORTRAN cuprindea 32 de
instrucțiuni și era destul de rudimentară după standardele actuale, funcționând doar pe
bază de cartele perforate. Versiunile FORTRAN II și FORTRAN III apar ambele în
1958 însă nu aduc mari modificari. Prima versiune care schimba tiparul celor anterioa-
re este FORTRAN IV, aparută în 1962, în care sunt eliminate funcțiile dependente de
arhitectura mașinii de calcul. Aparut in 1966, FORTRAN 66 schimbă radical tradiția
versiunilor precedente și devine standardul în industria de profil grație numeroaselor
standardizări. Versiunea FORTRAN 77, aparută în 1978, repară o parte din minusurile
versiunii precedente, prin completarea multora din instrucțiunile existente. Sunt intro-
duse, de pilda, brațele else si else if pentru instrucțiunea if, inexistente până atunci, se
flexibilizează buclele de tip do, se pot defini constante, etc. Versiunile urmatoare apar
sub numele Fortran, cele mai importante fiind Fortran 95 și Fortran 2003.
Apariția FORTRAN a fost urmată de apariția BASIC (Begginers All purpose Symbolic
Instruction Code) în 1964, limbajul fiind dezvoltat în 1964 de către John G. Kemeny
și Thomas E. Kurtz pentru studenții de la Darmouth College, SUA. Țelul acestora era
acela de a oferi tuturor studenților, nu doar celor care studiau fizica și matematica, un
limbaj prietenos, ușor de utilizat.
Dezvoltat în perioada 1969-1973 la AT&T Bell Labs de către Dennis Ritchie, limbajul
de programare C nu poate fi desprins de dezvoltarea sistemului de operare UNIX, al carui
kernel (nucleu, motor) a fost scris (după câteva încercări nereușite) în acest limbaj.
Codurile numerice dezvoltate în C sunt foarte numeroase iar varietatea algorimtilor
disponibili e cel mai bine ilustrată de o carte devenită acum clasică: Numerical recipes
in C. The art of scientific computing, W. H. Press, S. A. Teukolsky, W. T. Vetterling,
B. P. Flannery, Cambridge University Press, prima ediție 1988. Bjarne Stroustrup a
modificat C prin introducerea noțiunii de clase, de funcții virtuale, prin posibilitatea
de suprascriere a operatorilor, de mostenire multiplă, șamd, dând astfel naștere C++.
Apărut pentru prima oară în 1983 limbajul a fost folosit pentru crearea suitei Microsoft
Office, a Acrobat Reader și a browserului Firefox.
Python este un limbaj de programare creat în 1989 de Guido van Rossum folosit
atat pentru programarea aplicațiilor web (de către companii ca Google sau Yahoo!),
cât și pentru aplicații științifice. Python suportă mai multe paradigme de programare,
permițând programare imperativă, funcțională, procedurală sau object-oriented.
Menționăm în cele din urma Java, un limbaj de programare orientat pe obiecte, pu-
ternic tipizat, conceput de către James Gosling (cu contribuții ulterioare din partea lui
Mike Sheridan și Patrick Naughton) la începutul anilor ’90. Aplicațiile Java sunt compi-
late astfel încât să poată rula pe orice mașină virtuala Java, independent de arhitectura
hardware a mașinii fizice.
Istoria metodelor numerice este, în bună măsură, și istoria limbajelor de programare,
iar cine trece în revistă literatura de specialitate va observa cum majoritatea volumelor

13
1 Introducere

dedicate calculului numeric în anii 80 ofereau implementari FORTRAN (și conțineau

adesea un mic breviar de funcții FORTRAN alături de multe seturi de note explicative),
în timp ce volumele scrise în anii 90 arată deja o preferință pentru C (și, oarecum mai rar,
C++). Volumele ultimului deceniu arată deja creșterea interesului către noile limbaje,
cum sunt Python și Java, pentru care există deja librării extrem de performante. Tomuri
masive cum sunt A numerical library in Java for scientists and engineers, H. T. Lau,
Chapman & Hall/CRC (2004), A primer on scientific programming with Python, H. P.
Langtangen, Springer (2010), Numerical methods, algorithms and tools in C#, W. Dos
Passos, CRC Press (2010), și multe, multe altele reflectă extrem clar tendința de a muta
calculul științific către noile limbaje de programare.
Pe lângă limbajele de programare prezentate anterior, un rol aparte l-au avut așa-
numitele API-uri, Application Programming Interfaces, care permit împărțirea memoriei
pe mai multe procesoare, deschizând astfel calea spre calcul paralel și distribuit. Cel
mai cunoscut API este OpenMP (Open Multi-Processing) care lucrează cu C, C++,
și Fortran, și rulează pe marea majoritate a arhitecturilor de procesoare și sisteme de
operare. Pe un palier complementar se situează MPI, un protocol de comunicare extrem
de util în programarea paralelă și calcul de înaltă performanță.
Paragrafele precedente surprind doar o parte din dinamica calculului științific, cel mai
important aspect încă nediscutat fiind evoluția sistemelor de calcul. O discuție detaliată
depășește cu mult cadrul acestor note de curs2 , dar putem observa, de pildă, impactul
fenomenal pe care dezvoltarea așa-numitelor Graphical Processing Units - GPUs l-a avut
asupra calcului științific, oferind unei categorii extrem de largi de utilizatori posibilitatea
de a experimenta calculul științific puternic paralelizat. Unitățile de procesare grafică
sunt, de fapt, cele care au transformat calculul științific de înaltă performanță dintr-un
domeniu de interes mai degraba teoretic/academic, într-un subiect de maxim interest
pentru comunitatea științifică. O serie întreagă de metode extrem de cronofage (legate
de tranformările integrale, algebră liniară, procesarea de imagini, șamd) au cunoscut
dezvoltări uimitoare, fiind acum implementate în mod curent pe platforme de calcul
științific gen MATLAB și MATHEMATICA.

1.2 Bibliografie
1. Bashforth, F., Adams, J. C., An attempt to test the theories of capillary action by
comparing the theoretical and measured forms of drops of fluids, with an explanation
of the method of integration employed in constructing the tables which give the
theoretical forms of such drops, Cambridge at the University Press (1883)

2. Moulton, F. R., New methods in exterior ballistics, Chicago University Press (1926)

3. Goldstine, H. H, A history of numerical analysis from the 16th through the 19th
century, Springer-Verlag (1977)
2
Recomandăm cititorilor interesați remarcabilul volum scris de D. R. Hartree asupra Calculating instru-
ments and machines în care descrie începuturile mașinilor digitale (calculatoarelor) și implementările
hardware ale câtorva metode numerice (astăzi) elementare.

14
 1.2 Bibliografie

4. Bergin, T. J., și Gibson, R. G., History of programming languages, 3 volume,

Addison-Wesley (1996)

5. Brezinski, C., și Wuytack, L., editori, Numerical analysis: Historical developments

in the 20th century, Elsevier (2001)

15
2 Metode multipas
THE NUMERICAL SOLUTION OF DIFFERENTIAL EQUATIONS
177. Theory of the Method The best method of integrating
differential equations numerically is one devised by J. C. Adams; it is
applicable to equations of any order, but for simplicity we shall
describe its application to equations of the first order. [...]
178. Bibliographical Note Of the other methods which have
been proposed for integrating differential equations, the best known is
that of Runge, Math. Ann. 46 (1895), p. 167, improved and
extended by Kutta, Zeits. f. Math. u. Phys. 46 (1901), p. 435.
(The calculus of observations. A treatise on numerical mathematics,
E. T. Whittaker și G. Robinson, 50 Old Bailey, London (1924))

Considerate multă vreme a fi începutul și sfârșitul metodelor numerice pentru ecuații

diferențiale ordinare (așa cum foarte bine reiese din citatul ales ca moto al capitolului),
metodele multipas oferă o imagine extrem de intuitivă asupra modului în care putem
construi soluții aproximative pentru o ecuație diferențială dată. Să presupunem, de
pildă, că pentru ecuația (1.1) am putut construi o soluție aproximativă până la un punct
xn și că am pastrat, de asemenea, toate informațiile aferente punctului xn−1 . Știind,
așadar, punctele (xn , fn ) și (xn−1 , fn−1 ) putem construi fără probleme un polinom de
interpolare pentru f care ia forma

x − xn x − xn−1
F (x) = fn−1 + fn . (2.1)
xn−1 − xn xn − xn−1
Cu el putem determina imediat soluția aproximativă a ecuației diferențiale în punctul
xn+1 , calculând o integrală extrem de simplă, i.e.,
ˆ xn+1
yn+1 ≈ yn + dxF (x),
xn
ˆ xn+1 ( )
x − xn x − xn−1
= yn + dx fn−1 + fn ,
xn xn−1 − xn xn − xn−1
( )
3 1
= yn + h fn − fn−1 . (2.2)
2 2
Pentru a spori acuratețea rezultatului anterior, putem crește ordinul polinomului de in-
terpolare (în calculul precedent, doar o dreaptă), prin adăugarea punctului (xn−2 , fn−2 ).
Noul polinom de interpolare (acum de ordinul doi) este dat de

17
2 Metode multipas

x − xn x − xn−1 x − xn x − xn−2
F (x) = fn−2 + fn−1
xn−2 − xn xn−2 − xn−1 xn−1 − xn xn−1 − xn−2
x − xn−2 x − xn−1
+ fn , (2.3)
xn − xn−2 xn − xn−1
iar noua soluție aproximativă este dată de

ˆ xn+1
yn+1 ≈ yn + dxF (x),
xn
ˆ xn+1 [
1 1
= yn + 2 dx (x − xn ) (x − xn−1 ) − (x − xn ) (x − xn−2 )
h xn 2
]
1
+ (x − xn−2 ) (x − xn−1 ) ,
2
( )
23 16 5
= yn + h fn − fn−1 + fn−2 (2.4)
12 12 12
Soluția aproximativă obținută mai sus poate fi determinată mai simplu, fără a aproxima
f (x), scriind yn+1 ca o combinație liniară de fn , fn−1 și fn−2 , și dezvoltând în serie în
jurul lui xn , anume

yn+1 ≈ yn + hβ1 fn + hβ2 fn−1 + hβ3 fn−2 , (2.5)

yn+1 ≈ yn + hβ1 yn′ + hβ2 yn−1
′
+ hβ3 yn−2,′
(2.6)
h2 h3 ( )
yn + hyn′ + yn′′ + yn′′′ + O h4 ≈ yn + hβ1 yn′ + hβ2 yn′ − h2 β2 yn′′
2 6
h3
+ β2 yn′′′ + hβ3 yn′ − 2h2 β3 yn′′
2 ( )
+2h3 β3 yn′′′ + O h4 . (2.7)

Egalând derivatele de ordinul 1, 2 și 3 ale lui yn obținem sistemul

β1 + β2 + β3 = 1, (2.8)
1
−β2 − 2β3 = , (2.9)
2
β2 1
+2β3 = , (2.10)
2 6
a cărui soluție este dată de
23 16 5
β1 = , β2 = − , β3 = , (2.11)
12 12 12
așadar tocmai formula (2.4). Metode similare de determinare a coeficienților sunt folosite
și pentru metodele unipas de tip Runge-Kutta de ordin scăzut. Metoda nu funcționează

18
însă pentru metode (fie ele uni- sau multipas) de ordin înalt căci ecuațiile care rezultă
necesită condiții suplimentare pentru determinarea unică a coeficienților β.
Calculul polinomului de interpolare al lui f este relativ ușor atunci când sunt două
sau trei puncte, însă devine din ce în ce mai complicat odată cu creșterea numărului de
puncte. O soluție simplă și ingenioasă de calcul se bazează pe diferențele divizate din
tabelul de mai jos.

xn fn

fn − fn−1
xn−1 fn−1
h
fn−1 − fn−2 fn − 2fn−1 + fn−2
xn−2 fn−2
h 2h2
fn−2 − fn−3 fn−1 − 2fn−2 + fn−3 fn − 3fn−1 + 3fn−2 − fn−3
xn−3 fn−3
h 2h2 6h3

Cu ajutorul lor putem scrie polinomul de interpolare Newton

F0 (x) = fn , (2.12)
fn − fn−1
F1 (x) = fn + (x − xn ), (2.13)
h
fn − fn−1
F2 (x) = fn + (x − xn )
h
fn − 2fn−1 + fn−2
+ (x − xn ) (x − xn−1 ), (2.14)
2h2
fn − fn−1 fn − 2fn−1 + fn−2
F3 (x) = fn + (x − xn ) + (x − xn ) (x − xn−1 )
h 2h2
fn − 3fn−1 + 3fn−2 − fn−3
+ (x − xn ) (x − xn−1 ) (x − xn−2 ). (2.15)
6h3

Să observăm că polinoamele de interpolare folosite anterior (a se vedea ecuațiile (2.1)

și (2.3)) erau de tip Lagrange, adăugarea unui nou punct necesitând calculul de la zero
al polinomului de interpolare și nu doar simpla adăugare a unui noi termen, așa cum
se întâmplă în cazul polinoamelor de interpolare de tip Newton. Integralele necesare
calculului soluțiilor aproximative se pot scrie sub forma

19
2 Metode multipas

ˆ xn+1
dxF0 (x) = hfn , (2.16)
ˆxnxn+1
1
dxF1 (x) = hfn + (fn − fn−1 ) ,
xn 2
( )
3 1
= h fn − fn−1 , (2.17)
2 2
ˆ xn+1 ( )
3 1 5
dxF2 (x) = h fn − fn−1 + (fn − 2fn−1 + fn−2 ) ,
xn 2 2 12
( )
23 16 5
= h fn − fn−1 + fn−2 , (2.18)
12 12 12
ˆ xn+1 ( )
23 16 5 3
dxF3 (x) = h fn − fn−1 + fn−2 + (fn − 3fn−1 + 3fn−2 − fn−3 ) ,
xn 12 12 12 8
( )
55 59 37 9
= h fn − fn−1 + fn−2 − fn−3 . (2.19)
24 24 24 24

Am subliniat cu roșu ultimul termen al polinomului de interpolare și al soluției apro-

ximative pentru a sublinia ca am obținut o metodă iterativă cu ajutorul căreia putem
construi o soluție de ordin n + 1 pornind de la o soluție (cunoscută) de ordin n. Mai
clar, modul acesta de a construi metodele multipas subliniază ideea ca o metode de ordin
n + 1 corecteaza (în sensul ameliorării preciziei) o metodă de ordin n, iar termenul de
corecție este cel marcat cu roșu. Practic, dacă știm, de exemplu, aproximarea integralei

ˆ xn+1
dxf (x)
xn

pentru un polinom de grad doi, adica f (x) ≈ F1 (x), obținem aproximarea integralei
pentru un polinom de gradul trei calculând doar corecția, anume

ˆ xn+1 ˆ xn+1
dxF2 (x) = dxF1 (x)
xn xn
ˆ
xn+1
fn − 2fn−1 + fn−2
+ dx (x − xn ) (x − xn−1 ), (2.20)
xn 2h2
( )
3 1 5
= h fn − fn−1 + (fn − 2fn−1 + fn−2 ), (2.21)
2 2 12
( )
23 16 5
= h fn − fn−1 + fn−2 . (2.22)
12 12 12

La modul general putem scrie relația de recurență

20
 2.1 Metode Adams-Bashforth

 
ˆ xn+1 ˆ xn+1 ∑
s  s 
 
dxFs (x) = dxFs−1 (x) + γs (−1)j fn−j   , (2.23)
xn xn  
j=0
j

unde

s 1 2 3 4 5 6 7 8
.
1 5 3 251 95 19087 5257 1070017
γs
2 12 8 720 288 60480 17280 3628800

Pentru o discuție detaliată asupra formulei (2.23) recomandăm Ref. [3], capitolul III.1
Classical linear multistep formulas.

2.1 Metode Adams-Bashforth

After some trials a satisfactory micrometrical instrument was
constructed for the measurement of the forms of drops of fluid, but
my attempts to calculate their forms as surfaces of double curvature
failed entirely, and my undertaking must have ended here, if I had
depended upon my own resources. But at this point Professor J. C.
Adams furnished me with a perfectly satisfactory method of
calculating by quadratures the exact theoretical forms of drops of
fluids from the Differential Equation of Laplace, an account of which
he has now had the kindness to prepare for publication.
(An attempt to test the theories of capillary action by comparing the
theoretical and measured forms of drops of fluids, with an explanation
of the method of integration employed in constructing the tables
which give the theoretical forms of such drops, F. Bashforth și J. C.
Adams, Cambridge at the University Press (1883))

Metodele derivate anterior sunt cunoscute în literatura de specialitate drep metode

Adams-Bashforth, de la numele autorilor primului articol in care ele sunt introduse (a se
vedea Ref. [1]). Istoria articolului (din care am extras moto-ul acestui capitol) ne arată
că metoda a fost dezvoltată din nevoia practică de a descrie geometria unei picături
aflate pe o suprafață solidă. Citatul ales arată că paternitatea metodei îi aparține doar
lui Adams, Bashforth fiind citat mai mult datorită conjuncturii.
Folosind formula iterativa (2.23) am calculat mai jos formule Adams-Bashforth de la
ordinu l, formula care se identifica cu metoda Euler, pana la ordinul 8.

21
2 Metode multipas

yn+1 = yn + hfn (2.24)

( )
3 1
yn+1 = yn + h fn − fn−1 (2.25)
2 2
( )
23 16 5
yn+1 = yn + h fn − fn−1 + fn−2 (2.26)
12 12 12
( )
55 59 37 9
yn+1 = yn + h fn − fn−1 + fn−2 − fn−3 (2.27)
24 24 24 24
( )
1901 1387 109 637 251
yn+1 = yn + h fn − fn−1 + fn−2 − fn−3 + fn−4 (2.28)
720 360 30 360 720
( )
4277 2641 4991 3649 959 95
yn+1 = yn + h fn − fn−1 + fn−2 − fn−3 + fn−4 − (2.29)
fn−5
1440 480 720 720 480 288
(
198721 18637 235183 10754 135713
yn+1 = yn + h fn − fn−1 + fn−2 − fn−3 + fn−4
60480 2520 20160 945 20160
)
5603 19087
− fn−5 + fn−6 (2.30)
2520 60480
(
16083 1152169 242653 296053 2102243
yn+1 = yn + h fn − fn−1 + fn−2 − fn−3 + fn−4
4480 120960 13440 13440 120960
)
115747 32863 5257
− fn−5 + fn−6 − fn−7 (2.31)
13440 13440 17280
Determinarea metodelor Adams-Bashforth prezentată aici este una riguroasă, lipsită
de artificii, dar există și derivări mai simple, lipsite însă de transparență (a se vedea în
acest sens Exercițiul 1). Din punct de vedere practic pentru implementarea numerică a
metodei este nevoie doar de coeficienții β. Subliniem, însă, o particularitate a acestor
metode: ele necesită de obicei condiții inițiale mult mai numeroase decât este necesar
pentru a avea o problemă numerică corect definită din punct de vedere matematic. Mai
clar, la momentul t = 0 cunoaștem doar f0 ceea ce ne împiedică să folosim metodele
(2.25)-(2.31) din lipsa de condiții inițiale, singura opțiune fiind metoda de ordin 1 definită
de (2.24). La pasul următor vom cunoaște atât f0 cât și f1 , ceea ce ne permite să
folosim o metodă de ordin doi, iar la pasul urmator putem folosi o metodă de ordin trei,
căci cunoaștem f2 , șamd. Inițializarea metodei este principala problemă a metodelor
multipas și din această cauză o bună parte din implementările existente (inclusiv cea
din MATLAB) au la baza așa-numitele metode de ordin variabil ce permit schimbarea
ordinului metodei de la un pas la altul.
Menționăm, de asemenea, că formulele de interpolare (2.12)-(2.15) și metodele definite
de ecuațiile (2.24)-(2.31) pornesc de la presupunerea că punctele xj sunt distribuite
echidistant, i.e., x2 − x1 = x3 − x2 = x4 − x3 = .... Practica integrării numerice arată
însă că această distribuție a punctelor este deseori ineficientă, căci sunt regiuni în care
dinamica sistemului e minimală (e.g., oscilații de frecvență foarte mică) și se poate folosi
un pas de integrare relativ mare, în timp ce în alte regiuni dinamica este una complicată
(e.g., oscilații de frecvență foarte mare) și este necesar un pas de integrare relativ mic.

22
 2.2 Metode Adams-Moulton

2.2 Metode Adams-Moulton

The statement has often been made that I, being familiar with the
mathematical methods used by astronomers, naturally introduced this
method of computing trajectories. There seems to be a widespread
belief that astronomers use it in determining the orbits of comets, in
the lunar theory, and in perturbation theories. Nothing could be more
remote from the facts. Astronomers as a class are as innocent of any
acquaintanceship with this method of solving differential equations as
most mathematicians were before the war. Probably not one out of
twenty of them ever heard that such a method exists; and though the
underlying principles of it have been used by a few mathematicians in
computing special perturbations, the general procedure they employed
has been so different in form from that used in ballistics that one
would not easily detect the connection. So far as I am concerned, I did
not learn the method from astronomy, but from a general survey of
the various processes that are known for solving differential equations.
(New methods in exterior ballistics, F. R. Moulton, The American
Mathematical Monthly 35, 246 (1928))

În determinarea metodelor Adams-Bashforth prezentate anterior am folosit doar apro-

ximări explicite ale integralelor (2.16)-(2.19) care au avantul că generează, la rândul lor,
metode explicite, i.e., yn+1 poate fi determinat imediat pe baza valorii fn și a unei se-
rii fn−1 , fn−2 , etc. O alta abordare des întâlnită în literatura de specialitate pornește
de la aproximări implicite, i.e., integralele (2.16)-(2.19) sunt aproximate folosind inclu-
siv fn+1 . Calculele fiind similare cu cele prezentate anterior nu le mai repetăm aici ci
prezentăm mai jos varianta implicită a ecuațiilor (2.24)-(2.31). Aceste metode sunt cu-
noscute în literatură (a se vedea Ref. [2-4]) sub titulatura de metode Adams-Moulton si
reprezintă o îmbunătățire semnificativă. Diferența dintre metodele Adams-Bashforth și
metodele Adams-Moulton nu ține, așa cum ne-am putea imagina, de un plus de acura-
tețe, ci de un plus de stabilitate. În cuvinte simple, stabilitatea se traduce prin nevoia
unui pas de integrare foarte mic pentru a face schema numerică convergentă, iar me-
todele introduse F. R. Moulton1 au avantajul de a face schemele numerice convergente
pentru pași de integrare relativi mari, pași pentru care metodele Adams-Bashforth sunt
1
Impactul contribuției lui F. R. Moulton la studiul dinamicii proiectilelor reiese cel mai bine din recen-
ziile vremii, așa cum se poate vedea din urmatorul paragraf extras din Moulton on exterior ballistics,
J. E. Rowe, Bulletin of the American Mathematical Society 34, 229 (1928): „[...] The fundamental
differential equations of motion of a projectile in flight have been known since the days of Euler. Up to
the World War the progress in ballistics had consisted in devising approximate algebraic expressions
whereby to wrench these equations into soluble form. [...] These approximations had been improved
and improved, but still could not keep pace with the development of modern artillery. So finally
Professor F. R. Moulton of the University of Chicago, while serving as a Major in the U. S. Army
during the World War, cut the Gordian knot by going back to the Eulerian equations and solving
them in their original and exact form by numerical integration. Thus he laid the cornerstone for an
entirely new science of ballistics.”

23
2 Metode multipas

puternic divergente. Ca o observație de ordin general metodele implicite sunt mult mai
stabile numeric decât cele explicite în implementarea lor necesită un efort computațional
sporit.
Înainte de a prezenta formulele propriu-zise, subliniem că acestea nu au legătură cu
dinamica corpurilor cerești, acesta fiind și sensul moto-ului acestui capitol. Parcursul
biografic al lui Moulton este în mod vădit unul exotic, însă formulele care urmează se
datorează aproape exclusiv studiilor balistice făcute în perioada primului război mondia
(a se vedea în acest sens A. Gluchoff, Artillerymen and mathematicians: Forest Ray
Moulton and changes in American exterior ballistics, 1885–1934, Historia Mathematica
38, 506 (2011)).

yn+1 = yn + hfn+1 (2.32)

( )
1 1
yn+1 = yn + h fn+1 + fn (2.33)
2 2
( )
5 8 1
yn+1 = yn + h fn+1 + fn − fn−1 (2.34)
12 12 12
( )
9 19 5 1
yn+1 = yn + h fn+1 + fn − fn−1 + fn−2 (2.35)
24 24 24 24
( )
251 323 11 53 19
yn+1 = yn + h fn+1 + fn − fn−1 + fn−2 − fn−3 (2.36)
720 360 30 360 720
( )
95 1427 133 241 173 3
yn+1 = yn + h fn+1 + fn − fn−1 + fn−2 − fn−3 + (2.37)
fn−4
288 1440 240 720 1440 160
(
19087 2713 15487 586 6737
yn+1 = yn + h fn+1 + fn − fn−1 + fn−2 − fn−3
60480 2520 20160 945 20160
)
263 863
+ fn−4 − fn−5 (2.38)
2520 60480
(
5257 139849 4511 123133 88547
yn+1 = yn + h fn+1 + fn − fn−1 + fn−2 − fn−3
17280 120960 4480 120960 120960
)
1537 11351 275
+ fn−4 − fn−5 + fn−6 (2.39)
4480 120960 24192

Pentru a întelege de ce metodele implicite necesită un efort computațional sporit vom

scrie explicit metoda Adams-Moulton de ordin 2 definită în ecuația (2.33)

( )
1 1
yn+1 = yn + h f (xn+1 , yn+1 ) + f (xn , yn )
2 2

care definește acum o ecuație algebrică implicită în yn+1 ce trebuie rezolvată pentru
fiecare xn ! Daca dependența lui f de y este, de pildă, una polinomială, soluția numerică a
ecuației precedente nu este foarte dificilă, însă dacă dependența este una mai complicată,
sa zicem de tip sin (yn ), ecuația care rezultă este una transcendentă iar soluția numerică

24
 2.3 Implementări software disponibile

poate fi puternic netrivială. Ca o paranteză tematică menționăm următoarele: soluția

numerică a unei ecuații diferențiale nu este, în mod standard, un subiect ce ține de
calculul de înaltă performanță, în engleză High-Performance Computing (HPC), pentru
ca nu lasă loc paralelizării. Metodele implicite însă sunt potrivite pentru calculul de
înaltă performanță deoarece atunci când vor fi folosite pentru sisteme mari de ecuații
diferențiale ordinare puternic neliniare, să zicem, de ordinul miilor și zecilor de mii,
soluția ecuației algebrice ce trebuie să fie rezolvată la fiecare pas e o problemă în sine,
una care se pretează la calcul paralel atât pe sisteme cu mai multe core-uri cât și pe
sisteme de tip Graphical Processing Unit (GPU).
Mai menționăm în final un aspect deseori neglijat, dar extrem de important: prin creș-
terea exagerată a ordinului de marime al polinomului de interpolare apare așa-numitul
fenomen Runge, i.e., oscilații puternice (complet nefizice) la capetele domeniului de in-
terpolare atunci când sunt folosite interpolări polinomiale de ordin ridicat. Fenomenul
este puternic înrudit cu fenomenul Gibbs care se manifestă prin apariția unor oscilații
puternice în spectrul Fourier al unei unde neliniare periodice formată din funcții treaptă.
În mod deloc surprinzător metodele implicite atenuează oscilațiile puternice înrudite cu
fenomenul Gibbs, dar concluzia generală rămâne validă: metodele de ordin înalt nu sunt
în mod obligatoriu mai precise numeric.

2.3 Implementări software disponibile

Integratorii numerici de tip Adams-Bashforth și Adams-Moulton sunt implementați în
MATLAB, dar nu și în MATHEMATICA. În MATLAB integratorul de tip Adams este
ode113, care are la bază o metoda Adams-Bashforth-Moulton de ordin variabil. Inte-
gratorul acceptă un set standardizat de parametri care includ eroarea relativă (RelTol)
și eroarea absolută (AbsTol), pasul inițial de integrare (InitialStep), numărul maxim de
pași de integrare (MaxStep) și, foarte important, matricea Jacobi a sistemului de ecuații
diferențiale (Jacobian și JPattern). Această din urmă opțiune este extrem de utilă în
cazul ecuațiilor rigide și sporește considerabil eficiența integratorilor numerici.
În afara implementării din MATLAB cele mai cunoscute coduri dedicate metodelor
de tip Adams sunt acelea dezvoltate de Numerical Algorithms Group - NAG. Cea mai
importantă rutină a NAG este D02GFF care implementează o metodă Adams de or-
din variabil, de pas variabil, dedicată cu precăderelor sistemelor nerigide. Metoda este
disponibila în Fortran, varianta în C purtând eticheta D02QFC.

2.4 Exerciții

Exercițiul 1: (Adaptare după An approach to the Adams-Bashforth formula, M. Ritter, The Mathe-
matical Gazette 63, 48 (1979)) Considerând că F (x) este un polinom de gradul n iar F (x + h) se poate
scrie ca
( )
F (x + h) − F (x) = h B0 F ′ (x) + B1 F ′ (x − h) + ... + BN −1 F ′ (x − (n − 1)h) , (E.2.1)
să se arate că dezvoltarea în serie Taylor a funției F și a derivatelor sale și conduce la urmatorul sistem
de ecuații

25
2 Metode multipas





 B0 + B1 + ... + Bn = 1



     



 1

  1 2 3 ... ... n−1  B1   − 

     2 

     

     1 

  12 22 32 ... ... (n − 1)2  B2   

     3 

     
    
 13 23 33 (n − 1)3    −
1  (E.2.2)

  ... ...  B3   

    =  4 

     

     

  ... ... ... ... ... ...  ...   ... 

     

     

     

     

  ... ... ... ... ... ...  ...   ... 

     

     

 (−1)n−1
1n−1 2n−1 3n−1 ... ... (n − 1)n−1 Bn−1 −
n
Să se determine soluția numerică a ecuațiilor din sistemul (E.2.2) pentru n = 1, 2, 3, 4 și 5, și să se arate
că ecuația (E.2.1) descrie metodele Adams-Bashforth de ordinul 1, 2, 3, 4 și 5.

2.5 Bibliografie
1. Bashforth, F., și Adams, J. C., An attempt to test the theories of capillary action by comparing the
theoretical and measured forms of drops of fluids, with an explanation of the method of integration
employed in constructing the tables which give the theoretical forms of such drops, Cambridge at
the University Press (1883)
2. Butcher, J. C., Numerical methods for ordinary differential equations, John Wiley & Sons (2008)
3. Hairer, E., Nørsett, S.P., Wanner, G., Solving ordinary differential equations I. Nonstiff Problems,
Springer (2008)
4. Dormand, J.R., Numerical methods for differential equations. A computational approach, CRC
Press (1996)
5. Shampine, L.F., Gladwell, I., și Thompson, S., Solving ODEs with Matlab, Cambridge University
Press (2003)

26
3 Metode unipas
Metodele prezentate în capitolul precedent se bazează pe determinarea soluției numerice
folosind un polinom de interpolare global cu ajutorul căruia determinăm soluția unei
ecuații diferențiale la xn utilizând informațiile de la xn−1 , xn−2 , xn−3 , etc. Deoarece
interpolarea este globală, ar putea părea firesc ca aceste metode să fie prezentate după
metodele de Runge-Kutta care determină soluția la xn folosind exclusiv informația la
xn−1 . Numeroasele volume care preferă această abordare au în vedere faptul că meto-
dele unipas sunt metode fără memorie, în timp ce metodele de tip Adams sunt metode
cu memorie care necesită, de obicei, atenție sporită la implementarea numerică. Pre-
zentând metodele multipas înaintea celor unipas vrem să subliniem diferențele calitative
dintre polinoamele de interpolare. Astfel, chiar dacă metodele multipas se bazează pe
o interpolare globală, aceasta este una foarte simplă, polinomială, în timp ce interpola-
rea locala a metodelor unipas este, cum vom vedea, mult mai complicată, aspect care
motivează prezentarea acestora după metodele multipas.
În cele ce urmează vom prezenta metoda Euler și diferitele clase de metode Euler.
Discuția asupra metodelor Runge-Kutta pornește de la metodele Runge-Kutta explicite,
folosite în general pentru majoritatea problemelor care nu au probleme deosebite de
stabilitate, și continua apoi cu metodele implicite. Acestea sunt folosite cu precădere
pentru așa-numitele probleme rigide, i.e., pentru determinarea numerică a soluțiilor
ecuațiilor diferențiale ale căror soluții prezintă atât oscilații de frecvență foarte mică, cât
și oscilații de frecvență foarte mare. Capitolul se încheie cu o discuție detaliată asupra
metodelor Runge-Kutta împerecheate cu ajutorul cărora obținem controlul optim al
erorii.
Lăsând la o parte aspectele tehnice mesajul principal al acestui capitol este că metodele
de tip Runge-Kutta cunosc o diversitate mult mai mare decât metodele de tip Adams,
fiind utilă pentru o clasă mult mai largă de probleme, numele unui articol classic al lui
J. B. Rosser din 1967 fiind chiar A Runge-Kutta for all seasons (SIAM Review 9, 417
(1967)). Acest aspect este reflectat și în literatura de profil, extrem de amplă în ceea
ce privește metodele. Opțiunile bibliografice privitoare la metodele Runge-Kutta sunt
extrem de variate, unii autori optând pentru vaste articole de sinteză, în timp ce alții
optează pentru cărți de referință. În ceea ce ne privește am optat pentru autori clasici,
cu contribuții recunoscute internațional la dezvoltarea domeniului, ale căror cărți sunt
unanim recunoscute. Numele cel mai des întâlnit în paginele ce urmează este al Prof.
John C. Butcher, Profesor de Matematică la University of Auckland, care introduce în
anii ’60 metodele Runge-Kutta implicite și tabelul care îi poartă numele cu ajutorul
căruia se pot transmite într-o formă concisă toate informațiile privitoare la o anumită
metodă. Istoria metodelor Runge-Kutta scrisă de Prof. Butcher, Ref. [6] a acestui
capitol, a devenit un text de referință, în timp ce tratatul Numerical methods for ordinary

27
3 Metode unipas

differential equations (a se vedea Ref. [7] a acestui capitol) a devenit o carte clasică.
Munca Prof. Butcher a fost continuată de E. Hairer, S. P. Nørsett și G. Wanner ale
căror două volume dedicate Solving ordinary differential equations (a se vedea Ref. [8]
și Ref. [9]) sunt, pentru foarte multe probleme, referința ultimă în domeniu.

3.1 Metoda Euler

Cea mai simple dintre metodele folosite pentru rezolvarea numerică a unei ecuații dife-
rențiale ordinare este metoda Euler, pe care am întâlnit-o deja în capitolul precedent (a
se vedea ecuația (2.24)). Ideea din spatele metodei este extrem de simpla: soluția ecua-
ției diferențiale este aproximată de o linie polinomială formată din tangentele la soluția
locală a ecuației. Metoda a fost introdusă de Leonhard Euler în volumul Institutionum
calculi integralis (publicat în perioada 1768–1770). Astfel, dată fiind ecuația diferențială

y ′ = f (x, y) y (x0 ) = y0 (3.1)

soluția este construită ca

y (x0 + h) = y (x0 ) + hf (x0 , y0 ) , (3.2)

y (x0 + 2h) = y (x0 + h) + hf (x0 + h, y0 (x + h)) , (3.3)
y (x0 + 3h) = y (x0 + 2h) + hf (x0 + 2h, y0 (x + 2h)) , (3.4)
... (3.5)
y (x0 + N h) = y (x0 + (N − 1)h) + hf (x0 + (N − 1)h, y0 (x + (N − 1)h)) . (3.6)

Efortul computațional al metodei Euler este minim, dar precizia este de asemenea mi-
nima.

3.2 Metode Runge-Kutta

Runge-Kutta used to be what you used when (i) you didn’t know any
better, or (ii) you had an intransigent problem where Bulirsch-Stoer
was failing, or (iii) you had a trivial problem where computational
efficiency was of no concern. However, advances in Runge-Kutta
methods, particularly the development of higher-order methods, have
made Runge-Kutta competitive with the other methods in many
cases. Runge-Kutta succeeds virtually always.
(Numerical recipes The art of scientific computing, W. H. Press, S. A.
Teukolsky, W. T. Vetterling și B. P. Flannery, Cambridge University
Press (2007))

Prima extensie a metodei Euler este așa-numită metodă a punctului de mijloc, pentru
care, data fiind ecuația diferențială autonoma

28
 3.2 Metode Runge-Kutta

y ′ = f (x), y (x0 ) = y0 , (3.7)

cu soluția formală
ˆ x
y(x) = y0 + f (x), (3.8)
x0
aproximăm soluția numerică pe baza seriei

( )
h
y (x0 + h) ≈ y0 + hf x + , (3.9)
2
( )
h
y (x1 + h) ≈ y1 + hf x1 + , (3.10)
2
... (3.11)
( )
h
y (x) ≈ yN −1 + hf xN −1 + . (3.12)
2
Pornind de la această aproximare a soluției pentru sisteme dinamice autonome, vom
considera mai jos ecuația neautonoma

ẏ = f (x, y), y (x0 ) = y0 , (3.13)

a carei soluții numerice este construită după schema

( ( ))
h h
y (x0 + h) ≈ y0 + hf x0 + , y x0 + .
2 2
Această schemă se poate scrie sub forma

k1 = f (x0 , y0 ) , (3.14)
( )
h h
k2 = f x 0 + , y 0 + k1 , (3.15)
2 2
y1 = y0 + hk2 , (3.16)

care este, după vom vedea în paginile următoare, forma standard de scriere a unei metode
Runge-Kutta. Pentru a ne convinge că această metodă este de ordinul 2, vom scrie întâi
serie Taylor a lui y1 ca funcție de h
( )
h h
y1 = y0 + hf x0 + , y0 + f (x0 , y0 ) , (3.17)
2 2
h2
= y0 + hf (x0 , y0 ) + (fx + fy f ) (x0 , y0 )
2
h3 ( ) ( )
+ fxx + 2fxy f + fyy f 2 (x0 , y0 ) + O h4 , (3.18)
8

29
3 Metode unipas

c2 a21

c3 a31 a32

.. .. .. ..
. . . .

cs as1 as2 ... as,s−1

b1 b2 ... bs−1 bs

Tabel 3.1: Metode Runge-Kutta explicite. Tabel Butcher general.

iar apoi serie Taylor a lui y ca funcție de x, anume

h2
y (x0 + h) = y0 + hf (x0 , y0 ) + (fx + fy f ) (x0 , y0 )
2!
h3 ( ) ( )
+ fxx + 2fxy f + fyy f 2 + fy fx + fy2 f + O h4 . (3.19)
3!
( )
Din diferența y (x0 + h)−y1 = O h3 reiese clar că metoda precedentă este precisă până
la termeni de ordinul lui h3 , fiind așadar o metodă de ordinul 2.
Forma generală sub care este scrisă o metode Runge-Kutta este

k1 = f (x0 , y0 ) , (3.20a)
k2 = f (x0 + c2 h, y0 + ha21 k1 ) , (3.20b)
k3 = f (x0 + c3 h, y0 + h (a31 k1 + a32 k2 )) , (3.20c)
... (3.20d)
ks = f (x0 + cs h, y0 + h (as1 k1 + as2 k2 + ... + as,s−1 ks−1 )) , (3.20e)
y1 = y0 + h (b1 k1 + b2 k2 + ... + bs ks ) , (3.20f)

ecuațiile precedente putând fi comprimate într-un tabel de forma de mai jos, numit în
general tabel Butcher.

3.2.1 Metode Runge-Kutta explicite

Forma( generală
) a metodelor Runge-Kutta de ordin 2, i.e., cu erori locale de trunchiere
de O h , și cu doi pași, i.e., cu două valori k, este prezentată în tabelul Butcher (3.2).
3

În general vorbim de „metoda de ordin 2”, dar așa cum se observă în tabel avem o

30
 3.2 Metode Runge-Kutta

c2 c2

1 1
1−
2c2 2c2
Tabel 3.2: Metoda Runge-Kutta explicită (generală) cu s = p = 2 și c2 parametru liber.

întreagă clasă de metode, fiecare metodă fiind definită de valoarea lui c2 . De exemplu,
metoda definita de ecuațiile (3.14)-(3.16) se poate obține din tabelul Butcher (3.2) pentru
c2 = 0.5.
Parte din arta dezvoltării metodelor Runge-Kutta ține de calculul tabelelor Butcher.
Pentru metodele de ordin scăzut se poate aplica un raționament similar celui din jurul
ecuațiilor (3.17)-(3.19). Ideea generală este aceea de a folosi seriile Taylor ale lui y ca
funcție de x și h pentru a identifica parametrii unei metode de ordin dat. Astfel, pentru
o metodă în trei pași de ordinul trei, i.e., s = p = 3, obținem (din egalarea seriilor
Taylor) urmatoarele patru ecuații algebrice

b1 + b2 + b3 = 1, (3.21)
1
b2 c2 + b3 c3 = , (3.22)
2
1
b2 c22 + b3 c23 = , (3.23)
3
1
b3 a32 c2 = . (3.24)
6

Rezolvarea analitică a acestor ecuații definește familia metodelor Runge-Kutta de ordin

și număr de pași cunoscuți. Astfel, pentru sistemul definit de ecuațiile (3.21)-(3.24) avem
trei mari clase de metode. Prima clasă de metode cuprinde soluțiile pentru care c2 ̸= c3 ,
c2 ̸= 32 , c2 ̸= c3 ̸= 0. Cea de-a doua clasă de metode cuprinde metodele pentru care
c2 = c3 = 23 , b3 ̸= 0, și e definită de tabelul Butcher , în timp ce ultima clasă de soluții
cuprinde metodele pentru care c2 = 32 , c3 = 0, b3 ̸= 0.
Problema determinării metodelor Runge-Kutta de ordin superior este că, spre deosebi-
re de metodele Adams-Bashforth si Adams-Moulton, acestea nu se pot obține inductiv,
folosind informațiile de la o metodă de ordin inferior, fiind necesară determinarea lor de
la zero. Astfel, calculele facute pentru clasa de metode Runge-Kutta cu s = p = 3 sunt
irelevante pentru metodele cu s = p = 4. Acestea din urmă sunt definite de sistemul de
ecuații algebrice

31
3 Metode unipas

c2 c2
( )
c3 3c2 − 3c22 − c3 c3 (c3 − c2 )
c3
c2 (2 − 3c2 ) c2 (2 − 3c2 )
−3c3 + 6c2 c3 + 2 − 3c2 3c3 − 2 2 − 3c2
6c2 c3 6c2 (c3 − c2 ) 6c3 (c3 − c2 )

Tabel 3.3: Metoda Runge-Kutta explicită (generală) cu s = p = 3. Clasa I de soluții

corespunde metodelor pentru care c2 ̸= 0 ̸= c3 ̸= c2 ̸= 23 .

2 2
3 3
2 2 1 1
−
3 3 4b3 4b3
1 3
− b3 b3
4 4
Tabel 3.4: Metoda Runge-Kutta explicită (generală) cu s = p = 3. Clasa II de soluții
corespunde metodelor pentru care c2 = c3 = 23 , b3 ̸= 0.

2 2
3 3
1 1
0 −
4b3 4b3
1 3
− b3 b3
4 4
Tabel 3.5: Metoda Runge-Kutta explicită (generală) cu s = p = 3. Clasa III de soluții
̸ 0.
corespunde metodelor pentru care c2 = 23 , c3 = 0, b3 =

32
 3.2 Metode Runge-Kutta

b1 + b2 + b3 + b4 = 1 (3.25)
1
b1 c1 + b2 c2 + b3 c3 + b4 c4 = (3.26)
2
2 2 2 2 1
b1 c1 + b2 c2 + b3 c3 + b4 c4 = (3.27)
3
1
b3 a32 c2 + b4 (a42 c2 + a43 c3 ) = (3.28)
6
3 3 3 1
b2 c2 + b3 c3 + b4 c4 = (3.29)
4
1
b3 c3 a32 c2 + b4 c4 (a42 c2 + a43 c3 ) = (3.30)
8
2
( 2 2
) 1
b3 a32 c2 + b4 a42 c2 + a43 c3 = (3.31)
12
1
b4 a43 a32 c2 = (3.32)
24
ale cărui soluții sunt mult mai dificile decât cele ale ecuațiilor (3.21)-(3.24). Soluțiile
ecuațiilor (3.25)-(3.32) au fost determinate pentru prima oară de M. W. Kutta care le-a
împărțit în cinci clase distincte, fiecare având mai mulți parametrii liberi. Doua dintre
metodele determinate de Kutta au devenit extrem de populare în literatura de specia-
litate, anume metoda Runge-Kutta (cunoscută în engleză ca the Runge-Kutta method,
prezentată în tabelul Butcher (3.6)) și „regula 3/8” (prezentată în tabelul Butcher (3.7)).
Pentru o discuție detaliată asupra metodelor Runge-Kutta de ordin 4 recomandăm Ref.
[7], secțiunea 3.2 Low order explicit methods, subsecțiunea 3.2.2 Methods of order 4.
În alegerea unei metode contează povara computațională a evaluărilor funcției f . Din
acest motiv metoda din tabelul Butcher (3.6) a devenit metoda Runge-Kutta: având trei
elemente nule în tabelul Butcher necesită cu 30% mai puține evaluări decât metoda (3.7)
la aceeași acuratețe numerică.
În articolele de pionierat asupra metodele Runge-Kutta de la începutul secolului XX
ecuațiile algebrice care descriu o anumită metodă erau determinate manual, prin cal-
cule laborioase care duceau, cum vom vedea, la unele erori. Astăzi însă determinarea
ecuațiilor se face automat, grație așa-numiților arbori Runge-Kutta (a se vedea Ref. [7],
secțiunea 3.1 Order conditions). Deoarece ecuațiile algebrice care definesc metodele de
ordin înalt sunt extrem de complicate prezentăm mai jos doar tabelele Butcher pentru
două din metodele de ordin 5 introduse de Kutta, (3.8) și (3.9). A doua din aceste
metode este însă greșită, ecuațiile algebrice asociate fiind însă rezolvate incorect, varian-
ta corectă a tabelului Butcher fiind publicată de Nyström în 1925, tabel pe care noi îl
reproducem în (3.10).

3.2.2 Metode Runge-Kutta implicite

Cea mai simplă definiție a metodelor Runge-Kutta implicite este că acestea sunt acele
metode Runge-Kutta care au tabele Butcher cu elemente ne-nule deasupra diagonalei.

33
3 Metode unipas

1 1
2 2
1 1
0
2 2

1 0 0 1

1 1 1 1
6 3 3 6
Tabel 3.6: Metoda Runge-Kutta de ordin 4.

1 1
3 3
2 1
− 1
3 3

1 1 −1 1

1 3 3 1
8 8 8 8
Tabel 3.7: Metoda Runge-Kutta de ordin 4, cunoscută drept „regula 3/8”.

34
 3.2 Metode Runge-Kutta

1 1
5 5
2 2
0
5 5
9 15
1 −5
4 4
3 63 9 13 2
− −
5 100 5 20 25
4 6 4 2 8
− 0
5 25 5 15 75
17 25 1 25 25
0 −
144 36 72 72 48
Tabel 3.8: Metodă Runge-Kutta de ordin 5 determinată de Kutta.

1 1
3 3
2 4 6
5 25 25
1 15
1 −3
4 4
2 6 90 50 8
−
3 81 81 81 81
4 7 18 5 4
− 0
5 30 30 30 30
48 125 81 100
0 0 −
192 192 192 192
Tabel 3.9: Metodă Runge-Kutta de ordin 5 determinată de Kutta. Metoda este incorec-
tă, variantă corectă a coeficienților fiind determintă de Nyström în 1925 (a se
vedea tabelul Butcher (3.10)).

35
3 Metode unipas

1 1
3 3
2 4 6
5 25 25
1 15
1 −3
4 4
2 2 10 50 8
−
3 27 9 81 81
4 2 12 2 8
0
5 25 25 15 75
23 125 27 125
0 0 −
192 192 64 192
Tabel 3.10: Metodă Runge-Kutta de ordin 5 corectată de Nyström.

Evident, definiția este prea puțin utilă așa că vom încerca mai jos o abordare diferită.
Ideea fundamentală este următoarea: sistemele algebrice care descriu metodele Runge-
Kutta atunci când niciun element aij nu este în mod obligatoriu zero sunt imposibil de
rezolvat fără ipoteze simplificatoare. Altfel spus, fără a veni cu ipoteze sau presupuneri
cu privire la valorile numerice ale unor elemente din tabelul Butcher, nu putem obține
soluții de genul celor sistematizate în tabelele Butcher (3.3)-(3.5).
În paginile care urmează vom prezenta cea mai importantă clasă de metode Runge-
Kutta implicite ce are la baza cuadraturile gaussiene.

Cuadraturi gaussiene Această metodă are la bază polinoamele Legendre definite pe

intervalul [0, 1] definite ca

P0∗ = 1, (3.33)
P1∗ = 2x − 1, (3.34)
P2∗ = 6x − 6x + 1,
2
(3.35)
P3∗ = 20x − 30x + 12x − 1.
3 2
(3.36)

Cu ajutorul lor putem construi o metode Runge-Kutta de ordin p = 2 folosind doar un

singur pas, i.e., s = 1, folosind pentru c rădăcina ecuației P1∗ = 0, i.e., c = 12 . Metoda
este guvernată de tabelul Butcher
Observăm imediat că metoda este foarte similară metodei implicite Adams-Moulton de
ordin 2 (2.33) și că nu necesita algorimti de calcul dedicați. Lucrurile se schimbă radical
atunci când considerăm o metodă Runge-Kutta de ordin p = 4 folosind doar doi pași,
i.e., s = 2. Să observăm mai întâi că un asemenea nivel de precizie într-un număr așa de

36
 3.2 Metode Runge-Kutta

1 1
2 2

Tabel 3.11: Metodă Runge-Kutta bazată pe cuadraturi gaussiene cu s = 1 și p = 2.

√ √
1 3 1 1 3
− −
2 6 4 4 6
√ √
1 3 1 3 1
+ +
2 6 4 6 4
1 1
2 2
Tabel 3.12: Metodă Runge-Kutta bazată pe cuadraturi gaussiene cu s = 2 și p = 4.

mic de pași nu poate fi atins cu metode explicite. Pentru √ cele două valori √
numerice ale
1 3 1 3
lui c vom utiliza radăcinile ecuației P2∗ = 0, i.e., c1 = − și c2 = + . Calcule
2 6 2 6
detaliate arată ca tabelul Butcher al metodei este cel din (3.12).
Menționăm
√ în cele din urmă că√folosind pentru c rădăcinile ecuației P3∗ = 0, i.e., c1 =
1 15 1 1 15
− , c2 = și c3 = + , putem obține metoda Runge-Kutta implicită cu
2 10 2 2 10
s = 3 și p = 4 definită în tabelul Butcher (3.13).
Se impun câteva precizări privitoare la metodele Runge-Kutta implicite:

1. Faptul că aceste metode oferă un ordin mare cu un număr relativ mic de pași
nu implică și un plus de eficiență a metodelor, căci costul rezolvării numerice, la
fiecare pas, a ecuațiilor (3.20a)-(3.20f), atenție (!), acum implicite, poate fi extrem
de mare pentru sistemele dinamice în care funcția f este puternic neliniară.

2. Aceste metode devin eficiente pentru așa-numitele sisteme rigide, pentru care se

√ √ √
1 15 5 2 15 5 15
− − −
2 10 36 9 15 36 30
√ √
1 5 15 2 5 15
+ −
2 36 24 9 36 24
√ √ √
1 15 5 15 2 15 5
+ + +
2 10 36 30 9 15 36
5 4 5
18 9 18
Tabel 3.13: Metodă Runge-Kutta bazată pe cuadraturi gaussiene cu s = 3 și p = 4.

37
3 Metode unipas

observă suprapunerea a două tipuri de oscilații: unele pe o frecvență foarte joasă

și unele pe o frecvență foarte ridicată. Altfel spus, metodele sunt utile pentru
sistemele de ecuații diferențiale ale caror matrice Jacobi au valori proprii care
acoperă mai multe ordine de marime. Cu cât valoarea proprie cea mai mică e mai
de depărtată de valoarea proprie cea mai mare, cu atât metodele devin mai utile.
Să remarcă că pentru acest tip de probleme metodele Runge-Kutta explicite nu
converg decât pentru valori extremi de mici ale pașilor de integrare, ceea ce le face
ineficiente.

3. Tipul de probleme pentru care metodele Runge-Kutta implicite sunt eficiente este
excelent descris intr-un articol clasic scris de L. F. Shampine și C. W. Gear, A
user’s view of solving stiff ordinary differential equations, SIAM Review 21, 1
(1979): „The problems called ”stiff” are too important to ignore, and are too
expensive to overpower. They are too important to ignore because they occur
in many physically important situations. They are too expensive to overpower
because of their size and the inherent difficulty they present to classical methods,
no matter how great an improvement in computer capacity becomes available.
Even if one can bear the expense, classical methods of solution require so many
steps that roundoff errors may invalidate the solution. It is all the more frustrating
that the solutions of stiff problems look like they should be particularly easy to
compute. [...] By a stiff problem we mean one for which no solution component is
unstable (no eigenvalue has a real part which is at all large and positive) and at
least some component is very stable (at least one eigenvalue has a real part which
is large and negative). Further, we will not call a problem stiff unless its solution
is slowly varying with respect to the most negative real part of the eigenvalues.
(Roughly, we mean that the derivatives of the solution are small compared to
the corresponding derivatives [...]) Consequently, a problem may be stiff for some
intervals of the independent variable and not for others. [...] The essence of the
matter is that for most problems the accuracy requirement dictates the choice of
step size, but for some, the stiff problems, the stability requirement does. [...] One
worry should be dispelled at once. When implemented properly, the instability on
encountering stiffness of classical methods such as Euler’s is automatically detected
and handled by reducing the step size. Computer programs suitable for nonstiff
problems do not ”blow up” in the presence of stiffness, they just become inefficient.”

3.2.3 Metode Runge-Kutta împerecheate. Controlul erorii

Problema metodelor Runge-Kutta împerecheate este una veche și a pornit din nevoia
practică de a putea estima în mod eficient eroarea locală de trunchiere. Am subliniat
în mod eficient pentru că evaluarea în sine a erorii locale de trunchiere nu ridică, de
fapt, nicio problemă. Știind, de(pildă,) că avem o metodă de ordin p care e corectă, preț
de un pas, până la termeni O hp+1 , putem estima eroarea împărțind pasul în 10 și
calculând diferența dintre cele două rezultate. Pentru o estimare mai precisă a erorii
putem împărți pasul în 100 sau 1000 de pași mai mici, însă această abordare este extrem

38
 3.2 Metode Runge-Kutta

1 1
2 2
3 3
0
4 4
2 1 4
1
9 3 9
2 1 4
0
9 3 9
1 1 1 1
8 4 3 8
Tabel 3.14: Metode Runge-Kutta împerecheate. Prima metodă, cea superioară, este de
ordin p = 2 în timp ce a doua, cea inferioară, este de ordin p = 3.

de ineficientă sub raportul timpului de calcul.

Un mod de a rezolva problema are la baza așa-numitele metode Runge-Kutta împere-
cheate. Filozofia acestei abordări este următoarea: dată fiind o metodă Runge-Kutta de
ordin p și o metodă Runge-Kutta de ordin p + 1 care folosește aproape toate evaluările
funcției f ale metodei de ordin p se poate calcula eroarea locala de trunchiere calculând
diferența dintre cele doua estimări. În mod evident, marea problemă este determinarea
acestor două metode Runge-Kutta cuplate care folosesc aceleași estimări ale funcției f .
Una dintre cele mai populare perechi împerecheate este cea determinată de Bogacki și
Shampine, prezentată în tabelul Butcher (3.14). Metoda a fost introdusă de P. Bogacki
și L. F. Shampine în A 3(2) pair of Runge-Kutta formulas, Applied Mathematics Letters
2, 321 (1989), și a devenit folosită pentru ecuații diferențiale care nu necesitau un control
al erorii foarte precis și care nu sunt rigide. Această metodă este folosită de către funția
ode23 din MATLAB.
O alta metode extrem de populară este metoda Merson care se bazează pe împere-
cherea unei metode de ordin 3 cu una de ordin 4. Metoda este ilustrată în tabelul
Butcher (3.15). Metoda este extrem de populară în cadrul comunității utilizatorilor de
FORTRAN fiind implementată în două din rutinele standard ale Numerical Algorithms
Group - NAG, i.e., D02BGF și D02BHF. Metoda nu este implementată în MATLAB și
MATHEMATICA. Metoda Merson are particularitatea că, atunci când e aplicată unui
set de ecuații liniare, metodă inferioară, care este în mod curent o metodă de ordin p = 4
devine o metodă de ordin p = 5.
Menționăm, în final, metoda Zonneveld ca un exemplu de împerechere simplă a două
metode Runge-Kutta de nivel înalt. Tabelul Butcher al metodei este prezentat în tabelul
(3.16). Asemeni metodei lui Merson, și această metodă a aparut într-o lucrare obscură,
recte, J.A. Zonneveld, Automatic integration of ordinary differential equations, Report
R743, Mathematisch Centrum, Postbus 4079, 1009AB Amsterdam, 1963, cunoscând

39
3 Metode unipas

1 1
3 3
1 1 1
3 6 6
1 1 3
0
2 8 8
1 3
1 0 − 2
2 2
1 2 1
0 0
6 3 6
1 3 2 1
0
10 10 5 5
Tabel 3.15: Metode Runge-Kutta-Merson împerecheate. Prima metodă, cea superioară,
este de ordin p = 3 în timp ce a doua, cea inferioară, este de ordin p = 4.

apoi o larga utilizare.

Lista metodelor simplu împerecheate este una lungă, iar foarte scurta noastră trecere
în revistă a omis contribuții remarcabile din partea lui Ceschino, Fehlberg1 , Dormand și
Prince, etc. Pentru o tratare amănunțită a subiectului recomandă cititorilor interesați
Ref. [7], Ref. [8], Ref. [9] și Ref. [11].
În final dorim să facem următoarele comentarii: metodele Runge-Kutta explicite de
ordin 3 prezentate în tabelele Butcher (3.3), (3.4) și (3.5) pot fi transformate în tabe-
le Butcher pentru metode triplu împerecheate. Astfel, tabelul (3.3) se transformă cu
ușurință în (3.17), tabelul (3.4) în (3.18), iar tabelul (3.5) în (3.19).

1
Erwin Fehlberg este considerat a fi unul dintre părinții metodelor Runge-Kutta împerecheate, lucrările
sale de pionierat din perioada când lucra la NASA, toate cu un singur autor, fiind puncte de referință
în domeniu. Deoarece o parte din ele au fost avut o bună perioadă un caracter secret, abia recent
fiind disponibile pe serverul de documente al NASA, menționăm aici pe cele mai importante dintre
ele: Technical Note D-1834, On the error propagation of some interpolation formulas for second-order
differential equations, July 1963; Nasa Technical Report - 248, New one-step integration methods of
high-order accuracy applied to some problems in celestial mechanics, October 1966; Nasa Technical
Report - 287, Classical fifth-, sixth-, seventh-, and eighth-order Runge-Kutta formulas with stepsize
control, October 1968; Nasa Technical Report - 315, Low order classical Runge-Kutta formulas with
stepsize control and their application to some heat transfer problems, July 1969; Nasa Technical Report
- 352, Some experimental results concerning the error propagation in Runge-Kutta type integration
formulas, October 1970; Nasa Technical Report - 381, Classical eightth- and lower-order Runge-
Kutta-Nystrom formulas with stepsize control for special second-order differential equations, March
1972; Nasa Technical Report - 410, Classical eightth- and lower-order Runge-Kutta-Nyström formulas
with a new stepsize control procedure for special second-order differential equations, June 1973; Nasa
Technical Report - 432, Classical seventh-, sixth-, and fifty-order Runge-Kutta-Nyström formulas with
stepsize control for general second-order differential equations, October 1974.

40
 3.2 Metode Runge-Kutta

1 1
2 2
1 1
0
2 2

1 0 0 1

3 5 7 13 1
−
4 32 32 32 32
1 1 1 1
0
6 3 3 6
1 7 7 13 16
− −
2 3 3 6 3
Tabel 3.16: Metode Runge-Kutta-Zonneveld împerecheate. Prima metodă, cea superi-
oară, este de ordin p = 4 în timp ce a doua, cea inferioară, este de ordin
p = 5.

c2 c2

c3 (3c2 (1 − c2 ) − c3 ) c3 (c3 − c2 )
c3
c2 (2 − 3c2 ) c2 (2 − 3c2 )
2 − 3 (c2 + c3 ) + 6c2 c3 2 − 3c3 2 − 3c2
ordin 3
6c2 c3 6c2 (c2 − c3 ) 6c3 (c3 − c2 )
1 1
1− ordin 2
2c2 2c2

1 ordin 1

Tabel 3.17: Metode Runge-Kutta triplu împerecheate făcute pe baza tabelului Butcher
(3.3).

41
3 Metode unipas

2 2
3 3
1 1
0 −
4b3 4b3
1 3
− b3 b3 ordin 3
4 4
1 3
ordin 2
4 4

1 ordin 1

Tabel 3.18: Metode Runge-Kutta triplu împerecheate făcute pe baza tabelului Butcher
(3.4).

2 2
3 3
2 2 1 1
−
3 3 4b3 4b3
1 3
− b3 b3 ordin 3
4 4
1 3
ordin 2
4 4

1 ordin 1

Tabel 3.19: Metode Runge-Kutta triplu împerecheate făcute pe baza tabelului Butcher
(3.5).

42
 3.3 Implementări software disponibile

3.3 Implementări software disponibile

Cea mai bună implementarea a metodelor Runge-Kutta este cea oferită de MATHEMA-
TICA, care permite rezolvarea ecuațiilor diferențiale ordinare folosind atît metode expli-
cite, cât și metode implicite. Pentru metodele explicite MATHEMATICA are o serie de
metode Runge-Kutta împerecheate (ExplicitRungeKutta), cea mai puțin eficientă fiind
de ordin 2(1), iar cea mai eficientă fiind de ordin 9(8). Pentru metodele Runge-Kutta im-
plicite MATHEMATICA acoperă (prin opțiunea ImplicitRungeKutta) clasa metodelor
bazate pe cuadraturi gaussiene (prezentată de noi mai sus), alături de numeroase alte me-
tode implicite cum sunt Runge-Kutta-Lobatto-IIIA, Runge-Kutta-Lobatto-IIIB, Runge-
Kutta-Lobatto-IIIC, Runge-Kutta-Radau-IA, Runge-Kutta-Radau-IIA. Mai mult, MA-
THEMATICA are un pachet dedicat integratorilor simplectici (SymplecticPartitione-
dRungeKutta), care descriu precis sistemele hamiltoniene și sunt extrem de utili în mode-
larea numerică a sistemelor fizice. Principala funcție folosită pentru rezolvarea numerică
a ecuațiilor diferențiale ordinare este NDSolve care acceptă ca parametri tipul metodei
(Method), ordinul metodei (DifferenceOrder), tabelul de coeficienți (Coefficients, în ca-
zul în care metoda este una specială), eroarea (AccuracyGoal și PrecisionGoal), marimea
pasului inițial folosit pentru inițializarea integrării (StartingStepSize), marimea maximă
a pasului (MaxStepSize), numărul maxim de pași de integrare (MaxSteps), etc. Integra-
torul numeric poate verifica numeric dacă ecuația diferențială este rigidă sau nu, și, mai
important, poate schimba tipul și ordinul metodei pentru a obține un anumit nivel de
eroare.
Filosofia din spatele integratorului numeric folosit de MATHEMATICA a fost aceea
de a strange cele mai eficiente metode existente și a le oferi spre utilizare într-un mod
unificat. La polul opus, MATLAB are o selecție destul de săracă de integratori numerici,
care reflectă în bună măsură temele de cercetare ale creatorilor suitei, L. F. Shampine
și M. W. Reicheit (a se vedea în acest sens articolul The MATLAB ODE suite, SIAM
Journal of Scientific Computing 18, 1 (1997)). Suita MATLAB nu cuprinde integratori
Runge-Kutta impliciți de ordin înalt (ci doar trei integratori de ordin 3(2) care descriu
sisteme nu foarte rigide, i.e., ode23s, ode23t și ode23tb), iar la capitolul metode Runge-
Kutta explicite sunt implementate două metode împerecheată, recte aceea dezvoltată de
Bogacki și Shampine, prezentată de noi în tabelul Butcher (3.14), și o metodă de ordin
5(4) dezvoltată de Dormand și Prince (neprezentată de noi anterior datorită coeficienților
destul de complicați, dar care poate fi găsită în Ref. [11]). Toți integratorii numerici din
MATLAB acceptă un set standardizat de parametri care includ eroarea relativă (RelTol)
și eroarea absolută (AbsTol), pasul inițial de integrare (InitialStep), numărul maxim de
pași de integrare (MaxStep) și, foarte important, matricea Jacobi a sistemului de ecuații
diferențiale (Jacobian și JPattern). Această din urmă opțiune este extrem de utilă în
cazul ecuațiilor rigide și sporește considerabil eficiența integratorilor numerici.
Chiar daca codurile MATLAB dedicate ecuațiilor diferențiale ordinare sunt neimpre-
sionante, există o librărie - GniCodes - dezvoltată de E. Hairer și M. Hairer în care
sunt implementate metode Runge-Kutta implicite dezvoltate anume pentru ecuații di-
ferențiale rigide și cu invarianți (i.e., ecuații de tip Hamilton). Librăria este descrisă
pe larg în GniCodes - Matlab progras for geometric numerical integration, E. Hai-

43
3 Metode unipas

rer și M. Hairer, Frontiers in Numerical Analysis, Universitext, Springer (2003), co-

durile propriu-zise fiind prezentate în detaliu pe site-ul personal al lui Ernst Hairer,
[Link] O altă librărie populară în rândul
utilizatorilor este DiffMan, introduă în DiffMan - an object oriented MATLAB toolbox for
solving differential equations on manifolds, K. Engo, A. Marthinsen, și H. Munthe-Kaas,
raport al Departamentului de Informatică al Universității din Bergen, Norvegia, dispo-
nibil la adresa [Link]

3.4 Exerciții și probleme

Exercițiul 1: Arătați că metoda Runge-Kutta descrisă de tabelul Butcher

γ γ 0

1−γ 1 − 2γ γ ,

1 1
2 2
( √ )
unde γ = 3 ± 3 /6, este de ordinul 3. Comentați influența elementului nul din tabelul Butcher asupra
schemei numerice.
Exercițiul 2: Să se arate că tabelul Butcher

2 2
3 3
2 2
0
3 3
1 3
4 4
1 3 3
4 8 8
descrie o metoda Runge-Kutta cuplată de ordinul 2(3).
Problema 1: Să se rezolve numeric sistemul Lorentz

dx
= σ (y − x) (P.3.1a)
dt
dy
= x (ρ − z) − y (P.3.1b)
dt
dz
= xy − βz (P.3.1c)
dt
pentru ρ = 28, σ = 10, β = 8/3 și condițiile inițiale x(0) =, y(0) =, z(0) = și să se determine traiectoriile
din spațiul fazelor folosind metode Runge-Kutta implicite și explicite (de exemplu, Zonneveld și Gauss-
Legendre cu s = 2, p = 4). Se justifică utilizarea metodelor Runge-Kutta implicite? Să se determine
traiectoriile din spațiul fazelor pentru ρ = 99.96, păstrând ceilalți parametri constanți. Ce diferențe se
observă?

44
 3.5 Bibliografie

Problema 2: Să se rezolve numeric sistemul Rössler

dx
= −y − z (P.3.2a)
dt
dy
= x + ay (P.3.2b)
dt
dz
= b + z (x − c) (P.3.2c)
dt
și să se determine traiectoriile soluțiilor în planul x − y pentru a = b = 1 și c = 4, 6, 8.5, 8.7, 9, 12, 13
și 18. Se vor folosi condițiile inițiale x(0) =, y(0) =, z(0) = Să se determine perioada traiectoriilor
determinate numeric.

3.5 Bibliografie
1. Gear, C. W., Numerical initial value problems in ordinary differential equations, Prentice-Hall
(1971)
2. Iserles, A., A first course in the numerical analysis of differential equations, Cambridge University
Press (2009)
3. Ascher, U. M., Numerical methods for evolutionary differential equations, Society for Industrial
and Applied Mathematics (2008)
4. Griffiths, D. F., și Higham, D. J., Numerical methods for ordinary differential equations. Initial
value problems, Springer (2010)
5. Atkinson, K. E., Han, W., și Stewart, D. E., Numerical solution of ordinary differential equations,
John Wiley & Sons (2009)
6. Butcher, J. C., A history of Runge-Kutta methods, Applied Numerical Mathematics 20, 247
(1996).
7. Butcher, J. C., Numerical methods for ordinary differential equations, John Wiley & Sons (2008)
8. Hairer, E., Nørsett, S.P., și Wanner, G., Solving ordinary differential equations I. Nonstiff Pro-
blems, Springer (2008)
9. Hairer, E., și Wanner, G., Solving ordinary differential equations II. Stiff and differential-algebraic
problems, Springer (2010)
10. Hairer, E., Lubich, C., și Wanner, G., Geometric numerical integration. Structure-preserving
algorithms for ordinary differential equations, Springer (2000)
11. Dormand, J. R., Numerical methods for differential equations. A computational approach, CRC
Press (1996)
12. Leimkuhler, B., și Reich, S., Simulating Hamiltonian dynamics, Cambridge University Press
(2004)
13. Shampine, L. F., Gladwell, I., și Thompson, S., Solving ODEs with Matlab, Cambridge University
Press (2003)
14. Verner, J. H., Families of imbedded Runge-Kutta methods, SIAM Journal of Numerical Analysis
16, 857 (1979)

45