Programul de studii: Matematic a si Informatic a Aplicat a n Inginerie [605263]

UNIVERSITATEA POLITEHNICA BUCURES TI
FACULTATEA DE S TIINT  E APLICATE
Programul de studii: Matematic a  si Informatic a Aplicat a ^ n Inginerie
PROIECT DE DIPLOM A
Conduc ator  stiint ic , Absolvent: [anonimizat] sor COSTEA Florina-Adriana BOSIE
Bucure sti
2019

UNIVERSITATEA POLITEHNICA DIN BUCUREȘTI
FACULTATEA DE ȘTIINȚE APLICATE

Programul de studii: Matematică și Informatică Aplicată în Inginerie

Aprobat Decan,
Prof. dr. Emil PETRESCU

PROIECT DE DIPLOM Ă

Metode numerice pentru probleme
cu condiție inițială

CONDUCĂTOR ȘTIINȚIFIC, ABSOLVENT: [anonimizat] -Adriana BOSIE

București

2019

Introducere
^In lucrarea de fat  a sunt prezentate diverse metode analitice  si numerice de rezolvare
a unor ecuat ii diferent iale cu condit ie init ial a. Ecuat iile diferent iale de ordin npot
descrise sub forma
y(n)(x) =f(x;y(x);y0(x);:::;y(n1)(x)):
Numim problem a Cauchy o ecuat ie diferent ial a de ordinul ^ nt^ ai a c arei solut ie
trebuie s a verice o anumit a condit ie init ial a. Astfel, o problem a Cauchy ata sat a unei
ecuat ii diferent iale ordinare are urm atoarea form a
y0=f(x;y); x2R
y(x0) =y0:
Pe parcursul lucr arii, vom presupune c a x2[a;b]  si vom considera o diviziune
echidistant a a intervalului [ a;b], adic a o diviziune care are proprietatea c a, ^ ntre oricare
dou a noduri consecutive distant a este aceea si. A sadar, avem
a=x0<x 1<x 2<


<xN=b;
cuxn+1=xn+h,n=0;N1, iar pasul de discretizare h=ba
N.
Lucrarea este structurat a ^ n trei capitole, care includ prezentarea teoretic a a me-
todelor numerice, c^ at  si aplicat ii relevante.
^In capitolul 1 sunt descrise c^ ateva metode analitice de aproximare a solut iilor unor
probleme Cauchy. Mai ^ nt^ ai este prezentat a teorema lui Picard, care ofer a condit ii
suciente pentru existent a  si unicitatea solut iei pe un interval centrat ^ n punctul oferit
de condit ia init ial a. ^In general, solut ia unor astfel de probleme nu poate  determinat a
explicit, de aceea este necesar a implementarea unor metode numerice care s a permit a
aproximarea solut iei cu o eroare sucient de mic a.
Metodele prezentate ^ n primul capitol sunt: Metoda lui Picard, Metoda lui Taylor,
Metoda lui Euler, Metoda predictor-corector (Metoda lui Euler ^ mbun at at it a), care
i

Florina-Adriana Bosie Metode numerice pentru probleme cu condit ie init ial a
presupune determinarea solut iei prin doi pa si, unul predictor, reprezentat de metoda
lui Euler  si altul corector, reprezentat de metoda Runge-Kutta de ordinul II, care va
 descris a ^ n capitolul urm ator  si Metoda coecient ilor nedeterminat i. La nalul capi-
tolului sunt prezentate c^ ateva aplicat ii ale acestor metode, programele Maple utilizate
pentru rezolvarea lor  si grace care ilustreaz a apropierea solut iilor aproximative de
solut ia exact a pe m asur a ce pasul de discretizare scade sau num arul de iterat ii cre ste.
^In cel de-al doilea capitol sunt abordate metodele numerice unipas (directe) de tip
Runge-Kutta de ordinul I, II, III  si IV. Aceste metode ^ ncearc a s a obt in a o precizie
rezonabil a, f ar a a avea inconvenientul c a pasul de discretizare trebuie s a aib a o dimen-
siune foarte mic a, ca ^ n cazul metodei lui Euler  si f ar a a necesita calculul derivatelor
de ordin superior, ca ^ n cazul metodei lui Taylor. Metodele de tip Runge-Kutta au
importantul avantaj de a  metode cu autostart, adic a metode ^ n care este necesar a
cunoa sterea valorii lui ydoar ^ n punctul xnpentru a
area solut iei ^ n punctul urm ator,
xn+1. La fel ca primul capitol,  si acesta se ^ ncheie cu aplicat ii ale metodelor des-
crise teoretic, programe Maple utilizate ^ n rezolvarea lor  si grace care arat a cre sterea
preciziei solut iei aproximative odat a cu mic sorarea pasului de discretizare.
Ultimul capitol al lucr arii cont ine metodele numerice multipas (indirecte) pentru
probleme cu condit ie init ial a. ^In timp ce metodele numerice unipas (directe), prezen-
tate anterior, necesitau cunoa sterea valorii lui ydoar ^ ntr-un singur punct, xn, pentru
a
area solut iei ^ n punctul urm ator, xn+1, aceste metode necesit a valoarea ^ n punctul
xn, c^ at  si un num ar prestabilit de valori ^ n punctele anterioare xn1;xn2;:::;xnk
din discretizare.
Metodele prezentate ^ n acest capitol sunt: Metoda Adams-Bashforth  si Metoda
Adams-Moulton. Capitolul se ^ ncheie cu aplicat ii ale acestor metode, precum  si grace
care descriu apropierea solut iilor aproximative de solut ia exact a a problemei conside-
rate, pe masur a ce pasul de discretizare scade.
ii

Cuprins
Introducere 1
1 Metode analitice  si numerice pentru ecuat ii diferent iale 2
1.1 Metoda lui Picard (teorema de existent  a  si unicitate) . . . . . . . . . . 2
1.2 Metoda lui Taylor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.3 Metoda lui Euler . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.4 Metoda predictor-corector (metoda lui Euler ^ mbun at at it a) . . . . . . . 9
1.5 Metoda coecient ilor nedeterminat i . . . . . . . . . . . . . . . . . . . . 12
1.6 Aplicat ii . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
2 Metode de tip Runge-Kutta 22
2.1 Metoda Runge-Kutta de ordinul I . . . . . . . . . . . . . . . . . . . . . 23
2.2 Metoda Runge-Kutta de ordinul II . . . . . . . . . . . . . . . . . . . . 24
2.3 Metoda Runge-Kutta de ordinul III . . . . . . . . . . . . . . . . . . . . 25
2.4 Metoda Runge-Kutta de ordinul IV . . . . . . . . . . . . . . . . . . . . 26
2.5 Aplicat ii . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
3 Metode numerice multipas 39
3.1 Metoda Adams-Bashforth . . . . . . . . . . . . . . . . . . . . . . . . . 39
3.2 Metoda Adams-Moulton . . . . . . . . . . . . . . . . . . . . . . . . . . 46
3.3 Aplicat ii . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
Bibliograe 55
1

Capitolul 1
Metode analitice  si numerice
pentru ecuat ii diferent iale
1.1 Metoda lui Picard (teorema de existent  a  si uni-
citate)
^In acest capitol vom prezenta metode prin care se poate aproxima solut ia unei probleme
Cauchy de tipul:y0=f(x;y);
y(x0) =y0:(1.1)
Prezent am ^ n continuare condit ii suciente pentru existent a  si unicitatea solut iei
pentru o astfel de problem a.
Denit ia 1.1. Spunem c a o funct ie f:DR2!Rsatisface condit ia Lipschitz ^ n
raport cu a doua variabil a dac a exist a o constant a L>0astfel ^ nc^ at
jf(x;y 1)f(x;y 2)jLjy1y2j;8(x;y 1);(x;y 2)2D:
ConstantaLeste numit a constant a Lipschitz pentruf.
Teorema 1.2. (S uli & Mayers [3]) Fie f:D!Ro funct ie continu a pe dome-
niul dreptunghiular D=f(x;y)2R2:x0xXM; y0Cyy0+Cg. Pre-
supunem c a fsatisface condit ia Lipschitz ^ n raport cu y,jf(x;y 0)j K, pentru
x0xXM si
CK
L(eL(XMx0)1): (1.2)
2

Florina-Adriana Bosie Metode numerice pentru probleme cu condit ie init ial a
Atunci, exist a o unic a funct ie y: [x0;XM]!Rastfel ^ nc^ at y(x0) =y0 siy0=f(x;y),
pentrux2[x0;XM].^In plus,
jy(x)y0jC; x 0xXM:
Demonstrat ie. Denim  sirul de funct ii ( yn)n0prin:
y0(x)y0;
yn(x) =y0+Zx
x0f(t;yn(t))dt; n = 1;2;::: : (1.3)
Deoarecefeste continu a pe D, este evident c a ecare funct ie yneste continu a pe
[x0;XM]. Mai mult, deoarece
yn+1=y0+Zx
x0f(t;yn(t))dt;
rezult a, prin sc adere, c a:
yn+1yn=Zx
x0[f(t;yn(t))f(t;yn1(t))]dt: (1.4)
Acum, continu am prin induct ie  si presupunem c a pentru o valoare pozitiv a a lui n,
jyn(x)yn1(x)jK
L[L(xx0)]n
n!; x 0xXM (1.5)
 si c a
jyk(x)y0jK
LkX
j=1[L(xx0)]j
j!; x 0xXM; k= 1;:::;n: (1.6)
Trivial, ipoteza teoremei  si (1.3) implic a c a (1.5)  si (1.6) sunt adev arate pentru
n= 1. Acum, (1.6)  si (1.2) conduc la:
jyk(x)y0jK
L(eL(XMx0)1)C; x 0xXM; k= 1;:::;n:
Prin urmare ( x;yn1(x))2D si (x;yn(x))2D, pentru oricare x2[x0;XM]. A sadar,
folosind (1.4), condit ia Lipschitz  si (1.5),
jyn1(x)yn(x)jLZx
x0K
L[L(tx0)]n
n!dt
=K
L[L(xx0)]n+1
(n+ 1)!; (1.7)
3

Florina-Adriana Bosie Metode numerice pentru probleme cu condit ie init ial a
oricare ar  x2[x0;XM].^In plus, folosind (1.7)  si (1.6),
jyn+1y0jjyn+1(x)yn(x)j+jyn(x)y0j
K
L[L(xx0)]n+1
(n+ 1)!+K
LnX
j=1[L(xx0)]j
j!
=K
Ln+1X
j=1[L(xx0)]j+1
(j+ 1)!; (1.8)
oricare ar  x2[x0;XM]. Prin urmare, (1.5)  si (1.6) sunt adev arate pentru n^ nlocuit
cun+1, rezult a c a prin induct ie, sunt adev arate pentru toate numerele ^ ntregi pozitive
n.
Deoarece seria innit a1X
j=1cj
j!converge (c atre ec1) pentru orice valoare a lui
c2R, ^ n particular pentru c=L(XMx0), rezult a din (1.5) c a seria innit a
1X
j=1[yj(x)yj1(x)] converge absolut  si uniform pentru x2[x0;XM].
^In orice caz,
y0+nX
j=1[yj(x)yj1(x)] =yn(x);
arat a c a  sirul de funct ii continue ( yn) converge c atre o limit a, uniform pe [ x0;XM]  si,
prin urmare, c a limita respectiv a este o funct ie continu a.
Not^ and aceast a limit a cu y, deducem din (1.3) c a:
y(x) = lim
n!1yn+1(x)
=y0+ lim
n!1Zx
x0f(t;yn(t))dt
=y0+Zx
x0lim
n!1f(t;yn(t))dt
=y0+Zx
x0f(t;y(t))dt; (1.9)
unde am folosit convergent a uniform a a  sirului de funct ii ( yn) ^ n trecerea de la linia doi
la linia trei, pentru a schimba ordinea limitei  si a integralei  si continuitatea funct iei f
^ n trecerea de la linia trei la linia patru.
Cumt!f(t;y(t)) este o funct ie continu a ^ n tpe intervalul [ x0;XM], integrala
sa pe intervalul [ x0;x] este o funct ie diferent iabil a continu a ^ n x. Prin urmare, din
4

