INDRUMATOR DE LUCRARI DELABORATOR DE CALCUL NUMERIC
ALEXANDRU MIHAI BICA VIORICA PURJE
2
Cuprins
Introducere 5
1 Metode de interpolare a functiilor continue 71.1 Diferente divizate . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71.2 Polinomul de interpolare al lui Lagrange . . . . . . . . . . . . . . . . . . . 9
1.2.1 Algoritmul lui Aitken . . . . . . . . . . . . . . . . . . . . . . . . . . 101.2.2 Algoritmul lui Newton . . . . . . . . . . . . . . . . . . . . . . . . . 12
1.3 Polinomul de interpolare a lui Taylor . . . . . . . . . . . . . . . . . . . . . 151.4 Interpolarea Hermite neteda pe portiuni . . . . . . . . . . . . . . . . . . . 161.5 Interpolare spline liniara . . . . . . . . . . . . . . . . . . . . . . . . . . . . 191.6 Interpolare spline parabolica . . . . . . . . . . . . . . . . . . . . . . . . . . 201.7 Interpolare spline cubica . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
1.7.1 Spline cubic complet general de conditii bilocale . . . . . . . . . . . 221.7.2 Spline cubic natural generat de conditii initiale . . . . . . . . . . . 25
2 Ajustarea in medie patratica 292.1 Metoda celor mai mici patrate: dreapta de regresie . . . . . . . . . . . . . 292.2 Metoda celor mai mici patrate: parabola de regresie . . . . . . . . . . . . . 31
3 Formule de cuadratura 353.1 Exemple de functii continue cu primitiva in forma de serie in�nita . . . . . 353.2 Formula trapezului . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 373.3 Formula perturbata a trapezului . . . . . . . . . . . . . . . . . . . . . . . . 423.4 Formula lui Simpson . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 443.5 Formula lui Newton . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
4 Rezolvarea sistemelor liniare 514.1 Metoda directa a lui Gauss . . . . . . . . . . . . . . . . . . . . . . . . . . . 514.2 Metoda iterativa a lui Iacobi . . . . . . . . . . . . . . . . . . . . . . . . . . 53
5 Metode numerice pentru ecuatii neliniare 575.1 Metoda injumatatirii intervalului . . . . . . . . . . . . . . . . . . . . . . . 575.2 Metoda aproximatiilor succesive pentru ecuatii scalare . . . . . . . . . . . . 595.3 Metoda aproximatiilor succesive pentru sisteme neliniare . . . . . . . . . . 615.4 Metoda lui Newton pentru ecuatii scalare . . . . . . . . . . . . . . . . . . . 645.5 Metoda lui Newton pentru sisteme neliniare . . . . . . . . . . . . . . . . . 66
3
4 CUPRINS
5.6 Metoda coardei . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
6 Metode numerice pentru ecuatii diferentiale 756.1 Metoda retelelor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 756.2 Metoda lui Euler . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 776.3 Metoda lui Heun . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 796.4 Metoda Runge-Kutta . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 816.5 Metoda aproximatiilor succesive . . . . . . . . . . . . . . . . . . . . . . . . 83
Bibliogra�e 87
Introducere
Acest indrumator de laborator urmeaza cursul de calcul numeric [4] si este adresatstudentilor din anul II ai specializarilor Matematica si Informatica, din cadrul Facultatiide Stiinte. Desigur, acest indrumator poate � folosit si de catre studentii specializariloringineresti ce au in planul de invatamint un curs de metode numerice.Pe langa continutul clasic al unui indrumator de lucrari de laborator de analiza numer-
ica sau calcul numeric, prezentul indrumator contine multe elemente de noutate adaptatela evolutia actuala a metodelor numerice, obtinute prin consultarea unor elemente bib-liogra�ce noi (in acest sens putem mentiona lucrarile [2], [3], [5], [6], [9], [21]) inclusivrezultate proprii. Astfel, apar ca o noutate paragrafele 1.4, 1.5, 1.6, 1.7.2, 2.1, 2.2, 3.1,3.3, 6.5, tratand: procedura de interpolare neteda pe portiuni cubica de tip Hermitecu noduri duble folosind algoritmul lui Akima, interpolarea spline poligonala, interpo-larea spline parabolica generata de conditii initiale, interpolarea spline cubica naturalagenerata de conditii initiale, metoda celor mai mici patrate pentru obtinerea dreptei deregresie si a parabolei de regresie (utila la ajustarea datelor experimentale si a�ata catehnica matematica de cercetare in economie, psihologie, sociologie, statistica aplicata),cuadratura perturbata a trapezului (cu estimari recente ale restului formulei de cuadraturaatat pentru cuadratura trapezului, cat si pentru cea perturbata a trapezului) si metodaaproximatiilor succesive pentru problema Cauchy asociata ecuatiilor diferentiale ordinarede ordinul intai. De asemenea, la inceputul capitolului 3 sunt prezentate doua exemple deintegrale de�nite (provenite din geometrie si teoria probabilitatilor) ce nu se pot calculacu exactitate, motivand astfel utilizarea formulelor de cuadratura.Capitolul intai este dedicat procedeelor de interpolare a functiilor continue, prezentandu-
se polinomul de interpolare a lui Lagrange cu cele doua variante iterative de calcul gener-ate de algoritmul lui Aitken, respectiv algoritmul lui Newton, apoi polinomul lui Taylorsi procedeele de interpolare neteda pe portiuni (cubica de tip Hermite cu noduri duble,spline poligonala, spline parabolica generata de conditii initiale, spline cubica naturalagenerata de conditii initiale si spline cubica completa generata de conditii bilocale). Aldoilea capitol este dedicat aproximarii polinomiale in medie patratica utilizand metodacelor mai mici patrate exempli�cata pe polinoamele de gradul intai si doi si obtinanddreapta de regresie si respectiv, parabola de regresie.Al treilea capitol prezinta formulele de cuadratura de tip interpolator precum, for-
mula trapezului, formula lui Simpson, formula lui Newton (toate trei �ind formule de tipNewton-Cotes) si formula perturbata a trapezului (o formula de tip Euler-Mac Laurin cudoa noduri duble). Capitolele 4 si 5 prezinta metode numerice pentru ecuatii si sistemede ecuatii (liniare in capitolul 4, neliniare in capitolul 5). Astfel, in capitolul 3 aparemetoda exacta a lui Gauss bazata pe eliminarea ce utilizeaza tehnica pivotului maxim si
5
metoda iterativa a lui Iacobi. Rezolvarea aproximativa a sistemelor neliniare este abordataprin metoda aproximatiilor succesive si metoda lui Newton (cu algoritmi elaborati pentrusisteme de doua ecuatii cu doua necunoscute), iar a ecuatiilor neliniare prin metoda inju-matatirii intervalului, metoda tangentei (a lui Newton), metoda coardei, metoda aproxi-matiilor succesive, precedate de tehnica sirului lui Rolle de separare a radacinilor.Ultimul capitol este dedicat metodelor discrete de aproximare a solutiilor ecuatiilor
diferentiale ordinare cu conditii la limita. Intai este prezentata metoda retelelor pentruproblema bilocala asociata ecuatiei diferentiale liniare de ordinul doi si apoi sunt prezen-tate metodele unipas de tip Euler, Heun si Runge-Kutta. Capitolul se incheie cu metodaaproximatiilor succesive bazata pe tehnica punctului �x si cuadratura trapezului. Cei maimulti dintre algoritmii prezentati sunt insotiti si de programe corespunzatoare in limbajulC.
Capitolul 1
Metode de interpolare a functiilorcontinue
1.1 Diferente divizate
Consideram intr-un interval I � R, multimea de noduri
A = fak : ak 2 R; k = 0;m; g � I
si �e functia f : A ! R.De�nitia 1: Se numeste diferenta divizata de ordinul intai cu nodurile ar si ar+1
relativa la functia f , expresia
Df (ar)notatie= [ar; ar+1; f ] =
f (ar+1)� f (ar)
ar+1 � ar; r = 0;m� 1:
De�nitia 2: Se numeste diferenta divizata de ordinul k cu nodurile ar; ar+1; :::; ar+k;relativa la functia f; expresia
Dkf (ar)notatie= [ar; ar+1; :::; ar+k; f ] =
Dk�1f (ar+1)�Dk�1f (ar)
ar+k � ar; (1.1)
r = 0;m� 1; k = 1;m� r:
Prin conventie, se noteaza D0f (ar) = f (ar) ; r = 0;m: Deasemenea, este cunoscutfaptul ca operatorul Dk este liniar.Cand nodurile din A sunt echidistante atunci intre diferentele divizate si diferentele
�nite exista relatia
Dmf (a) =�mh f (a)
m!hm:
In general, pentru diferentele divizate exista exprimarea :
[ar; ar+1; :::; ar+k; f ] =kXi=0
f (ar+i)
u0 (ar+i)
7
8 CAPITOLUL 1. METODE DE INTERPOLARE A FUNCTIILOR CONTINUE
undeu (x) = (x� ar) � (x� ar+1) � ::: � (x� ar+k) :
Observatie: Pentru r = 0 si k = m
[a0; a1; :::; am; f ] =mXk=0
f (ak)
u0 (ak)=
mXk=0
f (ak)
uk (ak)
undeu (x) = (x� a0) � (x� a1) � ::: � (x� am)
si
uk (x) =u (x)
x� ak; 8k = 0;m:
In calcule, in special la constructia formulei de interpolare a lui Newton pentru polinomullui Lagrange, este util tabloul de diferente divizate.
x f Df D2f D3f ....... Dm�1f Dmfa0 f0 Df0 D2f0 D3f0 ....... Dm�1f0 Dmf0a1 f1 Df1 D2f1 D3f1 ....... Dm�1f1a2 f2 Df2 D2f2 D3f2 .......:::: :::: :::::: :::::: :::::am�2 fm�2 Dfm�2 D2fm�2am�1 fm�1 Dfm�1am fm
unde fi = f (ai) ; i = 0;m:La constructia tabelului s-a folosit formula de recurenta
Dkf (ai) =Dk�1f (ai+1)�Dk�1f (ai)
ai+k � ai; k = 1;m; i = 0;m� k:
Aceeasi formula de recurenta o foloseste algoritmul de mai jos :
I. Date de intrare :m ordinul maxim al diferentei divizateak; k = 0;m nodurile diferentelor divizatef (ak) ; k = 0;m valorile functiei pe noduri.II. Date de iesire : Tabloul diferentelor divizate pina la ordinul m, adica Dkf (ai) ;
k = 0;m; i = 0;m� k; sau numai valoarea Dmf (a0) :III. Pasii algoritmului :Pasul 1. Pentru i = 0;m ataseaza
D0fi := f (ai) :
Pasul 2. Pentru k = 1;m;Pentru i = 0;m� k; calculeaza
Dkfi :=Dk�1fi+1 �Dk�1fi
ai+k � ai
1.2. POLINOMUL DE INTERPOLARE AL LUI LAGRANGE 9
Pasul 3. A�seaza Dkfi; k = 1;m; i = 0;m� k; sau a�seaza numai Dmf0: STOP.
Probleme propuse :1. Sa se construiasca un program de calcul al tabloului de diferente divizate in limbajul
C.2. Considerând functia f : R ! R, f (x) = 1
1+x2sa se construiasca tabloul cu difer-
entele divizate relativ la nodurile �3;�1; 0; 1; 1:5 pina la ordinul 4.3. Pentru functia densitate de probabilitate a lui Gauss
f : R! R; f (x) =1p2�� e� 1
2x2
sa se construiasca tabloul de diferente divizate pina la ordinul 5 cu nodurile echidistanteak = a+ kh; k = 0; 5; h = 0:4:
1.2 Polinomul de interpolare al lui Lagrange
Sunt prezentate formula de interpolare a lui Newton ce foloseste recurent diferentele di-vizate (diagonala tabloului de diferente divizate) si procedura de interpolare a lui Aitkenpentru polinomul lui Lagrange.Se considera cunoscute valorile functiei f : [a; b]! R; pe nodurile distincte ak 2 [a; b];
k = 0;m; yk = f (ak) :De�nitia 1: Se numeste polinomul de interpolare al lui Lagrange relativ la functia
f si la nodurile ak; k = 0;m; polinomul de grad minim Lmf care satisface conditiile deinterpolare
Lmf (ak) = f (ak) = yk; 8k = 0;m:
Observatie: Gradul polinomului Lmf este m, iar conditiile de interpolare de maisus conduc la un sistem liniar de m+ 1 ecuatii cu cei m+ 1 coe�cienti ca necunoscute siastfel, polinomul de interpolare al lui Lagrange este unic determinat de aceste conditii deinterpolare. Se obtine expresia lui Lmf ,
Lmf (x) =
mXk=0
u (x)
(x� ak) � u0 (ak)� f (ak) ; 8x 2 [a; b]; (1.2)
undeu (x) = (x� a0) � (x� a1) � ::: � (x� am) :
Este cunoscut faptul ca operatorul de interpolare al lui Lagrange Lm este liniar. Deaseme-nea, in formula de interpolare a lui Lagrange
f (x) = Lmf (x) +Rmf (x) ; x 2 [a; b] (1.3)
are loc estimarea erorii :
jRmf (x)j ��
ju (x)j � [x; a0; :::; am; f ]; 8x 2 [a; b]; daca f 2 C[a; b]1
(m+1)!� ju (x)j �
f (m+1) C
8x 2 [a; b]; daca f 2 Cm+1[a; b]: (1.4)
10 CAPITOLUL 1. METODE DE INTERPOLARE A FUNCTIILOR CONTINUE
Din punct de vedere practic, dându-se o fonctie f : [a; b] ! R; cunoscuta doar penodurile ak 2 [a; b]; k = 0;m; si precizându-se o eroare admisa " > 0, se cere sa sedetermine m 2 N si sa se aproximeze valoarea f (�) prin Lmf (�) cu o eroare ce nu depas-este ": In acest caz, presupunandu-se cunoscut un majorant Mm+1 (f) ; pentru normelederivatelor
f (k) C; k = 0;m+ 1 se determina m astfel incat
ju (�)j(m+ 1)!
�Mm+1 (f) < ":
si atunci, jRmf (�)j < ": Dar aici, cerinta cunoasterii majorantului Mm+1 (f) nu poate �ocolita. Acest aspect restrange sfera de aplicabilitate a procedeului.
1.2.1 Algoritmul lui Aitken
Ometoda mult mai practica este stabilita prin algoritmul lui Aitken, bazat pe urmatoareaproprietate de recurenta a polinomului lui Lagrange:Proprietatea de recurenta: Notand L0f (x; aj) = f (aj) ; j = 0;m se obtine relatia
de recurenta
Lmf(x; a0; a1; :::; am) =x� aiak � ai
� Lm�1f(x; a0; :::; ai�1; ai+1; :::; am)�
� x� akak � ai
� Lm�1f(x; a0; :::; ak�1; ak+1; :::; am); k; i = 0;m; 8x 2 [a; b]:
Metoda lui Aitken consta in generarea tablouluia0 f00a1 f10 f11a2 f20 f21 f22a3 f30 f31 f32 f33...
......
......
am fm0 fm1 fm2 fm3 ..... fmmunde fi0 = f (ai) ; i = 0;m, iar elementele celorlalte coloane sunt calculate prin
fi;j+1 =1
ai � aj����� fjj aj � xfij ai � x
����pentru i = 1;m si j = 0; i� 1:Observam ca
f11 =x� a1a0 � a1
� f (a0) +x� a0a1 � a0
� f (a1) = L1f (x)
si analog se obtine pentru i = 2;m
fii = Lif (x)
adica polinomul de interpolare Lagrange relativ la nodurile a0; a1; :::; ai:
1.2. POLINOMUL DE INTERPOLARE AL LUI LAGRANGE 11
Astfel, sirul elementelor de pe diagonala tabloului f11; f22; :::; fii; :::; fmm devine un sirde aproximatii ale lui f (x) : Acest sir va converge catre f (x) atunci cand
limm!1
Lmf (x) = f (x) ;
ceea ce in baza criteriului lui Cauchy este echivalent cu
limi!1
jfii � fi�1;i�1j = 0:
De aici rezulta un criteriu de oprire a algoritmului lui Aitken dat de conditia
jfii � fi�1;i�1j < "
pentru " > 0 dat initial. Prin urmare, algoritmul lui Aitken genereaza primele i + 1coloane ale tabloului oprindu-se la prima coloana pentru care este indeplinita conditiade mai sus. Asadar, gradul maxim al polinomului de interpolare este m, dar algoritmulse poate opri la un grad mai mic (in urma atingerii preciziei dorite). Deoarece eroareaformulei de interpolare depinde de ju (x)j ; rezulta ca pentru a opri acest algoritm cat mairepede este indicat sa se ordoneze nodurile a0; :::; am dupa distanta lor fata de �; punctulin care dorim sa aproximam pe f (�) : Astfel, se ordoneaza nodurile pina este indeplinitaconditia
jai � xj � jaj � xj ; daca i < j; i = 0;m; j = 1;m:
Dupa cum se observa, metoda lui Aitken si criteriul sau de oprire necesita doar cunoast-erea valorilor f (ai) ; i = 0;m. Acest avantaj justi�ca prezenta acestui algoritm inbibliotecile matematice ale softurilor de calcul numeric (gen Matlab, Mathcad sau Math-ematica).
Algoritmul lui Aitken:I. Date de intrare:m gradul maxim al polinomului de interpolareai; i = 0;m; nodurile de interpolaref expresia functiei, sau valorile ei pe nodurile de interpolare f (ai) ; i = 0;m" eroarea absoluta maxima admisibila� 2 [a; b] punctul in care se aproximeaza functia f:II. Date de iesire:Lif (�) valoarea calculata cu precizia ceruta (i � m este gradul polinomului pentru
care s-a atins aceasta precizie)sau un mesaj de eroare in care se speci�ca faptul ca in conditiile date nu s-a putut
calcula o aproximatie cu precizia dorita.III. Pasii algoritmului:Pasul 1. Se ordoneaza nodurile in functie de distanta lor fata de punctul �; adicaPentru i = 0;m
Pentru j = i+ 1;mDaca jai � �j > jaj � �j atunci u := a[i]; a[i] := a[j]; a[j] := u;
Pasul 2. Pentru i = 0;m executa fi0 := f (ai) ;
12 CAPITOLUL 1. METODE DE INTERPOLARE A FUNCTIILOR CONTINUE
Pasul 3. Pentru i = 1;m si j = 0; i� 1 calculeaza
fi;j+1 =1
ai � aj� [fjj � (ai � �)� fij � (aj � �)];
Daca jfii � fi�1;i�1j < " atunci tipareste Lif (�) = fii: STOPAltfel, i := i+ 1;
Pasul 4. Daca jfmm � fm�1;m�1j < " atunci tipareste Lmf (�) = fmm:Altfel, tipareste " Nu s-a putut calcula o aproximatie cu precizia dorita ".
STOP.
Probleme propuse :1. Sa se construiasca un program in limbajul C pentru algoritmul lui Aitken.2. Valorile functiei f (x) = 3
px sunt date pe noduri in tabelul
x 1:0 1:1 1:3 1:5 1:6f (x) 1 1:032 1:091 1:145 1:17Folosind algoritmul lui Aitken, sa se calculeze f (1:15) cu o eroare absoluta maxima
de 10�3:3. Valorile unei functii f sunt date pe noduri prin intermediul tabeluluix 1:00 1:08 1:13 1:20 1:27 1:31 1:38f (x) 1:17520 1:30254 1:38631 1:50946 1:21730 1:22361 1:23470
Folosind algoritmul lui Aitken, sa se calculeze cu o precizie de 10�5 valorile functiei inpunctele 1:134; 1:151; 1:185:
1.2.2 Algoritmul lui Newton
Newton a dat o formula utila, din punct de vedere al calcului si al usurintei de progra-mare (in virtutea modului recurent de constructie), pentru polinomul de interpolare al luiLagrange.Pentru a obtine aceasta formula consideram diferentele divizate
[x; a0; f ]; [x; a0; a1; f ]; :::; [x; a0; a1; :::; am; f ]
exprimate cu ajutorul diferentelor de ordin inferior,
[x; a0; a1; :::; am; f ] =[x; a0; a1; :::; am�1; f ]� [a0; a1; :::; am; f ]
x� am
[x; a0; a1; :::; am�1; f ] =[x; a0; a1; :::; am�2; f ]� [a0; a1; :::; am�1; f ]
x� am�1
.....................................................................................
[x; a0; a1; f ] =[x; a0; f ]� [a0; a1; f ]
x� a1
[x; a0; f ] =f (x)� f (a0)
x� a0:
1.2. POLINOMUL DE INTERPOLARE AL LUI LAGRANGE 13
Eliminând din aceste m + 1 relatii diferentele divizate [x; a0; :::; ak; f ]; k = 0;m� 1 ,adica inmultind a doua egalitate cu 1
x�am , pe a treia egalitate cu1
(x�am�1)(x�am)si asa maideparte pina ultima egalitate se inmulteste cu
1
(x� a1) � ::: � (x� am)
si adunând egalitatile obtinute membru cu membru si apoi exprimând f (x) se obtine
f (x) = f (a0) +mXi=1
(x� a0) � ::: � (x� ai�1) � [a0; a1; :::; ai; f ] + u (x) � [x; a0; :::; am; f ]:
Astfel, s-a obtinut tocmai formula de interpolare a lui Lagrange
f (x) = Nmf (x) +Rmf (x) ; x 2 [a; b] (1.5)
cu polinomul de interpolare dat prin formula lui Newton
Nmf (x) = f (a0) +mXi=1
(x� a0) � ::: � (x� ai�1) � [a0; a1; :::; ai; f ] = (1.6)
= f (a0) +mXi=1
(x� a0) � ::: � (x� ai�1) �Dif0
iar restul dat cu ajutorul diferentelor divizate
Rmf (x) = u (x) � [x; a0; :::; am; f ]:
Se observa ca in formula (1.6) s-au folosit doar diferentele divizate de pe latura su-perioara a tabloului de diferente. Deasemenea, din (1.3) si (1.5) rezulta Lmf (x) =Nmf (x) ; 8x 2 [a; b]:Observatie: Din formula (1.6) rezulta ca polinomul de interpolare al lui Newton
veri�ca relatia de recurenta
Nkf (x; a0; :::; ak) = Nk�1f (x; a0; :::; ak�1) + (x� a0) � ::: � (x� ak�1) � [a0; :::; ak; f ];
8k = 1;m; unde N0f (x) = f (a0) : Aceasta proprietate permite construirea recurenta aacestui polinom (deci programare mai facila) precum si oprirea algoritmului odata atinsaprecizia dorita, dupa ce in prealabil s-au reordonat nodurile in functie de distanta lor pinala punctul in care dorim sa aproximam functia f: Astfel se obtin iterativ polinoamelelui Newton de gradele 1, 2, ... intr-un mod asemanator cu generarea polinoamelor luiLagrange prin metoda lui Aitken. Criteriul de oprire este asemanator. Daca
jNkf (�)�Nk�1f (�)j < "
atunci algoritmul se opreste la ordinul k, adica la polinomul Nkf (x; a0; :::; ak) : Insa gradulmaxim al polinomului generat de algoritmul de mai jos este m.
14 CAPITOLUL 1. METODE DE INTERPOLARE A FUNCTIILOR CONTINUE
Algoritmul polinomului lui Newton:I. Date de intrare:m gradul maxim al polinomului de interpolareak; k = 0;m; nodurile de interpolaref expresia functiei, sau valorile ei pe nodurile de interpolare f (ak) ; k = 0;m" eroarea absoluta maxima admisibila� 2 [a; b] punctul in care se aproximeaza valoarea functiei f:II. Date de iesire:Nkf (�) valoarea calculata cu precizia ceruta (k � m este gradul polinomului pentru
care s-a atins aceasta precizie)sau un mesaj de eroare in care se speci�ca faptul ca in conditiile date nu s-a putut
calcula o aproximatie cu precizia dorita.III. Pasii algoritmului:Pasul 1. Se ordoneaza nodurile in functie de distanta lor fata de punctul �; adicaPentru k = 0;m
Pentru j = k + 1;mDaca jak � �j > jaj � �j atunci u := a[k]; a[k] := a[j]; a[j] := u;
Pasul 2. Se calculeaza diferentele divizate [a0; :::; ak; f ] = Dkf0; k = 0;m utilizandprocedura din paragraful 1.2 sau obtinând aici iterativ,
D0fk := f (ak) ; pentru k = 0;m
si
Djfk :=Dj�1fk+1 �Dj�1fk
ak+j � akpentru j = 1;m si pentru k = 0;m� j:
Pasul 3. k := 0; N0f (�) = f (a0)Pasul 4. Pentru k � 0 cat timp k < m executa
Nkf (�) := Nk�1f (�) + (�� a0) � ::: � (�� ak�1) � [a0; :::; ak; f ]Daca jNkf (�)�Nk�1f (�)j < " atunci tipareste Nkf (�) : STOP.
Pasul 5. Daca jNmf (�)�Nm�1f (�)j < " atunci tipareste Nmf (�). STOP.Altfel, tipareste " Nu s-a putut calcula o aproximatie cu precizia dorita ".
STOP.
Probleme propuse :1.Sa se construiasca un program in limbajul C pentru algoritmul lui Newton.2. Valorile functiei f (x) = 3
px sunt date pe noduri in tabelul
x 1:0 1:1 1:3 1:5 1:6f (x) 1 1:032 1:091 1:145 1:17Folosind algoritmul lui Newton, sa se calculeze f (1:15) cu o eroare absoluta maxima
de 10�5:3. Valorile unei functii f sunt date pe noduri prin intermediul tabeluluix 1:00 1:08 1:13 1:20 1:27 1:31 1:38f (x) 1:17520 1:30254 1:38631 1:50946 1:21730 1:22361 1:23470
Folosind algoritmul lui Newton, sa se calculeze cu o precizie de 10�5 valorile functieiin punctele 1:134; 1:151; 1:185: Sa se compare cu valorile obtinute prin algoritmul luiAitken. Apoi, sa se calculeze cu o precizie de 10�5 valorile functiei in punctele 1.140,1.175, 1.195, cu ajutorul algoritmului lui Newton.
1.3. POLINOMUL DE INTERPOLARE A LUI TAYLOR 15
1.3 Polinomul de interpolare a lui Taylor
Consideram functia f : [a; b] ! R; pentru care sunt cunoscute pe nodurile ak � [a; b];k = 0;m; valorile f (j) (ak) ; j = 0; rk; k = 0;m; unde rk 2 N si �e n = m+r0+ :::+rm:De�nitia 1: Se numeste polinomul de interpolare al lui Hermite relativ la functia f
, la nodurile ak; k = 0;m; si la valorile f (j) (ak) ; j = 0; rk; k = 0;m; polinomul degrad minim Hnf care satisface conditiile de interpolare :
Hnf(j) (ak) = f (j) (ak) ; 8j = 0; rk; 8k = 0;m: (1.7)
Fieu (x) = (x� a0)
r0+1 � (x� a1)r1+1 � ::: � (x� am)
rm+1
si
uk (x) =u (x)
(x� ak)rk+1
8k = 0;m:
Conditiile (1.7) determina in mod unic polinomul de interpolare a lui Hermite, obtinandu-se
Hnf (x) =mXk=0
uk (x) �
rkXj=0
"(x� ak)
j
j!��f (t)
uk (t)
�(j)jt = ak
#!; 8x 2 [a; b] (1.8)
iar in formula de interpolare a lui Hermite
f (x) = Hnf (x) +Rnf (x) ; x 2 [a; b]
estimarea restului este
jRnf (x)j �1
(n+ 1)!� ju (x)j �
f (n+1) C
8x 2 [a; b]; daca f 2 Cn+1[a; b]:
Pentru m = 0 si n = r0 in formula (1.8) se obtine polinomul lui Taylor ca un polinomde interpolare cu un singur nod multiplu de ordinul n+1,
Tnf (x) =nXk=0
(x� a0)k
k!� f (k) (a0) (1.9)
iar estimarea restului este
jRnf (x)j �1
(n+ 1)!� jx� a0jn+1 �
f (n+1) C
daca f 2 Cn+1[a0 � h; a0 + h]:
Din formula (1.9) se observa ca polinomul lui Taylor poate � obtinut recurent prin
T0f (x) = f (a0)
Tkf (x) = Tk�1f (x) +(x� a0)
k
k!� f (k) (a0) ; pentru k = 1; n:
Astfel, se obtine algoritmul urmator :
16 CAPITOLUL 1. METODE DE INTERPOLARE A FUNCTIILOR CONTINUE
Algoritm pentru polinomul lui TaylorI. Date de intrare:n gradul polinomului de interpolarea nodul de interpolaref (k) (a) ; k = 0; n; valorile pe nodul au 2 [a; b] punctul in care se aproximeaza valoarea functiei fII. Date de iesire:Tnf (u) aproximatia lui f (u) prin polinomul lui TaylorIII. Pasii algoritmului:Pasul 1. �e f0 := f (a) , p (0) := 1 si pentru k = 1; n
fk := f (k) (a)
p (k) := p (k � 1) � k
Pasul 2. �e T0f (u) = f (a) si pentru k = 1; n
Tkf (u) = Tk�1f (u) +1
p (k)� (u� a)k � fk
Pasul 3. Tipareste Tnf (u) : STOP.
Probleme propuse :1.Sa se construiasca un program in limbajul C pentru algoritmul polinomului lui Tay-
lor.2. Pentru n = 7; a = 0; u = 1; f (k) (a) = 1; k = 0; n sa se calculeze T7f (u) si sa se
compare rezultatul cu e: Se obtine aici o aproximare a numarului e folosind polinomul luiTaylor pentru interpolarea functiei f (x) = ex:
3. Pentru n = 5, u = �4; a = 0; f (0) (a) = 0; f (1) (a) = 1; f (2) (a) = 0; f (3) (a) =
�1; f (4) (a) = 0; f (5) (a) = 1 sa se calculeze T5f (u) si sa se compare rezultatul cusin �
4=
p22: Se obtine o aproximare a numarului
p2 folosind polinomul lui Taylor pentru
interpolarea functiei f (x) = sinx: Sa se aplice apoi algoritmul pentru u = �18= 100
aproximand astfel numarul sin (100) :4. Sa se construiasca o problema similara pentru f (x) = cosx , n = 6 si u = �
6(se va
obtine o aproximatie pentrup3):
1.4 Interpolarea neteda pe portiuni prin polinomullui Hermite cu noduri duble. Metoda lui Akima
Pentru f : [a; b]! R; f 2 C1[a; b]; polinomul de interpolare al lui Hermite cu doua noduriduble este unic determinat de conditiile�
P (a) = f (a) ; P 0 (a) = f 0 (a)P (b) = f (b) ; P 0 (b) = f 0 (b)
1.4. INTERPOLAREA HERMITE NETEDA PE PORTIUNI 17
avand pe [a; b] expresia
P (x) =(b� x)2 (x� a)
(b� a)2� f 0 (a)� (x� a)2 (b� x)
(b� a)2� f 0 (b)+
+(b� x)2 [2 (x� a) + (b� a)]
(b� a)3� f (a) + (x� a)2 [2 (b� x) + (b� a)]
(b� a)3� f (b) :
Daca f 2 C4[a; b] atunci pentru estimarea erorii formulei de interpolare
f (x) = P (x) +R (x) (1.10)
avem
jR (x)j � (b� a)4
3 � 27 � f (4)
C8x 2 [a; b]:
Considerand o diviziune a intervalului [a; b]
� : a = x0 < x1 < ::: < xn�1 < xn = b
si valorile pe aceste noduri ale unei functii f 2 C1[a; b], f (xi) ; f 0 (xi) ; i = 0; n; se poateconstrui o functie neteda F 2 C1[a; b] pentru care
F (xi) = f (xi) ; F 0 (xi) = f 0 (xi) ; 8i = 0; n
cerând ca restrictia functiei F la �ecare subinterval [xi�1; xi]; i = 1; n sa �e Fi dat prin
Fi (x) =(xi � x)2 (x� xi�1)
(xi � xi�1)2 � f 0 (xi�1)�
(x� xi�1)2 (xi � x)
(xi � xi�1)2 � f 0 (xi)+ (1.11)
+(xi � x)2 [2 (x� xi�1) + (xi � xi�1)]
(xi � xi�1)3 �f (xi�1)+
(x� xi�1)2 [2 (xi � x) + (xi � xi�1)]
(xi � xi�1)3 �f (xi) :
Se observa ca F 00i exista pe intervalele deschise (xi�1; xi) pentru orice i = 1; n si este unpolinom de gradul intai. Atunci F 00 este continua pe portiuni si marginita pe [a; b]: Prinurmare, F 2 C1[a; b] \ L2[a; b]: In cele ce urmeaza vom nota yi = f (xi) ; 8i = 0; n:Algoritmul de mai jos realizeaza aproximarea valorii f (�) a functiei f intr-un punct
� 2 [a; b]; � 6= xi; 8i = 0; n prin valoarea F (�) : De cele mai multe ori, sunt cunos-cute (prin masuratori experimentale) doar valorile f (xi), i = 0; n; urmand ca valorilederivatei pe noduri sa �e aproximate. O metoda de aproximare, generata prin rationa-mente de natura geometrica, a fost propusa de H. Akima in 1970 si este expusa in [14]. Pescurt, dandu-se nodurile xi; i = 0; n si valorile yi; i = 0; n, se calculeaza intai pantele
mi =xi+1 � xiti+1 � ti
; i = 0; n� 1
iar apoi, se aproximeaza valorile derivatelor pe unele dintre nodurile interioare prin,
f 0 (xi) =jmi+2 �mi+1j �mi�1 + jmi�1 �mi�2j �mi+1
jmi+2 �mi+1j+ jmi�1 �mi�2j; i = 2; n� 3: (1.12)
18 CAPITOLUL 1. METODE DE INTERPOLARE A FUNCTIILOR CONTINUE
Formula (1.12) constituie un punct tare al metodei lui Akima, dar deoarece nu poatepropune o formula si pentru primele doua si ultimele trei noduri conduce la necesitateaintroducerii arti�ciale a unor pante corespunzatoare extremitatilor intervalului [a; b] pen-tru a extinde formula (1.12), ceea ce constituie un punct slab al metodei. In baza metodeilui Akima se introduc atunci doar pantele suplimentare
m�2 = 3m0 � 2m1; m�1 = 2m0 �m1
mn = 2mn�1 �mn�2; mn+1 = 3mn�1 � 2mn�2; mn+2 = 3mn � 2mn�1
si se extinde formula (1.12) pentru i = 0; n: In continuare functia F este construita inbaza formulei (1.11).
Algoritmul lui AkimaI. Date de intrare:nodurile xi; i = 0; nvalorile yi = f (xi) ; i = 0; n� punctul in care se realizeaza aproximareaII. Date de iesire:F (�) valoare ce aproximeaza pe f (�) :III. Pasii algoritmului:Pasul 1. Pentru i = 0; n� 1 calculeaza pantele
mi :=xi+1 � xiti+1 � ti
Pasul 2. Calculeaza
m�2 = 3m0 � 2m1; m�1 = 2m0 �m1
mn = 2mn�1 �mn�2; mn+1 = 3mn�1 � 2mn�2; mn+2 = 3mn � 2mn�1
Pasul 3. Pentru i = 0; n calculeaza
f 0 (xi) :=jmi+2 �mi+1j �mi�1 + jmi�1 �mi�2j �mi+1
jmi+2 �mi+1j+ jmi�1 �mi�2j
Pasul 4. Se determina indicele j 2 f1; :::; ng pentru care
xj�1 < � < xj:
Pasul 5. Calculeaza
F (�) :=(xj � �)2 (�� xj�1)
(xj � xj�1)2 � f 0 (xj�1)�
(�� xj�1)2 (xj � �)
(xj � xj�1)2 � f 0 (xj)+
+(xj � �)2 [2 (�� xj�1) + (xj � xj�1)]
(xj � xj�1)3 � f (xj�1) :+
1.5. INTERPOLARE SPLINE LINIARA 19
+(�� xj�1)
2 [2 (xj � �) + (xj � xj�1)]
(xj � xj�1)3 � f (xj)
Pasul. 6 Tipareste F (�) : STOP.
Probleme propuse:1. Sa se construiasca un program in limbajul C pentru algoritmul lui Akima.2. La momentele 7.5 (adica ora 730), 10.5 (ora 1030), 13, 15.5 (ora 1530), 18, 21, 24 si
27 (adica ora 3 A. M. a doua zi) s-au masurat valorile glicemiei, obtinandu-se 130, 121,128, 96, 122, 138, 114, 90 (masurate in mg/dl). Sa se aproximeze glicemia acestui pacientde la orele 12, 14, 23 si 26 (a doua zi A. M. la ora 2) folosind algoritmul lui Akima.
1.5 Interpolare spline liniara
Se prezinta functia spline poligonala de interpolare a datelor experimentale (utilizata inmulte domenii ca metoda empirica de trasare a gra�cului prin puncte, spre exemplu inspitale la curba temperaturii pacientilor internati) si doua tipuri de functii spline cubicegenerate de conditii bilocale si respectiv, conditii initiale. Pentru functia spline cubicagenerata de conditii bilocale s-a ales varianta mai usor programabila (fara expresii deforma (x� ai)+) din [8].Consideram un interval [a; b] si o diviziune a acestui interval
� : a = x0 < x1 < ::: < xn�1 < xn = b:
Notam hi = xi � xi�1 si Ii = [xi�1; xi]; 8i = 1; n:De�nitie: Functia s : [a; b] �! R; se numeste spline polinomial de gradul m de
interpolare a valorilor y0; y1; :::; yn pe nodurile xi; i = 0; n daca:(i) s 2 Cm�1[a; b](ii) restrictia lui s la �ecare subinterval Ii , i = 1; n; este polinom de gradul m(iii) s (xi) = yi; 8i = 0; n:Daca y0; y1; :::; yn sunt valorile pe nodurile xi; i = 0; n ale unei functii f : [a; b] �! R;
atunci functia spline este de interpolare a functiei f pe nodurile diviziunii �:Fie
�n : a = x0 < x1 < ::: < xn�1 < xn = b
o diviziune a intervaluluii [a; b]: Not¼am hi = xi�xi�1 si Ii = [xi�1; xi]; 8i = 1; n si sistemulde n + 1 numere reale, y = (y0; y1; :::; yn) 2 Rn+1 asociat diviziunii �n. Pentru �ecarei = 1; n se de�neste si : Ii �! R;
si(x) = yi�1 +yi � yi�1
hi� (x� xi�1); 8x 2 Ii:
Gra�cul functiei si este segmentul ce uneste punctele (xi�1; yi�1) si (xi; yi):Dac¼a construim s : [a; b] �! R; astfel încât s jIi= si; 8i = 1; n, atunci s este o functie
continu¼a pe [a; b], al c¼arei gra�c este linia poligonal¼a ce uneste punctele (x0; y0); (x1; y1); :::; (xn; yn):Functia s este o functie spline polinomiala de interpolare de gradul intai, numit splinepoligonal.
20 CAPITOLUL 1. METODE DE INTERPOLARE A FUNCTIILOR CONTINUE
Teorema: (a se vedea [21]) Daca exista f : [a; b] �! R; astfel incat f (xi) = yi; 8i =0; n si f 2 C2[a; b] atunci pentru functia spline poligonala de interpolare a functiei f avemestimarea erorii de aproximare:
jf (x)� s (x)j � kf 00kC8
� h2; 8x 2 [a; b]
unde h = maxfhi : i = 1; ng:
Probleme propuse:1. Sa se construiasca un algoritm pentru functia spline poligonala de interpolare al
c¼arei gra�c este linia poligonal¼a ce uneste punctele (x0; y0); (x1; y1):::; (xn; yn) si sa seprogrameze in limbajul C.2. La momentele 7.5 (adica ora 730), 10.5 (ora 1030), 13, 15.5 (ora 1530), 18, 21, 24 si 27
(adica ora 3 A. M. a doua zi) s-au masurat valorile glicemiei, obtinandu-se 130, 121, 128,96, 122, 138, 114, 90 (masurate in mg/100 ml). Sa se aproximeze glicemia acestui pacientde la orele 12, 14, 23 si 26 (a doua zi A. M. la ora 2) folosind functia spline poligonala deinterpolare. Sa se compare rezultatele cu cele obtinute prin algoritmul lui Akima.
1.6 Interpolare spline parabolica
Consideram o functie f : [a; b]! R; f 2 C1[a; b]; si diviziunea intervalului [a; b]
� : a = x0 < x1 < ::: < xn�1 < xn = b:
Deasemenea, �e y0; y1; :::; yn astfel incat yi = f (xi) ; 8i = 0; n si notam hi = xi � xi�1;8i = 1; n:Un spline parabolic de interpolare a functiei f : [a; b] ! R; este o functie spline
polinomiala de gradul 2 de interpolare a functiei f pe nodurile diviziunii�; s : [a; b] �! R;s 2 C1[a; b]: Vom nota prin si restrictiile functiei s la subintervalele [xi�1; xi]; i = 1; n alediviziunii �:Pentru functiile parabolice cubice se folosesc notatiile yi = s (xi) ; mi = s0 (xi) ; 8i =
0; n:Restrictiile functiei s la subintervalele [xi�1; xi]; i = 1; n sunt obtinute rezolvând prob-
lemele Cauchy:�s0i (x) = mi�1 +
mi�mi�1hi
� (x� xi�1) ; x 2 [xi�1; xi]si (xi�1) = yi�1
; i = 1; n: (1.13)
Se obtine,
si (x) =mi �mi�1
2hi� (x� xi�1)
2 +mi�1 � (x� xi�1) + yi�1; 8x 2 [xi�1; xi];8i = 1; n:(1.14)
Cerintele de netezime s 2 C[a; b]; s 2 C1[a; b] conduc la conditiile: si (xi) = yi;8i = 1; n:Se observ¼a c¼a s0i (xi�1) = mi�1; s
0i (xi) = mi,8i = 1; n: De aici, rezult¼a relatiile:
yi � yi�1 = (mi �mi�1) �hi2+mi�1hi
1.6. INTERPOLARE SPLINE PARABOLICA 21
care pot � puse în form¼a recurent¼a:
mi =2
hi� (yi � yi�1)�mi�1; 8i = 1; n: (1.15)
Din relatiile de recurenta (1.15) rezulta:Teorema: (a se vedea [21]) Fiind date valorile y0; y1; :::; yn sim0 exista o unica functie
spline parabolica pentru care,
s (xi) = yi; 8i = 0; n
si s0 (x0) = m0:In general valorile y0; y1; :::; yn sunt obtinute prin masuratori, iar m0 ramane liber.
Daca se cunoaste f 00 = f 0 (x0) atunci se poate lua m0 = f 00:
Algoritmul functiei spline parabolice generata de conditii initialeI. Date de intrare: a; b capetele intervalului
n numarul de subintervale ale diviziuniix[i]; i = 0; n niodurile de interpolarey[i]; i = 0; n valorile pe noduri� 2 (a; b) punctul in care se calculeaza functia splinef 00 valoarea derivatei intai pe primul nod
II. Date de iesire: s (�)III. Pasii algoritmului1. Pentru i = 1; n calculeaza h[i] = x[i]� x[i� 1]2. Fie m[0] := f 003. Pentru i = 1; n calculeaza
m[i] =2
h[i]� (y[i]� y[i� 1])�m[i� 1]
4. Se determina indicele j 2 f1; :::; ng pentru care
x[j � 1] < � < x[j]
5. Calculeaza
s (�) =m[j]�m[j � 1]
2h[j]� (�� x[j � 1])2 +m[j � 1] � (�� x[j � 1]) + y[j � 1]
6. Tipareste s (�) : STOP.
Probleme propuse1. Se se construiasca un program in limbajul C pentru algoritmul functiei spline
parabolice generata de conditii initiale.2. La momentele 7.5 (adica ora 730), 10.5 (ora 1030), 13, 15.5 (ora 1530), 18, 21, 24 si
27 (adica ora 3 A. M. a doua zi) s-au masurat valorile glicemiei, obtinandu-se 130, 121,128, 96, 122, 138, 114, 90 (masurate in mg/dl). Sa se aproximeze glicemia acestui pacientde la orele 12, 14, 23 si 26 (a doua zi A. M. la ora 2) folosind algoritmul de la problemaprecedenta. Sa se compare rezultatele cu cele obtinute prin algoritmul lui Akima si prininterpolare spline poligonala. Se da f 00 = �2:
22 CAPITOLUL 1. METODE DE INTERPOLARE A FUNCTIILOR CONTINUE
1.7 Interpolare spline cubica
1.7.1 Spline cubic complet general de conditii bilocale
Consideram o functie f : [a; b]! R; f 2 C2[a; b]; si diviziunea intervalului [a; b]
� : a = x0 < x1 < ::: < xn�1 < xn = b:
Deasemenea, �e y0; y1; :::; yn astfel incat yi = f (xi) ; 8i = 0; n si notam hi = xi � xi�1;8i = 1; n:Un spline cubic pentru f : [a; b]! R; este o functie spline polinomiala de gradul 3 de
interpolare a functiei f pe nodurile diviziunii �; s : [a; b] �! R; s 2 C2[a; b]: Vom notaprin si restrictiile functiei s la subintervalele [xi�1; xi]; i = 1; n ale diviziunii �:Pentru functiile spline cubice se folosesc notatiile yi = s (xi) ; mi = s0 (xi) ; Mi =
s00 (xi) ; 8i = 0; n:Deoarece s00i ; i = 1; n sunt functii polinomiale de gradul intai, in baza notatiilor de mai
sus se vor rezolva problemele bilocale8<: s00i (x) =Mi�(x�xi�1)+Mi�1�(xi�x)
hi; x 2 [xi�1; xi]
si (xi�1) = yi�1si (xi) = yi; 8i = 1; n
si se obtine
si (x) =
"(x� xi�1)
3
6hi� hi (x� xi�1)
6
#�Mi +
"(xi � x)3
6hi� hi (xi � x)
6
#�Mi�1+ (1.16)
+xi � x
hi� yi�1 +
x� xi�1hi
� yi; 8x 2 [xi�1; xi]; 8i = 1; n:
Pentru determinarea coe�cientilor Mi; i = 0; n; din conditia s 2 C1[a; b] rezulta cerinteles0i (xi) = s0i (xi+1) ; 8i = 1; n� 1: Acestea conduc la un sistem liniar de n-1 ecuatii cu n+1necunoscute Mi; i = 0; n
hiMi�1
6+hi + hi+1
3�Mi +
hi+1Mi
6=yi+1 � yihi+1
� yi � yi�1hi
; 8i = 1; n� 1:
La acestea se adauga inca doua conditii M0 =Mn = 0; sau s0 (x0) = f0; s0 (xn) = fn (inipoteza ca sunt cunoscute f 0 (x0) = f0 si f 0 (xn) = fn). Adaugand de exemplu, ultimeledoua conditii se obtine sistemul8><>:
hiMi�16
+ hi+hi+13
�Mi +hi+1Mi
6= yi+1�yi
hi+1� yi�yi�1
hi; 8i = 1; n� 1
h1M0
3+ h1M1
6= y1�y0
h1� f0
hnMn�16
+ hnMn
3= fn � yn�yn�1
hn
(1.17)
din care se determina in mod unic M0; M1; :::; Mn: Observam ca matricea acestui sistemeste tridiagonala diagonal dominanta si de aceea este nesingulara.In [14] se arata ca functia spline cubica (1.16) ce indeplineste oricare din cele doua
conditii suplimentare de mai sus are proprietatile :
1.7. INTERPOLARE SPLINE CUBICA 23
(i) Minimizarea functionalei
J (g) =
bZa
[g00 (x)]2dx
adicaJ (s) = minfJ (g) : g 2 C2[a; b]; g (xi) = yi; 8i = 0; ng:
(ii) Aproximare uniforma pentru orice functie f 2 C2[a; b]; cu estimarea globala aerorii
kf 0 � s0k �pb� a � kf 00kC �
ph
kf � sk �pb� a � kf 00kC � h
ph;
unde h = maxfhi : i = 1; ng:Pentru functia spline cubica (1.16) veri�cand conditiile s0 (x0) = f0; s0 (xn) = fn se
construieste urmatorul algoritm.
Algoritmul functiei spline cubice complete generata de conditii bilocaleI. Date de intrare:nodurile xi; i = 0; nvalorile yi = f (xi) ; i = 0; nvalorile derivatei pe nodurile extreme, f0; fn� 2 [a; b]; � 6= xi; 8i = 0; n; punctul in care se aproximeaza f (�) prin s (�)II. Date de iesire: s (�) :III. Pasii algoritmuluiPasul 1. Se determina indicele j 2 f1; :::; ng pentru care
xj�1 < � < xj:
Pasul 2. Executaa0 := 2; c0 := 1; bn := 1
Pentru i = 1; n calculeaza
hi := xi � xi�1:
ai := 2
Pasul 3. Executa
d0 :=6
h1��y1 � y0h1
� f0
�; dn =
6
hn��fn �
yn � yn�1hn
�Pentru i = 1; n� 1 calculeaza
bi :=hi
hi + hi+1
ci = 1� bi
24 CAPITOLUL 1. METODE DE INTERPOLARE A FUNCTIILOR CONTINUE
di :=6
hi + hi+1��yi+1 � yihi+1
� yi � yi�1hi
�Pasul 4. Executa
�0 :=c0a0
Pentru i = 1; n� 1 calculeaza
!i := ai � bi � �i�1
�i :=ci!i
Executa!n := an � bn � �n�1
Pasul 5. Executa
z0 :=d02
Pentru i = 1; n calculeaza
zi :=di � bi � zi�1
!i
Pasul 6. Determinarea solutiei M0; M1; :::; Mn a sistemului (1.17)Executa
Mn := zn
Pentru i = n� 1; 0 calculeaza
Mi := zi � �i �Mi+1:
Pasul 7. Calculeaza
s (�) :=
"(�� xj�1)
3
6hj� hj (�� xj�1)
6
#�Mj +
"(xj � �)3
6hj� hi (xj � �)
6
#�Mj�1+
+xj � �
hj� yj�1 +
�� xj�1hj
� yj
Pasul 8. Tipareste s (�) : STOP.
Probleme propuse1. Se se construiasca un program in limbajul C pentru algoritmul functiei spline cubice
generata de conditii bilocale.2. La momentele 7.5 (adica ora 730), 10.5 (ora 1030), 13, 15.5 (ora 1530), 18, 21, 24 si
27 (adica ora 3 A. M. a doua zi) s-au masurat valorile glicemiei, obtinandu-se 130, 121,128, 96, 122, 138, 114, 90 (masurate in mg/dl). Sa se aproximeze glicemia acestui pacientde la orele 12, 14, 23 si 26 (a doua zi A. M. la ora 2) folosind algoritmul de la problemaprecedenta. Sa se compare rezultatele cu cele obtinute prin algoritmul lui Akima si prininterpolare spline poligonala. Se dau f0 = �2; f7 = �3:
1.7. INTERPOLARE SPLINE CUBICA 25
1.7.2 Spline cubic natural generat de conditii initiale
Consideram o functie f : [a; b]! R; f 2 C2[a; b]; si diviziunea intervalului [a; b]
� : a = x0 < x1 < ::: < xn�1 < xn = b:
Deasemenea, �e y0; y1; :::; yn astfel incat yi = f (xi) ; 8i = 0; n si notam hi = xi � xi�1;8i = 1; n:Fie functia spline cubica de interpolare a functiei f pe nodurile diviziunii �; s :
[a; b] �! R; s 2 C2[a; b]: Vom nota prin si restrictiile functiei s la subintervalele [xi�1; xi];i = 1; n ale diviziunii �: In [15] se construieste o functie spline cubica de interpolaregenerata de conditii initiale. Astfel, restrictia sa la subintervalul [xi�1; xi]; i = 1; n estedeterminata prin rezolvarea problemei Cauchy8<:
s00i (x) =Mi�1 +1hi(Mi �Mi�1)(x� xi�1)
si(xi�1) = yi�1s0i(xi�1) = mi�1;
unde Mi = s00i (xi) si si = s j[xi�1;xi] :Se obtine
si(x) =1
6hi(Mi�Mi�1)(x�xi�1)3+
Mi�1
2(x�xi�1)2+mi�1(x�xi�1)+yi�1; 8x 2 [xi�1; xi]:
(1.18)Din conditia s 2 C2[a; b] si conditiile pe noduri : s(xi) = yi; s
0(xi) = mi; i = 1; n; seobtin relatiile, (
Mi + 2Mi�1 = 6 � yi�yi�1�mi�1�hih2i
Mi +Mi�1 =2(mi�mi�1)
hi; i = 1; n:
(1.19)
care sunt echivalente cu oricare din urm¼atoarele dou¼a sisteme de relatii :�mi�1 =
1hi� (yi � yi�1)� 1
6hi(Mi + 2Mi�1)
mi =1hi� (yi � yi�1) +
16hi(4Mi +Mi�1)
; i = 1; n (1.20)
si (Mi =
6h2i� (yi � yi�1)� 6mi�1
hi� 2Mi�1
mi =3hi� (yi � yi�1)� 2mi�1 � Mi�1hi
2
; i = 1; n; (1.21)
unde hi = xi � xi�1; 8i = 1; n:Prin utilizarea relatiilor de recurent¼a (1.21) se obtine in mod unic functia spline cubica
de interpolare pornind de la valorile y0; y1; :::; yn si de la valorile primelor doua derivate peprimul nod, m0 si M0: Estimarea erorii in formula corespunzatoare de interpolare splinepoate �gasita in [16] pentru cazurile f 2 Lip[a; b]; respectiv f 2 C1[a; b] cu f 0 2 Lip[a; b]:Observatie: Valorile m0 si M0 sunt in general libere (experimental, de obicei se
determina y0; y1; :::; yn) si pot � date astfel incat functia spline corespunzatoare sa aibaproprietati suplimentare. Astfel, daca M0 = Mn = 0 atunci functia spline cubica cu
26 CAPITOLUL 1. METODE DE INTERPOLARE A FUNCTIILOR CONTINUE
restrictiile la subintervalele diviziunii � date in (1.18) va avea proprietatea de minimizarea functionalei
J(y) =
bZa
[y00(x)]2dx (1.22)
pe multimea functiilor de clas¼a C2[a; b] ce interpoleaz¼a f pe nodurile diviziunii �n: Adica,se obtine o functie spline cubica naturala. Deasemenea, s si s0 aproximeaz¼a pe f sirespectiv, f 0 , cu o eroare estimat¼a global astfel :
kf � sk �pb� a � kf 00kC � h
32
kf 0 � s0k �pb� a � kf 00kC �
ph: (1.23)
Pentru a avea Mn = 0 se poate obtine m0 astfel:Se de�nesc a1 = �2; c1 = � 6
h1; b1 =
3h1� (y1 � y0); d1 =
6h21� (y1 � y0) si recurent se
calculeaza pentru i = 2; n
ai = �2ai�1 �hici�12
bi =3(yi � yi�1)
hi� 2bi�1 �
hidi�12
(1.24)
si
ci = �2ci�1 �6ai�1hi
di =6(yi � yi�1)
h2i� 6bi�1
hi� 2di�1: (1.25)
Punand
M0 = 0 si m0 = �dncn
(1.26)
se va obtine Mn = 0:
Algoritmul functiei spline cubice naturale generata de conditii initialeI. Date de intrare:nodurile xi; i = 0; nvalorile yi; i = 0; nm0 si M0 , valorile primelor doua derivate pe primul nod� 2 [a; b]; � 6= xi; 8i = 0; n; punctul in care se aproximeaza f (�) prin s (�)II. Date de iesire: s (�) :III. Pasii algoritmuluiPasul 1. Se determina indicele j 2 f1; :::; ng pentru care
xj�1 < � < xj:
Pasul 2. Pentru i = 1; n calculeaza
hi := xi � xi�1
1.7. INTERPOLARE SPLINE CUBICA 27
Pasul 3.Pentru i = 1; n calculeaza
mi :=3
hi� (yi � yi�1)� 2mi�1 �
Mi�1hi2
Mi :=6
h2i� (yi � yi�1)�
6mi�1
hi� 2Mi�1
Pasul 4. Calculeaza
s (�) :=1
6hj(Mj �Mj�1)(�� xj�1)
3 +Mj
2(�� xj�1)
2 +mj(�� xj�1) + yj�1:
Pasul 5. Tipareste s (�) : STOP.Observatie: Daca valorile m0 si M0 nu sunt date atunci se poate intercala intre pasii
2 si 3 subrutina urmatoare:Pasul 2.1. Executa
M0 := 0; a1 := �2; b1 :=3
h1� (y1 � y0)
c1 := �6
h1; d1 :=
6
h21� (y1 � y0):
Pasul 2.2. Pentru i = 2; n calculeaza
ai := �2ai�1 �hici�12
; ci := �2ci�1 �6ai�1hi
bi :=3(yi � yi�1)
hi� 2bi�1 �
hidi�12
di :=6(yi � yi�1)
h2i� 6bi�1
hi� 2di�1:
Pasul 2.3. Executa
m0 := �dncn:
Prin acesta subrutina se vor determina m0 si M0 astfel incat functia spline cubica core-spunzatoare sa minimizeze functionala (1.22) si sa aiba estimarea erorii (1.23).
Probleme propuse1. Sa se construiasca un program in limbajul C pentru algoritmul functiei spline cubice
naturale generata de conditii initiale.2. La momentele 7.5 (adica ora 730), 10.5 (ora 1030), 13, 15.5 (ora 1530), 18, 21, 24 si
27 (adica ora 3 A. M. a doua zi) s-au masurat valorile glicemiei, obtinandu-se 130, 121,128, 96, 122, 138, 114, 90 (masurate in mg/dl). Sa se aproximeze glicemia acestui pacientde la orele 12, 14, 23 si 26 (a doua zi A. M. la ora 2) folosind algoritmul de interpolarespline cubic generat de conditii initiale cu m0 siM0 date in (1.26) prin folosirea subrutineide mai sus. Sa se compare rezultatele cu cele obtinute prin algoritmul lui Akima si prinalgoritmul de interpolare spline cubic generat de conditii bilocale.
28 CAPITOLUL 1. METODE DE INTERPOLARE A FUNCTIILOR CONTINUE
Capitolul 2
Ajustarea in medie patratica adatelor experimentale
Atunci cand valorile unei functii pe nodurile xi; i = 0; n pot � obtinute doar cu aproxi-matie, problema interpolarii acestor valori îsi pierde justi�carea. O alternativa la aceastasituatie este obtinerea polinomului, de un anumit grad �xat, situat cel mai aproape deaceste valori. Un asemenea polinom se poate determina utilizand principiul celor mai micipatrate descoperit de K. F. Gauss la inceputul secolului XIX.Concret, �ind date nodurile xi; i = 0; n si valorile unei functii obtinute experimental
(prin masurare) yi; i = 0; n se cere sa se determine un polinom de grad p, p < n,
g (x) = apxp + :::+ a1x+ a0
pentru care reziduul (de variabilele a0; :::; ap),
R (a0; :::; ap) =nXi=0
[yi � apxpi � :::� a1xi � a0]
2
se minimizeaza. Determinarea acestui polinom va � astfel o aplicatie a metodei luiSylvester pentru punctele de extrem local ale functiilor de mai multe variabile (a se vedea[25]). Polinomul obtinut se numeste polinomul de regresie. Cel mai des folosite in de-ducerea empirica experimentala a legilor naturii sunt dreapta (utilizata si in psihologie,sociologie, economie, business) si parabola patratica de regresie (utilizata si in economie,business). Deasemenea, principiul celor mai mici patrate poate � folosit si la constructiafunctiilor spline de ajustare a datelor experimentale.
2.1 Metoda celor mai mici patrate: dreapta de re-gresie
Dorind sa stabilim legea de dependenta y = f (x) dintre doua marimi x si y se efectueazamasuratorile exprimate in tabelul
x x0 x1 ...... xny y0 y1 ...... yn:
29
30 CAPITOLUL 2. AJUSTAREA IN MEDIE PATRATICA
Punctele (x0; y0); (x1; y1):::; (xn; yn) se reprezinta intr-un sistem de axe XOY, iar imag-inea lor se numeste nor statistic. Daca forma norului statistic este asemanatoare uneidrepte atunci se poate presupune ca legea de dependenta dintre x si y este liniara. Astfel,se va cauta dreapta cea mai apropiata de punctele norului statistic. Fie ecuatia cartezianaa acestei drepte,
y = ax+ b
cu coe�cientii a si b necunoscuti. Consideram functia g (x) = ax + b si vom determinacoe�cientii ce minimizeaza reziduul (abatere patratica medie)
R (a; b) =
nXi=0
[f(xi)� g (xi)]2 =
nXi=0
[yi � axi � b]2:
In acest scop se considera sistemul8>>>>><>>>>>:@R@a= �2
nXi=0
(yi � axi � b) � xi = 0
@R@b= �2
nXi=0
(yi � axi � b) = 0
()
8>>><>>>:a �
nXi=0
x2i + b �nXi=0
xi =nXi=0
xiyi
a �nXi=0
xi + (n+ 1) b =nXi=0
yi:
(2.1)
Totodata,@2R
@a2= 2
nXi=0
x2i > 0
si
� =
���� @2R@a2
@2R@a@b
@2R@a@b
@2R@b2
���� = 4 (n+ 1)224 nXi=0
1
n+ 1� x2i �
nXi=0
1
n+ 1� xi
!235 > 0;in virtutea inegalitatii lui Jensen. Atunci, determinantul sistemului (2.1) este strict poz-itiv si solutia sistemului (2.1) este unicul punct de minim local pentru reziduul R (a; b) :Rezolvand acest sistem se determina coe�cientii din ecuatia dreptei de regresie.
Algoritmul dreptei de regresieI. Date de intrare:n numarul datelorxi; i = 0; n, valorile marimii xyi; i = 0; n, valorile marimii yII. Date de iesire:a; b coe�cientii din ecuatia carteziana a dreptei de regresieIII. Pasii algoritmului:
2.2. METODA CELOR MAI MICI PATRATE: PARABOLA DE REGRESIE 31
Pasul 1. Calculeaza
d := (n+ 1) �nXi=0
x2i �
nXi=0
xi
!2
d1 := (n+ 1) �nXi=0
xiyi �nXi=0
xi �nXi=0
yi
d2 :=
nXi=0
x2i �nXi=0
yi �nXi=0
xiyi �nXi=0
xi
Pasul 2. Calculeaza
a :=d1d
b :=d2d
Pasul 3. Tipareste a; b: STOP.
Probleme propuse:1. Sa se construiasca un program in limbajele C si Pascal pentru algoritmul dreptei
de regresie.2. (din [25]) Fie tabelul de valori ale functiei y = f (x)x 1 3 4 6 8 9y 1 2 4 4 5 3Sa se determine dreapta de regresie de ecuatie y = ax+b ce ajusteaza datele din tabel.
( Se gasesc a = 95281si b = 399
281).
2.2 Metoda celor mai mici patrate: parabola de re-gresie
Dorind sa stabilim legea de dependenta y = f (x) dintre doua marimi x si y se efectueazamasuratorile exprimate in tabelul
x x0 x1 ...... xny y0 y1 ...... yn:
Punctele (x0; y0); (x1; y1):::; (xn; yn) se reprezinta intr-un sistem de axe XOY, iar imag-inea lor se numeste nor statistic. Daca forma norului statistic este asemanatoare uneiparabole patratice atunci se poate presupune ca legea de dependenta dintre x si y este deforma
y = ax2 + bx+ c
cu coe�cientii a; b si c necunoscuti. Consideram functia g (x) = ax2 + bx + c si vomdetermina coe�cientii ce minimizeaza reziduul (abatere patratica medie)
R (a; b; c) =nXi=0
[f(xi)� g (xi)]2 =
nXi=0
[yi � ax2i � bxi + c]2:
32 CAPITOLUL 2. AJUSTAREA IN MEDIE PATRATICA
In acest scop se considera sistemul8>>>>>>>>>>>><>>>>>>>>>>>>:
@R@a= �2
nXi=0
[yi � ax2i � bxi + c] � x2i = 0
@R@b= �2
nXi=0
[yi � ax2i � bxi + c] � xi = 0
@R@c= �2
nXi=0
[yi � ax2i � bxi + c] = 0
()
8>>>>>>>><>>>>>>>>:
a
nXi=0
x4i + b
nXi=0
x3i + c
nXi=0
x2i =
nXi=0
x2i � yi
a
nXi=0
x3i + b
nXi=0
x2i + c
nXi=0
xi =
nXi=0
xiyi
a
nXi=0
x2i + bnXi=0
xi + c (n+ 1) =nXi=0
yi:
(2.2)
Totodata,@2R
@a2= 2
nXi=0
x4i > 0;
� =
���� @2R@a2
@2R@a@b
@2R@a@b
@2R@b2
���� = 424 nXi=0
x2i �nXi=0
x4i �
nXi=0
x3i
!235 > 0si
HR (a; b; c) =
������@2R@a2
@2R@a@b
@2R@a@c
@2R@a@b
@2R@b2
@2R@b@c
@2R@a@c
@2R@b@c
@2R@c2
������ = 8nXi=0
x4i �
24(n+ 1) � nXi=0
x2i �
nXi=0
xi
!235+
+8
nXi=0
x2i �
24 nXi=0
xi
!�
nXi=0
x3i
!�
nXi=0
x2i
!235++8
nXi=0
x3i �"
nXi=0
xi
!�
nXi=0
x2i
!� (n+ 1) �
nXi=0
x3i
#> 0:
Atunci, determinantul sistemului (2.2) este strict pozitiv si solutia sistemului (2.2) esteunicul punct de minim local pentru reziduul R (a; b; c) : Rezolvand acest sistem se deter-mina coe�cientii din ecuatia carteziana a parabolei de regresie,
y = ax2 + bx+ c:
Algoritmul parabolei de regresie
2.2. METODA CELOR MAI MICI PATRATE: PARABOLA DE REGRESIE 33
I. Date de intrare:n numarul datelorxi; i = 0; n, valorile marimii xyi; i = 0; n, valorile marimii yII. Date de iesire:a; b; c coe�cientii din ecuatia carteziana a parabolei de regresieIII. Pasii algoritmului:Pasul 1. Calculeaza
d := (n+ 1)
nXi=0
x2i
!�
nXi=0
x4i
!+ 2
nXi=0
xi
!�
nXi=0
x2i
!�
nXi=0
x3i
!�
�
nXi=0
x2i
!3�
nXi=0
xi
!2�
nXi=0
x4i
!� (n+ 1) �
nXi=0
x3i
!2
d1 := (n+ 1)
nXi=0
x2i
!�
nXi=0
x2i � yi
!+
nXi=0
xi
!�
nXi=0
x2i
!�
nXi=0
xiyi
!+
+
nXi=0
x3i
!�
nXi=0
xi
!�
nXi=0
yi
!�
nXi=0
x2i
!2�
nXi=0
yi
!�
�
nXi=0
xi
!2�
nXi=0
x2i � yi
!� (n+ 1)
nXi=0
x3i
!�
nXi=0
xiyi
!
d2 := (n+ 1) �
nXi=0
xiyi
!�
nXi=0
x4i
!+
nXi=0
x2i
!�
nXi=0
x3i
!�
nXi=0
yi
!+
+
nXi=0
xi
!�
nXi=0
x2i
!�
nXi=0
x2i � yi
!�
nXi=0
x2i
!2�
nXi=0
xiyi
!�
�
nXi=0
x4i
!�
nXi=0
xi
!�
nXi=0
yi
!� (n+ 1)
nXi=0
x2i � yi
!�
nXi=0
x3i
!
d3 :=
nXi=0
x4i
!�
nXi=0
x2i
!�
nXi=0
yi
!+
nXi=0
xi
!�
nXi=0
x2i � yi
!�
nXi=0
x3i
!+
+
nXi=0
x2i
!�
nXi=0
x3i
!�
nXi=0
xiyi
!�
nXi=0
x2i
!2�
nXi=0
x2i � yi
!�
�
nXi=0
xi
!�
nXi=0
xiyi
!�
nXi=0
x4i
!�
nXi=0
x3i
!2�
nXi=0
yi
!:
Pasul 2. Calculeaza
a :=d1d
34 CAPITOLUL 2. AJUSTAREA IN MEDIE PATRATICA
b :=d2d
c :=d3d
Pasul 3. Tipareste a; b; c: STOP.Probleme propuse:1. Sa se scrie un program in limbajul C pentru parabola de regresie.2. Fie tabelul de valori ale functiei y = f (x)x -2 -1 1 1.5 2 2.5 3 3.5 4 5y 3 1 6 8 9 13 16 21 25 36Sa se determine parabola de regresie de ecuatie y = ax2 + bx + c ce ajusteaza datele
din tabel.
Capitolul 3
Formule de cuadratura
Pentru o functie continua f : [a; b] �! R, cand primitiva sa se poate calcula sub forma�nita, formula lui Leibniz-Newton permite calculul integralei de�nite
bZa
f (x) dx = F (b)� F (a):
Cand primitiva F nu se poate determina sub forma �nita, atunci integrala de�nita
bZa
f (x) dx
se poate calcula doar aproximativ, motivand astfel studiul si utilizarea formulelor decuadratura numerica.In general, o formula de cuadratura este de forma
bZa
f (x) dx =nXk=1
Ak � f (xk) +R(f)
in care Ak 2 R se numesc coe�cientii, iar xk 2 [a; b]; k = 1; n nodurile formuleide cuadratura. Estimarea erorii comise se poate realiza stabilind o margine superioarapentru restul jR(f)j al formulei de cuadratura.Daca o formula de cuadratura se poate obtine prin integrarea pe intervalul [a,b] a
unei formule de interpolare
f(x) = P (x) +R(x); x 2 [a; b]atunci se numeste formula de cuadratura de tip interpolator.Formula trapezului, formula lui Simpson si in general formulele lui Newton-Cotes sunt
formule de cuadratura de tip interpolator( a se vedea [4], [11] si [24]).
3.1 Motivatia formulelor de cuadratura. Exemple defunctii cu primitiva in forma de serie in�nita
Printre exemplele de functii continue (chiar analitice) ce au primitiva inexprimabila subforma �nita, putem mentiona densitatea de probabilitate normala a lui Gauss pe un
35
36 CAPITOLUL 3. FORMULE DE CUADRATURA
interval compact si functiile ce conduc la integralele eliptice in sensul Legendre.Densitatea de probabilitate normala cu media zero si dispersia � = 1;
f : R �! R, f(x) =1p2�� e�x2
2
conduce la functia integrala a lui Laplace,
�(x) =1p2�
xZ0
e�t2
2 dt; x > 0
in care integrala se poate calcula doar aproximativ. Valorile functiei lui Laplace pentrux 2 [0; 4] sunt aproximate in 401 noduri echidistante din intervalul [0; 4] si se gasesctabelate in tratatele de calculul probabilitatilor.Consideram elipsa de ecuatie implicita
(E) :x2
a2+y2
b2� 1 = 0; a; b > 0; a > b
cu excentricitatea � =pa2�b2a
: Ecuatiile parametrice ale elipsei sunt
(E) :
�x = a cos ty = b sin t
; t 2 [0; 2�]:
Incercand sa calculam lungimea arcului de elipsa cuprins in cadranul I putem folosi atatecuatiile parametrice la calculul integralei
l
�E
4
�=
�2Z0
q[x0(t)]2 + [y0(t)]2dt
cat si ecuatia explicita y = ba
pa2 � x2; x 2 [0; a];
la calculul integralei
l
�E
4
�=
aZ0
q1 + [y0(x)]2dx:
Ambele integrale conduc la integrala eliptica
l
�E
4
�= a
�2Z0
p1� �2 sin2 tdt:
Spre exemplu,
l
�E
4
�=
aZ0
q1 + [y0(x)]2dx =
1
a
aZ0
ra4 + (b2 � a2)x2
a2 � x2dx;
3.2. FORMULA TRAPEZULUI 37
cu substitutia x = a sin t; t 2 [0; �2] se obtine
l
�E
4
�=1
a
�2Z0
ra2 + (b2 � a2) sin2 t
cos2 t� a cos tdt =
=
�2Z0
qa2 + (b2 � a2) sin2 tdt =
�2Z0
a
r1� (a
2 � b2)
a2sin2 tdt =
= a
�2Z0
p1� �2 sin2 tdt:
Integralele eliptice de genul I, II si III in sensul lui Legendre sunt respectiv,
F (k; ') =
'Z0
dtp1� k2 sin2 t
E(k; ') =
'Z0
p1� k2 sin2 tdt
'Z0
dt�1 + h sin2 t
�p1� k2 sin2 t
;
pentru k 2 (0; 1); h 2 R si ' > 0 constante.In secolul al XIX-lea J. Liouville a aratat ca pentru nici o integrala eliptica in sensul
lui Legendre primitiva nu se poate obtine sub forma �nita (este suma unei serii de put-eri). Deci, integralele eliptice nu se pot calcula exact. Spre exemplu, lungimea orbiteiPamantului in jurul Soarelui se determina prin calculul aproximativ al integralei
l(P ) = 4aE(�;�
2) = 4a
�2Z0
p1� �2 sin2 tdt;
unde � = 0; 016729 si a = 149; 6� 106 km.
3.2 Formula trapezului
Dup¼a cum se stie, pentru o functie continu¼a f : [a; b] �! R; are loc egalitatea,
bZa
f(x)dx =b� a
2� [f(a) + f(b)] +R(f):
38 CAPITOLUL 3. FORMULE DE CUADRATURA
Pentru valoarea absolut¼a a restului R(f) se pot stabili margini superioare. Corespunz¼atorcu propriet¼atile functiei f aceste margini superioare pot � mai apropiate sau mai înde-p¼artate de jR(f)j. Astfel, dac¼a f 2 C[a; b] atunci estimarea se face cu ajutorul modululuide continuitate,
jR(f)j � b� a
2�$(f ; b� a
2):
În [9] s-au obtinut recent urm¼atoarele estim¼ari :
jR(f)j �
8><>:L(b�a)2
4; daca f 2 Lip[a; b]
(b�a)2kf 0k4
; daca f 2 C1[a; b]:
Pentru f 2 C2[a; b] este cunoscuta estimarea clasic¼a,
jR(f)j � (b� a)3 kf 00k12
:
Pentru functii cu derivata Lipschitz are loc:Teorema: (in [6]) Dac¼a f 2 C1[a; b] cu f 0 2 Lip[a; b], iar L0 > 0 este constanta
Lipschitz pentru f 0, atunci
jR(f)j � (b� a)3L0
12:
Considerând o diviziune a intervalului [a; b];
� : a = x0 < x1 < ::: < xn�1 < xn = b
not¼am hi = xi+1 � xi; 8i = 0; n� 1; si
In(f;�) =1
2
n�1Xi=0
hi[f (xi) + f(xi+1)]:
Corespunz¼ator diviziunii de mai sus se obtine :Teorema: (in [6] si [9]) Formula trapezului corespunz¼atoare diviziunii � este
bZa
f(x)dx = In(f;�) +Rn(f;�);
iar pentru estimarea restului avem
jRn(f;�)j �
8>>>>>>>>>>>>>><>>>>>>>>>>>>>>:
12
n�1Pi=0
hi$(f ;12hi); dac¼a f 2 C[a; b]
L4
n�1Pi=0
h2i ; dac¼a f 2 Lip[a; b]
kf 0k4
n�1Pi=0
h2i ; dac¼a f 2 C1[a; b]
L0
12
n�1Pi=0
h3i ; dac¼a f 0 2 Lip[a; b]
kf 00k12
n�1Pi=0
h3i ; dac¼a f 2 C2[a; b]:
(3.1)
3.2. FORMULA TRAPEZULUI 39
Corolar: (in [6] si [9]) Dac¼a diviziunea � este echidistant¼a atunci,
bZa
f(x)dx =b� a
2n� [f(a) + 2
n�1Xi=1
f(a+i(b� a)
n) + f(b)] +Rn(f) (3.2)
si
jRn(f)j �
8>>>>><>>>>>:
b�a2�$(f ; b�a
2n); dac¼a f 2 C[a; b]
(b�a)2L4n
; dac¼a f 2 Lip[a; b](b�a)2�kf 0k
4n; dac¼a f 2 C1[a; b]
(b�a)3L012n2
; dac¼a f 0 2 Lip[a; b](b�a)3kf 00k
12n2; dac¼a f 2 C2[a; b]:
(3.3)
Observatie: Formula (3.2) se poate obtine si integrand formula de interpolare splinepoligonala cu noduri echidistante.Fie functia integrabila f : [a; b] �! R cu valori cunoscute pe noduri echidistante din
intervalul [a,b].De�nitia 1. Numim formula de cuadratura generalizata a trapezelor formula de
cuadratura obtinuta prin aplicarea formulei de cuadratura a trapezului pe �ecare subin-terval
[a+ (k � 1)h; a+ kh]; k = 1; n;
unde h = b�an:
1. Formula de cuadratura a trapezelor este
bZa
f (x) dx =b� a
2n
"f (a) + 2
m�1Xk=1
f (a+ kh) + f(b)
#+Rn (f) :
2. Daca f 2 C2[a; b] atunci restul in formula de cuadratura a trapezelor are expresia
Rn (f) = �(b� a)3
12n2f 00(�) = �(b� a)h2
12f 00(�); a < � < b;
de unde
jRn (f)j �(b� a)3
12n2M2(f) =
(b� a)h2
12M2(f):
3. Formula de cuadratura a trapezelor are gradul de exactitate 1.4. (Metoda lui Romberg) Daca se considera in formula de cuadratura a trapezelor
n = 2p; p = 0; 1; 2; :::
si se noteaza prin
hp =b� a
n=b� a
2p; p = 0; 1; 2; :::
Sp = f (a) + 2n�1Xk=1
f (ak) + f(b); ak = a+ khp; k = 0; n;
40 CAPITOLUL 3. FORMULE DE CUADRATURA
cuS0 = f (a) + f(b);
atunci urmatoarele formule de recurenta au loc
hp+1 =hp2;
Sp+1 = Sp + 22pXk=1
f (a+ (2k � 1)hp+1) ; p = 0; 1; :::
decibZa
f (x) dx � hp2Sp:
Algoritmul lui Romberg pentru cuadratura trapezului
Algoritmul calculeaza valoarea aproximativa a integralei
bZa
f (x) dx prin formula repetata
a trapezelor, cu o precizie data �:Intrari: a,b-limitele de integraref-expresia functiei��precizia dorita
Iesiri:
bZa
f (x) dx�valoarea aproximativa
Pasul 1: (Initializari)p := 0; hp := b� a; Sp :=
hp2[f (a) + f(b)] ; np := 1
Pasul 2: p := p+ 1Pasul 3: hp :=
hp�12
Pasul 4: sum:=0Pasul 5: Pentru i = 1; np�1 executasum:=sum+f (a+ (2i� 1)hp)Pasul 6: sum:=hp�sumPasul 7: Sp :=
Sp�12+sum
Pasul 8: np := 2np�1Pasul 9: Daca jSp � Sp�1j � � atuncisalt la pasul al doileaPasul 10: Tipareste Sp: STOP.
Algoritmul cuadraturii clasice a trapezului
I. Date de intrare:a, b: capetele intervaluluin: numarul de subintervale al diviziuniif : expresia functiei de sub integrala
3.2. FORMULA TRAPEZULUI 41
II. Date de iesire:
bZa
f (x) dx�valoarea aproximativa a acestei integrale.
III. Pasii algoritmului:1. Se citesc datele: a,b,n2. Se creaza procedura de introducere a functiei f3. Se construiesc nodurile diviziunii:
Pentru i = 0; n calculeaza
x[i] = a+i (b� a)
n
4. Calculul sumei integrale:S := 0Pentru i = 1; n calculeaza
S := S + (f (x[i]) + f (x[i� 1]))
5. Calculul aproximatiei integralei:
T :=(b� a)
2n� S
6. A�sarea rezultatului:Tipareste T . STOP
Probleme propuse:1. Sa se scrie un program in limbajul C pentru algoritmul cuadraturii clasice a trapezu-
lui.2. Sa se utilizeze formula de cuadratura a trapezului pentru n=10 la calculul integraleiZ 1
0
1
x+ 1dx
si sa se compare rezultatul cu valoarea exacta a acestei integrale ln 2:3. Sa se utilizeze formula de cuadratura a trapezului pentru n=10 la calculul integraleiZ 1
0
1
x2 + 1dx
si sa se compare rezultatul cu valoarea exacta a acestei integrale �4:
4. Sa se utilizeze formula de cuadratura a trapezului pentru n=10 la calculul integraleiZ 1:5
0
e�t2
2 dt:
si a integralei (constanta lui Catalan)Z 1
0
arctanx
xdx
42 CAPITOLUL 3. FORMULE DE CUADRATURA
folosind faptul ca
limx!0
arctanx
x= 1:
5. Sa se utilizeze formula de cuadratura a trapezului pentru n=20 la calculul integralei
�2Z0
p1� �2 sin2 tdt
cu � = 0; 016729 si apoi sa se calculeze o aproximatie pentru lungimea orbitei Pamantuluiin jurul Soarelui
L (P ) = 4a
�2Z0
p1� �2 sin2 tdt;
unde a = 149; 6� 106 km.
3.3 Formule de tip Euler-Mac Laurin. Formula per-turbata a trapezului
Uneori se utilizeaz¼a urm¼atoarea formul¼a de cuadratur¼a de tip Euler-Mac Laurin, numitasi formula perturbat¼a a trapezului (în care mai apar pe noduri doar derivate de ordinulîntâi) :
bZa
f(x)dx =b� a
2[f(a) + f(b)]� (b� a)2
12[f 0(b)� f 0(a)] +R(f): (3.4)
Estimarea restului acestei formule de cuadratur¼a a fost obtinut¼a de K.Petr în 1915 (pentrualte moduri de obtinere a acestei formule si a restului, a se vedea si [17], [24]),
jR(f)j �(b� a)5 �
f IV 720
(3.5)
pentru functii f 2 C4[a; b]: Recent, în 2001, N.Barnett si S.S.Dragomir au obtinut în [2]o estimare ce utilizeaz¼a derivata de ordinul trei,
jR(f)j � (b� a)4 � kf 000k160
(3.6)
pentru f 2 C3[a; b]: Ulterior, în 2002, au fost condusi la o estimare valabil¼a pentru functiide 3 ori derivabile cu derivata de ordinul 3 Lipschitz (având constanta Lipschitz L000, a sevedea [3]) :
jR(f)j � (b� a)5 � L000720
: (3.7)
Considerând o diviziune a intervalului [a; b] :
In : a = x0 < x1 < ::: < xn�1 < xn = b;
3.3. FORMULA PERTURBATA A TRAPEZULUI 43
hi = xi+1 � xi; 8i = 0; n� 1 se obtine în [2] urm¼atoarea formul¼a de cuadratur¼a ( numit¼aformul¼a de cuadratur¼a trapezoidal¼a perturbat¼a ) :
bZa
f(x)dx =1
2
n�1Xi=0
hi[f(xi) + f (xi+1)]�1
12
n�1Xi=0
h2i [f0(xi+1)� f 0 (xi)] +Rn(f; In) (3.8)
pentru care restul satisface estimarea
jRn(f; In)j �kf 000k160
n�1Xi=0
h4i : (3.9)
Dac¼a diviziunea In este echidistant¼a atunci hi = h = b�an; 8i = 0; n� 1; iar formula se
simpli�ca (3.8) devenind,
bZa
f(x)dx =b� a
2n[f(a) + 2
n�1Xi=1
f(xi) + f(b)]� (b� a)2
12n2[f 0(b)� f 0(a)] +Rn(f): (3.10)
Estimarea restului este în acest caz ( conform cu [2]),
jRn(f)j �(b� a)4 � kf 000k
160n3: (3.11)
iar inegalitatile (3.7) si (3.5) ne conduc la estimarile
jRn(f)j �
8><>:(b�a)5�L000720n4
(b�a)5�kfIV k720n4
:
(3.12)
Observatie: Formula (3.4)-(3.5) se poate obtine si prin integrarea pe [a; b] a formuleide interpolare Hermite cu doua noduri duble (1.10), iar formula (3.8) se poate obtineprin integrarea pe [a; b] a formulei de interpolare neteda pe portiuni generata pe �ecaresubinterval de formula (1.11). Formula (3.10) conduce la urmatorul algoritm :
Algoritmul cuadraturii perturbate a trapezuluiI. Date de intrare:a,b capetele intervalului de integraren, numarul nodurilorf functia ce se integreaza
II. Date de iesire : S, valoarea aproximativa a integraleibRa
f(x)dx
III. Pasii algoritmuluiPasul 1. Pentru i = 0; n calculeaza
xi := a+ i � b� a
n
44 CAPITOLUL 3. FORMULE DE CUADRATURA
Pasul 2. CalculeazaT := 0
Pentru i = 1; n calculeaza
T := T + [f (xi�1) + f (xi)]
Pasul 3. Calculeaza
S :=(b� a)
2n� T � (b� a)2
12n2� [f 0(b)� f 0(a)]:
Pasul 4. Tipareste S. STOP.
Probleme propuse :1. Sa se scrie un program in limbajul C pentru algoritmul cuadraturii perturbate a
trapezului.2. Sa se utilizeze formula de cuadratura perturbata a trapezului pentru n=10 la
calculul integralei Z 1
0
1
x+ 1dx
si sa se compare rezultatul cu valoarea exacta a acestei integrale ln 2:3. Sa se utilizeze formula de cuadratura a perturbata a trapezului pentru n=10 la
calculul integralei Z 1:5
0
e�t2
2 dt;
si sa se compare cu rezultatul de la paragraful precedent si cu cel dat in tabelele de valoriale functiei integrale a lui Laplace.4. Sa se utilizeze formula de cuadratura perturbata a trapezului pentru n=20 la
calculul integralei�2Z0
p1� �2 sin2 tdt
cu � = 0; 016729 si apoi sa se calculeze o aproximatie pentru lungimea orbitei Pamantuluiin jurul Soarelui
L (P ) = 4a
�2Z0
p1� �2 sin2 tdt;
unde a = 149; 6� 106 km.
3.4 Formula lui Simpson
Formula de cuadratura a lui Simpson este
bZa
f (x) dx =b� a
6
�f (a) + 4f
�a+ b
2
�+ f(b)
�+R (f) ;
3.4. FORMULA LUI SIMPSON 45
unde restul formulei, cand f 2 C4[a; b]; este
R (f) = �(b� a)5
2880f (iv) (�) ; a < � < b;
astfel pentru evaluarea erorii (a se vedea [4], [17]) avem
jR (f)j � (b� a)5
2880M4(f);
undeM4(f) = maxf
��f IV (x)�� : x 2 [a; b]g:De�nitia 1. Numim formula de cuadratura a lui Simpson repetata formula de cuadratura
obtinuta prin aplicarea formulei de cuadratura a lui Simpson pe �ecare subinterval[a+ 2 (k � 1)h; 2kh] ; k = 1; n; unde h = b�a
2n:
1. Formula de cuadratura repetata a lui Simpson este
bZa
f (x) dx =b� a
6n
"f (a) + 4
nXk=1
f (a+ (2k � 1)h) + 2n�1Xk=1
f (a+ 2kh) + f(b)
#+Rn (f) :
2. Daca f 2 C4[a; b] atunci restul in formula de cuadratura repetata a lui Simpsonare expresia
Rn (f) = �(b� a)5
2880n4f (iv) (�) = �(b� a)h4
180f (iv) (�) ; a < � < b;
de unde
jRn (f)j �(b� a)5
2880n4M4(f) =
(b� a)h4
180M4(f):
3. Formula de cuadratura repetata a lui Simpson are gradul de exactitate 3.4. Daca se considera in formula de cuadratura repetata a lui Simpson
n = 2p; p = 0; 1; 2; :::
si se noteaza prin
hp =b� a
2n=b� a
2p+1; p = 0; 1; 2; :::
cuS0 = f (a) + f(b);
S(1)p =
nXk=1
f (a+ (2k � 1)hp) ; p = 0; 1; 2; :::
S(2)p =n�1Xk=1
f (a+ 2khp) ; p = 0; 1; 2; :::
46 CAPITOLUL 3. FORMULE DE CUADRATURA
cu S(2)0 = 0; atunci au loc relatiile de recurenta
hp+1 =1
2hp;
S(2)p+1 = S(2)p + S(1)p ;
decibZa
f (x) dx � hp3
�S0 + 2S
(2)p + 4S(1)p
�:
Algoritmul lui Romberg pentru formula de cuadratura a lui Simpson
Algoritmul calculeaza valoarea aproximativa a integralei
bZa
f (x) dx prin formula repetata
a lui Simpson, cu o precizie data �:Intrari: a,b-limitele de integraref-expresia functiei��precizia dorita
Iesiri:
bZa
f (x) dx-valoarea aproximativa
Pasul 1: (Initializari)p := 0; hp :=
b�a2; S0 := f (a) + f(b); n0 := 1;
S(2)0 := 0; S
(1)0 := f(a+ h0); I0 :=
h3
�S0 + 2S
(2)p + 4S
(1)p
�Pasul 2: p := p+ 1Pasul 3: hp := 1
2hp�1
Pasul 4: np := 2np�1Pasul 5: S(2)p := S
(2)p�1 + S
(1)p�1
Pasul 6: S(1)p := 0Pasul 7: Pentru k := 1; np executaS(1)p := S
(1)p + f (a+ (2k � 1)hp)
Pasul 8: Ip :=hp3
�S0 + 2S
(2)p + 4S
(1)p
�Pasul 9: Daca jIp � Ip�1j � � atuncisalt la pasul al doileaPasul 10: Tipareste Ip: STOP.
Observatie: Formula de cuadratura repetata a lui Simpson se poate scrie si sub forma
bZa
f (x) dx =b� a
6n[f (a) + 2
nXk=1
f
�a+ i � b� a
n
�+
3.4. FORMULA LUI SIMPSON 47
+4nXk=1
[f
�a+ i � b� a
n
�+ f
�a+ (i� 1) � b� a
n
�] + f(b)] +Rn (f)
cu aceeasi estimare a erorii.Algoritmul cuadraturii clasice a lui SimpsonI. Date de intrare:a,b capetele intervalului de integraren, numarul nodurilorf functia ce se integreaza
II. Date de iesire : T , valoarea aproximativa a integraleibRa
f(x)dx
III. Pasii algoritmului1. Se citesc datele: a,b,n2. Se creaza procedura de introducere a functiei f3. Se construiesc nodurile diviziunii:
Pentru i = 0; n calculeaza
x[i] = a+i (b� a)
n
4. Calculul sumei integrale:S := 0; M := 0Pentru i = 1; n calculeaza
S := S + (f (x[i]) + f (x[i� 1]))
M :=M + f
�x[i� 1] + x[i]
2
�5. Calculul aproximatiei integralei:
T :=(b� a)
6n� (S + 4M)
6. A�sarea rezultatului:Tipareste T . STOP
Probleme propuse
1. Sa se scrie un program in limbajul C pentru algoritmul cuadraturii clasice a luiSimpson.2. Se considera functia f (x) = 1
1+x: Sa se calculeze valoarea aproximativa a integralei
1Z0
f (x) dx; cu o precizie � = 10�1; � = 10�2; � = 10�3; folosind formulele repetate ale
trapezului, respectiv clasica a lui Simpson.3. Sa se utilizeze formula de cuadratura clasica a lui Simpson pentru n=20 la calculul
integralei�2Z0
p1� �2 sin2 tdt
48 CAPITOLUL 3. FORMULE DE CUADRATURA
cu � = 0; 016729 si apoi sa se calculeze o aproximatie pentru lungimea orbitei Pamantuluiin jurul Soarelui
L (P ) = 4a
�2Z0
p1� �2 sin2 tdt;
unde a = 149; 6�106 km. Sa se compare rezultatele obtinute folosind cuadratura trapezu-lui, cuadratura perturbata a trapezului si cuadratura lui Simpson.
3.5 Formula lui Newton
Formula de cuadratura a lui Simpson este
bZa
f (x) dx =b� a
8
�f (a) + 3f
�a+
b� a
3
�+ 3f
�a+
2(b� a)
3
�+ f(b)
�+R (f) ;
unde restul formulei, cand f 2 C4[a; b]; este
R (f) = �(b� a)5
6480f (iv) (�) ; a < � < b;
astfel pentru evaluarea erorii (a se vedea [4], [17]) avem
jR (f)j � (b� a)5
6480M4(f);
undeM4(f) = maxf
��f IV (x)�� : x 2 [a; b]g:Pentru formula repetata a lui Newton se considera o diviziune echidistanta a interval-
ului [a; b];� : a = x0 < x1 < ::: < xn�1 < xn = b
cu
xi = a+i (b� a)
n; i = 0; n
si punctele care impart in trei segmente de lungimi egale �ecare subinterval [xi�1; xi]; i =1; n :
yi = xi�1 +xi � xi�1
3; zi = xi�1 +
2(xi � xi�1)
3; i = 1; n
si se obtine formula
bZa
f (x) dx =b� a
8
"nXk=1
(f (xi) + f (xi�1)) + 3
nXk=1
f (yi) + 3
nXk=1
f (zi)
#+Rn (f)
cu estimarea erorii
jRn (f)j �(b� a)5
6480n4M4(f):
3.5. FORMULA LUI NEWTON 49
Algoritmul formulei de cuadratura a lui Newton:I. Date de intrare:a,b capetele intervalului de integraren, numarul nodurilorf functia ce se integreaza
II. Date de iesire : T , valoarea aproximativa a integraleibRa
f(x)dx
III. Pasii algoritmului1. Se citesc datele: a,b,n2. Se creaza procedura de introducere a functiei f3. Se construiesc nodurile diviziunii:
Pentru i = 0; n calculeaza
x[i] = a+i (b� a)
n
4. Se construiesc nodurile din interiorul subintervalelorPentru i = 1; n calculeaza
y[i] := x[i� 1] + x[i]� x[i� 1]3
z[i] := x[i� 1] + 2(x[i]� x[i� 1])3
5. Calculul sumelor integraleFie S := 0; U := 0; V := 0Pentru i = 1; n calculeaza
S := S + (f (x[i]) + f (x[i� 1]))
U := U + f (y[i])
V := V + f (z[i])
6. Calculul aproximatiei integralei:
T =(b� a)
8n� (S + 3U + 3V )
7. A�sarea rezultatului:Tipareste T . STOP.
1. Sa se scrie un program in limbajul C pentru algoritmul cuadraturii lui Newton.2. Se considera functia f (x) = 1
1+x: Sa se calculeze valoarea aproximativa a integralei
1Z0
f (x) dx; cu o precizie � = 10�1; � = 10�2; � = 10�3; folosind formulele repetate ale
lui Simpson, respectiv Newton..
50 CAPITOLUL 3. FORMULE DE CUADRATURA
Capitolul 4
Rezolvarea sistemelor algebrice deecuatii liniare
4.1 Metoda directa a lui Gauss
Metodele de rezolvare a sistemelor de ecuatii liniare se impart in doua clase: metodedirecte care furnizeaza solutia exacta a sistemului si metode iterative care furnizeaza unsir recurent de aproximatii a acestei solutii.Fie sistemul Ax = b; A = (aij)i;j=1;n ; b = (b1; :::; bn)
T ; cu det (A) 6= 0:Metoda lui Gauss consta din doua etape (a se vedea [4] si [12]):1) transformarea sistemului dat intr-un sistem triunghiular echivalent:
nXi=p
a(p)pj xj = b(p); p = 1; n;
elementele matricei acestui sistem �ind date de relatiile:
a(1)ij := aij; b
(1)i := bi; i; j = 1; n;
a(p+1)ij := a
(p)ij � a
(p)pj
�a(p)ip =a
(p)pp
�;
b(p+1)i := b
(p)i � b(p)p
�a(p)ip =a
(p)pp
�;
pentru p = 1; n� 1; iar i; j = p+ 1; n; in ipoteza ca a(p)pp 6= 0; p = 1; n:
Observatia 1. (a se vedea [7]) Asigurarea conditiilor a(p)pp 6= 0; p = 1; n; se poate
face prin procedeul pivotului maxim. Astfel, pentru a asigura ca a(p)pp 6= 0; se determina,
dintre elementele a(p)ij ; i; j = p; n; elementul cu valoarea absoluta maxima. Fie acesta
a(p)ipjp
: Schimband intre ele liniile p si ip respectiv coloanele p si jp, elementul a(p)ipjp
vaocupa pozitia (p,p). Matricea A �ind nesingulara, acest element este diferit de zero.Observatia 2. Pentru a evita unele schimbari de coloane, se poate folosi o pivotare
partiala. Astfel, daca a(p)pp = 0, se schimba linia p cu linia ip pentru care avem��aip;p�� = maxn���a(p)ip ��� ; i = p+ 1; n
o:
51
52 CAPITOLUL 4. REZOLVAREA SISTEMELOR LINIARE
In cazul in care a(p)ip = 0; i = p; n; se va face uz si de schimbarea coloanelor.
2) rezolvarea sistemului triunghiular folosind metoda substitutiei inverse:
xn = b(n)n =a(n)nn ;
xp =1
a(p)pp
"b(p)p �
nXj=p+1
a(p)pj xj
#; p = n� 1; 1:
Algoritmul (Gauss)
Algoritmul rezolva sistemul liniar Ax = b prin metoda lui Gauss, alegandu-se dreptvarianta de implementare aceea a pivotului maxim.Intrari: n-ordinul matriceiA-matricea sistemuluib- vectorul termenilor liberiIesiri: x-vectorul solutiilor, sauun mesaj de eroare, in cazul in care sistemul nu e compatibil determinatPasul 1: (crearea unui vector auxiliar p) Pentru k = 1; n; executa: pk := kPasul 2: Pentru k = 1; n� 1; executaPasul 2.1: determina i; j 2 fk; k + 1; :::; ng astfel incat:jaijj = max
�jar;qj ; r; q = k; n
Pasul 2.2: Daca aij = 0 atunciTipareste:�Sistemul nu este compatibil determinat�. STOP.Pasul 2.3: Schimba linia i cu linia k din matricea APasul 2.4: Schimba componentele i si k in vectorul bPasul 2.5: Schimba coloana j cu coloana k in matricea APasul 2.6: Schimba componentele i si k ale vectorului p.Pasul 2.7: Pentru i = k + 1; n executaPasul 2.7.1:
temp :=aikakk
Pasul 2.7.2: Pentru j = k + 1; n executa
aij := aij � temp � akj
Pasul 2.7.3: bi := bi � bk � tempPasul 3: Daca ann = 0 atunciTipareste:�Sistemul nu este compatibil determinat�. STOP.Pasul 4: Executa
xn =bnann
Pentru p = n� 1; 1 executa
xp =1
app
bp �
nXj=p+1
apj � xj
!
4.2. METODA ITERATIVA A LUI IACOBI 53
Pasul 5: Ordoneaza in mod crescator componentele vectorului p si corespunzator ele-mentele vectorului x .Pasul 6: Tipareste x. STOP.
Probleme propuse1. Sa se scrie un program in limbajul C pentru algoritmul lui Gauss..2. Sa se rezolve sistemele:
(a)
8>><>>:2x1 + x2 � 4x3 + x4 = 0
�4x1 + 3x2 + 5x3 � 2x4 = 10x1 � x2 + x3 � x4 = 3
x1 + 3x2 � 3x3 + 3x4 = �1
;
(solutia exacta este�72; 72; 52;�1
2
�)
(b)
8<:2x1 + 2x2 + 2x3 = 6
x1 � x2 = 23x1 + x2 + 2x3 = 8
:
(solutia exacta este (2; 0; 1))
(c)
8>><>>:x1 + 2x2 + 3x3 � 2x4 = 62x1 � x2 � 2x3 � 3x4 = 83x1 + 2x2 � x3 + 2x4 = 42x1 � 3x2 + 2x3 + x4 = �8
(solutia exacta este (1; 2;�1;�2))folosind algoritmul lui Gauss si sa se compare rezultatul obtinut cu solutia exacta
punand in evidenta singura eroare care poate apare in acest caz, eroarea de trunchiere sirotunjire a masinii de calcul la efectuarea operatiilor aritmetice.
4.2 Metoda iterativa a lui Iacobi
Consider¼am sistemul liniar 8<:a11 � x1 + :::+ a1n � xn = b1
::::::::::::::::::::::::::an1 � x1 + :::+ ann � xn = bn
(4.1)
ce poate �scris sub forma vectorial¼a A�x = b; undeA =
0@ a11 ::: a1n::: ::: :::an1 ::: ann
1A ; x =
0@ x1:::xn
1Asi b =
0@ b1:::bn
1A : Presupunem c¼a detA 6= 0:
Putem scrie sistemul (4.1) sub forma echivalent¼a x = � � x+ �; unde
� = (�ij)i;j=1;n ; � = (�1; :::; �n)
54 CAPITOLUL 4. REZOLVAREA SISTEMELOR LINIARE
se pot exprima în functie de A si b: Aceasta se poate realiza întotdeauna astfel :Dac¼a aii 6= 0; 8i = 1; n; atunci împ¼artind �ecare ecuatie de ordin i cu aii; obtinemsistemul echivalent
fxi =biaii� ai1aii� x1 � :::� ai;i�1
aii� xi�1 �
ai;i+1aii
� xi+1 � :::� ainaii� xn + bi; i = 1; n:
Notând
� =
0BB@0 �a12
a11::: �a1n
a11
�a21a22
0 ::: �a2na22
::: ::: ::: :::� an1ann
� an2ann
::: 0
1CCA si � =
0@ b1a11
:::::bnann
1Aobtinem forma sistemului,
x = � � x+ �: (4.2)
Dac¼a exist¼a un aii = 0; atunci se rearanjeaz¼a ecuatiile si se renumeroteaz¼a necunoscuteleastfel încât sistemul s¼a poat¼a � adus la forma în care aii 6= 0; 8i = 1; n: Astfel, odat¼aajunsi la forma (4.2) a sistemului (4.1), se va putea aplica sistemului liniar metoda aprox-imatiilor succesive, obtinând iteratiile
x(k+1) = � � x(k) + �; k = 0; 1; 2; ::: (4.3)
Dac¼a not¼am x(k+1) =
0@ x(k+1)1
:::
x(k+1)n
1A ; atunci pe componente (4.3) devine
fx(k+1)i =nXj=1
�ij � x(k)j + �i , i = 1; n: (4.4)
Iteratiile (4.3) vor porni de la un vector initial arbitrar, x(0) =
0@ x(0)1
:::
x(0)n
1A : Astfel, este
de�nit un proces iterativ Picard. Aplicând principiul de punct �x al lui Banach vomobtine conditii su�ciente de convergent¼a a acestui proces. :Teorema: (in [4]) Dac¼a detA 6= 0 si k�k1 < 1; atunci sistemul (4.1) are o unic¼a
solutie x� 2 Rn, ce este limita sirului de�nit prin recurent¼a în (4.3) indiferent de vectorulinitial x(0): Dac¼a x(0) = � atunci are loc urm¼atoarea estimare a erorii (apriori, respectiva posteriori): x� � x(k)
1 �
k�kk11� k�k1
� k�k1 ; 8k 2 N�: (4.5)
x� � x(k) 1 �
k�k11� k�k1
� x(k) � x(k�1)
1
Observatia 1: O conditie su�cienta pentru a avea k�k1 < 1 se obtine tinand contde expresia normei pentru matrici,
k�k1 = maxi=1;n
nXj=1
j�i;jj : (4.6)
4.2. METODA ITERATIVA A LUI IACOBI 55
Astfel,
maxi=1;n
nXj=1
j�i;jj = maxi=1;n
nXj=1
�����ai;jai:i���� < 1, nX
j=1;j 6=i
jai;jj < jai;ij ; 8i = 1; n
si atunci o astfel de conditie (usor de veri�cat) este:
nXj=1;j 6=i
jai;jj < jai;ij ; 8i = 1; n:
Totodata, in conformitate cu [20], aceasta conditie este su�cienta pentru ca matricea Asa �e nesingulara, adica detA 6= 0:
Corolar 4.2.1 DacanX
j=1;j 6=i
jai;jj < jai;ij ; 8i = 1; n
atunci sirul de�nit de procesul iterativ (4.3) este convergent si sistemul (4.1) are o unic¼asolutie x� 2 Rn, obtinuta ca limita acestui sir.
In conditiile acestui corolar, sirul�x(k)�k2N este convergent, adica fundamental in R
n
si atunci un criteriu de oprire al algoritmului lui Iacobi este: pentru " > 0 dat, sa sedetermine primul numar natural k astfel incat,���x(k)i � x
(k�1)i
��� < "; 8i = 1; n:
Acest k va reprezenta ultima iteratie si x(k) va �aproximatia ceruta a solutiei, cu o eroarek�k1"1�k�k1
:
Algoritmul lui Iacobi:I. Date de intrare: n numarul de ecuatii si de necunoscute
matricea sistemului A = (aij)i;j=1;nvectorul termenilor liberi b = (b1; :::; bn)" > 0 eroarea admisa
II. Date de iesire: k, ultima iteratie, x(k; i) ; i = 1; n, aproximatia solutiei la ultimaiteratie
sau un mesaj de eroareIII. Pasii algoritmului1) Se citesc n si "2) Se citeste matricea A si vectorul b3) Se veri�ca conditia de convergenta
Pentru i = 1; n calculeaza v(i) := 0Pentru j = 1; n calculeaza v(i) := v (i) + jai;jj
Calculeaza s(i) := v (i)� jai;ijDaca s(i) � jai;ij atunci tipareste mesaj de eroare " Solutia sistemului nu poate
� aproximata folosind metoda lui Iacobi ".
56 CAPITOLUL 4. REZOLVAREA SISTEMELOR LINIARE
4) Pentru i = 1; n calculeaza beta (i) := biai;i
5) Pentru i = 1; n calculeaza x(0; i) := beta (i)u(1; i) := 0
Pentru j = 1; n calculeaza u(1; i) := u (1; i)+��ai;jai;i
��x (0; i)
x(1; i) := u (1; i) + x (0; i) + beta (i)6) Pentru k � 0; cat timp
maxi=1;n
jx (k + 1; i)� x (k; i)j � "
Pentru i = 1; ncalculeaza u(k + 1:i) := 0
Pentru j = 1; n calculeaza u(k + 1; i) := u (k + 1; i)+��ai;jai;i
��x (k; j)
calculeza x(k + 1; i) := u (k + 1; i) + x (k; i) + beta (i)7) Pentru k pentru care se realizeaza
maxi=1;n
jx (k; i)� x (k � 1; i)j < "
Tipareste k si tipareste x(k; i) ; i = 1; n: Stop.Probleme propuse
1. Sa se scrie un program in limbajul C pentru algoritmul lui Jacobi..2. Se da urmatorul sistem liniar, Ax = b :
A =
0BBBBBB@4 �1 0 �1 0 0�1 4 �1 0 �1 00 �1 4 0 0 �1�1 0 0 4 �1 00 �1 0 �1 4 �10 0 �1 0 �1 4
1CCCCCCA ; b =
0BBBBBB@212212
1CCCCCCA :
Acest sistem are solutia exacta x = (1; 1; 1; 1; 1; 1)T :Sa se rezolve sistemul folosind metoda iterativa a lui Jacobi. Ca si aproximatie initiala seva folosi x(0) = b: Precizia cu care se cer solutiile este � = 0; 0001: Comparati rezultatele.
3. Fie sistemul Ax = b; A =
0@ 2 �1 0�1 3 �10 �1 2
1A si b =
0@ 111
1A :
(i) Sa se determine o solutie aproximativa a sistemului folosind: metoda lui Jacobi cu3 pasi;
(ii) Stiind ca valoarea exacta a solutiei este x =
0@ 111
1A ; studiati erorile ce apar.
(iii) Sa se determine o solutie aproximativa a sistemului folosind: metoda lui Jacobiastfel incat diferenta dintre ultimii doi pasi sa �e mai mica decat " = 10�6:
Capitolul 5
Metode numerice pentru ecuatiineliniare
5.1 Separarea radacinilor si metoda injumatatirii in-tervalului
Considerand ecuatia f (x) = 0, cu f : I ! R; pentru determinarea solutiilor reale aleacestei ecuatii intai se izoleaza �ecare dintre aceste radacini si apoi se aplica �ecareiradacini un procedeu de aproximare. De obicei se presupune ca ecuatia are numai radacinisimple (in cazul radacinilor reale multiple de ordin k suntem condusi si la rezolvareaecuatiilor f (i) (x) = 0; i = 1; k � 1): Izolarea radacinilor se obtine urmand metoda deseparare generata de sirul lui Rolle. Concret, daca f 2 C1 (I) si ecuatia f 0 (x) = 0 aresolutiile distincte x1; :::; xn 2 I; iar c = inf I; d = sup I; atunci sirul de elemente nenulelim
x!c;x>cf (x) ; f (x1) ; :::; f (xn) ; lim
x!d;x<df (d) se numeste sirul lui Rolle. Intrucat conform
teoremei lui Rolle, intre doua radacini consecutive ale derivatei exista cel mult o radacina afunctiei, va rezulta ca oriunde in sirul lui Rolle doi termeni consecutivi au remne contrare,vom avea in subintervalul corespunzator exact o radacina a functiei.Dupa izolarea radacinilor se va aplica �ecarei radacini un procedeu iterativ de aproxi-
mare. Un astfel de procedeu este metoda injumatatirii intervalului de care ne vom ocupacin continuare in acest paragraf.Metoda injumatatirii pentru aproximarea solutiei izolate x� 2 (a; b) a ecuatiei f (x) =
0; f 2 C[a; b]; consta in construirea unui sir de intervale
[a0; b0] � [a1; b1] � ::: � [an; bn] � :::
astfel incat x� 2 [an; bn] :Pentru aceasta se considera [a0; b0] = [a; b]: In continuare se ia
[a1; b1] =
�[a; a+b
2]; daca f(a)f
�a+b2
�< 0;
[a+b2; b]; daca f(a)f
�a+b2
�> 0:
In general avem
[an+1; bn+1] =
�[an;
an+bn2]; daca f(an)f
�an+bn2
�< 0;
[an+bn2
; bn]; daca f(an)f�an+bn2
�> 0:
57
58 CAPITOLUL 5. METODE NUMERICE PENTRU ECUATII NELINIARE
Observatia 1. Se poate intampla ca in construirea sirului de intervale sa aveman = x� sau bn = x�; caz in care procedeul se incheie.Observatia 2. Sirurile (an)n�0 si (bn)n�0 reprezinta siruri de aproximatie prin
lipsa, respectiv prin adaos, pentru solutia x�.Observatia 3. Deoarece bn � an =
b�a2n; avem ca
jxn � x�j � b� a
2n;
daca xn = an sau xn = bn:
Algoritmul (Metoda injumatatirii)Algoritmul aproximeaza solutia izolata x� 2 [a; b] a ecuatiei f (x) = 0; f 2 C[a; b];prin metoda injumatatirii intervalului.Intrari: functia fcapetele intervalului: a, bprecizia dorita: epsIesiri: Aproximarea radacinii ecuatiei: SOLPasul 1: Fie c := a+b
2
Pasul 2: Daca b� c � eps atunci SOL:=c. STOP.Pasul 3: Daca f (b) f (c) � 0 atunci a := caltfel b := cPasul 4: Mergi la Pasul 1.
Observatia 4. Inainte de aplicarea acestui algoritm trebuie sa ne asiguram ca inintervalul [a,b] ecuatia are doar o singura solutie (una din conditii �ind f (a) f (b) < 0).Aceasta se poate face folosind una din metodele de izolare a radacinilor unei ecuatii: sirullui Rolle (prezentat mai sus), Sturm, etc.Observatia 5. Chiar daca nu ne-am asigurat ca ecuatia f(x) = 0 are o singura
solutie in [a,b], daca f (a) f (b) < 0; algoritmul de mai sus va converge spre una dinradacinile ce se a�a in [a,b].Observatia 6. Chiar si la aplicarea altor metode iterative de aproximare a unei
radacini (de exemplu metoda aproximatiilor succesive, metoda tangentei, metoda coardei,ce se vor prezenta in paragrafele urmatoare) se izoleaza radacinile prin metoda sirului luiRolle urmata de cativa pasi ai metodei injumatatirii intervalului pentru a micsora �ecareinterval in care se gaseste cate o radacina a ecuatiei. Prin micsorarea intervalului seasigura de asemenea si indeplinirea conditiilor de convergenta speci�ce �ecarei metode.Probleme propuse
1. Sa se scrie un program in limbajul C pentru algoritmul prezentat.2. Folosind procedura de la 1., sa se aproximeze radacina, situata in intervalul (0; 1) ;
a ecuatieix4 + 2x3 � x� 1 = 0:
3. Folosind procedura de la 1., sa se aproximeze cu 4 zecimale exacte, cea mai mareradacina a ecuatiei:
x6 � x� 1 = 0:
5.2. METODA APROXIMATIILOR SUCCESIVE PENTRU ECUATII SCALARE 59
5.2 Metoda aproximatiilor succesive pentru ecuatiineliniare scalare
Consider¼am ecuatia' (x) = x (5.1)
cu ' : R! R; continu¼a. Pentru x0 2 J ( unde J � R; este interval ) �xat se construiesteprin recurent¼a sirul aproximatiilor succesive
xn+1 = ' (xn) ; n 2 N: (5.2)
Utilizand principiul de punct �x al lui Banach se obtine:Teorema 1: Fie J � R; interval compact si ' : J ! R; astfel încât ' (x) 2 J; 8x 2
J si ' 2 C1 (J) : Dac¼a j'0 (x)j � q < 1; 8x 2 J; atunci ecuatia (5.1) are în J o unic¼asolutie x�, care este limita sirului (5.2). În plus, are loc estimarea erorii :
jxn � x�j � qn
1� q� jx1 � x0j ; 8n 2 N�: (5.3)
Folosind teorema lui Lagrange a cresterilor �nite se poate demonstra:Teorema 2: Fie x� un punct �x al lui ', iar (xn )n2N� sirul aproximatiilor succesive
generat de ' pentru un x0 dat, care converge la x�: Dac¼a ' 2 C1 (V (x�)) si j'0 (x)j <q; 8x 2 V (x�) ; atunci
jxn � x�j � q
1� q� jxn � xn�1j ; 8n 2 N�: (5.4)
Observatia 1: Propozitia precedent¼a precizeaz¼a o estimare a posteriori si d¼a posibil-itatea stabilirii unui criteriu de oprire a procedurii iterative Picard. Astfel, notând cu "eroarea absolut¼a maxim¼a admis¼a, adic¼a jxn � x�j � "; din (5.4) se obtine
jxn � xn�1j �1� q
q� ";
conditie care constituie un criteriu de oprire a calculului.Observatia 2: Conditia de convergent¼a a procedurii iterative Picard este dat¼a de
inegalitatea k'0k1 < 1: R¼asturnarea acestei inegalit¼ati poate conduce la divergenta it-eratiilor Picard. Dac¼a k'0k1 < 1 si '0 (x) > 0; 8x 2 J; atunci sirul (5.2) convergemonoton c¼atre x�; descresc¼ator dac¼a x0 > x� si cresc¼ator dac¼a x0 < x�: Dac¼a k'0k1 < 1si '0 (x) < 0; 8x 2 J; atunci sirul (5.2) converge oscilant c¼atre x�:Observatia 3: Ecuatia f (x) = 0 poate �pus¼a sub forma (5.1) si se poate g¼asi functia
' astfel încât s¼a �e îndeplinit¼a conditia k'0k1 < 1: Astfel, dac¼a f 2 C1 (J) si presupunândc¼a exist¼a m1;M1 > 0 astfel încât
0 < m1 � f 0 (x) �M1; 8x 2 J
( acest lucru se poate realiza dac¼a f 0 (x) 6= 0 si atunci când f 0 (x) < 0; c¼aci ecuatiaf (x) = 0 poate � pus¼a sub forma �f (x) = 0 ). Apoi, ecuatia f (x) = 0 se pune subforma echivalent¼a x = x� � � f (x) ; unde � > 0 se alege astfel încât s¼a avem k'0k1 < 1;
60 CAPITOLUL 5. METODE NUMERICE PENTRU ECUATII NELINIARE
si se va pune ' (x) = x � � � f (x) : Pentru a determina constanta � > 0 pornim de lainegalitatea
0 � '0 (x) = 1� � � f 0 (x) � q < 1
ceea ce va conduce la a alege � = 1M1si atunci vom obtine q = 1� m1
M1< 1:
Observatia 4: Dac¼a intervalul compact J nu contine originea si functia ' este bi-jectiv¼a, atunci ecuatia (5.1) poate � rezolvat¼a si când j'0 (x)j � k > 1: În acest cazprocesul iterativ ar deveni divergent. De aceea se pune ecuatia (5.1) sub forma echiva-lent¼a x = '�1 (x)
not= (x) : Vom �condusi la solutia ecuatiei (5.1), iar procedura iterativa
Picard corespunz¼atoare este convergent¼a, deoarece
j 0 (x)j =���� 1
'0 ( (x))
���� � 1
k= q < 1:
Algoritmul metodei aproximatiilor succesiveI. Date de intrare:x0 iteratia initiala' functia contractie" > 0 eroarea admisaII. Date de iesire:xn ultima iteratie ce aproximeaza solutia cu eroarea admisaIII. Pasii algoritmuluiPasul 1. Calculeaza
x1 = ' (x0)
Pasul 2. Pentru n � 0; cat timp
jxn � xn+1j � "
Calculeaza
xn+1 = ' (xn)
Pasul 3. Pentru n 2 N; pentru care
jxn � xn�1j < "
Tipareste xn: STOP.
Probleme propuse :1. Sa se construiasca un program in limbajul C pentru metoda aproximatiilor succe-
sive.2. Folosind metoda aproximatiilor succesive sa se aproximeze solutia din intervalul
(1,2) a ecuatieix4 � x� 1 = 0:
3. Folosind metoda aproximatiilor succesive sa se aproximeze radacina reala a ecuatiei
x3 � x� 1 = 0:
5.3. METODA APROXIMATIILOR SUCCESIVE PENTRU SISTEME NELINIARE61
4. Folosind metoda aproximatiilor succesive sa se aproximeze radacina reala a ecuatiei
x5 � 5x+ 1 = 0
situata in intervalul (�1; 1) :5. Folosind metoda aproximatiilor succesive sa se aproximeze radacina reala a ecuatiei
lui Keplerx = sinx+ 0:25
situata in intervalul��4; �2
�: Radacina se va izola, folosind metoda sirului lui Rolle si
injumatatirea intervalului, ajungand la intervalul (1:1; 1:3) si se va lua x0 = 1:2: Atunciq = 0:62 = cos 520 si aplicand procedeul recurent (5.2) se va obtine pentru x5 = 1:172eroare mai mica decat 0:0016:
5.3 Metoda aproximatiilor succesive pentru sistemede ecuatii neliniare
Consider¼am sistemul de ecuatii :8<:x1 = '1(x1; :::; xn)
::::::::::::::::xn = 'n(x1; :::; xn)
; (x1; :::; xn) 2 D; (5.5)
unde D � Rn este un domeniu, 'i : D ! R sunt functii continue pe D astfel încât
('1(x1; :::; xn); :::; 'n(x1; :::; xn)) 2 D; 8(x1; :::; xn) 2 D:
Acest sistem poate � pus sub forma vectorial¼a
x = ' (x) ; x 2 D; (5.6)
unde x = (x1; :::; xn) si ' = ('1; :::; 'n) :În conformitate cu metoda aproximatiilor succesive expus¼a în capitolul precedent se
obtine sirul de iteratii Picard :
x(0); x(1); :::; x(m); :::
ce porneste de la o valoare aleas¼a x(0) si se construieste prin recurenta :
x(m+1) = '�x(m)
�; m 2 N: (5.7)
Convergenta acestui sir este asigurat¼a de calitatea de contractie a functiei ': Pentru aobtine calitatea de contractie pentru ' va �su�cient s¼a presupunem c¼a 'i 2 C1 (D) ; 8i =1; n si s¼a g¼asim constanta de contractie cu ajutorul normei matricii Jacobi a lui ';'0 =
�@ 'i@ xj
�; i; j = 1; n: Astfel, în spatiul de matrici putem considera norma Cebîsev,
k'0k1 = �@'i@xj
�; i; j = 1; n
= max1�i�n
nXj=1
maxx2D
����@'i(x)@xj
���� :
62 CAPITOLUL 5. METODE NUMERICE PENTRU ECUATII NELINIARE
Aplicând principiul de punct �x al lui Banach si varianta pe Rn a teoremei cresterilor�nite a lui Lagrange se obtine :Teorema 1: Dac¼a 'i 2 C1 (D) ; 8i = 1; n si k'0k1 � q < 1; atunci ecuatia (5.6)
are jjn D o solutie unic¼a x� astfel încât sirul de iteratii Picard (5.7) este convergent sipentru orice x(0) 2 D;
limm!1
x(m) = x�:
În plus, are loc estimarea apriori a erorii de aproximare : x(m) � x� 1 �
qm
1� q� x(0) � x(1)
1 ; 8m 2 N�; (5.8)
unde k�k1 este norma Cebîsev din Rn:Observam c¼a o estimare a posteriori este x(m) � x�
1 �
q
1� q� x(m) � x(m�1)
1
ceea ce va furniza si un criteriu de oprire a procesului iterativ Picard.Observatia 1: Dac¼a aplicatia ' este de�nit¼a pe o submultime compact¼a a lui Rn; ca
de exempluS�x(0); r
�= fx 2 Rn :
x� x(0) � r; r > 0g
atunci pentru convergenta procesului iterativ trebuie îndeplinit¼a si conditia x(0) � '
�x(0)� �
(1� �) r ce va asigura c¼a '�x(m)
�2 S
�x(0); r
�cccond x(m) 2 S
�x(0); r
�:
Observatia 2: Metoda aproximatiilor succesive poate � aplicat¼a si ecuatiilor vec-toriale de forma f(x) = 0; deoarece acesta poate � adus¼a la forma x = x + Mf(x);unde M este o matrice nesingular¼a. Notând ' (x) = x +Mf(x) aceast¼a ecuatie devinex = ' (x) :MatriceaM se va determina astfel încât metoda generat¼a de operatorul I+Mfs¼a convearg¼a. În acest scop, în baza conditiei k'0k1 � q < 1, când f 2 C1 (D) si det(f0 �x(0)�) 6= 0; vom avea '0 (x) = I +M � f 0 (x) ; iar atunci matricea M se va determinaastfel încoot '0
�x(0)�= 0; adic¼a M = �(f 0 �x(0)�)�1: Dac¼a det(f 0 �x(0)�) = 0, atunci se
va alege o alt¼a valoare de pornire x(0):Observatia 3: Fie sistemul de dou¼a ecuatii scris sub forma�
x = F (x; y)y = G(x; y);
unde F;G : D �! R, domeniul D � R2 continand solutia izolata (x�; y�) a sistemului,consta in construirea sirului de puncte (xn; yn)n�0 dupa metoda recurenta,�
xn+1 = F (xn; yn)yn+1 = G(xn; yn);
(5.9)
cu (x0; y0) 2 D:Daca sirul (xn; yn)n�0 este convergent si F;G sunt continue pe D, atunci limita
sirului este solutia (x�; y�) :
5.3. METODA APROXIMATIILOR SUCCESIVE PENTRU SISTEME NELINIARE63
Proprietatea 1. Daca F;G : D �! D sunt continue, cu derivatele partiale deordinul intai continue pe D si
kWk =M < 1; W (x; y) =
� @F@x
@F@y
@G@x
@G@y
�;
atunci (xn; yn) �! (x�; y�) :Observatia 4. Un exemplu de norma ce poate � considerata este
kWk =M = max
�maxD
����@F@x����+maxD
����@F@y���� ;maxD
����@G@x����+maxD
����@G@y����� : (5.10)
Proprietatea 2. Daca sunt satisfacute conditiile dinainte cu norma k�km ; atunciau loc evaluarile
jxn � x�j � Mn
1�Mjx1 � x0j si jyn � y�j � Mn
1�Mjy1 � y0j :
Algoritmul metodei aproximatiilor succesive pentru sistem neliniarI. Date de intrare:x0 si y0 iteratia initialaF si G functiile din sistem" > 0 eroarea admisaII. Date de iesire:xn si yn ultima iteratie ce aproximeaza solutia cu eroarea admisaIII. Pasii algoritmuluiPasul 1. Calculeaza �
x1 := F (x0; y0)y1 := G(x0; y0);
Pasul 2. Pentru n � 0; cat timp
jxn � xn+1j � " sau jyn � yn+1j � "
Calculeaza�xn+1 := F (xn; yn)yn+1 := G(xn; yn);
Pasul 3. Pentru n 2 N�; pentru care
jxn � xn�1j < " si jyn � yn�1j < "
Tipareste xn; yn: STOP.Probleme propuse :1. Sa se construiasca un program in limbajul C pentru metoda aproximatiilor succesive
asociata sistemelor neliniare de doua ecuatii cu doua necunoscute.
64 CAPITOLUL 5. METODE NUMERICE PENTRU ECUATII NELINIARE
2. Sa se aproximeze solutia, cu ambele coordonate pozitive, a sistemului�f (x; y) = 2x2 � xy � 5x+ 1 = 0g (x; y) = x+ 3 lg x� y2 = 0:
Se va considera
F (x; y) =
rx (y + 5)� 1
2
si G (x; y) =px+ 3 lg x si atunci sistemul devine(
x =q
x(y+5)�12
y =px+ 3 lg x
luandu-se ca prima aproximatie punctul (x0; y0) = (3:5; 2:2) : Se va considera vecinatatealui (x0; y0) ca �ind V = f(x; y) : jx� x0j � 0:1; jy � y0j � 0:1g si pe acest domeniu,in conformitate cu (5.10) vom avea M = maxf0:54 + 0:27; 0:42 + 0g = 0:81 < 1: Apli-cand procedeul iterativ (5.9) vom obtine aproximatia cu patru cifre semni�cative exacte(x6; y6) = (3:487; 2:262) :
5.4 Metoda tangentei a lui Newton pentru ecuatiiscalare
Fie ecuatia algebrica sau transcendenta f(x) = 0; unde f : (�; �) �! R. Presupunemca ecuatia are o solutie unica x� 2 (a; b) � (�; �); adica x� este izolata.De�nitia 1. Daca se construieste sirul de aproximante (xn)n�0 ; pentru solutia
x� 2 (a; b) ; dupa regula
xn+1 = xn �f (xn)
f 0 (xn)(5.11)
cu x0 = a sau x0 = b; spunem ca se utilizeaza metoda tangentei (Newton) pentruaproximarea solutiei x�:Observatia 1. Daca sirul (xn)n�0 este convergent si f 2 C1[a; b]; atunci limita
sirului este x�:Observatia 2. Aproximanta initiala x0 2 [a; b] astfel incat f (x0) f
00 (x0) > 0:Proprietatea 1. Daca f 2 C2[a; b] si f 0; f 00 au semn constant pe (a; b) ; atunci
xn �! x�:Proprietatea 2. (a se vedea [22]) In conditiile proprietatii precedente, eroarea
absoluta comisa se poate evalua prin
jxn � x�j � jf (xn)jm1(f)
;
sau
jxn � x�j � M2(f)
2m1(f)(xn � xn�1)
2 ;
5.4. METODA LUI NEWTON PENTRU ECUATII SCALARE 65
undem1(f) = inf
[a;b]jf 0(x)j si M2(f) = sup
[a;b]
jf 00(x)j
de unde se vede ca ordinul de convergenta al metodei este 2. Totodata, din aceastainegalitate rezulta un criteriu de oprire a algoritmului: pentru " > 0 dat sa se determineprimul numar natural n pentru care jxn � xn�1j < ":Observatia 3. Daca M2(f) � 2m1(f); atunci are loc evaluarea
jxn � x�j � (xn � xn�1)2 :
Observatia 4. In ipoteza ca f 00 (x0) 6= 0; o modi�care a metodei tangentei este ceain care derivata se calculeaza numai in primul termen al sirului, obtinandu-se recurenta:
xn+1 = xn �f (xn)
f 0 (x0):
Algoritmul metodei tangentei (Newton)
Algoritmul calculeaza o aproximatie a solutiei izolate x� 2 [a; b] a ecuatiei f (x) = 0;f 2 C[a; b]; prin metoda lui Newton.Intrari: a,b-capetele intervaluluif-expresia functiei fdf-expresia derivatei functiei f (f 0)eps-precizia doritavaloarea derivatei a doua in a (f 00(a))Iesiri: SOL-aproximatia dorita a solutiei ecuatieiPasul 1: Daca f(a)f 00(a) > 0 atunci x0 := aaltfel x0 := bPasul 2: n := 0Pasul 3: n := n+ 1; xn := xn�1 � f(xn�1)=df(xn�1)Pasul 4: Daca jxn � xn�1j < eps atunci SOL:=xn; STOPaltfel mergi la Pasul 3.
Observatia 4. Algoritmul a fost dat pentru cazul in care solutia ecuatiei a fost izolatain intervalul [a,b], iar functiile f 0 si f 00 au semn constant pe acest interval. Alegerealui x0 la pasul 1 este evidenta.
x0 =
�a; f (a) f 00 (a) > 0b; f (b) f 00 (b) > 0
deoarece acopera toate situatiile posibile. Aceasta alegere va conduce la convergentametodei. De asemenea, se observa ca sirul dat in (5.11) este monoton convergent.In cazul general, cand nu avem informatii despre functia f se ia o valoare oarecare, iar
daca dupa un numar �xat de iteratii nu ajungem la precizia dorita se reia procedeul cu oalta valoare de pornire.
66 CAPITOLUL 5. METODE NUMERICE PENTRU ECUATII NELINIARE
Probleme propuse
1. Sa se construiasca algoritmul metodei modi�cate a lui Newton.2. Sa se scrie un program in C pentru algoritmul metodei tangentei si pentru modi�-
carea acesteia prezentata in Observatia 4.3. Sa se aproximeze cu 4 zecimale corecte solutia reala a ecuatiei
f(x) := x3 + x� 0; 1 = 0:
4. Folosind metoda lui Newton sa se aproximeze solutia reala a ecuatiei
x3 � x� 1 = 0:
5. Folosind metoda tangentei sa se calculeze cea mai mica radacina pozitiva a ecuatiei
x� tgx = 0:
6. Folosind metoda lui Newton sa se aproximeze radacina negativa a ecuatiei
x4 � 3x2 + 75x� 10000 = 0:
5.5 Metoda lui Newton pentru sisteme neliniare
Consideram sistemul 8<:f1 (x1; :::; xn) = 0
::::::::::::::::fn (x1; :::; xn) = 0
(5.12)
undef = (f1; :::; fn) este astfel incat fi 2 C1 (D) ; 8i = 1; n, cu D � Rn domeniu.Presupunem ca pe o multime deschisa V � D avem
Jf (a) =D (f1; :::; fn)
D (x1; :::; xn)(a) =
0@ @f1@x1(a) ::: @f1
@xn(a)
::::::: ::: ::::::@fn@x1(a) ::: @fn
@xn(a)
1A 6= 0; 8a 2 V:
Metoda lui Newton de aproximare a solutiilor sistemului (5.12) se bazeaza pe formula luiTaylor si furnizeaza un sir recurent de aproximatii a �ecarei solutii pornind de la un punctinitial x0 = (x01; :::; x
0n) 2 V � D: Mai precis, se construieste sirul de�nit recurent prin,
xk+1 = xk � J�1f�xk�� f�xk�; k 2 N: (5.13)
Asupra convergentei acestui sir are loc:Teorema 1. (a se vedea [12]): Presupunem ca pentru sistemul (5.12) avem fi 2
C2 (D) ; 8i = 1; n. Fie r > 0; x0 2 D si V = fx 2 D : kx� x0k < rg � D: Presupunemca:(i) matricea Jf (x0) este inversabila (adica det (Jf (x0)) 6= 0) si exista A0 > 0 astfel incatin norma matriceala de�nita in (4.6) avem
J�1f (x0) � A0;
5.5. METODA LUI NEWTON PENTRU SISTEME NELINIARE 67
(ii) exista B0 > 0 astfel incat J�1f (x0) � f (x0)
� B0 � r2;
(iii) exista C > 0 astfel incat
nXk=1
���� @2fi@xj@xk
(x)
���� � C; 8i; j = 1; n; 8x 2 V
(iv) � = 2nA0B0C � 1: Atunci procesul iterativ al lui Newton descris de sirul (5.13) esteconvergent,
limk!1
xk = x�
iar x� 2 V este unica solutie in V a sistemului (5.12). In plus, viteza de convergenta ametodei este data de inegalitatea
x� � xk � �1
2
�p�1� �2p�1B0; 8k 2 N�;
iar metoda este stabila in raport cu alegerea punctului initial intr-o bila centrata in x0 curaza mai mica decat 1��
2�B0:
Observatia 1: Metoda modi�cata a lui Newton este generata de sirul recurent:
xk+1 = xk � J�1f�x0�� f�xk�; k 2 N;
in care se evalueaza matricea Iacobi Jf (x0) la �ecare pas doar pe punctul initial.Observatia 2: Considerand sistemul de doua ecuatii neliniare cu doua necunoscute�
F (x; y) = 0G (x; y) = 0
(5.14)
cu F;G 2 C2 (D) ; D � R2; metoda lui Newton se obtine din formula lui Taylor neglijandtermenii ce contin derivatele partiale de ordinul doi:
F (x; y) = F (xn; yn)+@F
@x(xn; yn) (x� xn)+
@F
@y(xn; yn) (y � yn)+
1
2[@2F
@x2(�; �) (x� xn)
2+
+2@2F
@x@y(�; �) (x� xn) (y � yn) +
@2F
@y2(�; �) (y � yn)
2]
G (x; y) = G (xn; yn)+@G
@x(xn; yn) (x� xn)+
@G
@y(xn; yn) (y � yn)+
1
2[@2G
@x2(�0; �0) (x� xn)
2+
+2@2G
@x@y(�0; �0) (x� xn) (y � yn) +
@2G
@y2(�0; �0) (y � yn)
2];
�; �0 �ind situate intre x si xn; iar �; �0 �ind situate intre y si yn: Astfel,
F (x; y) � F (xn; yn) +@F
@x(xn; yn) (x� xn) +
@F
@y(xn; yn) (y � yn)
G (x; y) � G (xn; yn) +@G
@x(xn; yn) (x� xn) +
@G
@y(xn; yn) (y � yn) :
68 CAPITOLUL 5. METODE NUMERICE PENTRU ECUATII NELINIARE
Punand xn+1 = xn+hn, yn+1 = yn+kn si F (xn+1; yn+1) = 0; G (xn+1; yn+1) = 0 obtinem:�F (xn; yn) +
@F@x(xn; yn)hn +
@F@y(xn; yn) kn = 0
G (xn; yn) +@G@x(xn; yn)hn +
@G@y(xn; yn) kn = 0
sistem in necunoscutele (hn; kn). Determinantul acestuia este:
J (xn; yn) =
���� @F@x(xn; yn)
@F@y(xn; yn)
@G@x(xn; yn)
@G@y(xn; yn)
���� = @F
@x(xn; yn)�
@G
@y(xn; yn)�
@F
@y(xn; yn)�
@G
@x(xn; yn) 6= 0
(5.15)si rezulta,
hn = �1
J (xn; yn)����� F (xn; yn) @F
@y(xn; yn)
G (xn; yn)@G@y(xn; yn)
����kn = �
1
J (xn; yn)����� @F@x(xn; yn) F (xn; yn)
@G@x(xn; yn) G (xn; yn)
���� :Se obtine astfel, procedeul iterativ al metodei lui Newton, dat de relatiile recurente,
xn+1 = xn �1
J (xn; yn)� [F (xn; yn) �
@G
@y(xn; yn)�G (xn; yn) �
@F
@y(xn; yn)] (5.16)
yn+1 = yn�1
J (xn; yn)� [@F@x
(xn; yn) �G (xn; yn)�@G
@x(xn; yn) �F (xn; yn)]; n 2 N; (5.17)
pornind de un punct initial (x0; y0) 2 D: In conditiile teoremei 1, adaptate aici pentrun = 2 se obtine proprietatea (5.15) si convergenta procedeului iterativ (5.16)+(5.17)catreunica solutie din D a sistemului (5.14), (x�; y�) 2 D:Observatia 3: Se poate construi si metoda modi�cata a lui Newton pentru sisteme
neliniare de doua ecuatii cu doua necunoscute pentru care cerinta (5.15) se veri�ca doarpe punctul initial. Procedeul iterativ este in acest caz:
xn+1 = xn �1
J (x0; y0)� [F (xn; yn) �
@G
@y(x0; y0)�G (xn; yn) �
@F
@y(x0; y0)]
yn+1 = yn �1
J (x0; y0)� [@F@x
(x0; y0) �G (xn; yn)�@G
@x(x0; y0) � F (xn; yn)]; n 2 N:
Algoritmul metodei lui Newton pentru sisteme neliniareI. Date de intrare: coordonatele punctului initial x0; y0
expresia functiei F, expresia functiei Gexpresiile functiilor @F
@x; @F@y; @G@x, @G@y
" > 0 eroarea admisaII. Date de iesire: n, numar natural, ultima iteratie
xn; yn aproximatia solutieiIII. Pasii algoritmului1. Se citesc x0; y0; "2. Se de�nesc functiile F; G; @F
@x; @F@y; @G@x, @G@y
3. Calculeaza J (x0; y0) := @F@x(x0; y0) � @G@y (x0; y0)�
@F@y(x0; y0) � @G@x (x0; y0)
5.5. METODA LUI NEWTON PENTRU SISTEME NELINIARE 69
Calculeaza
x1 := x0 ��
1
J (x0; y0)
�� [F (x0; y0) �
@G
@y(x0; y0)�G (x0; y0) �
@F
@y(x0; y0)]
y1 := y0 ��
1
J (x0; y0)
�� [@F@x
(x0; y0) �G (x0; y0)�@G
@x(x0; y0) � F (x0; y0)]
4. Pentru n � 0; cat timp
jxn � xn+1j � " sau jyn � yn+1j � "
Calculeaza
J (xn; yn) :=@F
@x(xn; yn) �
@G
@y(xn; yn)�
@F
@y(xn; yn) �
@G
@x(xn; yn)
xn+1 := xn ��
1
J (xn; yn)
�� [F (xn; yn) �
@G
@y(xn; yn)�G (xn; yn) �
@F
@y(xn; yn)]
yn+1 := yn ��
1
J (xn; yn)
�� [@F@x
(x0; y0) �G (xn; yn)�@G
@x(x0; y0) � F (xn; yn)]
5. Pentru n 2 N�; pentru care
jxn � xn�1j < " si jyn � yn�1j < "
Tipareste n; xn; yn: STOP.Probleme propuse:1. Sa se construiasca algoritmul metodei modi�cate a lui Newton.2. Sa se scrie un program in C pentru algoritmul metodei tangentei si pentru modi�-
carea acesteia prezentata in Observatia 3.3. Cu metoda lui Newton sa se gaseasca solutia aproximativa a sistemului�
2x3 � y2 � 1 = 0xy3 � y � 4 = 0
pornind de la punctul initial (x0; y0) = (1:2; 1:7) cu " = 10�6:4. Cu metoda lui Newton sa se gaseasca solutia aproximativa a sistemului�
x2 + y2 = 10px+ y = 2
:
pornind odata de la punctul initial (x0; y0) = (0:9; 3:1) si apoi de la punctul initial(x0; y0) = (2:9; 1:1) cu " = 10�6 si sa se compare rezultatul obtinut cu solutia exacta(1; 3) ; respectiv (3; 1) :
70 CAPITOLUL 5. METODE NUMERICE PENTRU ECUATII NELINIARE
5.6 Metoda coardei
Fie ecuatia algebrica sau transcendenta f (x) = 0; unde f : (�; �) �! R. Presupunemca ecuatia are o solutie unica x� 2 (a; b) � (�; �); adica radacina x� este izolata.De�nitia 1. Daca se construieste sirul de aproximante (xn)n�0 ; pentru solutia
x� 2 (a; b) ; dupa una din formulele
xn+1 = xn �f (xn)
f(b)� f (xn)(b� xn) ; x0 = a; (5.18)
sau
xn+1 = xn �f (xn)
f(a)� f (xn)(a� xn) ; x0 = b; (5.19)
spunem ca se utilizeaza metoda coardei pentru aproximarea solutiei x�:Observatia 1. Daca sirul (xn)n�0 este convergent si f 2 C[a; b], atunci limita
sirului este x�:Observatia 2. Aproximarea initiala x0 se alege dupa cum urmeaza
x0 =
�a; daca f(a)f 00(a) < 0b; daca f(b)f 00(b) < 0
(5.20)
iar aceste doua posibilitati acopera toate cazurile in care putem obtine convergentametodei coardei. Se observa ca ambele siruri date de procesele iterative (5.18) sau (5.19)sunt monoton convergente.Proprietatea 1. Daca f 2 C2[a; b] si f 0; f 00 au semn constant pe [a,b], atunci
xn �! x�:Proprietatea 2. (a se vedea [22]) In conditiile proprietatii 1., eroarea absoluta
comisa se evalueaza prin
jxn � x�j � jf (xn)jm1(f)
;
sau
jxn � x�j � M1(f)�m1(f)
m1(f)jxn � xn�1j ;
undem1(f) = inf
x2[a;b]jf 0 (x)j si M1(f) = sup
x2[a;b]jf 0 (x)j ;
de unde se vede ca ordinul de convergenta al metodei este 1. Totodata, de aici rezultasi un criteriu de oprire a algoritmului: pentru " > 0 dat sa se determine primul numarnatural n pentru care jxn � xn�1j < ":Observatia 3. Daca M1(f) � 2m1(f) atunci are loc evaluarea
jxn � x�j � jxn � xn�1j :
Algoritmul metodei coardei:
I. Date de intrare: a; b capetele intervalului" > 0 eroarea admisa
5.6. METODA COARDEI 71
expresia functiei fvalorile f 00 (a) ; f 00 (b) necesare pentru alegerea primului termen de
iteratieII. Date de iesire: n ultima iteratie
xn aproximatia solutiei la ultima iteratie, sau sirul de aproximatiixi; i = 0; nIII. Pasii algoritmului1. Se introduc datele a; b; "2. Se de�neste expresia functiei f3. Se introduc valorile f 00 (a) ; f 00 (b)4. Daca f(a) � f 00(a) < 0 atunci x (0) := a4.1. Fie
x1 := x0 �f (x0)
f(b)� f (x0)(b� x0)
4.2. Pentru n � 0; cat timp
jxn � xn+1j � "
Calculeaza
xn+1 := xn �f (xn)
f(b)� f (xn)(b� xn)
4.3. Pentru n 2 N�; pentru care
jxn � xn�1j < "
Tipareste n; xn: STOP.altfel treci la pasul 5.
5. Daca f(b) � f 00(b) < 0 atunci x (0) := b5.1. Fie
x1 = x0 �f (x0)
f(a)� f (x0)(a� x0)
5.2. Pentru n � 0; cat timp
jxn � xn+1j � "
Calculeaza
xn+1 = xn �f (xn)
f(a)� f (xn)(a� xn)
5.3. Pentru n 2 N�; pentru care
jxn � xn�1j < "
Tipareste n; xn: STOP.altfel treci la pasul 4.
72 CAPITOLUL 5. METODE NUMERICE PENTRU ECUATII NELINIARE
Deoarece in oricare din cazurile considerate in (5.20) sirul recurent dat de metoda luiNewton si sirul dat de metoda coardei sunt monoton convergente avand tipuri contrare demonotonie, rezulta ca solutia ecuatiei se va situa intotdeauna intre cele doua aproximatiide acelasi rang date de aceste metode. Din acest motiv, precizia se poate imbunataticombinand cele doua metode, iar aproximatia solutiei se va lua ca media aritmetica aaproximatiilor de acelasi rang date de cele doua metode. Astfel de metode combinate seprezinta in continuare.De�nitia 2. Spunem ca avem metoda combinata secanta-tangenta, daca se constru-
ieste sirul de intervale �xs0; x
t0
���xs1; x
t1
�� ::: �
�xsn; x
tn
�� ::: (5.21)
sau �xt0; x
s0
���xt1; x
s1
�� ::: �
�xtn; x
sn
�� ::: (5.22)
astfel incat x� 2 (xsn; xtn) ; respectiv x� 2 (xtn; xsn) ; iar
xtn+1 = xtn �f (xtn)
f 0 (xtn);
xsn+1 = xsn �f (xsn)
f (xtn)� f (xsn)
�xtn � xsn
�:
Observatia 4. Daca f(a)f 0(a) < 0 se ia xs0 = a si xt0 = b; caz in care se obtinesirul (5.21), iar daca f(b)f 0(b) < 0 se ia xs0 = b si xt0 = a; caz in care se obtine sirul(5.22).Observatia 5. Daca xn = xtn sau xn = xsn avem evaluarea
jxn � x�j ���xtn � xsn
�� :Algoritmul pentru : Metoda Combinata Secanta-Tangenta (conform cu [7])
Algoritmul determina un sir de aproximatii ale solutiei izolate x� 2 [a; b] a ecuatieif (x) = 0; f 2 C [a; b] ; f 00 are semn constant pe [a,b] folosind metoda combinatasecanta-tangenta.Intrari: a,b-extremitatile intervaluluif-expresia functiei fdf-expresia derivatei functiei f (f 0)sign(f 00)-semnul derivatei a doua a functiei feps�precizia doritaIesiri: SOL-aproximatie a solutiei ecuatiei, cu precizia cerutaPasul 1: Daca f(a)sign(f 00) < 0 atunci xs0 := a; xt0 := b; n := 0altfel xs0 := b; xt0 := a; n := 0Pasul 2: xtn+1 = xtn � f (xtn) =df (x
tn)
Pasul 3: xsn+1 = xsn �f(xsn)
f(xtn)�f(xsn)(xtn � xsn)
5.6. METODA COARDEI 73
Pasul 4: n := n+ 1Pasul 5: Daca jxtn � xsnj < eps atunci SOL:=xtn: STOPaltfel mergi la Pasul 2.
Observatia 6. Reamintim ca algoritmul a fost dat pentru cazul in care f(a)f(b) < 0(exista o solutie in intervalul [a,b]) si semnul functiei f 00 e constant pe intervalul [a,b].O alta posibilitate de alegere a xs0 si xt0 este cea data de Observatia 4., in care nu enecesara cunoasterea semnului derivatei a doua a functiei f.Observatia 7. Ca aproximatie a solutiei ecuatiei se poate considera xsn sau xtn, dar
exista si o alta posibilitate
SOL:=xsn + xtn2
:
Observatia 8. In aplicatii se cere aproximarea unei solutii pentru o ecuatie data, cueroare precizata. De obicei, eroarea este precizata prin distanta intre 2 aproximatii succe-sive jxn � xn�1j : In general insa, jxn � xn�1j � � nu implica intotdeauna jxn � x�j < �:De aceea, sunt foarte utile metodele combinate prin care solutia x� este aproximata dinambele parti, ea a�andu-se mereu intre 2 aproximante succesive:
xsn � x� � xdn; n = 1; 2; :::
In aceste cazuri, eroarea comisa in aproximarea x� ��xsn + xdn
�=2 nu depaseste marimea�
xdn � xsn�=2: Deci, eroarea de aproximare poate �direct controlata. Un astfel de exemplu
il constituie metoda combinata secanta-tangenta.
Probleme propuse
1. Sa se implementeze in C algoritmul pentru metodacoardei si pentru metoda com-binata secanta-tangenta.2. Folosind metoda secantei sa se aproximeze cu 4 zecimale exacte radacina reala a
ecuatiei:x3 � x� 0; 1 = 0:
3. Sa se aproximeze cu 4 zecimale exacte cea mai mare radacina reala a ecuatiei:
x3 � x� 1 = 0:
4. Folosind metoda coardei sa se aproximeze radacina situata in intervalul (0; 1) aecuatiei
x3 � 6x2 + 10x� 4 = 0cu precizia " = 10�2:
74 CAPITOLUL 5. METODE NUMERICE PENTRU ECUATII NELINIARE
Capitolul 6
Metode numerice pentru ecuatiidiferentiale ordinare
6.1 Metoda retelelor
Aceasta metoda furnizeaza aproximatii ale solutiei problemelor la limita asociate ecuatiilordiferentiale pe un numar �nit de puncte echidistante ale intervalului pe care este de�nitaecuatia.Metoda retelelor aplicata la ecuatii diferentiale liniare utilizeaza diferentele divizate
in aproximarea derivatelor succesive ale unei functii pe noduri si conduce la rezolvari desisteme algebrice de ecuatii liniare (se vedea [4]).Astfel, folosind forma data de Newton polinomului lui Lagrange, se obtin aproximari
ale derivatelor succesive ale unei functii pe un nod x0
f 0 (x0) = [x0; x0 + h; f ] =f (x0 + h)� f (x0)
h(6.1)
sau
f 0 (x0) =f (x0 + h)� f (x0 � h)
2h(6.2)
respectiv
f 000 (x0) = 2[x0 � h; x0; x0 + h; f ] =f (x0 + h)� 2f (x0) + f (x0 � h)
h2(6.3)
si in general, folosind formula de legatura dintre diferentele divizate cu noduri echidistantesi diferentele �nite, se obtine
f (n) (x0) =�nhf (x0)
hn(6.4)
unde diferenta �nita de ordinul intai in punctul a se de�neste prin
�hf (a) = f (a+ h)� f (a) ;
iar diferenta �nita de ordinul k cu pasul h pe nodul a relativa la functia f; este datarecurent prin,
�khf (a) = �h
��k�1h f (a)
�= �k�1
h f (a+ h)��k�1h f (a) :
75
76 CAPITOLUL 6. METODE NUMERICE PENTRU ECUATII DIFERENTIALE
Consideram ecuatia diferentiala liniara de ordinul II
y00 (x) + p (x) � y0 (x) + q (x) � y (x) = f (x) ; x 2 [a; b] (6.5)
cu conditiile la limitay (a) = A; y (b) = B (6.6)
si p; q; f 2 C[a; b]: Presupunem ca problema bilocala (6.5)+(6.6) are solutie unica pe [a; b]:Se considera o diviziune echidistanta a intervalului [a; b] prin nodurile
xi = a+ i � h; i = 0; n
cu h = 1n� (b� a) si se fac notatiile
pi = p (xi) ; qi = q (xi) ; fi = f (xi) ; i = 0; n:
Se scrie ecuatia (6.5) cu conditiile (6.6) pe nodurile xi i = 1; n� 1 si se obtine�y00 (xi) + p (xi) � y0 (xi) + q (xi) � y (xi) = f (xi) ; i = 1; n� 1
y (x0) = A; y (xn) = B:(6.7)
Folosind formulele (6.2) si (6.3) pe nodurile xi i = 1; n� 1 se obtine
y0 (xi) =y (xi+1)� y (xi�1)
2h(6.8)
y00 (xi) =y (xi+1)� 2y (xi) + y (xi�1)
h2: (6.9)
Notand yi = y (xi) ; i = 0; n, din (6.7), (6.8) si (6.9) se obtine sistemul liniar�1h2� (yi+1 � 2yi + yi�1) +
12h� pi � (yi+1 � yi�1) + qi � yi = fi; i = 1; n� 1
y0 = A; yn = B:(6.10)
Rezolvand sistemul (6.10) de n-1 ecuatii cu n-1 necunoscute y1; :::; yn�1 se obtine aproxi-matia solutiei problemei bilocale (6.5)+(6.6) pe nodurile diviziunii considerate.
Algoritmul metodei retelelorI. Date de intrare:n numarul de nodurih pasul diviziuniia, b capetele intervaluluip; q; f functiile coe�cienti ai ecuatieiA; B valorile bilocaleII. Date de iesire:y1; :::; yn�1 aproximatia solutiei pe nodurile interioareIII. Pasii algoritmuluiPasul 1. Pentru i = 0; n calculeaza
xi := a+ i � h
6.2. METODA LUI EULER 77
Pasul 2. Pentru i = 1; n� 1 executa
pi := p (xi) ; qi := q (xi) ; fi := f (xi)
Pasul 3. Executay0 := A; yn := B
Pasul 4. Rezolva sistemul�(2� h � pi) � yi�1 + (2h2 � qi � 4) � yi + (2 + h � pi) � yi+1 = 2h2 � fi; i = 1; n� 1
y0 = A; yn = B:
Pasul 5. Tipareste y1; :::; yn�1. STOP.
Probleme propuse1. Sa se construiasca un program in limbajul C pentru algoritmul metodei retelelor.2. Folosind metoda retelelor sa se rezolve problema bilocala�
y00 (x) + xy0 (x) + x2y (x) = 1; x 2 [0; 1]y (0) = y (1) = 1
folosind 6 puncte ale retelei.
6.2 Metoda lui Euler
Intai prezentam metoda directa a seriilor Taylor astfel incat metoda practica a lui Eulersa �e un caz particular al acestei metode.Se considera problema Cauchy�
y0 = f(x; y)y(x0) = y0
; x 2 [x0; x0 + T ] (6.11)
Presupunand ca solutia y a problemei considerate este unica si satisface conditiile dediferentiabilitate necesare, functia y poate � dezvoltata in serie Taylor in vecinatateapunctului x0 :
y(x) = y(x0) +x� x01!
y0(x0) +(x� x0)
2
2!y00(x0) + ::: (6.12)
Avand in vedere (6.11) si folosind pentru derivatele partiale ale lui f notatiile
fx =@f
@x; fxy =
@2f
@x@y; ...
se obtiney0 = f
y00 = fx + fyy0 = fx + fyf
y000 = fxx + 2fxyf + fyyf2 + fxfy + f 2y f:
78 CAPITOLUL 6. METODE NUMERICE PENTRU ECUATII DIFERENTIALE
Deci,
y(x) = y0 +x� x01!
f(x0; y0) +(x� x0)
2
2!f (1)(x0; y0) + :::
:::+(x� x0)
p�1
p!f (p�1)(x0; y0) + ::: (6.13)
unde f (j) este derivata totala de ordinul j a functiei f(x; y (x)) in raport cu x. Seobserva ca f (j) = y(j+1):Prin transcrierea seriei (6.13) si alegerea unui pas de lungime h, x�x0 = h; se poate
evalua solutia y in punctul x1 := x0 + h: Mai departe, calculand derivatele y0; y00; ::: inpunctul x1; se poate aproxima y(x2); x2 := x1+h: Continuand in acest mod se obtineo multime discreta de valori yk care aproximeaza solutia y in punctele xk := x0 + kh;si yk � y (xk) ; k = 1; 2; :::.In felul acesta, suntem condusi la o metoda iterativa de aproximare discreta a solutiei
y a problemei (6.11).
yk+1 = yk + hf (xk; yk) + :::+hp
p!f (xk; yk) :
Introducand polinomul lui Taylor
(Tkf) (x; y) = f(x; y) +h
2!f (1)(x; y) + :::+
hk�1
k!f (k�1)(x; y); k = 1; 2; :::
formula (6.13) se scrie sub forma:
y(x) = y0 + h (Tkf) (x; y): (6.14)
Luand k = 1 se obtine metoda unipas a lui Euler:
y(x) = y0 + hf(x0; y0)
ceea ce conduce la algoritmul recurent:
yk+1 = yk + hf (xk; yk) ; k = 0; n
in care h = Tn; xi = x0 + ih; i = 1; n:
Datorita interpretarii sale geometrice, metoda lui Euler se mai numeste si metodaliniiilor poligonale. Estimarea erorii comise (a se vedea [1]) este:
jy (xk)� ykj �h2
2� ky00kC ; 8k = 1; n: (6.15)
Algoritmul metodei lui EulerI. Date de intrare: x0; y0; T (de tip real) n (de tip intreg)
expresia functiei f (argumente reale, valori reale)II. Date de iesire: y[i]; i = 1; n aproximatiile solutiei pe noduriIII. Pasii algoritmului
6.3. METODA LUI HEUN 79
1. Se introduc datele si se de�neste functia f2. Calculeaza
h =T
n; x[0] = x0; y[0] = y0
3. Pentru i = 1; n calculeazax[i] = x0 + i � h
4. Pentru i = 1; n calculeaza
y[i] = y[i� 1] + h � f (x[i� 1]; y[i� 1])
5. Pentru i = 1; n tipareste y[i]: Stop.
Probleme propuse1. Sa se scrie un program in limbajul C pentru metoda lui Euler.2. Sa se rezolve problema lui Cauchy�
y0 = x2 � yy(0) = 1
; x 2 [0; 4]
luand pe rand pasul h = 14; 18; 116: Pentru �ecare valoare a lui h sa se tipareasca solutia
adevarata, cea aproximata, eroarea la nodurile luate in considerare. Solutia exacta aproblemei Cauchy este
y(x) = x2 � 2x+ 2� e�x:
Sa se analizeze rezultatele obtinute.3. Se considera problema Cauchy�
y0 = y + x; x 2 [0; 1]y(0) = 1
:
Sa se aproximeze valorile solutiei y in punctele xi; i = 1; 10 folosind metoda lui Euler cuh = 0:1 si sa se compare valorile obtinute cu valorile corespunzatoare ale solutiei exactepe aceste puncte. Solutia exacta este
y (x) = 2ex � x� 1:
6.3 Metoda lui Heun
Observatie: Din inegalitatea (6.15) se vede ca y (xk) = yk + O (h2) ; adica eroarea deaproximare este de ordinul 2 si lim
h!0jy (xk)� ykj = 0; 8k = 1; n. Pentru a creste ordinul
erorii de aproximare s-au propus modi�cari ale algoritmului lui Euler. Una dintre acestemodi�cari este data prin recurenta:
yi+1 = yi + h � f�xi +
h
2; yi +
h
2� f (xi; yi)
�; i = 0; n
80 CAPITOLUL 6. METODE NUMERICE PENTRU ECUATII DIFERENTIALE
care pentru o programare mai facila poate � scrisa sub forma
K1 = f (xi; yi)
K2 = f
�xi +
h
2; yi +
h
2�K1
�yi+1 = yi + hK2; i = 0; n:
O alta modi�care este metoda lui Heun data prin recurenta
K1 = f (xi; yi)
K2 = f (xi + h; yi + hK1)
yi+1 = yi +h
2(K1 +K2) ; i = 0; n:
Aceste doua modi�cari conduc la cresterea cu 1 a ordinulului erorii de aproximare (a sevedea [1]).Algoritmul metodei lui HeunDate de intrare: x0; y0; T (de tip real) n (de tip intreg)
expresia functiei f (argumente reale, valori reale)II. Date de iesire: y[i]; i = 1; n aproximatiile solutiei pe noduriIII. Pasii algoritmului1. Se introduc datele si se de�neste functia f2. Calculeaza
h =T
n; x[0] = x0; y[0] = y0
3. Pentru i = 1; n calculeazax[i] = x0 + i � h
4. Pentru i = 1; n calculeaza
K1 = f (x[i� 1]; y[i� 1])K2 = f (x[i� 1] + h; y[i� 1] + hK1)
y[i] = y[i� 1] + h
2� (K1 +K2)
5. Pentru i = 1; n tipareste y[i]: Stop.
Probleme propuse1. Sa se scrie un program in limbajul C pentru metoda lui Heun.2. Se considera problema Cauchy�
y0 = y + x; x 2 [0; 1]y(0) = 1
:
Sa se aproximeze valorile solutiei y in punctele xi; i = 1; 10 folosind metoda lui Heun cuh = 0:1 si sa se compare valorile obtinute cu valorile corespunzatoare ale solutiei exactepe aceste puncte. Solutia exacta este
y (x) = 2ex � x� 1:Sa se compare rezultatele obtinute cu metoda lui Euler si cele obtinute cu metoda luiHeun.
6.4. METODA RUNGE-KUTTA 81
6.4 Metoda Runge-Kutta
Polinomul lui Taylor
(Tkf) (x; y) = f(x; y) +h
2!f (1)(x; y) + :::+
hk�1
k!f (k�1)(x; y); k = 1; 2; :::
si formulay(x) = y0 + h (Tkf) (x; y)
genereaza metodele de tip Runge-Kutta pentru k = 1; 2; 3; 4:Pentru k = 1 se obtine metoda lui Euler (Runge-Kutta de ordinul I).Pentru k = 2 se obtine metoda Runge-Kutta de ordinul al II-lea data recurent prin
yi+1 = yi +1
2(K0 +K1) ; i = 0; n; n 2 N,
undeK0 = hf(xi; yi); K1 = hf(xi + h; yi +K0); y0 dat.
Pentru k = 3 se obtine metoda Runge-Kutta de ordinul al III-lea (numita si metoda luiRunge) generata recurent prin
yi+1 = yi +1
6(K0 + 4K1 +K2) ; i = 0; n; n 2 N, y0 dat, (6.16)
K0 = hf(xi; yi);
K1 = hf(xi +1
2h; yi +
1
2K0);
K2 = hf(xi + h; yi + 2K1 �K0):
Cea mai des folosita este metoda Runge-Kutta de ordinul al IV-lea (numita si metoda luiKutta �ind obtinuta pentru k = 4) generata recurent prin
yi+1 = yi +1
6(K0 + 2K1 + 2K2 +K3) ; i = 0; n; n 2 N, y0 dat, (6.17)
K0 = hf(xi; yi);
K1 = hf(xi +1
2h; yi +
1
2K0);
K2 = hf
�xi +
1
2h; yi +
1
2K1
�;
K3 = hf(xi + h; yi +K2):
Ordinul de marime al erorii de trunchiere in metodele Runge-Kutta este de hk+1:(a sevedea [8], [1]). Pentru estimarea erorii in metoda Runge-Kutta de ordinul IV are loc:Teorema: (a se vedea [19]) Considerand problema Cauchy�
y0 = f(x; y)y(x0) = y0
82 CAPITOLUL 6. METODE NUMERICE PENTRU ECUATII DIFERENTIALE
pe domeniul D = f(x; y) : jx� x0j � a; jy � y0j � bg cu a:b > 0 date, daca f 2 C4 (D)atunci eroarea de aproximare a solutiei acestei probleme pe noduri y (xi) ; i = 1; n prinvalorile yi; i = 1; n date de formula (6.17) este:
jy (xi)� yij �23163
4320�Kh5; 8i = 1; n
unde K > 0 este o constanta ce depinde de kfkC si @if
@xj@yk
C
; i = 1; 4; j; k = 0; i; j + k = i:
Metoda Runge-Kutta de ordinul IV este descrisa de urmatorul algoritm:Algoritmul metodei Runge-KuttaI. Date de intrare: x0; y0; T (de tip real) n (de tip intreg)
expresia functiei f (argumente reale, valori reale)II. Date de iesire: y[i]; i = 1; n aproximatiile solutiei pe noduriIII. Pasii algoritmului1. Se introduc datele si se de�neste functia f2. Calculeaza
h =T
n; x[0] = x0; y[0] = y0
3. Pentru i = 1; n calculeazax[i] = x0 + i � h
4. Pentru i = 1; n calculeaza
K1 = h � f (x[i� 1]; y[i� 1])
K2 = h � f�x[i� 1] + h
2; y[i� 1] + K1
2
�K3 = h � f
�x[i� 1] + h
2; y[i� 1] + K2
2
�K4 = h � f (x[i� 1] + h; y[i� 1] +K3)
y[i] = y[i� 1] + 16� (K1 + 2K2 + 2K3 +K4)
5. Pentru i = 1; n tipareste y[i]: Stop.Observatia 1. In cazul in care functia f nu depinde de y, procedeul lui Runge (6.16)
se reduce la regula de integrare a lui Simpson.
Probleme propuse1. Sa se scrie un program in limbajul C pentru metoda lui Runge-Kutta.2. Se considera problema Cauchy�
y0 = y + x; x 2 [0; 1]y(0) = 1
:
6.5. METODA APROXIMATIILOR SUCCESIVE 83
Sa se aproximeze valorile solutiei y in punctele xi; i = 1; 10 folosind metoda lui Runge-Kutta cu h = 0:1 si sa se compare valorile obtinute cu valorile corespunzatoare ale solutieiexacte pe aceste puncte. Solutia exacta este
y (x) = 2ex � x� 1:
Sa se compare rezultatele cu cele date de metoda lui Euler.3. Fie problema Cauchy �
y0 = 11+x2
� 2y2y0(0) = 0
; x 2 [0; 4]:
Sa se aproximeze valorile solutiei y in punctele xi = hi; i = 1; n unde h = 14; iar
n = 16; folosind metoda lui Runge-Kutta. Stiind ca solutia exacta a problemei luiCauchy este y = x
1+x2; sa se compare rezultatele obtinute cu valorile solutiei exacte pe
punctele considerate. Sa se observe ce se intampla odata cu departarea de punctul initial.
6.5 Metoda aproximatiilor succesive
Consider¼am problema Cauchy, �y0 = f(x; y)y(a) = y0
(6.18)
unde f : [a; b]� R! R , y0 2 R:Presupunem ca f 2 C ([a; b]� R) si ca exista M > 0; � > 0; > 0 astfel incat
jf (x; y)j �M; 8 (x; y) 2 [a; b]� R
jf (t; u)� f (s; v)j � � ju� vj+ js� tj ; 8 (t; u) ; (s; v) 2 [a; b]� R
si � (b� a) < 1: Atunci, aplicand principiul de punct �x al lui Banach, problema Cauchy(6.18) are o unica solutie y� 2 C[a; b] catre care converge uniform sirul aproximatiilorsuccesive dat in mod recurent prin
y0 (x) = y0; 8x 2 [a; b]
ym(x) = y0 +
xZa
f(s; ym�1(s))ds; 8x 2 [a; b];8m 2 N�: (6.19)
Considerând pentru n 2 N� o diviziune echidistanta a intervalului [a; b];
� : a = x0 < x1 < ::: < xn�1 < xn = b
si notam h = b�an: Atunci,
xi = a+i (b� a)
n; 8i = 0; n
84 CAPITOLUL 6. METODE NUMERICE PENTRU ECUATII DIFERENTIALE
si formula (6.19) se poate scrie pe nodurile acestei diviziuni prin
y0 (xi) = y0; 8i = 0; n
ym(xi) = y0 +
xiZa
f(s; ym�1(s))ds; 8i = 0; n; 8m 2 N�: (6.20)
Pentru calculul aproximatiilor succesive discrete ale solutiei date in (6.20) se folosesteformula de cuadratura a trapezului cu estimarea erorii data in [9] pentru functii Lipschitz,obtinandu-se algoritmul:
y0 (xi) = y0; 8i = 0; ny1(x0) = y0
y1(xi) = y0 +(b� a)
2n�
iXj=1
[f(xj�1; y0) + f(xj; y0)] +R1;i = y1(xi) +R1;i; 8i = 1; n
y2(xi) = y0 +(b� a)
2n�
iXj=1
[f(xj�1; y1(xj�1) +R1;j�1) + f(xj; y1(xj) +R1;j)] +R2;i =
y2(x0) = y0
= y0 +(b� a)
2n�
iXj=1
[f(xj�1; y1(xj�1)) + f(xj; y1(xj))] +R2;i = y2(xi) +R2;i; 8i = 1; n:
Inductiv, pentru m � 3 se obtine,ym(x0) = y0
ym(xi) = y0 +(b� a)
2n�
iXj=1
[f(xj�1; ym�1(xj�1) +Rm�1;j�1)+
+f(xj; ym�1(xj) +Rm�1;j)] +Rm;i = y0 +(b� a)
2n�
iXj=1
[f(xj�1; ym�1(xj�1))+
+f(xj; ym�1(xj))] +Rm;i = ym(xi) +Rm;i; 8i = 1; n:Valorile calculate efectiv sunt ym(xi); i = 1; n;m 2 N� si ele vor aproxima solutia penodurile diviziunii �; y (xi) : Asupra erorii de aproximare are loc:Teorema (a se vedea [5], [6]): Daca functia f este continua, marginita si Lipschitz in
�ecare argument, iar � (b� a) < 1, atunci:���y (xi)� ym(xi)��� � [� (b� a)]m
1� � (b� a)�M (b� a) +
( + �M) (b� a)2
4n[1� � (b� a)]; 8i = 1; n; 8m 2 N�:
(6.21)Observatie: Daca f 2 C2 ([a; b]� R) atunci estimarea erorii este mai buna decat in
(6.21) �ind:���y (xi)� ym(xi)��� � [� (b� a)]m
1� � (b� a)�M (b� a) +
M 00 (b� a)3
12n2[1� � (b� a)]; 8i = 1; n; 8m 2 N�
(6.22)
6.5. METODA APROXIMATIILOR SUCCESIVE 85
undeM 00 =M2(M + 1)2 +M2
1 (M + 1)
M2 = maxf @2f@x2
C
;
@2f@x@y
C
;
@2f@y2 C
g
M1 = maxf @f@x
C
;
@f@y C
g:
Observatie: Din inegalitatile (6.21) si (6.22) rezulta ca
limn!1;m!1
���y (xi)� ym(xi)��� = 0; 8i = 1; n
asigurand convergenta metodei. Pentru oprirea algoritmului se foloseste urmatorul cri-teriu: Fiind dat " > 0 sa se determine primul numar natural m 2 N� astfel incat���ym(xi)� ym�1(xi)
��� < "; 8i = 1; n:
Algoritmul metodei aproximatiilor succesiveI. Date de intrare: a; b; y0; " (de tip real), n (de tip intreg)
expresia functiei fII. Date de iesire: m; ultima iteratie
ym(xi); i = 1; n aproximatiile solutiei pe noduri obtinute la ultimaiteratieIII. Pasii algoritmului:1. Se citesc datele si se de�neste functia f2. Calculeaza x[0] = a; h = b�a
n
3. Pentru i = 1; n calculeaza x[i] = a+ i � h4. Pentru i = 0; n calculeaza y[0; i] = y05. Calculeaza y[1; 0] = y06. Pentru i = 1; n calculeaza
y[1; i] = y0 +(b� a)
2n�
iXj=1
(f (x[j � 1]; y0) + f (x[j]; y0)
7. Pentru m � 1 (intreg), cat timp
maxi=1;n
jy[m; i]� y[m� 1; i]j � "
calculeazay[m; 0] = y0
Pentru i = 1; n calculeaza
y[m; i] = y0 +(b� a)
2n�
iXj=1
(f (x[j � 1]; y[m� 1; j � 1]) + f (x[j]; y[m� 1; j]))
86 CAPITOLUL 6. METODE NUMERICE PENTRU ECUATII DIFERENTIALE
8. Pentru m pentru care
maxi=1;n
jy[m; i]� y[m� 1; i]j < "
Tipareste mesaj "ultima iteratie este m "Pentru i = 0; n tipareste y[m; i]: Stop.
Probleme propuse1. Sa se scrie un program in limbajul C pentru metoda aproximatiilor succesive.2. Se considera problema Cauchy�
y0 = y + x; x 2 [0; 1]y(0) = 1
:
Sa se aproximeze valorile solutiei y in punctele xi; i = 1; 10 folosind metoda aproximatiilorsuccesive cu h = 0:1; n = 10; si sa se compare valorile obtinute cu valorile corespunzatoareale solutiei exacte pe aceste puncte. Solutia exacta este
y (x) = 2ex � x� 1:
Se va lua " = 10�6: Sa se compare rezultatele cu cele date de metoda lui Euler si metodaRunge-Kutta. Sa se observe la ce iteratie s-a oprit algoritmul. Micsorand pasul princresterea numarului de puncte, n = 20; sa se observe eroarea obtinuta si sa se determineiteratia la care s-a oprit algoritmul.
Bibliogra�e
[1] O. Agratini, I. Chiorean, Gh. Coman, R. Trîmbitas, Analiz¼a numeric¼a si teoria aprox-im¼arii, vol.3, Presa Univ.Clujean¼a, Cluj-Napoca, 2002.
[2] N.S.Barnett, S.S.Dragomir, A perturbed trapezoid inequality in terms of the thirdderivative and applications, RGMIA Research Report Collection, vol.4, no.2 (2001),221-233.
[3] N.S.Barnett, S.S.Dragomir, A perturbed trapezoid inequality in terms of the fourthderivative, The Korean J. of Computational & Applied Math., vol.9, no.1, 2002,45-60.
[4] A. Bica, S. G. Gal, Analiza numerica. Metode numerice, partea a doua, Litogra�aUniv. Oradea 1998.
[5] A. M. Bica, Survey on the numerical methods for ODE�s using the sequence of succe-sive approximations, RGMIA Research Report Collection, vol.8, no.1, article 7, 1-16,2005.
[6] A.M. Bica, Metode numerice iterative pentru ecuatii operatoriale, Editura Universi-tatii din Oradea, 2006.
[7] P. Blaga, S. Pop, Gh. Coman, Analiza numerica : lucrari de laborator, UniversitateaBabes-Bolyai Cluj-Napoca, 1994.
[8] C. M. Bucur, C. A. Popeea, Gh. Gh. Simion, Matematici speciale. Calcul numeric,Ed. Didactic¼a si Ped., Bucuresti, 1983.
[9] P.Cerone, S.S. Dragomir, Trapezoidal and Midpoint-Type Rules from InequalitiesPoint of View, Handbook of Analytic Computational Methods in Applied Mathe-matics ( G.A. Anastassiou ed.), Chapman&Hall/CRC Press, Boca Raton, London,New York , Washington DC, 2000.
[10] Gh.Coman, Sever Groze, Rezolvarea numeric¼a a ecuatiilor operatoriale, Litogra�aUniv. �Babes-Bolyai �, 1981.
[11] Gh.Coman, Analiz¼a numeric¼a, Ed. Libris, Cluj-Napoca, 1995.
[12] B.Demidovitch, I.Maron, Elements de calcul numerique, Ed. Mir, Moscow, 1979.
87
88 BIBLIOGRAFIE
[13] S.G. Gal, Analiza numerica. Aproximarea functiilor, partea intai, Editura Universi-tatii din Oradea, 1994.
[14] C.Iacob (coord.), Matematici clasice si moderne, vol.4, Ed. Tehnic¼a, Bucuresti, 1983.
[15] C.Iancu, On the cubic spline of interpolation, Seminar of Functional Analysis andNumerical Methods, Preprint 4, Cluj-Napoca, (1981), 52-71.
[16] C.Iancu, C.Must¼ata, Error estimation in the approximation of functions by interpo-lation cubic splines, Mathematica, Tome 29 (52), No1, (1987),33-39.
[17] D.V.Ionescu, Cuadraturi numerice, Ed. Tehnica, Bucuresti, 1957.
[18] D.V.Ionescu, Aplicarea metodei aproximatiilor succesive in integrarea numerica aecuatiilor diferentiale, Studii si Cerc. Mat., Cluj-Napoca, 11, (1960), 273-286.
[19] D.V. Ionescu, Ecuatii diferentiale si integrale, editia a doua, Ed. Didactica si Peda-gogica, Bucuresti 1972.
[20] Gh. Marinescu, Analiza numerica, Editura Academiei R.S. Romania, Bucuresti 1974.
[21] R. Plato, Concise Numerical Mathematics, Graduate Studies in Mathematics vol.57,AMS Providence, Rhode Island 2003.
[22] Gh. Siretchi, Calcul diferential si integral, vol.1, Ed. Stiinti�ca si Enciclopedica, Bu-curesti 1985.
[23] D.D.Stancu, Gh.Coman, O.Agratini, R.Trîmbitas, Analiz¼a numeric¼a si teoria aprox-im¼arii, vol.1, Presa Univ. Clujean¼a, Cluj-Napoca, 2001.
[24] D.D.Stancu, Gh.Coman, P.Blaga, Analiz¼a numeric¼a si teoria aproxim¼arii, vol.2,Presa Univ.Clujean¼a, Cluj-Napoca, 2002.
[25] O. Stanasila, Analiza matematica, Editura Didactica si Pedagogica Bucuresti 1981.