Laborator Calcul Num

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)


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:


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) ;


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:


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.


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 :


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)


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.


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:


1.7 Interpolare spline cubica

1.7.1 Spline cubic complet general de conditii bilocale


� : 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


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


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.2 Spline cubic natural generat de conditii initiale


� : 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


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


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.


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:


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


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):


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)


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;


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


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:


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


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


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


Pasul 2. CalculeazaT := 0


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


L (P ) = 4a

�2Z0


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


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; :::


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:


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



L (P ) = 4a

�2Z0


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


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:


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..


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)


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 ".


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;


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

�� :


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.


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.


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) :


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) :


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.


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:


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:


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


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



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


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



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


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


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]))


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.

Date post:	09-Aug-2015
Category:	Documents
Upload:	iulia-ghelbere
View:	59 times
Download:	4 times

Laborator Calcul Num

Documents