Florina-Adriana Bosie Metode numerice pentru probleme cu condit ie init ial a
(1.9),yeste o funct ie diferent iabil a continu a ^ n xpe [x0;XM], adic ay2C1[x0;XM].
Diferent iind (1.9) deducem c a:
y0=f(x;y);
a sa cum este necesar. De asemenea y(x0) =y0.
Am v azut deja c a ( x;yn(x))2Dc^ andx0xXM. Deoarece Deste o mult ime
^ nchis a ^ n R2, c^ andn!1 , rezult a c a  si ( x;y(x))2Dpentrux0xXM.
Pentru a ar ata c a solut ia problemei cu valoare init ial a este unic a, presupunem, prin
metoda reducerii la absurd, c a exist a dou a solut ii diferite, y siz. Atunci, prin sc adere,
y(x)z(x) =Zx
x0[f(t;y(t))f(t;z(t))]dt; x2[x0;XM];
din care rezult a c a
jy(x)z(x)jLZx
x0jy(t)z(t)jdt; (1.10)
oricare ar  x2[x0;XM]. Presupunem c a meste valoarea maxim a a expresiei
jy(x)z(x)jpentrux0xXM sim> 0.
Atunci,
jy(x)z(x)jmL(xx0); x 0xXM:
^Inlocuind aceast a inegalitate ^ n membrul drept al ecuat iei (1.10), g asim
jy(x)z(x)jL2mZx
x0(tx0)dt=m[L(xx0)]2
2!:
Similar, prin induct ie, rezult a c a
jy(x)z(x)jm[L(xx0)]k
k!; k = 1;2;::: ;
oricare ar  x2[x0;XM].^In orice caz, membrul drept al ultimei inegalit at i este
m arginit superior de m[L(xx0)]k
k!oricare ar  x2[x0;XM], care poate  mic sorat,
aleg^ andksucient de mare. Prin urmare,
jy(x)z(x)j= 0;8×2[x0;XM]:
Rezult a c a solut iile y sizsunt identice.

Demonstrat ia Teoremei 1.2 sugereaz a urm atoarea abordare numeric a de aproxi-
mare a solut iei Problemei 1.1 pe domenii rectangulare, adic a f:DR2!R si
5

Florina-Adriana Bosie Metode numerice pentru probleme cu condit ie init ial a
D:=f(x;y)2R2jjxx0ja;jyy0jbg, numit a metoda lui Picard saumetoda
aproximat iilor succesive :
8
<
:yn+1(x) =y0+Zx
x0f(t;yn(t))dt; n0
y0dat:(1.11)
S irul de funct ii ( yn)n0converge uniform la solut ia exact a y(x). Dezavantajul major
al metodei const a ^ n faptul c a necesit a integrarea funct iilor. Principalul avantaj al me-
todei este faptul c a furnizeaz a o aproximare global a pentru solut ia problemei Cauchy.
De asemenea, metoda este rapid convergent a:
jyn(x)y(x)jCjxx0jn+1
(n+ 1)!;
undeCeste o constant a ce depinde de f.
Remarca 1.3. Metoda lui Picard poate  utilizat a  si pentru aproximarea solut iilor
unor sisteme de ecuat ii diferent iale cu condit ii init iale de tipul:
8
>>>><
>>>>:y0
1=f1(x;y 1;:::;yN);
y0
2=f2(x;y 1;:::;yN);
::: ::: ::: ::: ::: :::
y0
N=fN(x;y 1;:::;yN);
yj(x0) =y0
k; k21;N;(1.12)
dac a funct iile fk:DRN+1!R,D=f(x;y 1;:::;yN)2RN+1:jxx0j 
a;jyy0
kjbk; k21;Ngsunt Lipschitz ^ n raport cu variabila yk.
De exemplu, pentru un sistem de dou a ecuat ii diferent iale
8
>><
>>:y0=f(x;y;z );
z0=g(x;y;z );
y(x0) =y0;
z(x0) =z0;
metoda aproximat iilor succesive presupune generarea  sirurilor de funct ii (yn)n si(zn)n
astfel: 8
>>>>>><
>>>>>>:yn+1(x) =y0+Zx
x0f(t;yn(t);zn(t))dt;
zn+1(x) =z0+Zx
x0g(t;yn(t);zn(t))dt;
y0; z0date:(1.13)
6

Florina-Adriana Bosie Metode numerice pentru probleme cu condit ie init ial a
1.2 Metoda lui Taylor
Consider am o ecuat ie diferent ial a de ordinul ^ nt^ ai, cu valori init iale, de forma
y0=f(x;y); y(x0) =y0: (1.14)
Funct iafpoate  liniar a sau neliniar a, dar presupunem c a feste diferent iabil a ^ n
raport cux siy. S tim c a (1.14) are solut ie unic a dac a@f
@yeste continu a pe domeniul
de interes. Dac a y(x) este solut ia exact a a ecuat iei (1.14), putem c auta y(x) sub forma
unei serii Taylor ^ n jurul punctului x=x0, astfel:
y(x) =y0+ (xx0)y0(x0) +(xx0)2
2!y00(x0) +


: (1.15)
Derivatele din aceast a dezvoltare nu se cunosc ^ n mod explicit deoarece solut ia nu este
cunoscut a. ^Ins a, dac a feste diferent iabil a, ele pot  obt inute deriv^ and total (1.14)
^ n raport cu x si lu^ and ^ n considerare faptul c a yeste o funct ie de x. Pentru primele
c^ ateva derivate, obt inem:
y0=f(x;y);
y00=f0=fx+fyy0=fx+fyf;
y000=f00=fxx+fxyf+fyxf+fyyf2+fyfx+f2
yf
=fxx+ 2fxyf+fyyf2+fxfy+f2
yf: (1.16)
Continu^ and ^ n acest mod, putem exprima orice derivat a a lui y^ n funct ie de f(x;y)
 si derivatele sale part iale. ^Insa, except^ and cazul ^ n care f(x;y) este o funct ie foarte
simpl a, derivatele de ordin superior devin din ce ^ n ce mai complexe. Deci, din motive
practice, trebuie s a reducem num arul termenilor din dezvoltarea (1.15) la un num ar
rezonabil, iar aceast a restrict ie conduce la o restrict ie asupra valorilor lui x, pentru
care (1.15) este o aproximare rezonabil a. Dac a presupunem c a seria trunchiat a (1.15)
conduce la o bun a aproximare pentru un pas de lungime h, pentruxx0=h, putem
evaluaylax0+h. Apoi folosim (1.15) pentru a trece la pasul urm ator. Dac a continu am
^ n acest mod, obt inem o mult ime discret a a valorilor yncare sunt aproxim ari ale solut iei
exacte ^ n punctul xn=x0+nh(n= 0;1;2;:::).^In aceast a sect iune trebuie indicat a
^ ntotdeauna valoarea solut iei exacte ^ ntr-un punct xn, dat a dey(xn)  si cea a unei
solut ii aproximative dat a de yn.
Pentru a dezvolta aceast a procedur a, mai introducem operatorul
Tk(x;y) =f(x;y) +h
2!f0(x;y) +


+hk1
k!f(k1)(x;y); k = 1;2;:::; (1.17)
7

Florina-Adriana Bosie Metode numerice pentru probleme cu condit ie init ial a
unde presupunem c a este folosit un pas hde m arime xat a  si unde f(j)indic a aj-a
derivat a total a a funct iei f(x;y(x)) ^ n raport cu x. Atunci putem introduce urm atorul
algoritm:
Algoritm. (Conte & de Boor [1] )Algoritmul lui Taylor de ordin kpentru a g asi
o solut ie aproximativ a a ecuat iei diferent iale
y0=f(x;y);
y(a) =y0
pe un interval [ a;b]:
1.Alegem un pas h=(ba)
N. Punem
xn=a+nh; n = 0;1;:::;N:
2.Gener am aproximarea yna luiy(xn) din recursivitatea
yn+1=yn+hTk(xn;yn); n = 0;1;:::;N1;
undeTk(xn;yn) este denit de (1.17).
At^ at algoritmul lui Taylor, c^ at  si alte metode bazate pe acest algoritm, care calcu-
leaz aylax=xn+1folosind doar informat ii despre y siy0^ ntr-un singur punct, x=xn,
sunt numite metode unipas.
Eroarea local a a algoritmului lui Taylor de ordin keste:
"=hk+1f(k)(;y())
(k+ 1)!; xn< <xn+h
=hk+1
(k+ 1)!y(k+1)():
Spunem c a algoritmul lui Taylor este de ordinul kdac a eroarea local a ", denit a
mai sus, esteO(hk+1).
1.3 Metoda lui Euler
Fix^ andk= 1 ^ n algoritmul anterior, obt inem metoda lui Euler  si eroarea sa local a:
yn+1=yn+hf(xn;yn); " =h2
2y00(): (1.18)
8

Florina-Adriana Bosie Metode numerice pentru probleme cu condit ie init ial a
Pentru a obt ine o precizie bun a cu metoda lui Euler, trebuie s a lu am o valoare
relativ mic a a lui h. Din acest motiv, metoda lui Euler nu este folosit a, ^ n general,
pentru integrarea ecuat iilor diferent iale.
Desigur, putem aplica algoritmul lui Taylor pentru un ordin mai mare ^ n scopul
obt inerii unei precizii mai bune  si, ^ n general, cu c^ at ordinul algoritmului este mai
mare, cu at^ at precizia este mai bun a, pentru o m arime a pasului dat a. Dac a f(x;y)
este o funct ie relativ simpl a de x siy, atunci putem genera rapid derivatele cerute, pe
un calculator, folosind diferent ierea simbolic a sau altfel, t in^ and cont de propriet at ile
particulare pe care le poate avea funct ia f(x;y).
^In orice caz, necesitatea de a calcula derivate de ordin superior face algoritmul lui
Taylor neadecvat pentru calculatoarele de mare vitez a ^ n scopul general al integr arii.
Cu toate acestea, este de mare interes teoretic deoarece majoritatea metodelor practice
a steapt a s a obt in a aceea si precizie ca  si algoritmul lui Taylor pentru un ordin dat, f ar a
dezavantajul de a  necesar calculul derivatelor de ordin superior.
De si algoritmul lui Taylor nu este folosit, ^ n general, ^ n scopuri practice, cazul
special al metodei lui Euler este luat ^ n considerare datorit a implicat iilor sale teoretice.
1.4 Metoda predictor-corector (metoda lui Euler
^ mbun at at it a)
^In aceast a sect iune vom introduce denumirea de formul a de tip ^ nchis. Spunem c a
o formul a este de tip deschis dac a este obt inut a cu ajutorul polinoamelor prin care
s-a determinat valoarea funct iei ^ n punctul xnpe baza valorilor funct iei dinainte de
xn. Formulele de tip ^ nchis sunt obt inute pe baza interpol arii polinomiale ^ n punctul
xn+1; xn si punctele dinainte de xn. Cea mai simpl a formul a de acest tip se obt ine
dac a aproxim am
yn+1=yn+6X
i=1ciki
cu formula trapezului ( Conte & de Boor [1] ):
I(f)T=1
2(ba)[f(a) +f(b)]; eroare =f00()(ba)3
12; 2(a;b):
Aceasta conduce la formula
yn+1=yn+h
2[f(xn;yn) +f(xn+1;yn+1)]; n = 1;2;::: : (1.19)
9

Florina-Adriana Bosie Metode numerice pentru probleme cu condit ie init ial a
Eroarea^ n aceast a formul a este h3
12y000 si asta^ nseamn a o^ mbun at at ire fat  a de metoda
lui Euler. ^In orice caz, (1.19) este o ecuat ie implicit a de yn+1deoareceyn+1apare ca
 si argument ^ n membrul drept.
Dac af(x;y) este o funct ie neliniar a, ^ n general, nu o s a putem rezolva (1.19) pentru
o valoare exact a yn+1. Cu toate acestea, putem ^ ncerca s a obt inem yn+1prin mijloace
de iterat ie. Prin urmare, x^ and xn, obt inem aproximarea y(0)
n+1a luiyn+1cu formula
lui Euler
y(0)
n+1=yn+hf(xn;yn): (1.20)
Evalu amf(xn+1;y(0)
n+1)  si ^ nlocuim ^ n membrul drept al formulei (1.19) pentru a obt ine
aproximarea
y(1)
n+1=h
2[f(xn;yn) +f(xn+1;y(0)
n+1)]:
Apoi evalu am f(xn+1;y(1)
n+1)  si, folosind din nou (1.19), obt inem urm atoarea aproxi-
mare. ^In general, iterat ia este denit a de
y(k)
n+1=yn+h
2[f(xn;yn) +f(xn+1;y(k1)
n+1)]; k = 1;2;::: : (1.21)
Iterat ia ia sf^ ar sit c^ and dou a iterat ii succesive produc precizia dorit a. Aceast a pro-
cedur a pentru obt inerea unor valori ^ mbun at at ite yn+1, cuxn+1xat, este numit a,
uneori, iterat ie intern a pentru a  deosebit a de (1.19), care este folosit a pentru a
genera valorile lui yn; n= 0;1;::: . Vom rezuma aceast a procedur a ^ n urm atorul
algoritm:
Algoritm. (Conte & de Boor [1] )Metoda predictor-corector de ordinul II.
Pentru ecuat ia diferent ial a
y0=f(x;y); y(x0) =y0;
cuhdat  sixn=x0+nh, pentrunxat,n= 0;1;::::
1.Calcul amy(0)
n+1folosind (1.20).
2.Calcul amy(k)
n+1(k= 1;2;:::) folosind (1.21). Iter am p^ an a c^ and
jy(k)
n+1y(k1)
n+1j
jy(k)
n+1j<", pentru"convenabil ales.
Pentru stabilirea lui "^ n acest algoritm trebuie s a t inem cont c a precizia a steptat a
la ecare pas este limitat a de eroarea formulei de baz a (1.19)  si de dimensiunea pa-
suluih. Pentru a adapta acest algoritm la solut ia unei probleme specice, trebuie s a
ment ion am:
10

Florina-Adriana Bosie Metode numerice pentru probleme cu condit ie init ial a
(a)num arulNde pa si dorit i;
(b)num arul maxim Kde iterat ii interne;
(c)procedura care trebuie urmat a ^ n cazul ^ n care kdep a se steK.
^In mod uzual, numim o formul a explicit a, cum ar  formula lui Euler, formul a de
tip deschis, ^ n timp ce o formul a implicit a, cum ar  (1.19) este numit a formul a de tip
^ nchis.
C^ and sunt folosite ca o pereche de formule, formula de tip deschis este numit a
predictor , iar formula de tip ^ nchis este numit a corector . Formula corector are, ^ n
general, o precizie mai mare dec^ at formula predictor, chiar  si ^ n cazul ^ n care ambele au
eroarea de discretizare de acela si ordin, ^ n principal deoarece coecientul din termenul
erorii este mai mic.
Rescriem metoda predictor-corector utiliz^ and aceste noi notat ii. Pasul predictor ^ l
obt inem cu ajutorul formulei lui Euler, iar pasul corector este dat de metoda Runge-
Kutta de ordinul II (care este mai performant a, aceasta urm^ and s a e prezentat a ^ n
capitolul urm ator), dup a cum urmeaz a:
pas predictor:
yp
n+1=yn+hf(xn;yn);
pas corector:
yc
n+1=yn+1
2(k0+k1)
=yn+1
2[hf(xn;yn) +hf(xn+h;yn+k0)]
=yn+h
2[f(xn;yn) +f(xn+1;yp
n+1)]:
Dou a ^ ntreb ari apar, ^ n mod natural, ^ n privint a formulelor corector. Prima este
,,^In ce condit ii va converge iterat ia intern a pe k?\  si a doua "De c^ ate iterat ii va  nevoie
pentru a produce precizia cerut a?\. R aspunsul la cea de-a doua ^ ntrebare depinde de
mult i factori. Cu toate acestea, dac a formulele predictor  si corector sunt de acela si
ordin, experient a a ar atat c a doar una sau dou a aplic ari ale corectorului sunt suciente,
cu condit ia ca m arimea pasului hs a e aleas a corect. Dac a consider am c a una sau
dou a corect ii nu sunt suciente este de preferat s a reducem dimensiunea pasului h
dec^ at s a continu am iterarea. R aspunsul primei ^ ntreb ari este cont inut ^ n urm atoarea
teorem a.
Teorema 1.4. (Conte & de Boor [1]) Dac a f(x;y) si@f
@ysunt continue ^ n x siype
intervalul ^ nchis [a;b], iterat ia intern a denit a de (1.21) va converge, cu condit ia ca h
11

Florina-Adriana Bosie Metode numerice pentru probleme cu condit ie init ial a
s a e ales sucient de mic pentru a face acest lucru posibil pentru x=xn si oriceycu
jyyn+1jjy(0)
n+1yn+1j,jyyn+1jjyp
n+1yn+1j;
@f
@yh<2:
Demonstrat ie. Pentru a demonstra acest lucru, mai ^ nt^ ai observ am c a ^ n iterat ia
(1.21),xneste xat. Prin urmare, dac a not am y(k)
n+1cuY(k), putem scrie (1.21) sub
forma
Y(k)=F(Y(k1));
unde
F(Y) =h
2f(xn+1;Y) +C
 si undeCdepinde de n, dar nu depinde de Y. O astfel de iterat ie va converge dac a
F0(Y) este continu a  si satisface jF0(Y)j<1;oricare ar  YcujYyn+1jjY(0)yn+1j,
undeyn+1este un punct xat al lui F(Y).
DeoareceF0(Y) =h
2@f
@y si@f
@yeste m arginit a, iterat ia (1.21) va converge dac a
jF0(Y)j=h
2@f
@y<1;adic ah<2@f
@y:DeoareceF0(Y) =h
2@f
@y, teorema este demon-
strat a.

1.5 Metoda coecient ilor nedeterminat i
Metoda coecient ilor nedeterminat i este o metod a analitic a des utilizat a ^ n studiul
ecut iilor diferent iale.
Pentru expunerea acestei metode vom considera problema Cauchy:
8
<
:y00+f(x)y0+g(x)y=h(x)
y(x0) =y0
y0(x0) =y0
0:(1.22)
^In continuare, presupunem c a sunt ^ ndeplinite condit iile din teorema de existent  a
 si unicitate  si c a ecare coecient al ecuat iei considerate se poate dezvolta ^ n serie de
12

Florina-Adriana Bosie Metode numerice pentru probleme cu condit ie init ial a
puteri ale lui xx0, adic a:
f(x) =1X
n=0an(xx0)n;
g(x) =1X
n=0bn(xx0)n;
h(x) =1X
n=0cn(xx0)n:
C aut am solut ia ecuat iei (1.22) sub forma:
y(x) =1X
n=0kn(xx0)n:
Pentru a determina solut ia problemei Cauchy luat a ^ n considerare, trebuie s a identi-
c am valorile coecient ilor kn; n0.
Consider am x0= 0, iar ecuat ia liniar a devine:
8
<
:y00+f(x)y0+g(x)y=h(x)
y(0) =y0
y0(0) =y0
0:(1.23)
Solut ia ecuat iei (1.23) se caut a, deci, sub forma:
f(x) =1X
n=0anxn;
g(x) =1X
n=0bnxn;
h(x) =1X
n=0cnxn:
Diny(0) =y0)k0=y0.
Deriv^ andy(x) =1X
n=0knxn, obt inem:
y0(x) =1X
n=1nknxn1=1X
n=0(n+ 1)kn+1xn:
13

Florina-Adriana Bosie Metode numerice pentru probleme cu condit ie init ial a
Diny0(x) =y0
0)k1=y0
0.
Coecient ii k0 sik1s-au determinat din condit iile init iale.
Dac a deriv am y(x) de dou a ori, rezult a:
y00(x) =1X
n=1n(n+ 1)kn+1xn1
=1X
n=0(n+ 1)(n+ 2)kn+2xn:
^Inlocuind ^ n (1.23), avem:
8
>>>>>>>>><
>>>>>>>>>:1X
n=0(n+ 1)(n+ 2)kn+2xn+1X
n=0anxn
1X
n=0(n+ 1)kn+1xn
+1X
n=0bnxn
1X
n=0knxn=1X
n=0cnxn
k0=y0;
k1=y0
0;
din care, prin identicare, a
 am valorile coecient ilor kn. Prin identicarea coecient ilor,
obt inem: 8
>>><
>>>:2k2+a0k1+b0k0=c0
3
2k3+ 2a0k2+a1k1+b0k1+b1k0=c1
…
(n+ 2)
(n+ 1)kn+2+L(k0;k1;:::;kn+1) =cn;
undeLeste o funct ie liniar a ^ n k0;k1;:::;kn+1.
^In cazul ^ n care funct iile f; g sihdin (1.23) sunt periodice de perioad a 2 , atunci
putem c auta solut ia sub forma unei serii trigonometrice ( a se vedea Pitea & Postolache
[2]):
y(x) =1X
n=0(pncosnx+qnsinnx):
Coecient ii se determin a  si ^ n acest caz, prin identicare.
Dezavantajul major al acestei metode const a ^ n faptul c a se lucreaz a cu un num ar
mare de necunoscute.
14

Florina-Adriana Bosie Metode numerice pentru probleme cu condit ie init ial a
1.6 Aplicat ii
I. Consider am problema Cauchy:
y0=yx2+ 1
y(0) = 0:5:
1) S a se calculeze primii 7 termeni din  sirul aproximat iilor succesive.
Solut ie. Pentru rezolvarea acestei cerint e, vom folosi formula de recurent  a dat a de
metoda aproximat iilor succesive:
yn+1(x) =y0+Zx
x0f(t;yn(t))dt; n = 0;1;2;::: :
^In cazul nostru, metoda lui Picard conduce la recurent a:
yn+1(x) = 0:5 +Zx
0(yn(t)t2+ 1)dt; n = 0;1;2;::: :
Dup a calcule, obt inem:
y1(x) =1
3×3+3
2x+1
2;
y2(x) =1
12×41
3×3+3
4×2+3
2x+1
2;
y3(x) =1
60×51
12×41
12×3+3
4×2+3
2x+1
2;
y4(x) =1
360×61
60×51
48×41
12×3+3
4×2+3
2x+1
2;
y5(x) =1
2520×71
360×61
240×51
48×41
12×3+3
4×2+3
2x+1
2;
y6(x) =1
20160×81
2520×71
1440×61
240×51
48×41
12×3+3
4×2+3
2x+1
2;
y7(x) =1
181440×91
20160×81
10080×71
1440×61
240×51
48×41
12×3
+3
4×2+3
2x+1
2:
Dac a x am x= 1, g asim:
y1(1) = 1:6666667; y 5(1) = 2:6384921;
y2(1) = 2:3333334; y 6(1) = 2:6405259;
y3(1) = 2:5666667; y 7(1) = 2:6408180:
y4(1) = 2:6263889;
15

Florina-Adriana Bosie Metode numerice pentru probleme cu condit ie init ial a
Programul Maple pentru rezolvarea problemei Cauchy considerat a, cu metoda lui
Picard, este urm atorul:
restart;
Digits := 8;
ec := diff(y(x), x) = y(x)-x^2+1;
CI := y(0) = 1/2;
dsolve({ec, CI});
sol[e] := unapply(%, x);
evalf(sol[e](1));
a := 1;
f := (x, y) -> y-x^2+1;
x[0] := 0;
y[0] := 1/2;
N := 7;
for n from 0 to N-1 do
y[n+1] := unapply(y[0]+int(f(t, y[n](t)), t = x[0] .. x), x);
print(nprintf("Iteratia %d", n+1));
print(y[n+1](x));
print(nprintf("Valoarea aproximativa in %d este", a));
print(y[n+1](1.));
end do:
2) S a se apoximeze solut ia problemei considerate printr-un polinom Taylor de grad 5.
Solut ie. Vom c auta solut ia sub forma unei serii de puteri, astfel:
y(x)'y(x0) +y0(x0)
1!(xx0) +y00(x0)
2!(xx0)2+y000(x0)
3!(xx0)3+
+y(4)(x0)
4!(xx0)4+y(5)(x0)
5!(xx0)5:
^In cazul nostru, x0= 0, deci solut ia va  c autat a sub forma unei serii de puteri
centrat a ^ n origine:
y(x)'y(0) +y0(0)
1!x+y00(0)
2!x2+y000(0)
3!x3+y(4)(0)
4!x4+y(5)(0)
5!x5:
16

Florina-Adriana Bosie Metode numerice pentru probleme cu condit ie init ial a
Din condit ia init ial a, avem y(0) = 0:5 =1
2 siy0(x) =y(x)x2+ 1, deci
y0(0) =y(0)0 + 1 = 1:5 =3
2.
^In continuare, pentru a
area solut iei, vom deriva relat ia y0=yx2+ 1 de c^ ate
ori este nevoie. Deoarece, ^ n enunt , se cere a
area solut iei problemei Cauchy prin
aproximarea cu un polinom Taylor de grad 5, vom aplica derivata de ^ nc a 4 ori.
y00(x) =y0(x)2x)y00(0) =y0(0)2
0)y00(0) = 1:5 =3
2;
y000(x) =y00(x)2)y000(0) =y00(0)2)y000(0) =0:5 =1
2;
y(4)(x) =y000(x))y(4)(x) =y000(0))y(4)(0) =0:5 =1
2;
y(5)(x) =y(4)(x))y(5)(0) =y(4)(0))y(5)(0) =0:5 =1
2:
Deci, aproximarea solut iei este:
y(x)'1
2+3
2x+3
2
2!x21
2
3!x31
2
4!x41
2
5!x5;adic a
y(x)'1
2+3
2x+3
4×21
12×31
48×41
240×5:
Observat ie. Aproximarea solut iei unei probleme Cauchy cu un polinom Taylor
ofer a acela si rezultat ca ^ n cazul aproxim arii solut iei cu ajutorul metodei coecient ilor
nedeterminat i.
^In urm atorul tabel, vom xa x= 1  si vom calcula at^ at solut ia exact a, c^ at  si solut iile
aproximative, determinate cu ajutorul metodei aproximat iilor succesive  si metodei lui
Taylor, ^ n acest punct. Apoi, pe baza acestor rezultate, vom calcula eroarea ce apare
datorit a aproxim arii solut iei exacte prin aceste 2 metode. Vom lucra cu 4 zecimale
exacte.
y(1) Eroarea
y(x) = 1 + 2 x+x21=2ex2.6408 0
1=2 + 3=2x1=3×31.6666 0.9742
1=2 + 3=2x+ 3=4×21=12×41=3×32.3333 0.3075
1=2 + 3=2x+ 3=4×21=12×31=60×51=12×42.5666 0.0742
1=2 + 3=2x+ 3=4×21=12×31=48×41=60×51=360×62.6263 0.0145
1=2 + 3=2x+ 3=4×21=12×31=48×41=240×51=360×61=2520×72.6384 0.0024
1=2 + 3=2x+ 3=4×21=12×31=48×41=240×52.6416 0.0008
17

Florina-Adriana Bosie Metode numerice pentru probleme cu condit ie init ial a
Gracul solut iei exacte  si al solut iilor aproximative, determinate cu ajutorul meto-
delor Picard  si Taylor este prezentat ^ n Figura 1.1 de mai jos:
Figura 1.1: Solut ia exact a  si solut iile aproximative
II. Fie problema Cauchy:y0=xe3x2y
y(0) = 0:
S a se determine solut ia problemei propuse, folosind metoda lui Euler pe intervalul
[0;1], cu pasulh= 0:2.
Solut ie. Pentru rezolvarea acestei probleme, vom utiliza relat ia:
yn+1=yn+hf(xn;yn); n =0;n1:
^In cazul nostru, obt inem urm atoarea relat ie de recurent  a:
yn+1=yn+ 0:2(xne3xn2yn):
18

Florina-Adriana Bosie Metode numerice pentru probleme cu condit ie init ial a
Dup a calcule, obt inem urm atoarele rezultate:
Pentrun= 0, avem:
x0= 0;
y1=y0+hf(x0;y0) = 0;
Pentrun= 1, avem:
x1= 0:2;
y2=y1+hf(x1;y1) = 0:07288475200;
Pentrun= 2, avem:
x2= 0:4;
y3=y2+hf(x2;y2) = 0:3093402050;
Pentrun= 3, avem:
x3= 0:6;
y4=y3+hf(x3;y3) = 0:9115618186;
Pentrun= 4, avem:
x4= 0:8;
y5=y4+hf(x4;y4) = 2:310645312:
Urm atorul tabel cont ine rezultatele obt inute prin alegerea pasului h= 0:1.
nxny(xn)yn Eroarea
0 0 0 0 0
10.1 0.005752 0 0.005752053967
20.2 0.026813 0.01349858808 0.01331421376
30.3 0.071145 0.04724124646 0.02390328120
40.4 0.150778 0.1115810905 0.0391967450
50.5 0.283617 0.2220695493 0.0615469726
60.6 0.496020 0.4017400929 0.0942794728
70.7 0.826481 0.6843709221 0.1421099476
80.8 1.330857 1.119128632 0.211728394
90.9 2.089774 1.777157016 0.312617380
10 13.219099 2.760901468 0.458197851
Programul Maple care rezolv a aceast a problem a este:
restart;
f := (x, y) -> x*exp(3*x)-2*y;
19

Florina-Adriana Bosie Metode numerice pentru probleme cu condit ie init ial a
ec := diff(y(x), x) = f(x, y(x));
CI := y(0) = 0;
dsolve({ec, CI});
Sol := x -> ((-1/25+(1/5)*x)*exp(5*x)+1/25)*exp(-2*x);
a := 0;
b := 1;
h := 0.2;
N := (b-a)/h;
x[0] := 0;
y[0] := 0;
f := (x, y) -> x*exp(3*x)-2*y;
for n from 0 to N-1 do
x[n+1] := x[n]+h;
y[n+1] := y[n]+h*f(x[n], y[n]);
print(nprintf("Solutia exacta y(x[%d])=%f", n+1, Sol(x[n+1])));
print(nprintf("Aprox %d", n+1));
print(y[n+1]); print(nprintf("Eroarea in punctul x[%d]", n+1));
print(abs(Sol(x[n+1])-y[n+1]));
end do;
with(plots):
p1 := plot([seq([x[k], y[k]], k = 0 .. N)], color = blue):
p2 := plot(Sol(x), x = a .. b, color = red):
display(p2, p1);
^In urm atorul grac sunt reprezentate solut ia exact a  si solut ia obt inut a cu metoda
lui Euler, pentru problema Cauchy considerat a.
Figura 1.2: Euler h=0.2
20

Florina-Adriana Bosie Metode numerice pentru probleme cu condit ie init ial a
Pentru a ilustra apropierea solut iei aproximative de solut ia exact a pe m asura
mic sor arii pasului de discretizare, vom reprezenta solut iile obt inute cu metoda lui
Euler pentru diverse m arimi ale pasului.
Figura 1.3: Euler h=0.15
Figura 1.4: Euler h=0.1
Figura 1.5: Euler h=0.03
21

Capitolul 2
Metode de tip Runge-Kutta
Dup a cum am ment ionat anterior, metoda lui Euler nu este foarte util a ^ n probleme
practice pentru c a necesit a o dimensiune foarte mic a a pasului de discretizare ^ n scopul
obt inerii unei precizii rezonabile. Algoritmul lui Taylor pentru ordin superior nu este
acceptat ca  si procedur a cu scop general din cauza necesit at ii de a obt ine derivate
totale de ordin superior ale lui y(x).
Metodele de tip Runge-Kutta ^ ncearc a s a obt in a o precizie mai mare  si ,^ n acela si
timp, s a evite necesitatea calculului derivatelor de ordin superior, evalu^ and funct ia
f(x;y) ^ ntr-un anumit punct pe ecare subinterval.
Consider am problema Cauchy
y0=f(x;y); y(x0) =y0: (2.1)
Fiepun num ar natural  si x0<x 1<


<xm=b; xn=x0+nh; n =1;m, o diviziune
echidistant a a intervalului [ x0;b].
Din formula lui Taylor, avem formula de recurent  a
yn+1=yn+hf(xn;yn) +


+hp+1
(p+ 1)!f(p)(xn;yn); (2.2)
undey0este dat.
Metodele de tip Runge-Kutta presupun g asirea unor formule de recurent  a de tipul
yn+1=yn+h[0f(xn;yn) +1f(xn+1h;yn+1h) +2f(xn+2h;yn+2h)
+


+pf(xn+ph;yn+ph)];
cui; i; i; i=1;pdeterminate din condit ia ca dezvoltarea^ n serie Taylor a membru-
lui drept, ^ n funct ie de h, s a coincid a cu membrul drept al formulei Taylor de ordinul
(p+ 1) dat a de relat ia (2.2).
22

Florina-Adriana Bosie Metode numerice pentru probleme cu condit ie init ial a
Pentru a simplica calculele, vom c auta o formul a iterativ a de tipul:
yn+1=yn+0k0+1k1+2k2+


+pkp;
undek0; k1; k2;:::;kpse determin a sub forma
k0=hf(xn;yn);
k1=hf(xn+1h;yn+10k0);
k2=hf(xn+2h;yn+20k0+21k1);
…
kp=hf(xn+ph;yn+p0k0+p1k1+


+p;p1kp1);
iary0este dat.
^In funct ie de valorile lui pse dezvolt a diverse formule de tipul Runge-Kutta.
Vom ^ ncepe cu cea mai simpl a metod a de tip Runge-Kutta.
2.1 Metoda Runge-Kutta de ordinul I
Din formula lui Taylor, avem
yn+1=yn+hf(xn;yn);
iar formula Runge-Kutta devine
yn+1=yn+0hf(xn;yn);
y0dat:
Prin identicare ^ ntre aceste dou a formule, g asim 0= 1. Deci formula Runge-
Kutta de ordinul I este: 8
<
:yn+1=yn+k0
k0=hf(xn;yn)
y0dat;
adic a formula lui Euler.
Din capitolul 1,  stim c a eroarea, ^ n cazul metodei lui Euler, este de ordinul h2, ceea
ce ^ nseamn a c a  si ^ n cazul formulei Runge-Kutta de ordinul I, eroarea este de ordinul
h2.
23

Florina-Adriana Bosie Metode numerice pentru probleme cu condit ie init ial a
2.2 Metoda Runge-Kutta de ordinul II
Pentru a determina formula Runge-Kutta de ordinul II, c aut am o expresie de forma:
yn+1=yn+0k0+1k1; (2.3)
unde
k0=hf(xn;yn);
k1=hf(xn+h;yn+k0);
unde0;  1; ;  sunt constante care trebuie determinate astfel ^ nc^ at (2.3) s a se iden-
tice cu algoritmul lui Taylor de ordin superior. Dezvolt^ and y(xn+1) ^ n serie Taylor
p^ an a la termeni de ordinul h3, obt inem:
y(xn+1) =y(xn) +hy0(xn) +h2
2y00(xn) +h3
6y000(xn) +



=y(xn) +hf(xn;yn) +h2
2(fx+ffy)n (2.4)
+h3
6(fxx+ 2ffxy+fyyf2+fxfy+f2
yf)n+O(h4);
unde am folosit dezvoltarea (1.16)  si indicele n^ nseamn a c a toate funct iile implicate
vor  evaluate ^ n ( xn;yn).
Pe de alt a parte, folosind dezvoltarea ^ n serie Taylor pentru funct ii de dou a varia-
bile, g asim
k1
h=f(xn+h;yn+k0) =f(xn;yn) +hfx+k0fy
+2h2
2fxx+hk 0fxy+2k0
2fyy+O(h3);
unde toate derivatele sunt evaluate ^ n ( xn;yn).
Acum ^ nlocuim aceast a dezvoltare pentru k1^ n (2.3)  si t in^ and cont de faptul c a
k0=hf(xn;yn), g asim, dup a rearanjarea ^ n funct ie de puterile lui h, c a
yn+1=yn+ (0+1)hf+1h2(fx+ffy)
+1h32
2fxx+ffxy+2
2f2fyy
+O(h4): (2.5)
Compar^ and cu (2.4), observ am c a pentru a face puterile h sih2s a se potriveasc a,
trebuie s a avem
0+1= 1; (2.6)
1=1=1
2: (2.7)
24

Florina-Adriana Bosie Metode numerice pentru probleme cu condit ie init ial a
Cu toate c a avem patru necunoscute, avem doar trei ecuat ii, prin urmare, ^ nc a avem
un grad de libertate pentru solut ia lui (2.6). Putem spera la o utilizare a acestui grad
de libertate adit ional pentru a obt ine o potrivire a coecient ilor termenilor h3. Este
evident c a acest lucru nu este posibil pentru toate funct iile f(x;y). Sunt multe solut ii
pentru (2.6), cea mai simpl a ind
0=1=1
2;  == 1:
Am obt inut, astfel, formula de tip Runge-Kutta de ordinul II
8
>>><
>>>:yn+1=yn+1
2(k0+k1)
k0=hf(xn;yn)
k1=hf(xn+h;yn+k0)
y0dat:(2.8)
Formula Runge-Kutta de ordinul II mai poate  scris a astfel:
(
yn+1=yn+h
2[f(xn;yn) +f(xn+h;yn+hf(xn;yn))]
y0dat:
Eroarea local a pentru (2.8) este de forma
y(xn+1)yn+1=h3
12(fxx+ 2ffxy+f2fyy2fxfy2ff2
y) +O(h4):
Complexitatea coecientului din formula acestei erori este specic a tuturor meto-
delor de tip Runge-Kutta  si reprezint a una dintre cele mai indezirabile caracteristici
ale unor astfel de metode deoarece estim arile erorii locale sunt foarte greu de obt inut.
Eroarea local a a (2.8) este, totu si, de ordinul h3, pe c^ and cea obt inut a prin metoda
lui Euler este de ordinul h2. Prin urmare, ne a stept am s a putem folosi o m arime a
pasului mai mare cu (2.8). Inconvenientul este c a trebuie s a evalu am funct ia f(x;y)
de dou a ori pentru ecare pas al integr arii. Formulele de tip Runge-Kutta de orice
ordin pot  obt inute prin metoda folosit a mai sus, ^ ns a driv arile devin din ce ^ n ce mai
complicate.
2.3 Metoda Runge-Kutta de ordinul III
Pentru a determina formula Runge-Kutta de ordinul III proced am similar cazului ^ n
care am determinat formula Runge-Kutta de ordinul II. C aut am o formul a de tipul:
yn+1=yn+0k0+1k1+2k2; (2.9)
25

Florina-Adriana Bosie Metode numerice pentru probleme cu condit ie init ial a
^ n care
k0=hf(xn;yn);
k1=hf(xn+1h;yn+10k0);
k2=hf(xn+2h;yn+20k0+21k1);
undey0este dat, iar 0;  1; 1; 2; 10; 20; 21sunt constante ce trebuie determinate
prin identicarea cu dezvoltarea ^ n serie Taylor. ^In urma acestei identic ari, obt inem
formula Runge-Kutta de ordinul III:
8
>>>>>>><
>>>>>>>:yn+1=yn+1
6(k0+ 4k1+k2)
k0=hf(xn;yn)
k1=hf(xn+h
2;yn+k0
2)
k2=hf(xn+h;yn+ 2k1k0)
y0dat:
^In acest caz, eroarea este de ordinul h4.
2.4 Metoda Runge-Kutta de ordinul IV
Cea mai cunoscut a  si des folosit a formul a de tip Runge-Kutta este cea de ordinul IV,
pe care o vom prezenta ^ n cele ce urmeaz a.
Pentru ecuat ia
y0=f(x;y); y(x0) =y0
gener am aproxim arile ynpentruy(x0+nh), cuhxat  sin= 0;1;2;:::, folosind
formula recursiv a
yn+1=yn+1
6(k0+ 2k1+ 2k2+k3); (2.10)
unde
k0=hf(xn;yn);
k1=hf(xn+h
2;yn+k0
2);
k2=hf(xn+h
2;yn+k1
2);
k3=hf(xn+h;yn+k2):
Eroarea local a de discretizare a acestui algoritm este de ordinul h5. La fel ca  si
^ n cazul celorlalte metode de tip Runge-Kutta, inconvenientul metodei este c a sunt
26

Florina-Adriana Bosie Metode numerice pentru probleme cu condit ie init ial a
necesare evalu arile mai multor funct ii (^ n acest caz, 4) la ecare iterat ie pentru a
obt ine o eroare de discretizare favorabil a.
Formula (2.10) este folosit a pe scar a larg a ^ n practic a, cu un succes considerabil.
Are avantajul important de a  o metod a cu autostart, adic a este nevoie doar de
valoarea lui y^ ntr-un punct x=xnpentru a a
a y siy0^ n punctul x=xn+1.
^In cele ce urmeaz a, vom expune formulele de tip Runge-Kutta pentru sisteme.
Metoda Runge-Kutta de ordinul II pentru sisteme
Consider am sistemul de ecuat ii diferent iale
8
>><
>>:y0=f(x;y;z )
z0=g(x;y;z )
y(x0) =y0
z(x0) =z0:
Solut ia acestui sistem, obt inut a cu ajutorul metodei Runge-Kutta de ordinul II
este de forma
yn+1=yn+1
2(k0+k1); z n+1=zn+1
2(l0+l1);
k0=hf(xn;yn;zn); l 0=hg(xn;yn;zn);
k1=hf(xn+h;yn+k0;zn+l0); l 1=hg(xn+h;yn+k0;zn+l0):
Metoda Runge-Kutta de ordinul III pentru sisteme
Fie dat sistemul:
8
>><
>>:y0=f(x;y;z )
z0=g(x;y;z )
y(x0) =y0
z(x0) =z0
Solut ia acestui sistem se determin a astfel:
27

Florina-Adriana Bosie Metode numerice pentru probleme cu condit ie init ial a
yn+1=yn+1
6(k0+ 4k1+k2); z n+1=zn+1
6(l0+ 4l1+l2);
k0=hf(xn;yn;zn); l 0=hg(xn;yn;zn);
k1=hf(xn+h
2;yn+k0
2;zn+l0
2); l 1=hg(xn+h
2;yn+k0
2;zn+l0
2);
k2=hf(xn+h;yn+ 2k1k0;zn+ 2l1l0); l 2=hg(xn+h;yn+ 2k1k0;zn+ 2l1l0):
Metoda Runge-Kutta de ordinul IV pentru sisteme
Pentru sistemul
8
>><
>>:y0=f(x;y;z )
z0=g(x;y;z )
y(x0) =y0
z(x0) =z0;
c aut am solut ia sub forma:
yn+1=yn+1
6(k0+ 2k1+ 2k2+k3); zn+1=zn+1
6(l0+ 2l1= 2l2+l3);
k0=hf(xn;yn;zn); l 0=hg(xn;yn;zn);
k1=hf(xn+h
2;yn+k0
2;zn+l0
2); l 1=hg(xn+h
2;yn+k0
2;zn+l0
2);
k2=hf(xn+h
2;yn+k1
2;zn+l1
2); l 2=hg(xn+h
2;yn+k1
2;zn+l1
2);
k3=hf(xn+h;yn+k2;zn+l2); l 3=hg(xn+h;yn+k2;zn+l2):
2.5 Aplicat ii
I. Folosind formula Runge-Kutta de ordinul II,  stiind c a intervalul de interes este [0 ;1],
iar pasul folosit este 0.2, s a se rezolve urm atoarea problem a Cauchy:
y0=xe3x2y
y(0) = 0:
28

Florina-Adriana Bosie Metode numerice pentru probleme cu condit ie init ial a
Solut ie. Pentru rezolvarea acestei probleme, vom c auta solut ia sub forma:
yn+1=yn+k0+k1
2
y0dat;
unde
k0=hf(xn;yn);
k1=hf(xn+h;yn+k0):
Din datele problemei, avem: x0= 0,y0= 0,h= 0:2)x1= 0:2;x2= 0:4,
x3= 0:6;x4= 0:8;x5= 1,f(x;y) =xe3x2y.
Pentrun= 0, obt inem:
x0= 0;
k0= 0;
k1= 0:072884;
y1= 0:036442:
Pentrun= 1, obt inem:
x1= 0:2;
k0= 0:058308;
k1= 0:22770;
y2= 0:17945:
Pentrun= 2, obt inem:
x2= 0:4;
k0= 0:19382;
k1= 0:57666;
y3= 0:56469:
Pentrun= 3, obt inem:
x3= 0:6;
k0= 0:50008;
k1= 1:3378;
y4= 1:4836:
Pentrun= 4, obt inem:
x4= 0:8;
k0= 1:1702;
k1= 2:9556;
y5= 3:5465:
Dac a,^ n acest a problem a, alegem pasul de discretizare h= 0:1, obt inem urm atoarele
rezultate:
29

Florina-Adriana Bosie Metode numerice pentru probleme cu condit ie init ial a
nxny(xn)yn Eroarea
0 0 0 0 0
10.1 0.005752 0.0067495 0.0009971
20.2 0.026813 0.029155 0.002342
30.3 0.071142 0.075379 0.004237
40.4 0.150780 0.15773 0.00695
50.5 0.283610 0.29450 0.01089
60.6 0.496030 0.51261 0.01658
70.7 0.826480 0.85135 0.02487
80.8 1.330800 1.3676 0.0368
90.9 2.089700 2.1438 0.0541
10 13.219200 3.2979 0.0787
Programul Maple, utilizat pentru rezolvarea acestei probleme cu ajutorul metodei
Runge-Kutta de ordinul II, este urm atorul:
restart;
f := (x, y) -> x*exp(3*x)-2*y;
ec := diff(y(x), x) = f(x, y(x));
CI := y(0) = 0;
dsolve({ec, CI});
Sol := x -> ((1/25)*(5*x-1)*exp(5*x)+1/25)*exp(-2*x);
Digits := 5;
f := (x, y) -> x*exp(3*x)-2*y;
x[0] := 0;
y[0] := 0;
a := 0;
b := 1;
h := 0.2;
N := (b-a)/h;
for n from 0 to N-1 do
print(nprintf("Iteratia %d", n+1));
x[n+1] := x[n]+h;
k[0] := h*f(x[n], y[n]);
k[1] := h*f(x[n]+h, y[n]+k[0]);
y[n+1] := y[n]+(k[0]+k[1])*(1/2);
print(nprintf("Solutia exacta y(x[%d])=%f", n+1, Sol(x[n+1])));
30

Florina-Adriana Bosie Metode numerice pentru probleme cu condit ie init ial a
print(nprintf("Eroarea in punctul x[%d]=%g", n+1, x[n+1]));
print(abs(Sol(x[n+1])-y[n+1]));
end do;
with(plots);
p1 := plot([seq([x[k], y[k]], k = 0 .. N)], color = blue):
p2 := plot(Sol(x), x = a .. b, color = red):
display(p2, p1);
^In urm atoarele grace, vom reprezenta solut ia exact a  si solut ia obt inut a aplic^ and
metoda Runge-Kutta de ordinul II, pentru diverse valori ale pasului de discretizare.
Figura 2.1: Runge-Kutta II h=0.5
Figura 2.2: Runge-Kutta II h=0.2
31

Florina-Adriana Bosie Metode numerice pentru probleme cu condit ie init ial a
Figura 2.3: Runge-Kutta II h=0.1
II. Se consider a urm atoarea problem a Cauchy:
y0= 1 + (xy)2; x2[2;3]
y(2) = 1:
S a se rezolve aceast a problem a cu metodele Runge-Kutta de ordinul III  si IV, con-
sider^ and pasul de discretizare h= 0:2.
Solut ie.
(i) Vom ^ ncepe rezolvarea acestei probleme cu aplicarea metodei Runge-Kutta de
ordinul III.
Formula utilizat a este:
yn+1=yn+k0+ 4k1+k2
6
y0dat;
unde
k0=hf(xn;yn);
k1=hf(xn+h
2;yn+k0
2);
k2=hf(xn+h;yn+ 2k1k0):
Datele problemei sunt: x0= 2,y0= 1,h= 0:2)x1= 2:2;x2= 2:4;x3= 2:6,
x4= 2:8;x5= 3,f(x;y) = 1 + (xy)2.
La prima iterat ie obt inem:
x0= 2;
k0= 0:4;
k1= 0:36200;
k2= 0:35348;
y1= 1:3669:
32

Florina-Adriana Bosie Metode numerice pentru probleme cu condit ie init ial a
La a doua iterat ie obt inem:
x1= 2:2;
k0= 0:33882;
k1= 0:31664;
k2= 0:30910;
y2= 1:6860:
La a treia iterat ie obt inem:
x2= 2:4;
k0= 0:30196;
k1= 0:28792;
k2= 0:28198;
y3= 1:9752:
La a patra iterat ie obt inem:
x3= 2:6;
k0= 0:27808;
k1= 0:26864;
k2= 0:26398;
y4= 2:2446:
La ultima iterat ie obt inem:
x4= 2:8;
k0= 0:26170;
k1= 0:25504;
k2= 0:25140;
y5= 2:5001:
^In urm atorul tabel, avem rezultatele obt inute pentru h= 0:1:
nxny(xn)yn Eroarea
0 2 1 0
12.1 1.190900 1.1909 0
22.2 1.366700 1.3667 0
32.3 1.530800 1.5308 0
42.4 1.685700 1.6857 0
52.5 1.833300 1.8334 0.0001
62.6 1.975000 1.9750 0
72.7 2.111800 2.1118 0
82.8 2.244400 2.2444 0
92.9 2.373700 2.3736 0.0001
10 32.500000 2.5000 0
33

Florina-Adriana Bosie Metode numerice pentru probleme cu condit ie init ial a
Programul Maple pentru rezolvarea problemei, folosind metoda Runge-Kutta de
ordinul III, este:
restart;
f := (x, y) -> 1+(x-y)^2;
ec := diff(y(x), x) = f(x, y(x));
CI := y(2) = 1;
dsolve({ec, CI});
Sol := x -> (x^2-x-1)/(x-1);
Digits := 5;
f := (x, y) -> 1+(x-y)^2;
x[0] := 2;
y[0] := 1;
a := 2;
b := 3;
h := 0.2;
N := (b-a)/h;
for n from 0 to N-1 do
print(nprintf("Iteratia %d", n+1));
x[n+1] := x[n]+h;
k[0] := h*f(x[n], y[n]);
k[1] := h*f(x[n]+(1/2)*h, y[n]+(1/2)*k[0]);
k[2] := h*f(x[n]+h, y[n]+2*k[1]-k[0]);
y[n+1] := y[n]+(k[0]+4*k[1]+k[2])*(1/6);
print(nprintf("Solutia exacta y(x[%d])=%f", n+1, Sol(x[n+1])));
print(nprintf("Eroarea in punctul x[%d]=%g", n+1, x[n+1]));
print(abs(Sol(x[n+1])-y[n+1]));
end do;
with(plots);
p1 := plot([seq([x[k], y[k]], k = 0 .. N)], color = blue):
p2 := plot(Sol(x), x = a .. b, color = red):
display(p2, p1);
^In urm atorul grac este reprezentat a solut ia exact a  si solut ia obt inut a ^ n urma
aplic arii metodei Runge-Kutta de ordinul III, pentru problema Cauchy
y0= 1 + (xy)2; y(2) = 1.
34

Florina-Adriana Bosie Metode numerice pentru probleme cu condit ie init ial a
Figura 2.4: Solut ia exact a  si solut ia Runge-Kutta III h=0.5
Figura 2.5: Solut ia exact a  si solut ia Runge-Kutta III h=0.2
(ii)^In continuare, vom rezolva problema cu ajutorul metodei Runge-Kutta de or-
dinul IV.
Formula utilizat a este:
yn+1=yn+k0+ 2k1+ 2k2+k3
6
y0dat;
unde
k0=hf(xn;yn);
k1=hf(xn+h
2;yn+k0
2);
k2=hf(xn+h
2;yn+k1
2);
k3=hf(xn+h;yn+k2):
35

Florina-Adriana Bosie Metode numerice pentru probleme cu condit ie init ial a
Dup a calcule, obt inem urm atoarele rezultate:
Pentrun= 0:
x0= 2;
k0= 0:4;
k1= 0:36200;
k2= 0:36892;
k3= 0:33814;
y1= 1:3668:
Pentrun= 1:
x1= 2:2;
k0= 0:33884;
k1= 0:31668;
k2= 0:32010;
k3= 0:30170;
y2= 1:6859:
Pentrun= 2:
x2= 2:4;
k0= 0:30198;
k1= 0:28794;
k2= 0:28980;
k3= 0:27796;
y3= 1:9751:
Pentrun= 3:
x3= 2:6;
k0= 0:27810;
k1= 0:26864;
k2= 0:26976;
k3= 0:26162;
y4= 2:2444:
Pentrun= 4:
x4= 2:8;
k0= 0:26174;
k1= 0:25506;
k2= 0:25578;
k3= 0:24996;
y5= 2:5000:
Urm atorul tabel cont ine rezultatele obt inute pentru h= 0:1:
36

Florina-Adriana Bosie Metode numerice pentru probleme cu condit ie init ial a
nxny(xn)yn Eroarea
0 2 1 1 0
12.1 1.190900 1.1908 0.0001
22.2 1.366700 1.3664 0.0003
32.3 1.530800 1.5305 0.0003
42.4 1.685700 1.6855 0.0002
52.5 1.833300 1.8332 0.0001
62.6 1.975000 1.9749 0.0001
72.7 2.111800 2.1117 0.0001
82.8 2.244400 2.2443 0.0001
92.9 2.373700 2.3736 0.0001
10 32.500000 2.4999 0.0001
Programul Maple pentru rezolvarea problemei, folosind metoda Runge-Kutta de
ordinul IV, este:
restart;
f := (x, y) -> 1+(x-y)^2;
ec := diff(y(x), x) = f(x, y(x));
CI := y(2) = 1;
dsolve({ec, CI});
Sol := x -> (x^2-x-1)/(x-1);
Digits := 5;
f := (x, y) -> 1+(x-y)^2;
x[0] := 2;
y[0] := 1;
a := 2;
b := 3;
h := .2;
N := (b-a)/h;
for n from 0 to N-1 do
print(nprintf("Iteratia %d", n+1));
x[n+1] := x[n]+h;
k[0] := h*f(x[n], y[n]);
k[1] := h*f(x[n]+(1/2)*h, y[n]+(1/2)*k[0]);
k[2] := h*f(x[n]+(1/2)*h, y[n]+(1/2)*k[1]);
k[3] := h*f(x[n]+h, y[n]+k[2]);
37

Florina-Adriana Bosie Metode numerice pentru probleme cu condit ie init ial a
y[n+1] := y[n]+(k[0]+2*k[1]+2*k[2]+k[3])*(1/6);
print(nprintf("Solutia exacta y(x[%d])=%f", n+1, Sol(x[n+1])));
print(nprintf("Eroarea in punctul x[%d]=%g", n+1, x[n+1]));
print(abs(Sol(x[n+1])-y[n+1]));
end do;
with(plots);
p1 := plot([seq([x[k], y[k]], k = 0 .. N)], color = blue):
p2 := plot(Sol(x), x = a .. b, color = red):
display(p2, p1);
Gracul solut iei exacte  si al solut iei obt inut a cu metoda Runge-Kutta IV este:
Figura 2.6: Solut ia exact a  si solut ia Runge-Kutta IV h=0.5
Figura 2.7: Solut ia exact a  si solut ia Runge-Kutta IV h=0.2
38

Capitolul 3
Metode numerice multipas
Algoritmul lui Taylor de ordinul k si metodele de tip Runge-Kutta sunt ambele
exemple de metode unipas (directe). Acestea necesit a cunoa sterea solut iei ^ ntr-un sin-
gur punctx=xn, de la care metodele continu a pentru a
area lui y^ n punctul urm ator
x=xn+1. Metodele multipas (indirecte) folosesc informat ii despre toate valorile cu-
noscute p^ an a la punctul x=xn+1. Presupunem c a am obt inut deja aproxim arile lui
y0 siy^ n punctele x0; x1; x2;:::;xnale unei diviziuni echidistante. O clas a de metode
multipas se bazeaz a pe principiul integr arii numerice.
3.1 Metoda Adams-Bashforth
Consider am problema Cauchyy0=f(x;y)
y(x0) =y0:
Dac a integr am ecuat ia diferent ial a y0=f(x;y) de laxnlaxn+1, obt inem formula
\exact a": Zxn+1
xny0dx=Zxn+1
xnf(x;y(x))dx;
y(xn+1)y(xn) =Zxn+1
xnf(x;y(x))dx;
y(xn+1) =y(xn) +Zxn+1
xnf(x;y(x))dx (3.1)
sau
yn+1=yn+Zxn+1
xnf(x;y(x))dx: (3.2)
Pentru a efectua integrarea^ n formula (3.2), aproxim am funct ia f(x;y(x)) cu un po-
linom care interpoleaz a funct ia f(x;y(x)) ^ n (m+1) puncte, xn; xn1; xn2;:::;xnm.
39

Florina-Adriana Bosie Metode numerice pentru probleme cu condit ie init ial a
Dac a folosim notat ia
f(xk;y(xk)) =fk;
putem utiliza formula lui Newton de grad m, astfel
pm(x) =mX
k=0(1)ks
k
rkfn; s =xxn
h:
Introduc^ and aceast a formul a ^ n (3.2)  si not^ and dx=hds, obt inem
yn+1=yn+hZ1
0mX
k=0(1)ks
k
rkfnds
=yn+hmX
k=0(1)krkfnZ1
0s
k
ds
=yn+hmX
k=0(1)kZ1
0s
k
rkfnds: (3.3)
Introducem coecient ii auxiliari

k= (1)kZ1
0s
k
ds; (3.4)
iars
k
=s(s1)(s2)


(sk+ 1)
k!. Prin urmare, am obt inut formula
de tip Adams-Bashforth:
yn+1=yn+hmX
k=0
krkfn:
^In funct ie de valorile lui m, obt inem diferite clase de formule. ^In primul r^ and, vom
calcula primele c^ ateva valori ale coecient ilor
k.
40

Florina-Adriana Bosie Metode numerice pentru probleme cu condit ie init ial a

0= (1)0Z1
0s
0
ds= 1;

1= (1)1Z1
0s
1
ds=Z1
0s
1!ds=s2
21
0=1
2;

2= (1)2Z1
0s
2
ds=Z1
0s(s1)
2!ds=Z1
0s2+s
2ds=1
2s3
3+s2
21
0
=1
21
3+1
2
=5
12;

3= (1)3Z1
0s
3
ds=Z1
0s(s1)(s2)
3!ds=Z1
0s3+ 3s2+ 2s
6ds
=1
6s4
4+3s3
3+2s2
21
0=1
61
4+ 1 + 1
=3
8;

4= (1)4Z1
0s
4
ds=Z1
0s(s1)(s2)(s3)
4!ds
=Z1
0s4+ 6s3+ 11s2+ 6s
24ds=1
24s5
5+6s4
4+11s3
3+6s2
21
0=251
720:
Deci am obt inut primii 5 coecient i
k:

0= 1;
1=1
2;
2=5
12;
3=3
8;
4=251
720:
Formula (3.3) este cunoscut a ca  si metoda Adams-Bashforth. Cel mai simplu caz,
obt inut pentru m= 0 ^ n (3.3), conduce c atre metoda lui Euler.
yn+1=yn+h
0r0fn;
yn+1=yn+hfn;
yn+1=yn+hf(xn;yn):
^In general, utilizarea formulei (3.3) necesit a valorile lui y0=f^ n (m+ 1) puncte,
xn; xn1;:::;xnm. Pornind de la cunoa sterea acestora, putem forma diferent ele rkfn.
Din (3.3) putem calcula yn+1, iar din ecuat ia diferent ial a putem calcula
41

Florina-Adriana Bosie Metode numerice pentru probleme cu condit ie init ial a
fn+1=f(xn+1;yn+1). Acum renot am punctul xn+1cuxn, form am o nou a linie de
diferent e  si repet am procesul.
Pentrum= 1, formula de tip Adams-Bashforth devine:
yn+1=yn+h1X
k=0
krkfn
=yn+h(
0r0fn+
1rfn)
=yn+h
fn+1
2rfn
: (3.5)
Pentrum= 2, obt inem:
yn+1=yn+h2X
k=0
krkfn
=yn+h(
0r0fn+
1rfn+
2r2fn)
=yn+h
fn+1
2rfn+5
12r2fn
: (3.6)
Pentrum= 3, care este ^ n mod uzual folosit ^ n practic a, avem:
yn+1=yn+h3X
k=0
krkfn
=yn+h(
0r0fn+
1rfn+
2r2fn+
3r3fn)
=yn+h
fn+1
2rfn+5
12r2fn+3
8r3fn
: (3.7)
^In practic a este mai convenabil s a se lucreze cu coordonate dec^ at cu diferent e. Din
denit ia operatorului dieferent  a nit a la st^ anga, r, avem:
r0fn=fn;
r1fn=fnfn1;
r2fn=fn2fn1+fn2;
r3fn=fn3fn1+ 3fn2fn3:
^Inlocuind aceste valori ^ n (3.5), (3.6)  si (3.7)  si regrup^ and, obt inem:
42

Florina-Adriana Bosie Metode numerice pentru probleme cu condit ie init ial a
Pentrum= 1:
yn+1=yn+h
fn+1
2(fnfn1)
=yn+h3
2fn1
2fn1
=yn+h
2(3fnfn1):
Pentrum= 2:
yn+1=yn+h
fn+1
2(fnfn1) +5
12(fn2fn1+fn2)
=yn+h
fn+1
2fn1
2fn1+5
12fn10
12fn1+5
12fn2
=yn+h23
12fn16
12fn1+5
12fn2
=yn+h
12(23fn16fn1+ 5fn2):
Pentrum= 3:
yn+1=yn+h
fn+1
2(fnfn1) +5
12(fn2fn1+fn2)
+3
8(fn3fn1+ 3fn2fn3)
=yn+h
fn+1
2fn1
2fn1+5
12fn10
12fn1+5
12fn2+3
8fn
9
8fn1+9
8fn23
8fn3
=yn+h55
24fn59
24fn1+37
24fn29
24fn3
=yn+h
24(55fn59fn1+ 37fn29fn3):
Eroarea de trunchiere ^ n cazul metodei Adams-Bashforth se obt ine sc az^ and relat ia
(3.2) din (3.1), dup a cum urmeaz a ( Pitea & Postolache [2] ):
jy(xn+1)yn+1jjy(xn)ynj+Zxn+1
xnf(x;y(x))dxZxn+1
xnpm(x)dx:
43

Florina-Adriana Bosie Metode numerice pentru probleme cu condit ie init ial a
Observ am c a ^ n membrul drept avem suma erorilor obt inute la metodele directe (uni-
pas)  si la integrarea numeric a.
^In continuare, vom evalua erorile obt inute ^ n cazurile m= 1; m= 2  sim= 3.
Pentrum= 1:Zxn+1
xnf(x;y(x))dxZxn+1
xnp1(x)dx=Zxn+1
xn[f(x;y(x))p1(x)]dx
=f00(;y())
2Zxn+1
xn(xxn)(xxn+1)dx;
unde2[xn;xn+1]. Not^ andx=xn+ht, avem:
dx=hdt;
x=xn)t= 0;
x=xn+1)t= 1;
xxn=xn+htxn=ht;
xxn+1=xn+htxn+1=xn+htxnh=h(t1);
iarZxn+1
xnf(x;y(x))dxZxn+1
xnp1(x)dxM
2Z1
0jh3t(t1)jdt=h3
12M;
undeM= sup
abjf00(;y())j. Deci eroarea de integrare ^ n cazul metodei Adams-
Bashforth pentru m= 1 este de ordinul h3.
Pentrum= 2:Zxn+1
xnf(x;y(x))dxZxn+1
xnp2(x)dx=Zxn+1
xn[f(x;y(x))p2(x)]dx
=f000(;y())
3!Zxn+1
xn(xxn)(xxn+1)(xxn+2)dx:
Proced am ca  si ^ n cazul precedent  si not am x=xn+ht. Rezult a:
dx=hdt;
x=xn)t= 0;
x=xn+1)t= 1;
xxn=xn+htxn=ht;
xxn+1=xn+htxn+1=xn+htxnh=h(t1);
xxn+2=xn+htxn+2=xn+htxn2h=h(t2)
 si g asim
Zxn+1
xnf(x;y(x))dxZxn+1
xnp2(x)dxM
6Z1
0jh4t(t1)(t2)jdt
=M
6Z1
0jh4(t33t+ 2t)jdt=h4
24M:
44

Florina-Adriana Bosie Metode numerice pentru probleme cu condit ie init ial a
Similar cazului precedent, M= sup
abjf000(;y())j. Am ar atat c a eroarea de integrare
^ n cazul metodei Adams-Bashforth pentru m= 2 este de ordinul h4.
Pentrum= 3:
Zxn+1
xnf(x;y(x))dxZxn+1
xnp3(x)dx=Zxn+1
xn[f(x;y(x))p3(x)]dx
=f(4)(;y())
4!Zxn+1
xn(xxn)(xxn+1)(xxn+2)(xxn+3)dx:
Not amx=xn+ht:
dx=hdt
x=xn)t= 0
x=xn+1)t= 1
xxn=ht
xxn+1=h(t1)
xxn+2=h(t2)
xxn+3=xn+htxn+3=xn+htxn3h=h(t3)
Prin urmare
Zxn+1
xndxZxn+1
xnp3(x)dxM
24Z1
0jh5t(t1)(t2)(t3)jdt
=M
24Z1
0jh5(t46t3+ 11t26t)jdt=19h5
720M;
undeM= sup
abjf(4)(;y())j.^In nal, concluzion am c a eroarea de integrare ^ n cazul
metodei Adams-Bashforth pentru m= 3 este de ordinul h5.
Un dezavantaj major al formulelor multipas este c a nu sunt cu autostart. Prin
urmare, ^ n cazul metodelor Adams-Bashforth trebuie s a avem valorile f(x;y) prece-
dente pasului din care putem aplica metoda (de exemplu, ^ n metoda Adams-Bashforth
pentrum= 3 trebuie s a avem patru valori succesive ale funct iei f(x;y) ^ n punctele
diviziunii echidistante ^ nainte s a putem aplica formula). Aceste valori de ^ nceput tre-
buie s a e obt inute cu ajutorul unor metode independente. Am putea, de exemplu,
s a folosim algoritmul lui Taylor sau una dintre metodele Runge-Kutta pentru a obt ine
aceste valori de start. De asemenea, trebuie s a ne asigur am c a aceste valori de start
au precizia necesar a pentru a obt ine precizia cerut a ^ n ansamblu. Pe de alt a parte,
formulele multipas necesit a doar o evaluare la ecare pas, ^ n comparat ie cu metodele
Runge-Kutta, care necesitau mai multe evalu ari la ecare pas  si, deci, acestea sunt
considerabil mai rapide  si necesit a mai put in a munc a de calcul.
45

Florina-Adriana Bosie Metode numerice pentru probleme cu condit ie init ial a
3.2 Metoda Adams-Moulton
Formulele corectoare de ordin superior pot  obt inute folosind un polinom care inter-
poleaz a ^ nxn+1; xn;:::;xnmpentru un num ar ^ ntreg m> 0. Formula lui Newton cu
diferent e nite la st^ anga care interpoleaz a ^ n aceste m+ 2 puncte este
pm+1(x) =m+1X
k=0(1)ks
k
rkfn+1; (3.8)
undes=xxn+1
h.
Aceste diferent e se bazeaz a pe valorile fn+1; fn;:::;fnm. Metoda Adams-Moulton
const a ^ n ^ nlocuirea formulei \exacte" (3.1) cu formula \aproximativ a"
yn+1=yn+h(0r0fn+1+1r1fn+1+


+m+1rm+1fn+1); (3.9)
unde
k= (1)kZ0
1s
k
ds; k = 0;1;:::;m + 1:
Mai departe, vom calcula primii c^ at iva coecient i k.
0= (1)0Z0
1s
0
ds= 1;
1= (1)Z0
1s
1
ds=Z0
1(s)
1!ds=1
2;
2= (1)2Z0
1s
2
ds=Z0
1(s)(s1)
2!ds=1
12;
3= (1)3Z0
1s
3
ds=Z0
1(s)(s1)(s2)
3!ds=1
24;
4= (1)4Z0
1s
4
ds=Z0
1(s)(s1)(s2)(s3)
4!ds
=19
720:
46

Florina-Adriana Bosie Metode numerice pentru probleme cu condit ie init ial a
Prin urmare, am obt inut primele 5 valori ale coecient ilor k:
0= 1;  1=1
2;  2=1
12;  3=1
24;  4=19
720:
^Inlocuind aceste valori ^ n (3.9)  si d^ and diverse valori lui m, obt inem diferite clase
de formule.
Pentrum= 0:
yn+1=yn+h(0r0fn+1+1r1fn+1)
=yn+h
fn+11
2rfn+1
: (3.10)
Pentrum= 1:
yn+1=yn+h(0r0fn+1+1r1fn+1+2r2fn+1)
=yn+h
fn+11
2rfn+11
12r2fn+1
: (3.11)
Pentrum= 2:
yn+1=yn+h(0r0fn+1+1r1fn+1+2r2fn+1+3r3fn+1)
=yn+h
fn+11
2rfn+11
12r2fn+11
24r3fn+1
: (3.12)
Vom calcula valorile diferent elor cu ajutorul denit iei operatorului diferent  a nit a
la st^ anga, dup a cum urmeaz a:
r0fn+1=fn+1;
rfn+1=fn+1fn;
r2fn+1=r(rfn+1) =r(fn+1fn) =rfn+1rfn
=fn+1fn(fnfn1) =fn+12fn+fn1;
r3fn+1=r(r2fn+1) =r(fn+12fn+fn1) =rfn+12rfn+rfn1
=fn+1fn2(fnfn1) +fn1fn2
=fn+13fn+ 3fn1fn2:
47

Florina-Adriana Bosie Metode numerice pentru probleme cu condit ie init ial a
Vom^ nlocui aceste valori^ n (3.10), (3.11)  si (3.12) pentru a obt ine formulele Adams-
Moulton.
Pentrum= 0:
yn+1=yn+h
fn+11
2(fn+1fn)
=yn+h
2(fn+1+fn):
Pentrum= 1:
yn+1=yn+h
fn+11
2(fn+1fn)1
12(fn+12fn+fn1)
=yn+h
12(5fn+1+ 8fnfn1):
Pentrum= 2:
yn+1=yn+h
fn+11
2(fn+1fn)1
12(fn+12fn+fn1)
1
24(fn+13fn+ 3fn1fn2)
=yn+h
24(9fn+1+ 19fn5fn1+fn2): (3.13)
Eroarea de trunchiere a metodei Adams-Moulton este dat a de formula ( Pitea &
Postolache [2] )
jy(xn+1)yn+1jjy(xn)ynj+Zxn+1
xnf(x;y(x))Zxn+1
xnpm+1(x)dx
 si, la fel ca ^ n cazul erorii de trunchiere obt inut a la metoda Adams-Bashforth, aceasta
reprezint a suma dintre eroarea de trunchiere la metodele directe  si eroarea de trunchi-
ere la integrarea numeric a.
Cazul metodei Adams-Moulton pentru m= 2 este cel mai des utilizat. Formula
Adams-Moulton de ordinul IV, (3.13), este ^ n mod cert o formul a corector de tip ^ nchis
deoarecefn+1=f(xn+1;yn+1) implic a cantitatea necunoscut a yn+1. Prin urmare,
trebuie s a e rezolvat a prin iterat ie. Un predictor convenabil ce poate  utilizat cu
acest corector este formula Adams-Bashforth de ordinul IV. ^In acest caz, predictorul  si
corectorul au acela si ordin. Dac a heste ales corespunz ator, atunci o singur a aplicare
48

Florina-Adriana Bosie Metode numerice pentru probleme cu condit ie init ial a
a corectorului va conduce la o ^ mbun at at ire semnicativ a a preciziei. Dac a precizia
obt inut a astfel nu este sucient a, este mai bine s a reducem dimensiunea pasului dec^ at
s a aplic am formula corectoare din nou.
Metodele Adams descrise anterior sunt des folosite ^ n metodele cu ordin variabil  si
pas variabil. Obiectivul acestor metode este s a selecteze automat ordinul potrivit  si
pasul potrivit, fapt ce va minimiza cantitatea de munc a necesar a pentru obt inerea unei
precizii specicate pentru o problem a dat a. Un alt avantaj important al acestor metode
este c a sunt cu autostart deoarece o metod a de ordin inferior poate  folosit a la ^ nceput
 si ele pot  ajustate cu u surint  a pentru a furniza valorile lips a c^ and dimensiunea pasului
se schimb a.
3.3 Aplicat ii
Consider am problema Cauchy:
(
y0=1
x2y
xy2
y(1) =1:
S a se aproximeze solut ia aceastei probleme pe intervalul [1;2], cu pasulh= 0:2, folosind
metoda Adams-Bashforth, pentru m= 2, calcul^ and valoarea de start cu metoda Runge-
Kutta de ordinul II.
Solut ie. Primele 2 iterat ii le vom calcula cu ajutorul metodei Runge-Kutta de ordinul
II, astfel:
yn+1=yn+k0+k1
2
y0dat;
unde
k0=hf(xn;yn);
k1=hf(xn+h;yn+k0):
Pentru urm atoarele iterat ii, vom utiliza formula de recurent  a:
yn+1=yn+h
12(23fn16fn1+ 5fn2); n2:
Din datele problemei, avem: x0= 1,×1= 1:2,×2= 1:4,×3= 1:6,×4= 1:8,×5= 2,
y0=1,f(x;y) =1
x2y
xy2.
^Incepem prin a calcula valoarea funct iei ^ n punctul ( x0;y0).
f0=f(x0;y0) =f(1;1) = 1:
49

Florina-Adriana Bosie Metode numerice pentru probleme cu condit ie init ial a
Pentru prima iterat ie, calculat a folosind metoda Runge-Kutta de ordinul II, avem:
f0= 1;
x0= 1;
k0= 0:2;
k1= 0:14422;
y1=y0+k0+k1
2=0:82789:
La iterat ia a doua obt inem:
f1= 0:69900;
x1= 1:2;
k0= 0:13980;
k1= 0:10565;
y2=y1+k0+k1
2=0:70516:
Pentru iterat ia a treia, de la care ^ ncepem s a utiliz am metoda de tip Adams-
Bashforth, obt inem:
f2= 0:51665;
x2= 1:4;
y3=y2+h
12(23f216f1+ 5f0) =0:61018:
Pentru iterat ia a patra, g asim:
f3= 0:39966;
x3= 1:6;
y4=y3+h
12(23f316f2+ 5f1) =0:5365:
La ultima iterat ie obt inem:
f4= 0:31887;
x4= 1:8;
y5=y4+h
12(23f416f3+ 5f2) =0:47779:
Urm atorul tabel cont ine valorile obt inute pentru h= 0:1.
50

Florina-Adriana Bosie Metode numerice pentru probleme cu condit ie init ial a
nxny(xn)yn Eroarea
0 1 1 -1 0
11.1 -0.909090 -0.90827 0.00082
21.2 -0.833330 -0.83187 0.00146
31.3 -0.769230 -0.76716 0.00207
41.4 -0.714290 -0.71174 0.00255
51.5 -0.666670 -0.66371 0.00296
61.6 -0.625000 -0.6217 0.00330
71.7 -0.588240 -0.58462 0.00362
81.8 -0.555560 -0.55165 0.00391
91.9 -0.526320 -0.52214 0.00418
10 2-0.500000 -0.49556 0.00444
Programul Maple pentru rezolvarea acestei probleme, folosind metoda Adams-
Bashforth, pentru m= 2, este urm atorul:
restart;
f := (x, y) -> 1/x^2-y/x-y^2;
ec := diff(y(x), x) = f(x, y(x));
CI := y(1) = -1;
dsolve({ec, CI});
Sol := x -> -1/x;
Digits := 5;
m := 2; a := 1; b := 2; h := 0.2; N := (b-a)/h;
f := (x, y) -> 1/x^2-y/x-y^2;
x[0] := 1; y[0] := -1;
for n from 0 to N-1 do
print(nprintf("Iteratia %d:", n+1));
x[n+1] := x[n]+h;
phi[n] := f(x[n], y[n]);
if n+1 <= m then
k[0] := h*f(x[n], y[n]);
k[1] := h*f(x[n]+h, y[n]+k[0]);
y[n+1] := y[n]+(k[0]+k[1])*(1/2)
else
y[n+1] := y[n]+(1/12)*h*(23*phi[n]-16*phi[n-1]+5*phi[n-2])
end if;
print(nprintf("y[%d]=%g", n+1, y[n+1]));
51

Florina-Adriana Bosie Metode numerice pentru probleme cu condit ie init ial a
print(nprintf("Solutia exacta y(x[%d])=%f", n+1, Sol(x[n+1])));
print(nprintf("Eroarea in punctul x[%d]=%g", n+1, x[n+1]));
print(abs(Sol(x[n+1])-y[n+1]));
end do;
with(plots);
p1 := plot([seq([x[n], y[n]], n = 0 .. N)], color = blue):
p2 := plot(Sol(x), x = a .. b, color = red):
display(p2, p1);
^In urm atoarele grace este reprezentat a solut ia exact a  si solut ia obt inut a cu metoda
Adams-Bashforth, pentru m= 1,m= 2  sim= 3, cu diverse valori ale pasului de
discretizare.
Figura 3.1: Solut ia exact a  si solut ia Adams-Bashforth (m=1), autostart R-K II, h=0.5
Figura 3.2: Solut ia exact a  si solut ia Adams-Bashforth (m=2), autostart R-K II, h=0.5
52

Florina-Adriana Bosie Metode numerice pentru probleme cu condit ie init ial a
Figura 3.3: Solut ia exact a  si solut ia Adams-Bashforth (m=3), autostart R-K II, h=0.5
Figura 3.4: Solut ia exact a  si solut ia Adams-Bashforth (m=1), autostart R-K II, h=0.2
Figura 3.5: Solut ia exact a  si solut ia Adams-Bashforth (m=2), autostart R-K II, h=0.2
53

Florina-Adriana Bosie Metode numerice pentru probleme cu condit ie init ial a
Figura 3.6: Solut ia exact a  si solut ia Adams-Bashforth (m=3), autostart R-K II, h=0.2
54

Bibliograe
[1] S.D. Conte, C. de Boor, Elementary Numerical Analysis: An Algorithm Approach ,
McGraw-Hill Book Company, New York, 1980.
[2] A. Pitea, M. Postolache, Analiz a Numeric a Modele  si Tehnici pentru Ecuat ii
Diferent iale , Fair Partners, 2010.
[3] E. S uli, D. Mayers, An Introduction to Numerical Analysis , Cambridge University
Press, Cambridge, 2003.
55