Capitolul 1.

Rezolvarea sistemelor de ecuatii liniare

Dupa cum se stie, sistemele de ecuatii liniare care admit solutie unica pot fi
rezolvate analitic, prin metode elementare. De altfel, posibilitatea rezolvarii analitice
se intalneste si in multe alte situatii, uneori la probleme foarte complicate (de
exemplu la unele ecuatii diferentiale sau integrale), dar nu intotdeauna se cunoaste
apriori aceasta posibilitate. Programele de calcul simbolic (Mathematica, Matlab,
Mapple, etc.) pot identifica astfel de situatii, astfel ca inainte de abordarea numerica a
problemei este util sa se verifice daca se pot identifica solutii analitice, care ar
conduce la o mult mai mare viteza de calcul si la o precizie oricat de buna.
Motivul principal pentru care sistemele de ecuatii liniare au beneficiat de o
atentie deosebita in cadrul metodelor numerice il constituie faptul ca solutiile
analitice, desi exista, implica de regula calcule algebrice extensive si genereaza
expresii extrem de lungi. Abordarea numerica aducea o economie importanta de timp
chiar si inainte de aparitia calculatoarelor numerice (primele metode de acest tip au
aparut in secolele 18 si 19), dar in momentul de fata se poate considera ca rezolvarea
numerica a unor sisteme destul de mari de ecuatii liniare se face practic instantaneu.
Un alt punct de interes pentru metodele numerice de rezolvare a sistemelor de
ecuatii liniare il constituie faptul ca astfel de sisteme, de regula de dimensiuni mari,
apar in numeroase alte probleme care nu au o rezolvare analitica directa si fac
obligatorie abordarea numerica (de exemplu in probleme de aproximare si interpolare
a functiilor, in anumite metode de rezolvare a ecuatiilor diferentiale si integrale, etc.).
De aceea, cunoasterea cat mai completa a rezolvarilor numerice a sistemelor de
ecuatii liniare este absolut necesara pentru multe dintre capitolele de calcul numeric
care vor fi prezentate in cele ce urmeaza.
Exista si unele dificultati legate de abordarea numerica in acest domeniu. Se
va observa ca exista situatii in care apar erori mult mai mari decat cele datorate
preciziei finite a sistemului de calcul (care are in mod inerent erori de rotunjire
nenule) datorate faptului ca sistemul este rau conditionat. Aceasta inseamna ca una
sau mai multe ecuatii din sistem nu difera foarte mult de o combinatie liniara a
celorlalte ecuatii, ceea ce il apropie de situatia in care numarul de necunoscute este
mai mare decat numarul de ecuatii liniar independente, astfel ca sistemul tinde sa
devina nedeterminat. In legatura cu acest fenomen se va defini un asa numit numar de
conditionare al matricei coeficientilor care poate da indicatii legate de erorile
estimarii solutiilor, indiferent de algoritm si de sistemul de calcul.
Forma generala a unui sistem de n ecuatii cu n necunoscute este:
⎧a11 x1 + a12 x2 + ... + a1n xn = b1
⎪a x + a x + ... + a x = b
⎪ 21 1 22 2 2n n 2
⎨ (1.1)
⎪ ...........................................
⎪⎩an1 x1 + an 2 x2 + ... + ann xn = bn
sau concis:
n

∑a x
j =1
ij j = bi , i = 1, 2,..n (1.2)

Considerand trei matrici: A matricea coeficientilor variabilelor (patrata), X

matricea necunoscutelor (coloana) si B matricea termenilor liberi (coloana), avand
urmatoarele forme:
 ⎛ x1 ⎞ ⎛ b1 ⎞
⎛ a11 a12 … a1n ⎞ ⎜ ⎟ ⎜ ⎟
⎜ ⎟ ⎜ x2 ⎟ ⎜ b2 ⎟
⎜ a a ... a ⎟
21 22 2 n
⎜. ⎟ ⎜. ⎟
A=⎜ ⎟; X = ⎜ ⎟; B = ⎜ ⎟ (1.3)
⎜ ⎟ ⎜ . ⎟ ⎜ . ⎟
⎜ an1 an 2 ann ⎟
⎜. ⎟ ⎜. ⎟
⎜ ⎟ ⎜ ⎟ ⎜⎜ ⎟⎟
⎝ ⎠ ⎜x ⎟
⎝ n⎠ ⎝ bn ⎠
se observa ca in ecuatiile (1.2) membrul stang este echivalent cu inmultirea
0H

matricilor A si X, iar membrul drept contine elementele matricii B, deci sistemul se

poate scrie matriceal sub forma:
A.X = B (1.4)
Cea mai evidenta cale de rezolvare ar fi inmultirea la stanga a ambilor membri
ai acestei ecuatii cu A −1 , de unde rezulta direct solutiile:
X = A −1.B (1.5)
Exista mai multe metode de rezolvare care se bazeaza pe aflarea matricii
inverse a coeficientilor, dar exista si altele, indirecte, bazate pe iteratii, sau care
speculeaza anumite caracteristici speciale ale matricii in unele cazuri particulare. Un
astfel de caz va fi prezentat in continuare, pentru ca unele metode incearca mai intai
sa reduca problema la conditiile acestui caz particular, care este usor de rezolvat.

1.1 Matrici superior si inferior triunghiulare. Substitutie inversa si

directa.

Vom numi matrice superior triunghiulara o matrice care are toate elementele
de sub diagonala principala nule, deci avand urmatoarea proprietate:
⎪⎧aij = 0, pentru i > j
⎨ ; i, j = 1, 2,..n (1.6)
⎪⎩aij ≠ 0, pentru i ≤ j
si urmatoarea forma
⎛ a11 a12 … a1n ⎞
⎜ ⎟
⎜ 0 a22 ... a2 n ⎟
A=⎜ ⎟ (1.7)
⎜ ⎟
⎜ 0 0 ann ⎟
⎜ ⎟
⎝ ⎠
Un sistem caracterizat de o matrice a coeficientilor superior triunghiulara este
rezolvabil usor prin procedeul denumit substitutie inversa. Denumirea este data de
faptul ca se incepe de la ultima ecuatie si se calculeaza pe rand substituind variabilele
deja calculate catre primele ecuatii.
Intr-adevar, din ultima ecuatie, avand in vedere ca toti coeficientii primelor
n − 1 variabile sunt nuli, rezulta:
b
xn = n (1.8)
ann
In penultima ecuatie, numai ultimele doua necunoscute au coeficienti nenuli.
Dar deoarece xn este deja cunoscut, se poate determina xn −1 :
b −a x
xn −1 = n −1 n −1,n n (1.9)
an −1, n −1
 Continuind in acest mod, pentru ecuatia i , coeficientii primelor i − 1
necunoscute sunt nuli iar necunoscutele xi +1 , xi + 2 ,..., xn au fost deja determinate.
Ecuatia se scrie deci
aii xi + ai ,i +1 xi +1 + ... + ain xn = bi (1.10)
de unde se obtine necunoscuta xi :
n
bi − ∑a x
j = i +1
ij j

xi = (1.11)
aii
Procedeul se va repeta pana la obtinerea tuturor necunoscutelor folosind
formula generala (1.11).
1H

In alte cazuri este intilnita o matrice inferior triunghiulara, in care toate

elementele de deasupra diagonalei sunt nule. Rezolvarea este similara, dar in ordine
inversa: se porneste de la prima ecuatie si se fac substitutii pe rand, catre ultima
ecuatie. Acest procedeu se va numi in continuare substitutie directa.

1.2 Metoda de eliminare Gauss. Pivotarea

Daca matricea coeficientilor nu este triunghiulara, cum se intampla in marea

majoritate a cazurilor, se pot face asupra sa anumite transformari care sa o aduca la
forma triunghiulara. Aceste transformari se bazeaza pe urmatoarele operatii care se
pot efectua asupra sistemului (1.1) si care in mod evident pastreaza validitatea
2H

sistemului:
1. Inmultirea unei ecuatii cu o constanta, echivalenta cu inmultirea unei linii din
matricea coeficientilor si a elementului corespunzator din matricea termenilor
liberi cu constanta respectiva.
2. Adunarea sau scaderea a doua ecuatii, echivalenta cu adunarea sau scaderea a
doua linii din matricea coeficientilor si a elementelor corespunzatoare din
matricea termenilor liberi.
3. Schimbarea ordinii de sumare in membrul stang (echivalenta cu permutarea a
doua coloane din matricea coeficientilor ) sau a ordinii ecuatiilor (echivalenta
cu permutarea a doua linii din matricea coeficientilor si a elementelor
corespunzatoare din matricea termenilor liberi ).
Transformarile de acest tip se fac de regula in scopul obtinerii de elemente nule
sub diagonala principala, astfel ca matricea sa devina inferior triunghiulara iar
sistemul sa fie rezolvat apoi prin substitutie inversa.
Sa consideram ca, folosind transformari de tipul celor mentionate, dorim sa facem
0 toate elementele de sub diagonala din prima coloana. Pentru aceasta, va trebui sa
inmultim prima linie cu o constanta corespunzatoare m2 si sa o scadem din a doua,
apoi sa inmultim prima linie cu o alta constanta m3 si sa o scadem din a treia, s.a.m.d.
Pentru anularea elementului de pe prima coloana din lina i , va trebui satisfacuta
conditia:
ai1 − a11mi = 0 (1.12)
de unde rezulta constanta cu care trebuie inmultita prima linie pentru anularea
elementului din prima coloana a liniei i :
a
mi = 1i (1.13)
a11
 Elementul a11 folosit pentru determinarea constantei de multiplicare a primei
linii il vom numi in continuare element pivot, desi denumirea nu este in intregime
corecta, dupa cum se va vedea mai tarziu.
Facand operatia:
aij' = aij − a1 j mi ; i = 2,3..., n; j = 1, 2..., n
pentru liniile 2,3,..., n si coloanele 1, 2,..., n , se va obtine o matrice echivalenta, dar
cu elementele aflate in prima coloana sub diagonala nule. Prima coloana a unei
matrici superior triunghiulare a fost formata.
⎛ a11 a12 … a1n ⎞
⎜ ' ' ⎟
⎜ 0 a22 ... a2 n ⎟
⎜ 0 a ' ... a ' ⎟
⎜ 32 3n
⎟
A=⎜ ⎟ (1.14)
⎜ ' ' ⎟
⎜ 0 an 2 ann ⎟
⎜ ⎟
⎜ ⎟
⎜ ⎟
⎝ ⎠
In continuare, pentru a forma a doua coloana, trebuiesc anulate elementele
a32 , a42 ,..., an' 2 . Aceasta se va face in mod similar, dar folosindu-se ca pivot elementul
' '

'
de pe diagonala celei de a doua linii, a22 . Conditia de anulare a unui element din linia
i a celei de a doua coloane este:
ai' 2 − a22
'
mi = 0 (1.15)

Constanta de multiplicare a celei de a doua linii este deci:

a'
mi = 2'i (1.16)
a22
iar operatiile asupra celorlalte linii:
aij'' = aij' − a2' j mi ; i = 3, 4..., n; j = 2,3,..., n (1.17)
Este evident ca elementele nule din prima coloana raman nule deoarece din ele
se scade constanta inmultita cu a1' j care este deja nul. Din aceasta cauza nu mai este
nevoie sa le calculam, deci indexarea lui j poate sa inceapa de la 2.
S-au obtinut astfel zerouri sub diagonala si in a doua coloana.
⎛ a11 a12 … a1n ⎞
⎜ ' ' ⎟
⎜ 0 a22 ... a2 n ⎟
⎜ 0 0 ... a3' n ⎟
⎜ ⎟
A=⎜ ⎟ (1.18)
⎜ ' ⎟
⎜ 0 0 ann ⎟
⎜ ⎟
⎜ ⎟
⎜ ⎟
⎝ ⎠
Acelasi procedeu se aplica succesiv si pentru restul coloanelor (cu exceptia
ultimei care nu mai trebuie sa contina nici un 0), astfel ca in final matricea sistemului
ramane superior triunghiulara.
 Trebuie avut in vedere, ca pentru ca sistemul sa nu se modifice, aceleasi
operatii trebuiesc efectuate si asupra matricii coloana a termenilor liberi (ceea ce se
schimba in membrul drept al fiecarei ecuatii trebuie schimbat similar si in membrul
stang) .
Este mai avantajos, pentru simplificarea programului ca matricile A si B sa fie
reunite intr-o matrice extinsa, numita si matricea sistemului, cu n linii si n + 1
coloane, iar operatiile descrise sa se faca pentru j = 1, 2,..., n + 1 . Exista insa si
posibilitatea ca operatiile in cele doua matrici sa se efectueze separat.
Sintetizand, in etapa in care se anuleaza elementele de sub diagonala pe
coloana k se vor face urmatoarele operatii:
a'
mi = ik' ; k = 1, 2, n − 1 (1.19)
akk
aij'' = aij' − akj' mi ; k = 1, 2,..., n − 1; i = k + 1, k + 2..., n; j = k , k + 1..., n + 1 (1.20)
Programul 1.1 rezolva un sistem de ecuatii liniare prin metoda Gauss si
substitutie inversa. Intrarile in program sunt matricea sistemului si matricea
termenilor liberi, iar iesirea este o coloana cu solutiile, la care se aduaga si listele de
verificare.

Programul 1.1 Metoda Gauss pentru un sistem de ecuatii liniare

ai = 88−3, 1, 0, 4<, 82, 3, 0, −3<, 85, 0, 2, −5<, 81, 0, 0, 3<<; H∗ Matricea sistemului ∗L
b = 82, −1, 1, 3<; H∗ Matricea termenilor liberi ∗L
n = Length@bD; H∗Dimensiunea matricii ∗L
a = Table@0, 8n<, 8n + 1<D; H∗ Formarea matricii extinse ∗L
x = Table@0, 8n<D;
For@i = 1, i ≤ n, i++,
For@j = 1, j ≤ n, j++,
a@@i, jDD = ai@@i, jDD; a@@i, j + 1DD = b@@iDD DD;
H∗ Ultinul element al fiecarei linii este termenul liber ∗L
Print@ MatrixForm@aDD;
ForAk = 1, k ≤ n − 1, k++,
ForAi = k + 1, i ≤ n, i++,
a@@i, kDD
; H∗ Constanta pentru linia k ∗L
a@@k, kDD
m=

For@j = k, j ≤ n + 1, j++,
a@@i, jDD = a@@i, jDD − m a@@k, jDD; H∗ Noile elemente din linia k ∗L
D;
E;
E;

Print@ MatrixForm@aDD; H∗ Matricea extinsa, superior triunghiulara∗L

H∗ Rezolvarea prin substitutie inversa ∗L
a@@n, n + 1DD
x@@nDD =
a@@n, nDD
ForAj = n − 1, j ≥ 1, j−−,
a@@j, n + 1DD − ⁄nk=j+1 a@@j, kDD x@@kDD
x@@jDD =
a@@j, jDD
;

Print@x@@jDD D;
E H∗ Afisarea solutiilor ∗L

ai.x H∗ Verificarea solutiilor ∗L

b
 Trebuie insa subliniat ca de multe ori metoda Gauss prezentata simplificat mai sus
nu poate fi aplicata, sau da erori foarte mari. Aceasta se intampla cand in cadrul uneia
dintre etapele diagonalizarii, pivotul este nul sau repectiv foarte mic.
In cazul in care unul din pivoti este nul, expresia (1.19) va da rezultat infinit,
3H

semnalat ca o eroare (generand in procesor o intrerupere interna numita exceptie de

tip abort ). De asemenea, daca pivotul de la numitor are o valoare foarte mica,
rezultatul va fi foarte mare si va fi rotujit automat cu o eroare apreciabila sau poate
chiar sa nu incapa in locatia de memorie care i-a fost destinata, fenomen denumit
depasire de registru.

Programul 1.2 Metoda Gauss cu pivotare partiala (pe coloane)

ai = 880, 1, 0, 4<, 82, 3, 0, −3<, 85, 0, 2, −5<, 81, 0, 0, 3<<; H∗ Matricea sistemului ∗L
b = 82, −1, 1, 3<; H∗ Matricea termenilor liberi ∗L
n = Length@bD; H∗Dimensiunea matricii ∗L
a = Table@0, 8n<, 8n + 1<D;
x = Table@0, 8n<D;
H∗ Formarea directa a matricii extinse ∗L
For@i = 1, i ≤ n, i++, a@@iDD = Flatten@Join@ai@@iDD, 8b@@iDD<DD D;
Print@MatrixForm@aDD;
ForAk = 1, k ≤ n − 1, k++,
H∗ Pivotare partiala Hintre liniiL, cautare element maxim din coloana ∗L
For@r = k + 1, r <= n, r++, If@a@@k, kDD < a@@r, kDD, a@@8k, r<DD = a@@8r, k<DD; D;D;
ForAi = k + 1, i ≤ n, i++,
a@@i, kDD
; H∗ Constanta pentru linia k ∗L
a@@k, kDD
m=

For@j = k, j ≤ n + 1, j++,
a@@i, jDD = a@@i, jDD − m a@@k, jDD; H∗ Noile elemente din linia k ∗L
D;
E;
E;
Print@MatrixForm@aDD; H∗ Matricea extinsa, superior triunghiulara∗L
H∗ Rezolvarea prin substitutie inversa ∗L
a@@n, n + 1DD
x@@nDD =
a@@n, nDD
ForAj = n − 1, j ≥ 1, j−−,
a@@j, n + 1DD − ⁄kn=j+1 a@@j, kDD x@@kDD
x@@jDD =
a@@j, jDD
;

Print@x@@jDD D;
E H∗ Afisarea solutiilor ∗L
ai.x H∗ Verificarea solutiilor ∗L
b

Din aceste considerente, de cele mai multe ori programele care folosesc
metoda Gauss se completeaza cu instructiuni care fac alegerea cea mai buna a
pivotului inainte de prelucrarea fiecarei coloane, folosind posibilitatea de permutare a
liniilor si coloanelor matricii (extinse), procedeu denumit in general pivotare. Prin
permutare intre coloane si/sau linii se aduce in pozitia de pivot elementul cel mai
mare in valoare absoluta din restul liniei si/sau coloanei repective. In acest mod,
impartirea (1.19) se va face cu erorile de rotunjire cele mai mici posibile.
4H

Daca pivotarea se face numai intre linii sau numai intre coloane, avem o
pivotare partiala, iar daca se folosesc si pivotari de linii si de coloane avem pivotare
totala. Programul 1.2 este realizat prin completarea programului 1.1 cu unele
instructiuni de pivotare partiala si pot rezolva sisteme care nu ar putea fi rezolvate de
acesta, de exemplu daca primul element al primei linii este 0.
 Se observa ca in linia referitoare la pivotare s-a folosit sintaxa specifica
programului Mathematica pentru inversarea a doua elemente dintr-o lista
{x, y} = {y, x}
Aceasta este echivalenta cu secventa
t = a@@kDD; a@@kDD = a@@rDD; a@@rDD = t
in care a[[ k ]] desemneaza toate elementele liniei k.
Evident, intr-un limbaj obisnuit de programare (de exemplu C++) care nu
poate lucra simultan cu toate elementele unei linii se va folosi secventa mai detaliata:
For@s = 1, s ≤ n + 1, s++, t = a@@k, sDD; a@@k, sDD = a@@r, sDD; a@@r, sDD = tD

1.2 Metoda Gauss-Jordan

Analiza metodei Gauss de aducere a unei matrici oarecare la o forma superior

triunghiulare sugereaza ca, pentru cazurile in care acest lucru este posibil, exista o
matrice M care inmultita cu matricea coeficientilor A sa conduca la o matrice superior
trunghiulara U:
M.A = U (1.21)
Tot asa, pentru o matrice nesingulara A, exista o matrice, numita matricea
inversa, A −1 , care inmultita cu A sa conduca la o matrice unitara I:
A -1 .A = I (1.22)
unde
⎛ 1 0 ... 0 ⎞
⎜ ⎟
⎜ 0 1... 0 ⎟
I= (1.23)
⎜ ............. ⎟
⎜ ⎟
⎝ 0 0 ...1 ⎠
Aceasta echivalenta intre un set de transformari aplicate cu o matrice si
inmultirea unei alte matrici cu matricea data, sugereaza posibilitatea efectuarii de
calcule matriceale pentru identificarea unor transformari utile. Un exemplu in acest
sens il constituie metoda Gauss Jordan care, folosind un program destul de simplu
genereaza o cantitate surprinzator de mare de informatii.
Pentru rezolvarea problemei initiale data de ecuatia (1.4) se formeaza mai
5H

intai o matrice extinsa de dimensiuni n × (2n + 1) , prin reunirea coloanelor matricii

coeficientilor, a matricii termenilor liberi si a coloanelor unei matrici unitare,
Ae = A ∪ B ∪ I (1.24)
Daca efectuam anumite prelucrari ale acestei matrici care sa aduca primele n
coloane la forma unei matrici unitare, acest set de prelucrari ar fi echivalent cu
inmultirea la stanga a matricei extinse cu matricea inversa a coeficientilor. Intr-
adevar, prin aceasta inmultire, tinand cont ca A −1.A = I , A −1.B = X si A −1.I = A −1 ,
obtinem:
A −1.A e = A −1. ( A ∪ B ∪ I ) = I ∪ X ∪ A −1 (1.25)
Se observa ca se obtine in primele n coloane o matrice unitara si in plus, ca
“produse secundare” se obtin: in coloana n + 1 chiar solutia, iar in urmatoarele n
coloane matricea inversa a coeficientilor.
Prin urmare, efectuand transformari care sa aduca primele n coloane la o
matrice unitara vom obtine in final atat solutia cat si matricea inversa, care poate fi de
asemenea utila in problema respectiva.
 Programul 1.3 Metoda Gauss Jordan pentru rezolvarea sistemelor de ecuatii
liniare
A = 884, 3, 1, 3<, 82, 3, 0, 1<, 83, 1, 0, 3<, 81, 0, 1, 0<<;
b = 812, 24, 36, −6<;
n = Length@ bD;

H∗Formarea matricei extinse partiale ∗L

AE1 = Table@0, 8n<, 8n + 1<D;
For@i = 1, i ≤ n, i++, AE1@@iDD = Flatten@Join@ A@@iDD, 8 b@@iDD<DDD;

H∗Formarea matricei extinse totale ∗L

a = Transpose@Join@Transpose@AE1D, IdentityMatrix@nDDD;
Print@ MatrixForm@aDD;

H∗Permutarea liniilor matricei extinse partiale ∗L

p = 0;
For@r = 1, r ≤ n, r++,
For@k = r + 1, k <= n, k++, If@ Abs@a@@r, rDDD < Abs@a@@k, rDDD, p++;
For@c = 1, c ≤ 2 n + 1, c++, t = a@@r, cDD; a@@r, cDD = a@@k, cDD; a@@k, c DD = tD DD;D
Print@ pD;H∗ numarul de permutari de linii ∗L;
Print@ MatrixForm@aDD;

d = H−1Lp; H∗Valoare initiala a determinantului, daca este necesar∗L

ForAi = 1, i ≤ n, i++,
d = d a@@i, iDD;
a@@iDD
a@@iDD =
a@@i, iDD
;

Print@ MatrixForm@aDD;
For@s = 1, s ≤ n, s++,
If@s ≠ i, a@@sDD = a@@sDD − a@@iDD a@@s, iDDD;
D;E;
Print@ MatrixForm@aDD;
Print@"Det@aD=", dD
x = Table@a@@j, n + 1DD, 8j, n<D;
Print@"x=", MatrixForm@xDD;
Ai = Table@a@@i, jDD, 8i, n<, 8j, n + 2, 2 n + 1<D;
PrintA"A−1=", MatrixForm@ AiDE;

[Link] êê MatrixForm H∗ Verificarea matricii inverse obtinute ∗L

A.x H∗ Verificarea solutiilor ∗L

b

Det@ AD H∗ Verificarea determinantului ∗L

Aceste transformari trebuiesc facute in asa fel incat sa se obtina elemente nule
atat dedesubtul cat si deasupra diagonalei. In plus, fiecare element pivot trebuie sa se
imparta la el insusi pentru ca in final sa se obtina numai valori 1 pe diagonala
primelor n coloane.
 De mentionat ca se poate obtine chiar si determinantul matricei coeficientilor
prin inmultirea elementelor matricii diagonalizate (inainte de a fi facute 1) si tinand
seama ca fiecare permutare schimba semnul determinantului.
Practic, mai intai se imparte fiecare linie cu elementul sau de diagonala, astfel
ca el sa devina 1. In continuare, relatiile care vor alcatui algoritmul sunt aceleasi ca si
la metoda Gauss, doar ca se aplica la fiecare linie atat pentru liniile de dedesubtul ei
cat si pentru cele de deasupra, ceea ce face ca limitele de variatie ale indicilor sa se
modifice:
a'
mi = ik' ; k = 1, 2, n; k ≠ i (1.26)
akk
aij'' = aij' − akj' mi ; i = 1, 2,..., n; k = 1, 2..., n; k ≠ i; j = 1, 2..., 2n + 1 (1.27)
Un exemplu de implementare a acestui algoritm este ilustrat de programul 1.3.
S-a facut o mica simplificare fata de programele anterioare folosind posibilitatea
mediului Mathematica de a opera cu o intreaga linie a unei matrici, astfel ca se poate
elimina bucla care opereaza succesiv asupra elementelor unei linii. Astfel, a[[i,j]]
reprezinta elementul din linia i si coloana j, iar a[[i]] reprezinta intreaga linie i.
Programul furnizeaza matricea solutiilor, matricea inversa si valoarea determinantului
matricii coeficientilor.

1.4. Factorizarea Doolittle

Exista numeroase metode de rezolvare a unor sisteme de ecuatii liniare care

utilizeaza decompozitia (factorizarea) unei matrici, adica transformarea acesteia intr-
un produs de doua sau mai multe matrici cu proprietati avantajoase.
Un prim exemplu, pe care il tratam in continuare, il constituie factorizarea
Doolittle, care presupune exprimarea matricii coeficientilor printr-un produs intre o
matrice superior triunghiulara U (“Upper”) si una inferior triunghiulara L (“Lower”).
In plus se impune ca elementele diagonale ale matricii L sa fie unitare:
A = L.U (1.28)
⎛ u11 u12 … 1n ⎞
u ⎛1 0 … 0⎞
⎜ ⎟ ⎜ ⎟
⎜ 0 u22 ... u2 n ⎟ ⎜ l12 1 ... 0 ⎟
U=⎜ ⎟; L = ⎜ ⎟ (1.29)
⎜ ⎟ ⎜ ⎟
⎜ 0 0 unn ⎟ ⎜ ln1 ln 2 1⎟
⎜ ⎟ ⎜ ⎟
⎝ ⎠ ⎝ ⎠
Elementele acestora au deci proprietatile:
⎧lij = 0, pentru i > j
⎧⎪uij = 0, pentru i > j ⎪
⎨ ; ⎨lii = 1 ; i, j = 1, 2,..n (1.30)
⎩⎪uij ≠ 0, pentru i ≤ j ⎪l ≠ 0, pentru i < j
⎩ ij
Explicitand relatia (1.28) pe elemente sub forma:
6H

n
aij = ∑ lik ukj ; i, j = 1, 2,..., n (1.31)
k =1
si tinand cont de proprietatile (1.30), rezulta relatiile:
7H

j j −1
aij = ∑ lik ukj = ∑ lik ukj + lij u jj ; i, j = 1, 2,..., n (1.32)
k =1 k =1
 i i −1
aij = ∑ lik ukj = ∑ lik ukj + uij ; i, j = 1, 2,..., n (1.33)
k =1 k =1

In prima relatie s-a tinut cont ca primul indice al lui ukj , care este chiar indice
de sumare, poate fi doar mai mic sau egal cu al doilea. In a doua relatie s-a tinut cont
ca al doilea indice al lui lik , care este chiar indice de sumare, poate fi doar mai mic
sau egal cu primul si in plus lii = 1 .
Rezulta urmatoarele relatii care permit determinarea succesiva a tuturor
elementelor matricilor L si U:
j −1
aij − ∑ lik ukj
lij = k =1
; i, j = 1, 2,..., n (1.34)
u jj
i −1
uij = aij − ∑ lik ukj ; i, j = 1, 2,..., n (1.35)
k =1
Astfel, pentru i = 1 se obtin:
l11 = 1; u11 = a11
u12 = a12
(1.36)
.............
u1n = a1n
Pentru i = 2 obtinem:
a21
l21 = ; u22 = a22 − l21u12
u11
l22 = 1; u23 = a22 − l21u12 (1.37)
.............
u2 n = a2 n − l21u1n
Pentru i = 3 obtinem:
a
l31 = 31 ; u33 = a33 − l31u13 − l32u23
u11
a32 − l11u12
l22 = ; u34 = a34 − l31u13 − l32u24
u22
l33 = 1; (1.38)
.............
u3n = a3n − l31u1n − l32u2 n
Se continua astfel pana la n si se poate observa ca elementele necesare intr-o
etapa sunt deja cunoscute dintr-o etapa anterioara, astfel ca se pot determina toate
elementele celor doua matrici.
 Programul 1.4 Rezolvarea sistemului prin factorizare Doolittle

dim = 10;
A = Table@ Random@ Real, 10D, 8dim<, 8dim<D; H∗ Numere aleatoare intre 1 si 10∗L
B = Table@Random@Complex, 10 + 10 ID, 8dim<D;
Print@MatrixForm@ADD;
Print@MatrixForm@ BDD;
n = Length@BD;

L = Table@0, 8n<, 8n<D;

U = L; H∗ S−au initializat cu 0 matricile L si U ∗L
ForAk = 1, k ≤ n, k++,
LP k,kT = 1; H∗ Elementele de pe diagonala lui L sunt 1 ∗L
ForAj = k, j ≤ n, j++,

UP k,jT = AP k,jT − ‚ LP k,m T UP m,jT; H∗ Pentru linia k se calculeaza elemntele lui U ∗L

k−1

H∗ pentru k=1, rezulta U@@1,jDD=A@@1,jDD ∗L

m=1

H∗ U@@m,jDD se cunosc, deoarece suma este cu m<k deci s−au calculat anterior∗L
H∗ L@@k,mDD se calculasera anterior din bucla urmatoare ∗L
ForAi = k + 1, i ≤ n, i++,
i y
LPi,kT = j
j
j
z
z ì UP k,kT; H∗ Pentru linia i, se calculeaza elementele lui L ∗L
jAPi,kT − ‚ LPi,m T UP m,kTz
z
k−1

k m=1 {
H∗ pentru k=1, rezulta L@@2,1DD=A@@2,1DDêU@@1,1DD= A@@2,1DDêA@@1,1DD∗L
H∗ U@@m,kDD si L@@i,mDD sunt deja calculate, deoarece suma este pentru m<k ∗L
E;E;E

Print@MatrixForm@ LDD;
Print@MatrixForm@ UDD;

MatrixForm@A − [Link] H∗ Verificari factorizare ∗L

Max@Abs@A − [Link]

Y = Table@0, 8n<D;
H∗ Forward Substitution pentru a afla Y ∗L
BP nT
YP nT =
LP n,nT
;

ForAi = 1, i ≤ n, i++,
i
j y
z
YPiT = j
j z ì LP i,iT;E;
j BPiT − ‚ LPi,jT YPjTz
z
i−1

k j=1 {
H∗ Back substitution pentru a afla X ∗L
X = Table@0, 8n<D;
YP nT
XP nT =
UP n,nT
;

ForAi = n − 1, i ≥ 1, i−−,
i
j y
z
XPiT = j
j z ì UPi,iT;E;
j YPiT − ‚ UPi,jT XPjTz
z
n

k j=i+1 {
Print@MatrixForm@XDD;
A.X H∗ Verificari solutie ∗L
B
Max@Abs@A.X − BDD

Pentru rezolvarea unui sistem de ecuatii liniare de tipul (1.4) se scrie: 8H

A.X = B ⇒ L.(U.X) = B ⇒ L.Y = B (1.39)

unde s-a introdus o matrice auxiliara Y = U.X .
Mai intai se rezolva in raport cu Y ecuatia L.Y = B tinand cont ca L este
inferior triunghiulara. Aceasta permite utilizarea procedeului de substitutie directa
mentionat in paragraful 1.1: din prima linie se determina prima necunoscuta, din a
doua linie a doua necunoscuta, s.a.m.d.
Dupa determinarea lui Y, se poate determina X din ecuatia Y = U.X , unde se
tine seama ca U este superior triunghiulara si deci se poate aplica procedeul de
substitutie inversa.
Programul 1.4 rezolva un sistem de ecuatii liniare prin factorizare Doolittle,
urmata de substitutie directa si apoi substitutie [Link] sunt calculate aici
aleator, avand si unele valori complexe, avand in vedere ca programul Mathematica
poate trata direct astfel de valori. In unele limbaje de programare (Fortran) acest lucru
este de asemenea posibil, in timp ce in majoritatea celor uzuale (de exemplu C++)
trebuiesc incluse si modulele pentru calculul cu numere complexe.

1.5 Factorizarea Crout

Ca si in cazul factorizarii Doolitle, la factorizarea Crout se face exprimarea

matricii coeficientilor sistemului printr-un produs intre o matrice inferior triunghiulara
si una superior triunghilara. Deosebirea consta in faptul ca de aceasta data matricea
superior triunghiulara are elementele diagonale egale cu 1.
A=L.U (1.40)
⎛1 u12 … u1n ⎞ ⎛ l11 0 … 0⎞
⎜ ⎟ ⎜ ⎟
⎜ 0 1 ... u2 n ⎟ ⎜ l12 l22 ... 0 ⎟
U=⎜ ⎟; L = ⎜ ⎟ (1.41)
⎜ ⎟ ⎜ ⎟
⎜0 0 1 ⎟ ⎜ ln1 ln 2 lnn ⎟
⎜ ⎟ ⎜ ⎟
⎝ ⎠ ⎝ ⎠
Elementele celor doua matrici au deci proprietatile
⎧ukj , k < j
ukk = 1; ukj = ⎨ (1.42)
⎩0, k> j
⎧l , k < i
lik = ⎨ ik (1.43)
⎩0, k > i

Conform relatiilor (1.40)- (1.43), separand termenul diagonal u11 = 1 , se poate

9H 10H

scrie
i i −1
aii = ∑ lij u ji = ∑ lij u ji + lii uii (1.44)
j =1 j =1 1
de unde se poate obtine termenul diagonal al matricei U
i −1
lii = aii − ∑ lij u ji (1.45)
j =1

Pentru un element oarecare al matricii A, tinand cont de relatia (1.43) avem

11H

i i −1
aij = ∑ lik ukj =∑ lik ukj + lii uij (1.46)
k =1 k =1

de unde rezulta elementul omolog al matricii U

i −1
aij − ∑ lik ukj
uij = k =1
(1.47)
lii
 Pentru acelasi element al matricii A, tinand cont de relatia (1.42) 12H

j j −1
aij = ∑ lik ukj =∑ lik ukj + lij u jj (1.48)
k =1 k =1

Deoarece u jj = 1 , putem afla si elementul corespunzator al matricii L prin

relatia
j −1
lij = aij − ∑ lik ukj (1.49)
k =1
Programul 1.5 rezolva un sistem de ecuatii liniare prin metoda Crout.
Matricile sunt calculate aici aleator, avand si unele valori complexe, avand in vedere
ca programul Mathematica poate trata direct astfel de valori. Asa cum s-a specificat
anterior, in unele limbaje de programare (Fortran) acest lucru este de asemenea
posibil, in timp ce in majoritatea celor uzuale (de exemplu C++) trebuiesc incluse si
modulele pentru calculul cu numere complexe.

Programul 1.5 Rezolvarea unui sistem de ecuatii liniare prin metoda Crout
dim = 20;
A = Table@Random@Real, 10D, 8dim<, 8dim<D;
B = Table@ Random@Complex, 10 + 10 ID, 8dim<D;
Print@MatrixForm@ADD;
Print@MatrixForm@BDD;
n = dim;

L = Table@0, 8n<, 8n<D;

U = L;
ForAk = 1, k ≤ n, k++,

LP k,kT = AP k,kT − ‚ LP k,m T UP m,kT;

k−1

ForAj = k, j ≤ n, j++,
m=1

i y
UP k,jT = j
j
j
z
z ì LP k,kT;
jAP k,jT − ‚ LP k,m T UP m,jTz
z
k−1

k m=1 {
H∗ Se observa ca pentru j=k,
U@@k,kDD=1 deoarece numitorul L@@ k,kDD stabilit anterior este exact numaratorul∗L
ForAi = k + 1, i ≤ n, i++,
i y
LPi,kT = j
j
j
z
z;
jAPi,kT − ‚ LPi,m T UP m,kTz
z
k−1

k m=1 {
E;E;E

Print@MatrixForm@LDD;
Print@MatrixForm@UDD;

MatrixForm@A − [Link]
Max@Abs@A − [Link]

Y = Table@0, 8n<D;
BP nT
YP nT =
LP n,nT
;

ForAi = 1, i ≤ n, i++,
i
j y
z
YP iT = j z ì LPi,iT;E;
j BP iT − ‚ LPi,jT YPjTz
j z
i−1

k j=1 {
X = Table@0, 8n<D;
YP nT
XP nT =
UP n,nT
;

ForAi = n − 1, i ≥ 1, i−−,
i
j y
z
XP iT = j z ì UPi,iT;E
j YP iT − ‚ UP i,jT XPjTz
j z
n

k j=i+1 {
Print@MatrixForm@XDD;
Max@Abs@A.X − BDD
 1.6 Factorizarea Cholesky

In cazul in care matricea sistemului este simetrica, se poate folosi aceasta

particularitate pentru a se micsora volumul de calcule, prim metoda factorzarii
Choleski. Matricea coeficientilor se exprima de asemenea ca produs intre o matrice
inferior triunghiulara si una superior triunghiulara, dar cu particularitatea ca una este
transpusa celeilalte, deci elementele celor doua matrici sunt simetrice, deci trebuie
determinat doar un set de astfel de elemente:

A=U.L (1.50)
uij = l ji (1.51)
Din relatiile (1.50) si (1.51), pentru un element diagonal, prin separarea
13H 14H

ultimului termen se deduce

n k k −1
akk = ∑ lkmumk = ∑ lkm
2
= ∑ lkm
2
+ lkk2 (1.52)
m =1 m =1 m =1
de unde rezulta elementele de pe diagonale
k −1
lkk = akk − ∑ lkm
2
(1.53)
m =1

De asemenea pentru un element nediagonal se poate scrie

k k k −1
aik = ∑ limumk = ∑ lim ⋅ lkm =∑ limlkm + lik lkk (1.54)
m =1 m =1 m =1
de unde rezulta elementele nediagonale ale matricilor triunghiulare
k −1
aik − ∑ limlkm
lik = m =1
(1.55)
lkk
Programul 1.6 rezolva un sistem de ecuatii liniare cu matricea coeficientilor
simetrica prin metoda Cholesky. Matricile sunt calculate aici aleator, avand ca si in
cazurile anterioare si unele valori complexe.

Programul 1.6 Rezolvarea unui sistem de ecuatii liniare cu matricea

coeficientitor simetrica prin metoda Cholesky

dim = 10;
A = Table@0, 8dim<, 8dim<D; H∗ Se initializeaza matricea ∗L
H∗ Se construieste matricea ca matrice simetrica ∗L
For@i = 1, i ≤ dim, i++,
For@j = 1, j ≤ i, j++,
A@@i, jDD = Random@ Real, 10D;
A@@j, iDD = A@@i, jDD;D;D;
B = Table@ Random@Complex, 10 + 10 ID, 8dim<D;
Print@ MatrixForm@ADD;
Print@ MatrixForm@BDD;
n = dim;
 L = Table@0, 8n<, 8n<D;
U = L;
H∗ Aici este algoritmul Cholesky ∗L
ForAk = 1, k ≤ n, k++,

LP k,kT = $%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
AP k,kT − ⁄ HLP k,m TL2% ;
k−1
m=1

ForAi = k + 1, i ≤ n, i++,
i y
LPi,kT = j
j
j
z
z ì LP k,kT;
jAPi,kT − ‚ LPi,m T LP k,m Tz
z
k−1

k m=1 {
E;
E;
U = Transpose@ LD;
Print@ MatrixForm@LDD;
Print@ MatrixForm@UDD;

A − L.U êê MatrixForm H∗ Verificarea factorizarii ∗L

Y = Table@0, 8n<D;
H∗ Forward substitution ∗L
BP nT
YP nT = ;
LP n,nT
ForAi = 1, i ≤ n, i++,
i
j y
z
YPiT = j
j z ì LP i,iT;E;
jBP iT − ‚ LPi,jT YP jTz
z
i−1

k j=1 {
H∗ Back substitution ∗L
X = Table@0, 8n<D;
YP nT
XP nT = ;
UP n,nT
ForAi = n − 1, i ≥ 1, i−−,
i
j y
z
XPiT = j
j z ì UPi,iT;E
jYP iT − ‚ UPi,jT XPjTz
z
n

k j=i+1 {
Print@ MatrixForm@XDD;
Max@Abs@A.X − BDD
H∗ Verificarea solutiei ∗L

1.7 Metodele iterative Jacobi si Gauss-Siedel. Suprarelaxarea

In anumite probleme este preferata in mod natural o apropiere iterativa de

solutia exacta. In cazul in care iteratiile implica si sisteme de ecuatii liniare, se poate
aplica metoda Jacobi de rezolvare a sistemului matriceal
A.X=B (1.56)
Pentru aceasta se face o descompunere a matricii coeficientilor intr-o suma de trei
matric: una formata cu elementele de pe diagonala, alta cu elementele de sub
diagonala si una cu elementele de deasupra diagonalei
A=D+L+U (1.57)
Elementele acestor matrici sunt deci deja cunoscute, fiind urmatoarele
 dij = aijδ ij
⎧aij , i > j
lij = ⎨ (1.58)
⎩0, i < j
⎧aij , i < j
uij = ⎨
⎩0, i > j
Tinand seama de relatiile (1.56) si (1.57) se poate scrie urmatoarea ecuatie
15H 16H

matriceala
D.X=B-L.X-U.X (1.59)
Prin inmultire la stanga in ambii membrii cu inversa matricii diagonale se
poate exprima matricea necunoscutelor la iteratia k + 1 in functie de matricea
necunoscutelor la iteratia anterioara k
X(k +1) = - D-1 .(L+U).X(k ) +D-1 .B (1.60)
Relatia (1.59) scrisa pe componente este
17H

n i −1 n

∑ dij x(jk +1) = bi − ∑ aij x(jk ) −

j =1 j =1
∑a x
j =i +1
ij
(k )
j (1.61)

Tinand cont de definitia elementelor matricii diagonale

dij = aijδ ij (1.62)
rezulta ecuatia metodei iterative Jacobi (1.60) scrisa pe componente
18H

1 ⎛ i −1 n ⎞
xi( k +1) = ⎜ bi − ∑ aij x (jk ) − ∑ aij x (jk ) ⎟ (1.63)
aii ⎝ j =1 j = i +1 ⎠
Implementarea acestei metode este mai simpla decat la metodele anterioare
deoarece nu este necesara factorizarea matricii ci doar descompunerea ei. Din pacate
metoda este de multe ori slab convergenta, in special cand gradul de independenta
linaiara ala ecuatiilor este redus (numarul de conditionare este mare). Aceasta poate
conduce la erori si timp de calcul mai mari decat in cazul celorlalte metode.
Pentru accelerarea calculelor este preferata in practica metoda Gauss-Siedel
care este foarte apropiata de aceasta.
Tinand seama tot de relatiile (1.56) si (1.57) se poate scrie urmatoarea ecuatie
19H 20H

matriceala

(D+L).X=B-U.X (1.64)
Prin inmultire la stanga in ambii membrii cu inversa matricii (D+L) diagonale
se poate exprima matricea necunoscutelor la iteratia k + 1 in functie de matricea
necunoscutelor la iteratia anterioara k
X( k +1) = -(D+L)-1 .U.X(k ) + (D+L)-1 .B (1.65)
Relatia (1.64) scrisa pe componente este
21H

i n

∑ (D+L)ij x(jk +1) = bi −

j =1
∑ (U)
j = i +1
ij x (jk ) (1.66)

unde s-a tinut cont ca (D+L)ij = 0 pentru j > i si (U)ij =0 pentru j ≤ i . Aceasta se
poate scrie si sub forma
i −1 n
aii xi( k +1) + ∑ aij x (jk +1) = bi − ∑a x ij
(k )
j (1.67)
j =1 j = i +1

unde s-a tinut cont de definitiile (1.58) ale elementelor matricii diagonale si ale
22H

matricilor triunghilare.
 Rezulta ecuatia metodei iterative Gauss-Siedel scrisa pe componente
1 ⎛ i −1 n
⎞
xi( k +1) = ⎜ bi − ∑ aij x (jk +1) − ∑ aij x (jk ) ⎟ (1.68)
aii ⎝ j =1 j = i +1 ⎠

Se observa ca ecuatia metodei Gauss-Siedel difera de cea a metodei Jacobi

numai in ceea ce priveste termenul al doilea din membrul drept, care contine valorile
necunoscutelor la iteratia curenta x (jk +1) si nu pe cele de la iteratia anterioara x (jk ) .
Acest lucru este posibil, deoarece daca in fiecare iteratie valorile necunoscutelor se
calculeaza pe rand cu indice crescator, la calculul necunoscutei xi( k +1) sunt cunoscute
deja valorile x0( k +1) , x1( k +1) ,…, xi(−k1+1) . Evident, deoarece acestea sunt mai apropiate de
valoarea corecta decat omoloagele lor de la iteratia anterioara, procedeul are o
convergenta mai mare si este intotdeauna preferat in practica.
In programul 1.7 este exemplificata aplicarea metodei Gauss-Siedel pentru un
sistem de ecuatii liniare rau conditionat (prima ecuatie si a doua sunt destul de
aproape proportinale). Se constata ca, pentru un numar dat de iteratii, meoda conduce
la eroari mult mai mici decat metoda Jacobi.

Programul 1.7 Rezolvarea unui sistem de ecuatii liniare prin metoda Gauss-
Siedel
Clear@"`∗"D;
A = 886, 3, 2<, 83, 1.51, 1<, 82, −1, 3<<;
B = 81.2, −3.1, 2<;
Print@ MatrixForm@ADD;
Print@ MatrixForm@ BDD;
n = Length@BD;
Xo = Table@100., 8n<D; H∗ Valorile ghicite, initiale ∗L
X = Xo; H∗ Noile valori produse dupa iteratie ∗L
Nk = 2500; H∗ Numar maxim de iteratii ∗L
k = 0;
WhileAk ≤ Nk,
ForAi = 1, i ≤ n, i++,
i
j y
z
XPiT = j
j z ì APi,iT ;
j BPiT − ‚ APi,jT XPjT − ‚ APi,jT XoP jTz
z
i−1 n

k j=1 j=i+1 {
E; H∗ Gasirea noilor aproximatii ∗L
Xo = X; H∗ Pentru iteratia urmatoare ele devin valori vechi ∗L
k++E;
Print@ MatrixForm@ XDD;
Max@Abs@A.X − BDD

O imbunatatire in continuare a vitezei de convergenta se poate obtine folosind

metoda suprarelaxarii.
Daca in relatia (1.68) se aduna si se scade xi( k ) , incluzind termenul scazut in
23H

ultima suma se obtine relatia echivalenta xi( k )

 1 ⎛ i −1 n
⎞
⎜ i ∑ ij j
xi( k +1) = xi( k ) + ∑
( k +1)
b − a x − aij x (jk ) ⎟ = xi( k ) + Ri( k ) (1.69)
aii ⎝ j =1 j =i ⎠
Semnificatia acestei relatii este ca la valoarea din iteratia precedenta se aduna
un corector dat de restul membrului drept, Ri( k ) .
Uneori acest corector se ajusteaza cu un coeficient subunitar sau supraunitar
w , care asigura o convergenta marita, metoda numindu-se in acest caz suprarelaxare
si este descrisa de relatia
xi( k +1) = xi( k ) + wRi( k ) ; w ∈ [0, 2] (1.70)
sau desfasurat
w⎛ i −1 n ⎞
xi( k +1) = xi( k ) + ⎜ bi − ∑ aij x (jk +1) − ∑ aij x (jk ) ⎟ (1.71)
aii ⎝ j =1 j =i ⎠
Deoarece coeficientul de suprarelaxare w , optimal pentru un anumit sistem nu
poate fi cunoscut apriori, acest procedeu este util atunci cand trebuiesc rezolvate cat
mai rapid un numar mare de sisteme asemanatoare. Se experimenteaza pe un sistem
diferite valori ale coeficientului de suprarelaxare w cautandu-se pe cel care asigura
cea mai buna precizie dupa un numar dat de iteratii. Acest coeficient va fi considerat
in continuare optimal si pentru celelalte sisteme.
In programul 1.8 este exemplificata aplicarea suprarelaxarii pentru un sistem
de ecuatii liniare foarte rau conditionat (prima ecuatie si a doua sunt aproape
proportinale). Se constata ca, pentru un numar dat de iteratii, folosirea unui coeficient
w = 1.633 micsoreaza eroarea cu mai multe ordine de marime.

Programul 1.8 Rezolvarea unui sistem de ecuatii liniare prin iteratie Gauss-
Siedel cu suprarelaxare
Clear@"`∗"D;
A = 886, 3, 2<, 83, 1.5001, 1<, 82, −1, 3<<;
B = 81.2, −3.1, 2<;
Print@MatrixForm@ADD;
Print@MatrixForm@BDD;
n = Length@ BD;
Xo = Table@100., 8n<D; H∗ Valorile ghicite, initiale ∗L
X = Xo; H∗ Noile valori produse dupa iteratie ∗L
Nk = 2000; H∗ Numar maxim de iteratii ∗L
k = 0;
WhileAk ≤ Nk,
ForAi = 1, i ≤ n, i++,
i
j y
z
XPiT = XoPiT + 1.633 j
j z ì APi,iT ;
j BPiT − ‚ AP i,jT XPjT − ‚ APi,jT XoPjTz
z
i−1 n

k j=1 j=i {
E; H∗ Gasirea noilor aproximatii ∗L
Xo = X; H∗ Pentru iteratia urmatoare ele devin valori vechi ∗L
k++E;
Print@MatrixForm@XDD;
Max@Abs@A.X − BDD

1.9 Matrici tridiagonale. Eliminare Gauss

Exista numeroase cazuri in care sistemele de ecuatii sunt caracterizate de
matrici rare, adica matrici care au un numar important de elemente nule. Este de
asteptat ca astfel de sisteme sa poata fi rezolvate mai rapid, nefiind necesare
numeroase operatii aritmetice care implica aceste elemente nule. De la caza la caz,
astfel de situatii pot fi exploatate favorabil prin algoritmi specifici.
Exista si cazuri, cu un grad destul de mare de generalitate, pentru care s-au
dezvoltat deja algoritmi, in special daca elementele nenule prezinta anumitre simetrii
topologice. Cea mai raspandita situatie care implica matrici rare, atat in aplicatii
directe cat si provenita din alte metode numerice (se va vedea la fuctii spline si la
metoda diferentelor finite), este cea a matricilor tridiagonale, care va fi tratata in
continuare.
O matrice tridiagonala are elementele de pe diagonala principala si pe cele de
deasupra sau dedesubtul ei nenule iar restul sunt nule. Un sistem care implica o astfel
de matrice are deci forma
⎛ d1 c1 0............0 ⎞ ⎛ x1 ⎞ ⎛ b1 ⎞
⎜ ⎟ ⎜ ⎟ ⎜ ⎟
⎜ a1 d 2 c2 ............0 ⎟ ⋅ ⎜ x2 ⎟ = ⎜ b2 ⎟ (1.72)
⎜ ............................... ⎟ ⎜ ⎟ ⎜ ⎟
⎜ ⎟ ⎜ ⎟ ⎜ ⎟
⎝ 0 0 0...an −1 d n ⎠ ⎝ xn ⎠ ⎝ bn ⎠
Se observa ca matricea coeficientilor sistemului este
⎛ ............................................ ⎞
⎜ ⎟ ⎧di = aii ⎧di +1 = ai +1,i +1
0...... a d c ................0 ⎪ ⎪
A=⎜ ⎟
−1
ci = ai ,i +1 ⎨ai = ai +1,i
i i i
⎨ (1.73)
⎜ 0...............ai di +1 ci +1.....0 ⎟ ⎪ ⎪
⎜ ⎟ ⎩ai −1 = ai ,i −1 ⎩ci +1 = ai +1,i + 2
⎝ ............................................. ⎠
Avand in vedere ca majoritatea elementelor de sub diagonala sunt deja nule, se
poate aplica mult mai usor decat in cazul unei matrici dense (normale) metoda Gauss.
Punand conditia de anulare a unui element subdiagonal

ai' = 0 = ai − mi +1di (1.74)

rezulta multiplicator pentru linia i + 1
a
mi +1 = i (1.75)
di
unde di este elementul pivot.
Tinand cont ca operatia trebuie facuta pe intraga linie i + 1
di'+1 = di +1 − mi +1 ⋅ ci (1.76)
rezulta ca celelalte elemente care se modifica pe linia i + 1 sunt
(1.77)
a
di'+1 = di +1 − i ci (1.78)
di
a
bi'+1 = bi +1 − i bi (1.79)
di
In aceste relatii avem i = 1, 2,..., n − 1 .
Matricea devine si mai rara, si in plus superior triunghiulara
⎛ b1' ⎞
⎜ '⎟
⎛ ..................................... ⎞ ⎜ b2 ⎟
⎜ ⎟ ⎜ ⎟
............0 di' ci .......... ⎟
A' = ⎜ ; B '
= ⎜ ⎟ (1.80)
⎜ ............0 0 d ' c ⎟ ⎜ b '
⎟
⎜ i +1 i +1
⎟ ⎜ ⎟
i
⎜ ...................................... ⎟
⎝ ⎠ ⎜ ⎟
⎜ b' ⎟
⎝ n⎠
Prin urmare necunoscutele se pot determina din relatii cu numai trei termeni
 di' ⋅ xi + ci ⋅ xi +1 = bi' (1.81)
Deci si substitutia inversa se realizeaza mult mai rapid deoarece in relatiile
finale nu mai apar sume ca in cazul matricilor dense
b' − c x
xi = i 'i i +1 (1.82)
di
unde i = n − 1, n − 2,...,1 .
Evident, ultima necunoscuta este data de
b'
xn = n' (1.83)
dn
Programul 1.9 rezolva un sistem de ecuatii liniare cu matrice tridiagonala prin
metoda Gauss folosind relatiile (1.78), (1.79), (1.82) si (1.83), unde s-au notat cu
24H 25H 26H 27H

d [i ] elementele diagonale ale matricii superior triunghiulare obtinute.

Programul 1.9 Rezolvarea unui sistem de ecuatii liniare cu matrice

tridiagonala prin metoda Gauss
i2
j y
0z
j
j 0z
z
j z
5 0 0
j z
A= j
j
j 0z
z
z H∗ Matricea tridiagonala ∗L
3 4 2 0
j
j z
j
j 1z
z
z
0 1 5 2 ;
j0 z
k0 4{
0 3 3
0 0 3
b = 83, 2, 1, 4, 5<; H∗ Matricea termenilor liberi ∗L
B = b; H∗ Copie pentru verificare ∗L
n = Length@ bD;
Print@MatrixForm@ADD;
x = Table@0, 8n<D;
d@1D = A@@1, 1DD;
ForAi = 1, i ≤ n − 1, i++, H∗ Obtinerea matricii superior triunghiulare ∗L
A@@i + 1, iDD
d@i + 1D = A@@i + 1, i + 1DD − A@@i, i + 1DD;
d@iD
A@@i + 1, iDD
b@@i + 1DD = B@@i + 1DD − b@@iDD;
d@iD
E;

b@@nDD
x@@nDD =
d@nD
;

ForAi = n − 1, i ≥ 1, i−−, H∗ Substitutie inversa ∗L

b@@iDD − A@@i, i + 1DD x@@i + 1DD
x@@iDD =
d@iD
;

E;
x H∗ Afisarea solutiilor ∗L
A.x − B
H∗ Verificarea ∗L

1.10 Matrici tridiagonale. Factorizare Doolittle

Tot pentru cazul matricilor rare de tip tridiagonal prezentam algoritmul

Doolittle pentru care volumul de calcule se micsoreaza de asemenea foarte mult.
Pentru sistemul
A⋅X = B (1.84)
caracterizat de o matrice tridiagonala, se face factorizarea
 A = L⋅U (1.85)
unde matricile factori au formele
⎛ 1 0............0 0 ⎞ ⎛ u11 u12 0..............0 ⎞
⎜ ⎟ ⎜ ⎟
l21 1............0 0 ⎟ ; U = ⎜ 0 u22 u23 ............0 ⎟
L=⎜ (1.86)
⎜ ⎟ ⎜ ⎟
⎜ ⎟ ⎜ ⎟
⎝ 0 0............ln −1n −1 1 ⎠ ⎝0 0 0..............unn ⎠
si deci elementele lor sunt

⎧lij / .i = j , i = j + 1 ⎧uij / .i = j , i = j − 1
lij = ⎨ ; uij = ⎨ (1.87)
⎩0 / .i < j , i > j + 1 ⎩0 / .i > j , i < j − 1
Din relatia (1.85) rezulta pentru elementele de deasupra diagonalei
28H

n
ak ,k +1 = ∑ lkm , um ,k +1 = lkk uk ,k +1 (1.88)
m =1 1
unde s-a tinut cont de (1.87) din care rezulta ca m ∈ {k − 1, k} ∩ {k + 1, k} = {k} .
29H

Prin urmare, elementele nediagonale ale matricii U sunt chiar elementele

nediagonale ale matricii initiale A.
uk , k +1 = ak ,k +1 (1.89)
Tot din relatia (1.85) rezulta pentru elementele de dedesubtul diagonalei
30H

n
ak ,k −1 = ∑ lkmum ,k −1 = lk ,k −1uk −1,k −1 (1.90)
m =1
unde s-a tinut cont de (1.87) din care rezulta ca m ∈ {k − 1, k} ∩ {k − 2, k − 1} = {k − 1} .
31H

Prin urmare, elementele nediagonale ale matricii L sunt

a
lk ,k −1 = k , k −1 (1.91)
uk −1, k −1
Din relatia (1.85) se observa imediat ca primul element diagonal al matricii U
32H

este egal cu primul element diagonal al matricii A. a11 = u11 ⇒

u11 = a11 (1.92)

Din relatia pentru elementele diagonale

n
akk = ∑ lkmumk = lk ,k −1uk −1,k + lkk ukk (1.93)
m =1 1
unde s-a tinut cont de (1.87) din care rezulta m ∈ {k − 1, k} ∩ {k − 1, k} = {k − 1, k} ,
33H

rezulta celelalte elemente diagonale

ukk = akk − lk ,k −1 ⋅ uk −1, k (1.94)
ak −1,k

si deci
ukk = akk − lk , k −1 ⋅ ak −1, k (1.95)
Relatiile de factorizare (1.89), (1.91), (1.92) si (1.95) implica mult mai putine
34H 35H 36H 37H

calcule decat in cazul matricilor dense, deoarece bu implica sumari dupa n termeni.
De asemenea, ca si in cazul metodei Gauss pentru matrici tridiagonale, determinarea
necunoscutelor prin substitutie inversa si similar prin substitutie directa implica un
volum mult mai mic de calcule.
 Capitolul II
Interpolarea si aproximarea functiilor
Desi din punct de vedere semantic interpolarea are un sens diferit de
aproximare, in practica analizei numerice ele se bazeaza pe aceleasi procedee
matematice, astfel ca in cele ce urmeaza vor fi tratate unitar.
Dupa cum arata si nuantarea denumirilor, exista doua categorii importante de
aplicatii ale acestor procedee.
O prima categorie este intalnita in aplicatii de natura experimentala, atunci
cand se colecteaza o lista de date numerice, reprezentand valorile determinate
experimental ale unei functii Y = { f ( x1 ) , f ( x2 ) ,..., f ( xn )} , pentru anumite valori
(alese sau impuse) ale argumentului X = { x1 , x2 ,..., xn } , ca in figura 2.1. Daca dorim
sa determinam prin calcul valoarea functiei intr-un punct xi care nu apartine multimii
de valori experimentale, xi ∉ X , trebuie sa recurgem la un procedeu de interpolare.
Astfel de procedee sunt recomandate in special pentru cazul unor functii cu variatii
rapide si cu un numar redus de valori cunoscute, deoarece o simpla unire a punctelor
cunoscute prin segmente de dreapta poate furniza, dupa cum se observa in figura 2.1,
o valoare destul de departata de cea reala.

f(x)
Functia exacta
Functie aproximata
f(xi )
fa(x i)

xi x
x1 x2 x3 x4 x5 x6
Figura 2.1 Interpolarea unei functii intr-un punct

O a doua categorie o constituie aplicatiile cu caracter mai teoretic, in care o

functie, cunoscuta analitic dar foarte complicata, este folosita in calcule de mare
complexitate, si/sau cronofage, eventual fara solutii analitice in termeni de functii
elementare. In astfel de cazuri este utila aproximarea functiei de la bun inceput prin
functii mai simple care sa permita continuarea calculelor analitice sau cel putin
usurarea lor. O astfel de situatie este prezentata in figura 2.2, in care o functie dificil
de prelucrat analitic este aproximata printr-o functie mult mai comoda, de tip
polinomial.
 1.25

1.2

1.15

1.1

1.05

-0.75 -0.5 -0.25 0.25 0.5 0.75

0.95

0.9
Figura 2.2 Aproximarea functiei hypergeometrice 2 F1 (1/ 3,1/ 3, 2 / 3; x) (linie
continua) printr-un polinom de grad 9 (linie punctata)

Aproximarea unei functii complicate f ( x ) se poate face cu o combinatie

liniara de n functii elementare liniar independente ϕi
n
f ( x) = ∑ aiϕi ( x) (2.1)
i =1
Se pot folosi in acest scop monoame, functii trigonometrice elementare sau
diverse tipuri de polinoame ortogonale dupa cum se va prezenta mai departe.
Uneori este necesara o combinare a metodelor de rezolvare pentru cele doua
clase de probleme pentru aproximarea unei functii necunoscute ale carei valori sunt
determinate experimental intr-un numar de puncte dar sunt afectate de perturbatii si
erori. Cel mai raspandit procedeu folosit in astfel de cazuri este cel de minimizare
globala a erorilor in sensul celor mai mici patrate.
Mentionam de asemena ca, in afara de importanta practica, aceste metode au si
o importanta teoretica, fiind folosite de numeroase alte metode numerice mai
complexe (derivarea, rezolvarea ecuatiilor diferentiale si a celor integrale, etc.).

2.1 Interpolarea Newton

Cea mai raspandita alegere a functiilor elementare pentru formula de

aproximare sau interpolare (2.1) o constituie functiile de tip monom x i , astfel ca se
38H

obtine o functie polinomiala. Conform teoremei Weierstrass, orice functie continua pe

un interval finit [ a, b] poate fi aproximata cu o precizie oricat de buna printr-o functie
polinomiala.
O prima posibilitate de aproximare de acest tip o constituie interpolarea
Newton, care foloseste urmatoarea dezvoltare a unei functii f ( x ) definita pe un
interval [ a, b] :
f ( x) = a0 + a1 ( x − x1 ) + a2 ( x − x1 )( x − x2 ) +
(2.2)
+... + an −1 ( x − x1 )( x − x2 )...( x − xn −1 ) + Rn ( x) = Pn −1 ( x) + Rn ( x)
unde Pn −1 ( x) este polinomul aproximant, de grad n − 1 , iar Rn ( x) este restul
aproximarii, presupus foarte mic in modul si care este un polinom de grad n .
Coeficientii ai ai dezvoltarii sunt dati de urmatoarele formule care definesc asa
numitele diferente divizate DD ( xk , , xk −1 ,..., x1 ) :
a0 = f ( x1 );
f ( x2 ) − f ( x1 )
a1 = = DD( x2 , x1 );
x2 − x1
DD( x3 , x2 ) − DD( x2 , x1 )
a2 = = DD( x3 , x2 , x1 ); (2.3)
x3 − x1
...................................................................
DD( xn , xn −1 ,...x2 ) − DD( xn −1 , xn − 2 ,...x1 )
an −1 = = DD( xn , xn −1,..., x1 );
xn − x1
Prin urmare, formula de interpolare Newton poate fi scrisa sub forma:
n k −1 n
f ( x) = f ( x1 ) + ∑ DD ( xk , , xk −1 ,..., x1 )∏ ( x − xi ) + DD ( xn , xn −1 ,..., x1 , x)∏ ( x − xi ); (2.4)
k =2 i =1 i =1
Demonstratia se poate face prin inductie completa. Astfel pentru n = 1
formula este adevarata deoarece ea se poate scrie:
− f ( x1 ) + f ( x)
f ( x) = f ( x1 ) + ( x − x1 ) = f ( x); (2.5)
− x1 + x
Presupunand acum ca ea este valabila pentru n , putem sa o scriem pentru
n = n +1 :
n +1 k −1 n +1
f ( x) = f ( x1 ) + ∑ DD( xk , xk −1 ,...., x1 )∏ ( x − xi ) + DD( xn +1 , xn ,..., x1 , x)∏ ( x − xi ) (2.6)
k =2 i =1 i =1

Separand factorul cu indice i = n + 1 din produs si explicitand coeficientul

acestuia conform definitiilor diferentelor divizate (2.3), termenul final din (2.6)
39H 40H

devine:
n +1
DD( xn +1 , xn ,...., x1 , x)∏ ( x − xi )
i =1
(2.7)
DD( xn +1 , xn ,...., x1 ) − DD( xn , xn −1 ,...., x1 , x) n
= ( x − xn +1 )∏ ( x − xi )
xn +1 − x i =1

Separand si ultimul termenul cu indice k = n + 1 din suma din relatia (2.6) si 41H

schimband semnul si la numitorul si la numaratorul fractiei din (2.7) obtinem 42H

n k −1 n+ 1 − 1
f ( x) = f ( x1 ) + ∑ DD( xk , xk −1 ,..., x1 )∏ ( x − xi ) + DD( xn +1 , xn ,..., x1 ) ∏ (x − x ) i
k =2 i =1 i =1
(2.8)
− DD( xn +1 , xn ,...., x1 ) + DD( xn , xn −1 ,...., x1 , x) n
+ (− xn +1 + x) ∏ ( x − x1 )
− xn +1 + x i =1

Se observa ca aceasta relatie se reduce la relatia (2.4) pentru n , pe care am 43H

presupus-o adevarata, ceea ce incheie demonstratia.

Formula de interpolare Newton poate fi scrisa si ca formula de aproximare a
unei functii printr-un polinom:
n k −1
f ( x) ≈ f ( x1 ) + ∑ DD ( xk , , xk −1 ,..., x1 )∏ ( x − xi ) (2.9)
k =2 i =1
Restul interpolarii polinomiale este prin urmare ultimul termen din
dezvoltarea (2.4)
44H
 n
Rn ( x) = DD ( xn , xn −1 ,..., x1 , x)∏ ( x − xi ); (2.10)
i =1
Acesta se poate estima folosind functia ajutatoare:
n
Q (t ) = f (t ) − Pn −1 (t ) − DD ( xn , xn −1 ,..., x1 , x)∏ (t − xi ) (2.11)
i =1

Aceasta functie are pe intervalul [ a, b] un numar de n + 1 radacini deoarece

pentru t = x ea este nula conform relatiei (2.4), iar pentru cele n valori t = xi
45H

produsul se anuleaza si functia Q( xi ) se anuleaza de asemenea conform relatiei (2.4) 46H

t = x; Q (t ) = 0;
(2.12)
t = xi ; Q ( xi ) = 0; i = 1, 2,..., n
Daca Q (t ) are n + 1 radacini, inseamna ca derivata de ordinul intai Q '(t ) are
n radacini, derivata de ordinul doi Q ''(t ) are n − 1 radacini, si asa mai departe, iar
derivata de ordinul n are o radacina. Fie ξ ∈ [ a, b] aceasta radacina
Q ( n ) (ξ ) = 0 (2.13)
Luad t = ξ in relatia (2.11) si tiand cont ca derivata de ordinul n in raport cu
47H

t a produsului din (2.11) este n! , obtinem

48H

f ( n ) (ξ ) − P ( n ) n −1 (ξ ) − DD( xn , xn −1 ,..., x1 , x)n ! = 0 (2.14)

unde al doilea termen este nul deoarece este derivata de ordin n a unui polinom de
grad n − 1 .
Rezulta
f ( n ) (ξ )
DD( xn , xn −1 ,..., x1 , x) = (2.15)
n!
si inlocuind in relatia (2.10) se obtine formula
49H

f ( n ) (ξ ) n
Rn ( x) =
n ! i =1
∏ ( x − xi ) (2.16)

Aceasta este doar o estimare a restului, deoarece nu se cunoaste valoarea

exacta ξ pentru care se anuleaza derivata de ordin n a functiei Q (t ) , dar arata ca
valoarea absoluta a restului scade foarte rapid cu n datorita factorialului de la
numitor. Totusi, in zonele in care functia variaza foarte rapid, cum ar fi de exemplu la
capetele intervalelor (dincolo de care functia ar fi identic nula), este posibil ca si
numaratorul sa aiba valori mari iar restul sa devina important, contribuind la asa-
numitul fenomen Runge care va fi prezentat mai departe.
Programul 2.1 prezinta aproximarea unei functii printr-un polinom prin
metoda de interpolare Newton, iar programul 2.2 prezinta aproximarea unei functii
cunoscuta doar prin un numar finit de valori numerice in puncte de esantionare
neechidistante.
 Programul 2.1 Aproximarea unei functii printr-un polinom de interpolare
Newton de grad 15
f@x_D = 3 Cos@2 x + 1D − 2 SinA − 3E; H∗ functia de aproximat ∗L
x

g1 = Plot@f@xD, 8x, 0, 2 π<, PlotStyle → BlueD; H∗ graficul functiei ∗L

2

n = 16; H∗ numarul de termeni ai polinomului ∗L

ForAk = 1, k <= n, r@kD = k; k++E; H∗ punctele de esantionare∗L
2. π
n

ForAi = 1, i ≤ n,
f@r@i + 1DD − f@r@iDD
DD@i + 1, iD = ; H∗ diferentele finite directe ∗L
r@i + 1D − r@iD
i++E;

ForAj = 3, j <= n,
ForAi = 1, i < j − 1,
DD@j, j − iD − DD@j − 1, j − i − 1D
DD@j, j − i − 1D = ; H∗diferente finite recursive∗L
r@jD − r@j − i − 1D
i++E;
j++E;

p@x_D = f@r@1DD + „ DD@k, 1D ‰ Hx − r@lDL; H∗ Polinomul de interpolare ∗L

n k−1

k=2 l=1

g2 = Plot@p@xD, 8x, 0, 2 π<, PlotStyle → RedD;

Show@g1, g2D; H ∗ graficele ambelor functii ∗L
Print@"PHxL=", Expand@p@xDDD;

Cele doua functii afisate in final au grafice perfect suprapuse:

1 2 3 4 5 6

-1

Polinomul de grad 15 obtinut de program este urmatorul:

PHxL=1.88067 − 3.86244 x − 4.01909 x2 + 4.95796 x3 − 1.28767 x4 +
1.74569 x5 − 1.9508 x6 + 1.07521 x7 − 0.418778 x8 + 0.134868 x9 − 0.0341521 x10 +
0.00618301 x11 − 0.00075202 x12 + 0.0000580942 x13 − 2.57855 × 10−6 x14 + 5.02162 × 10−8 x15
 Programul 2.2 Aproximarea unei functii cunoscuta intr-un numar de puncte,
printr-un polinom de interpolare Newton

r = 81., 3, 4, 6, 8, 9, 11, 12, 14<; H∗ punctele de esantionare ∗L

y = 82, −1, 0, 3, 5, 7, 8, 7, 5<; H∗ valorile corespunzatoare ale functiei ∗L
n = Length@rD;
l = 8<;
For@k = 1, k ≤ n, k++, AppendTo@l, 8r@@kDD, y@@kDD<DD;
g1 = ListPlot@l, PlotJoined → True, PlotStyle → BlueD;
ForAi = 1, i <= n − 1,
y@@i + 1DD − y@@iDD
DD@i + 1, iD =
r@@i + 1DD − r@@iDD
;

i++E;
ForAj = 3, j <= n,
ForAi = 1, i < j − 1,
DD@j, j − iD − DD@j − 1, j − i − 1D
DD@j, j − i − 1D =
r@@jDD − r@@j − i − 1DD
;

H∗Print@"DD@",j,",",j−i−1,"D=",DD@j,j−i−1DD;∗L
i++E;
j++E;

pN@x_D = y@@1DD + „ DD@k, 1D ‰ Hx − r@@lDDL;

n k−1

k=2 l=1

g2 = Plot@ pN@xD, 8x, 1, Max@rD<, PlotStyle → RedD;

Show@g1, g2D;
Print@"PNHxL=", Expand@pN@xDDD;

Graficul obtinut in urma interpolarii polinomiale pentru punctele date este

urmatorul:

2 4 6 8 10 12 14

Polinomul de interpolare obtinut pentru punctele date este urmatorul:

PNHxL=− 37.5631 + 88.1934 x − 71.2328 x2 + 27.9811 x3 −

6.10102 x4 + 0.777721 x5 − 0.0575849 x6 + 0.00229166 x7 − 0.0000378749 x8
 2.2 Interpolarea Lagrange

O alta metoda de interpolare polinomiala, echivalenta cu cea prezentata

anterior dar beneficiind de formule mai simple este interpolarea Lagrange.
O functie f ( x ) continua pe intervalul [ a, b] si ale carei valori sunt cunoscute
in n puncte poate fi aproximata printr-un polinom de grad n − 1
f ( x) ≈ Pn −1 ( x) (2.17)
folosind formula
n n x−x
f ( x) ≈ Pn −1 ( x) = ∑ f ( xi )∏ j
(2.18)
i =1 j =1 xi − x j
j ≠i

Aceasta relatie se mai poate scrie si sub forma

n
f ( x) ≈ ∑ f ( xi )Li ( x) (2.19)
i =1

unde functiile Li ( x) sunt polinoamele Lagrange avand expresia

n x−x
Li ( x) = ∏ j
(2.20)
j =1 xi − x j
j ≠i

Ca si in cazul interpolarii Newton, trebuie sa existe egalitate stricta in relatia

(2.19) in punctele in care functia este cunoscuta, xi , iar in afara lor sa existe o
50H

diferenta cat mai mica, reprezentata de acelasi rest dat de relatia (2.16). 51H

n
f ( xi ) = ∑ f ( xi )Li ( xi ) (2.21)
i =1
n
f ( x) = ∑ f ( xi )Li ( x) + Rn ( x) (2.22)
i =1
Pentru demonstarea formulei, scriem mai intai forma desfasurata a
polinoamelor Lagrange
( x − x1 )( x − x2 )...( x − xk )...( x − xi −1 )( x − xi +1 )...( x − xn )
Li ( x) = (2.23)
( xi − x1 )( xi − x2 )...( xi − xk )...( xi − xi −1 )( x − xi +1 )...( xi − xn )
Se observa ca pentru x = xi rezulta Li ( xi ) = 1 , deoarece fiecare factor de la
numarator se simplifica cu unul de la numitor.
Pentru x = xk unde k ≠ i , exista la numarator un factor ( xk − xk ) = 0 , astfel ca
Li ( xk ) = 0 .
Rezulta urmatoarea proprietate a polinoamelor Lagrange:
Li ( xk ) = δ ik (2.24)
Calculand cu relatia (2.18) valoarea functiei polinomiale intr-un punct x j ,
52H

rezulta:
n
Pn −1 ( x j ) = ∑ f ( xi )δ ij = f ( x j ) (2.25)
i =1
astfel ca este demonstrata relatia (2.21). 53H

Aparent, formula de interpolare Lagrange este mai simpla decat cea de

interpolare Newton si este folosita mai des in calculele analitice. De remarcat insa ca
in calcule numerice interpolarea Newton poate sa fie mai eficienta, deoarece
coeficientii dezvoltarii se calculeaza recursiv, spre deosebire de cazul interpolarii
Lagrange.
Programul 2.3 ilustreaza metoda de aproximare si interpolare a unei functii
folosind polinoame Lagrange. In mediul Mathematica acest program este deosebit de
rapid datorita performantelor deosebite ale rutinelor de calculare a produselor.

Programul 2.3 Aproximarea unei functii folosind un polinom Lagrange de

gradul 20
è!!!!!!!!!!!
f@x_D = Cos@ Log@2 x − 3DD − SinA x2 + 1 E;
a = 3;
b = 9;
g1 = Plot@f@xD, 8x, a, b<, PlotStyle → BlueD;
n = 21;
r = Table@0, 8n<D;
AbsoluteTimingA ForAk = 1, k <= n, k++, r@@kDD = a + kE;
b−a
n
ii−1
j x − r@@jDD z yji n x − r@@jDD z y
p@x_D = „ f@r@@iDDD j
j
j‰
z
z
zj
j‰ z
z
z;E
n

@@ DD @@ DD j
i=1 kj=1 r i − r j { kj=i+1 r@@iDD − r@@jDD {
g2 = Plot@p@xD, 8x, a, b<, PlotStyle → RedD;
Show@g1, g2D;
Print@"PHxL=", Expand@ p@xDDD;

2.3. Algoritmul Neville

O schema practica mult mai eficienta de calcul a polinomului de interpolare
este data de algoritmul Neville, care va fii descris in continuare.
Notam mai intai cu P1 polinomul de gradul zero (o constanta) care trece prin
punctul { { f ( x1 ), x1} , deci P1 = f ( x1 ) , cu P2 polinomul de gradul zero (o constanta)
care trece prin punctul { { f ( x2 ), x2 } , deci P2 = f ( x2 ) , si asa mai departe pana la
Pn = f ( xn ) .
Din aceste polinoame se pot forma polinoame de grad unu (o dreapta) notate
astfel: cu P12 polinomul care trece prin punctele { f ( x1 ), x1} si { f ( x2 ), x2 } , cu P23
polinomul care trece prin punctele { f ( x2 ), x2 } si { f ( x3 ), x3} , si asa mai departe pana
la Pn −1,n polinomul care trece prin punctele { f ( xn −1 ), xn −1} si { f ( xn ), xn } .
Se continua acest proces prin formarea de polinoame din grad din ce in ce mai
mare, pana se ajunge la P1,2,...,n −1,n polinomul care trece prin punctele
{ f ( x1 ), x1} , { f ( x2 ), x2 } ,…, { f ( xn−1 ), xn−1} si { f ( xn ), xn } , adica tocmai polinomul
cautat.
Se poate alcatui deci o schema arborescenta in care fiecare coloana se
calculeaza pe rand pe baza coloanei precedente.
 x1 : f ( x1 ) = P1
P12
x2 : f ( x2 ) = P2 P123
P23 P1234 (2.26)
x3 : f ( x3 ) = P3 P234
P34
x4 : f ( x4 ) = P4

In aceasta schema, fiecare polinom de grad k se obtine din doua polinoame de

grad k − 1 dupa urmatoarea formula obtinuta pe baza diferentelor divizate:
( x − xi + k ) Pi (i +1)...(i + k −1) + ( xi − x) P( i +1)( i + 2)...(i + k )
Pi (i +1)...(i + k ) = (2.27)
xi − xi + k
Programul 2.4 ilustreaza calculul polinomului de interpolare prin algoritmul
Neville. Trebuie mentionat totusi ca acest algoritm este recomandat pentru limbaje de
programare de uz general, unde viteza va fi intr-adevar mai mare decat cea a
interpolarii Newton sau Lagrange. Folosind programul Mathematica, viteza obtinuta
cu formula de interpolare Lagrange este superioara celorlalte metode datorita
eficientei deosebite de calcul a produselor care apar in aceasta metoda.

Programul 2.4 Aproximarea unei functii prin polinomului de interpolare

obtinut prin algoritmul Neville.
f@x_D := Sin@x + 1D;
a = 0.;
b = 4 π;
g1 = Plot@f@xD, 8x, a, b<, PlotStyle → BlueD;
M = 21;
r = Table@0, 8 M<D;
AbsoluteTimingAForAk = 1, k ≤ M, k++,
r@@kDD = a + k; d@k, kD = f@r@@kDDDE;
b−a

ForA m = 1, m <= M, m ++,

M

ForAi = 1, i <= M − m, i++,

Hx − r@@i + mDDL d@i, i + m − 1D − Hx − r@@iDDL d@i + 1, i + mD
d@i, i + mD = ExpandA EE;EE
r@@iDD − r@@i + mDD
p = d@1, MD;
g3 = Plot@ Evaluate@pD, 8x, a, b<, PlotStyle → RedD;
Show@g1, g3D;
Print@"PHxL=", pD;

2.4 Fenomenul Runge

Atunci cand se incearca aproximarea unei functii oarecare printr-o functie

polinomiala unica pentru intregul domeniu de definitie [ a, b ] se constata de regula ca
erorile sunt mai mari inspre capete. Uneori este chiar posibil ca erorile spre capete sa
creasca odata cu cresterea gradului polinomului de interpolare, ceea ce face ca
metodele Newton sau Lagrange sa nu fie convergente, in sensul ca erorile nu pot fi
micsorate oricat de mult. Acesta este asa-numitul fenomen Runge si este intalnit in
special in cazul in care punctele in care egalitatea functiei cu polinomul aproximant
(punctele de esantionare) sunt echidistante.
Fenomenul Runge este ilustrat de programul 2.5, care evidentiaza (chiar
pentru functia data ca exemplu de Runge) erorile spre capetele intervalului. Rularea
programului cu un numar din ce in ce mai mare de puncte intensifica aceste erori.

Programul 2.5 Aparitia fenomenului Runge la o interpolare polinomiala

globala folosind puncte de esantionare echidistante
f@x_D = ; H∗ functia Runge ∗L
1
x2 + 1
a = −3.0;
b = 3.0;
g1 = Plot@f@xD, 8x, a, b<, PlotStyle → BlueD; H∗ graficul functiei ∗L
n = 10; H∗ numarul de termeni ai polinomului ∗L
ForAk = 1, k <= n, r@kD = a + k; k++E; H∗ punctele de esantionare∗L
b− a
n

ForAi = 1, i ≤ n,
f@r@i + 1DD − f@r@iDD
DD@i + 1, iD = ; H∗ diferentele finite directe ∗L
r@i + 1D − r@iD
i++E;

ForAj = 3, j <= n,
ForAi = 1, i < j − 1,
DD@j, j − iD − DD@j − 1, j − i − 1D
DD@j, j − i − 1D = ; H∗diferente finite recursive∗L
r@jD − r@j − i − 1D
i++E;
j++E;

p@x_D = f@r@1DD + „ DD@k, 1D ‰ Hx − r@lDL; H∗ Polinomul de interpolare ∗L

n k−1

k=2 l=1

g2 = Plot@p@xD, 8x, a, b<, PlotStyle → RedD;

Show@g1, g2D; H∗ graficele ambelor functii ∗L
Print@"PHxL=", Expand@p@xDDD;

Graficul de interpolare in puncte echidistante a functiei Runge (cea care a fost

exemplificata prima oara in legatura cu acest fenomen) printr-un polinom cu 10
termeni este urmatorul:

1.5

0.5

-3 -2 -1 1 2 3

Daca se creste numarul de termeni ai polinomului la 20, graficul de interpolare

in puncte echidistante a functiei Runge este urmatorul:
 1.5

0.5

-3 -2 -1 1 2 3
-0.5

-1

Se observa ca desi in portiunea centrala erorile scad cand se creste ordinul

polinomului de interpolare, spre capetele intervalului erorile tind sa creasca, relevand
fenomenul Runge.
Desigur, o solutie ar fi sa se foloseasca o functie polinomiala distincta pentru
fiecare subinterval, prin procedeul functiilor spline care va fi prezentat mai tarziu.
Totusi exista solutii pentru evitarea fenomenului Runge si printr-o interpolare
globala, printr-o alegere speciala a punctelor de esantionare a domeniului.

2.4 Utilizarea polinoamelor Chebyshev si teorema mini-max

Erorile la interpolare si aproximare sunt date de restul Rn ( x) dat de relatia

(2.16)
54H

f ( n ) (ξ ) n
Rn ( x) =
n ! i =1
∏ ( x − xi )
Deoarece, asa cum s-a aratat nu se cunoaste punctul ξ si valoarea derivatei de
ordinul n in acest punct, singura cale de a se minimiza restul este sa se aleaga
corespunzator punctele de esantionare xi astfel incat produsul
n
pn ( x) = ∏ ( x − xi ) (2.28)
i =1
sa fie minim.
Acesta este un polinom monic (polinom cu coeficientul termenului de gradul
cel mai inalt egal cu 1) de gradul n . In legatura cu acest tip de polinoame exista
urmatoarea teorema:
Teorema mini-max: dintre toate polinoamele monice de grad n ,
polinoamele monice Chebyshev au cea mai mica valoare absoluta maxima pe
intervalul [ −1,1] .
Rezulta ca daca punctele de esantionare xi sunt chiar radacinile unui polinom
monic Chebyshev de grad n , produsul (2.28) este un astfel de polinom si deci restul
55H

interpolarii polinomiale va avea cea mai mica valoare maxima posibila pe intervalul
[ −1,1] . Printr-o schimbare convenabila de variabila, acest procedeu poate fi extins
pentru orice interval finit [ a, b ] si teoretic chiar si pentru intervale cu una din limite
infinita.
Demonstratia acestei teoreme se va face dupa ce se trec in revista anumite
proprietati utile ale polinoamelor Chebyshev.
Polinoamele Chebyshev se definesc pe intervalul [ −1,1] prin urmatoarea
relatie:
 Tn ( z ) = cos(n arccos( z )) (2.29)
Notand
θ = arccos z; z ∈ [ −1,1] (2.30)
rezulta forma echivalenta:
Tn ( z ) = cos nθ (2.31)
unde
z = cos θ (2.32)
Pornind de la aceasta definitie rezulta ca primele doua polinoame Chebyshev
sunt:
T0 ( z ) = cos 0θ = 1
(2.33)
T1 ( z ) = cos θ = z
Restul polinoamelor Chebyshev se pot calcula si prin relatia de recurenta:
Tn ( z ) = 2 zTn −1 ( z ) − Tn − 2 ( z ); n ≥ 2 (2.34)
Relatia (2.34) se poate demonstra pe baza urmatoarei identitati
56H

trigonometrice:
1
cos(α ) cos( β ) = [cos(α + β ) + cos(α − β )] (2.35)
2
Luand in aceasta relatie α = θ si β = (n − 1)θ se obtine:
1
cos θ cos(n − 1)θ = [cos nθ + cos(n − 2)θ ] (2.36)
z
T (z)
2 Tn ( z ) T (z)
n−1 n− 2

de unde se deduce relatia (2.34). 57H

Se poate demonstra ca functiile Tn ( z ) sunt polinoame de grad n , avand

urmatoarea forma generala:
Tn ( z ) = 2n −1 z n + Rn −1 ( z ) (2.37)
unde Rn −1 ( z ) este un polinom de grad n − 1 . Vom folosi notatia Rn −1 ( z ) ∈ {1n −1 ( z )} ,
unde {1n −1 ( z )} este multimea polinoamelor de grad cel mult n − 1 .
Demonstratia se poate face prin inductie completa. Intr-adevar pentru n = 1
relatia (2.37) este evidenta. Presupunand acum ca ea este valabila pentru n , in cazul
58H

n → n + 1 , tinand cont de relatia (2.34) se obtine:

59H

Tn +1 ( z ) = 2 zTn ( z ) − Tn −1 ( z ) = 2 z 2n −1 z n + 2 zRn −1 ( z ) − 2n − 2 z n −1 − Rn − 2 ( z ) = 2n z n +1 + Rn ( z ) (2.38)

Rn ( z )

ceea ce arata ca relatia este valabila si in acest caz.

Cele n radacini ale polinomului Chebashev de grad n se pot determina
simplu din relatia de definitie
π
Tn ( zk ) = 0; ⇒ cos nθ k = 0; nθ k = (2k − 1) ; k = 1, 2,..., n (2.39)
2 n radacini

⎛ 1 ⎞π
θk = ⎜ k − ⎟ ; (2.40)
⎠ ⎝ 2 n
⎛ 1 ⎞π
zk = cos ⎜ k − ⎟ (2.41)
⎝ 2⎠ n
Extremele polinomului au valoarea 1 sau -1:
Tn ( z ) = cos nθ ∈[−1,1] (2.42)
Ele se pot obtine simplu punand conditia de extrem in definitia (2.31): 60H
 π
| cos nθ k |= 1 ⇒ nθ k = kπ ; θ k = k (2.43)
n
⎛ π⎞
zk = cos ⎜ k ⎟ ; k = 0,1, 2,..., n (2.44)
⎝ n⎠ n +1 extreme

De remarcat ca cele n + 1 valori ale lui zk date de acesta relatie sunt ordonate
descrescator, functia cosinus fiind descrescatoare de la 1 la -1 pe intervalul [ 0, π ] .
Astfel, pentru k = 0 avem z0 = 1 iar pentru k = n avem zn = −1 .
Inlocuind θ k in (2.42) rezulta valorile in extreme:
61H

Tn ( zk ) = (−1) k (2.45)
Conform celor precizate in alineatul precedent, rezulta ca la extrema dreapta a
domeniului de definitie avem extremul z0 = 1 , pentru care valoarea polinomului este
intotdeauna (−1)0 = 1 , in timp ce la extrema stanga avem extremul zn = −1 , pentru
care valoarea polinomului este 1 daca gradul este par si -1 daca gradul este impar.
Programul 2.6 calculeaza prin recurenta primele 7 polinoame Chebyshev si
afiseaza graficele si expresiile lor individuale, precum si graficele suprapuse. Din
aceste reprezentari grafice se pot observa pozitiile celor n radacinilor si ale celor
n + 1 extreme.

Programul 2.6 Calculul si reprezentarea polinoamelor Chebyshev

T@0D = 1;
T@1D = x;
T@n_D := 2 x T@n − 1D − T@n − 2D;
gr = 8<;
ForAn = 0, n ≤ 6, n++, Print@Expand@T@nDDD
AppendToAgr, PlotAT@nD, 8x, −1, 1<,
PlotStyle → RGBColorA , EEEE;
n n 10 − n
,
Show@grD;
10 10 10

Reprezentarea pe un singur grafic a primelor 7 polinoame Chebyshev,

generata de acest program este urmatoarea:
1

0.5

-1 -0.5 0.5 1

-0.5

-1

De asemenea programul genereaza urmatoarele reprezentari pe grafice

individuale pentru aceste polinoame:
1
2

1.5

0.5

-1 -0.5 0.5 1

x
1

0.5

-1 -0.5 0.5 1

-0.5

-1

−1 + 2 x2
1

0.5

-1 -0.5 0.5 1

-0.5

-1

3
−3 x + 4 x
1

0.5

-1 -0.5 0.5 1

-0.5

-1

1 − 8 x2 + 8 x4
1

0.5

-1 -0.5 0.5 1

-0.5

-1

5 x − 20 x3 + 16 x5
1

0.5

-1 -0.5 0.5 1

-0.5

-1

2 4 6
−1 + 18 x − 48 x + 32 x
1

0.5

-1 -0.5 0.5 1

-0.5

-1
 Numim polinoame monice acele polinoame care au coeficientul termenului de
gradul cel mai inalt egal cu 1. Din relatia (2.37) se observa ca se poate obtine simplu
62H

un polinom monic Chebyshev de gradul n , pe care il vom nota in continuare T n ( z ) ,

prin impartirea unui polinom Chebyshev cu 2n −1
T ( z)
Tn ( z ) = n n −1 = 21− n Tn ( z ) = z n + 21− n Rn −1 ( z ) (2.46)
2
Conform definitiei (2.31) se poate scrie
63H

Tn ( z ) ≤ 1 (2.47)
deci polinomul monic Chebyshev are proprietatea
1
Tn ( z ) ≤ n −1 = 21− n (2.48)
2
Demonstratia teoremei mini-max.
Aceasta demonstratie se va face prin reducere la absurd. Sa presupunem ca ar
exista un polinom monic f n ( z ) cu proprietatea
f n ( z ) < Tn ( z ) ≤ 21− n (2.49)
Daca f n ( z ) este un polinom monic se poate scrie
f n ( z ) = z n + Fn −1 ( z ) (2.50)
cu Fn −1 ( z ) ∈ {1n −1 ( z )}
Sa formam functia ajutatoare
Q( z ) = Tn ( z ) − f n ( z ) (2.51)
Tinand seama de relatiile (2.46) si (2.50) se observa ca avem
64H 65H

Q( z ) = 21− n Rn −1 ( z ) − Fn −1 ( z )
Q( z ) ∈ {1n −1 ( z )} (2.52)
deci functia Q( z ) este un polinom de grad cel mult n − 1 .
Sa studiem comportarea acestei functii intre doua extreme consecutive ale
polinomului monic Chebyshev. Evisent, acest polinom monic are tot n + 1 extreme ca
si polinomul Chebyshev normal.
Orice interval este intre doaua extreme consecutive are la un capat un extrem
cu numar par 2k si la celalalt un extrem cu numar impar 2k ± 1 .
In extremul cu numar par, conform relatiei (2.45), rezulta ca polinomul monic
66H

este pozitiv
(−1) 2 k
Tn (ξ 2 k ) = n −1 > 0 (2.53)
2
deci
Tn (ξ 2 k ) = Tn (ξ 2 k ) (2.54)
iar functia ajutatoare este pozitiva
Q(ξ 2 k ) = Tn (ξ 2 k ) − f n (ξ 2 k ) = Tn (ξ 2 k ) ∓ f n (ξ 2 k ) > 0 (2.55)
deoarece conform ipotezei modulul lui f n ( z ) este mai mic decat al polinomului
monic Chebyshev.
Pe de alta parte, in extremul cu numar impar, conform relatiei (2.45), rezulta 67H

ca polinomul monic este negativ

 (−1) 2 k ±1
Tn (ξ 2 k ±1 ) = <0 (2.56)
2n −1
deci
Tn (ξ 2 k ±1 ) = −Tn (ξ 2 k ±1 ) (2.57)
iar functia ajutatoare este negativa
Q(ξ 2 k ±1 ) = Tn (ξ 2 k ±1 ) − f n (ξ 2 k ±1 ) = − Tn (ξ 2 k ±1 ) ∓ f n (ξ 2 k ±1 ) < 0 (2.58)
deoarece conform ipotezei modulul lui f n ( z ) este mai mic decat al polinomului
monic Chebyshev.
Deoarece intr-un capat al intervalului functia Q( z ) este pozitiva iar in celalalt
capat este negativa, inseamna ca ea ar avea cel putin o radacina in acest interval.
Deoarece exista n + 1 extreme, intre acestea se afla n asemenea intervale si deci
functia Q( z ) ar avea cel putin n radacini.
Aceasta insa este imposibil deoarece, dupa cum s-a aratat, functia Q( z ) este
un polinom de grad cel mult n − 1 , ceea ce arata ca ipoteza este falsa si deci este
adevarata contrara ei:
f n ( z ) ≥ Tn ( z ) ≤ 21− n (2.59)
ceea ce demonstreaza teorema mini-max.
In concluzie, exista posibilitatea minimizarii erorilor de interpolare
polinomiala globala a unei functii pe intervalul [ −1,1] prin algoritmi de tip Lagrange
sau Newton prin alegerea punctelor de esantionare a abscisei (nodurile) in radacinile
unui polinom Chebyshev de ordin n data de relatia (2.41). Daca avem un interval
68H

oarecare [ a, b ] , se poate utiliza o schimbare de variabila de tipul

b−a b+a
x+ z= (2.60)
2 2
astfel ca nodurile de interpolare sunt date de relatia
b+a b−a ⎡⎛ 1 ⎞π ⎤
z= + cos ⎢⎜ k − ⎟ ⎥ , k = 1, 2,..., n (2.61)
2 2 ⎣⎝ 2⎠ n⎦
Programul 2.7 realizeaza interpolarea polinomiala Newton a unei functii (in
acest exemplu, functia Runge) prin esantionarea abscisei in radacinile unui polinom
Chebyshev. Singura deosebire fata de cazul esantionarii echidistante (programul 2.5)
este calculul nodurilor folosind relatia (2.61), conform teoremei mini-max.
69H
 Programul 2.7- Interpolare polinomiala globala a functiei Runge cu
esantionare in noduri Chebashev
f@x_D = ; H∗ functia Runge ∗L
1
x2 + 1
a = −3.0;
b = 3.0;
g1 = Plot@f@xD, 8x, a, b<, PlotStyle → BlueD; H∗ graficul functiei ∗L
n = 10; H∗ numarul de termeni ai polinomului ∗L
ForAk = 1, k <= n, r@kD = CosAJk − N E; k++E; H∗ punctele de esantionare∗L
a+ b b − a 1 π
+
2 2 2 n

ForAi = 1, i ≤ n,
f@r@i + 1DD − f@r@iDD
DD@i + 1, iD = ; H∗ diferentele finite directe ∗L
r@i + 1D − r@iD
i++E;

ForAj = 3, j <= n,
ForAi = 1, i < j − 1,
DD@j, j − iD − DD@j − 1, j − i − 1D
DD@j, j − i − 1D = ; H∗diferente finite recursive∗L
r@jD − r@j − i − 1D
i++E;
j++E;

p@x_D = f@r@1DD + „ DD@k, 1D ‰ Hx − r@lDL; H∗ Polinomul de interpolare ∗L

n k−1

k=2 l=1

g2 = Plot@p@xD, 8x, a, b<, PlotStyle → RedD;

Show@g1, g2D; H∗ graficele ambelor functii ∗L
Print@"PHxL=", Expand@ p@xDDD;

Graficul de interpolare a functiei Runge in puncte date de radacinile

polinomului Chebayshev de grad 10, printr-un polinom cu 10 termeni este urmatorul:

0.8

0.6

0.4

0.2

-3 -2 -1 1 2 3

Se observa ca in regiunea centrala erorile sunt mai mari decat in cazul

esantionarii echidistante (datorita densitatii mai mici a nodurilor), in schimb spre
capete erorile sunt substantial mai mici decat in cazul respectiv.
Daca se creste numarul de termeni ai polinomului la 20, graficul de interpolare
a functiei Runge in puncte date de radacinile polinomului Chebyshev de grad 20 este
urmatorul:
 1

0.8

0.6

0.4

0.2

-3 -2 -1 1 2 3

Se observa ca erorile scad foarte mult pe tot intervalul atunci cand se mareste
gradul polinomului de interpolare, ceea ce arata ca, spre deosebire de cazul
interpolarii echidistante, procesul este convergent.

2.5 Interpolarea cu functii spline

In cazul in care interpolarea cu un polinom pe intregul interval de interes

conduce sau la erori prea mari sau la un polinom de grad prea mare pentru a fi usor
manipulat, se poate recurge la impartirea in subintervale si aproximarea pe fiecare
subinteval cu un polinom distinct dar de grad mic.
Cea mai simpla varianta este interpolarea cu un polinom de grad 1, adica o
dreapta care uneste doua puncte consecutive cunoscute. Daca aceste puncte sunt
suficient de multe, erorile la aproximarea functiei necunoscute prin segmente de
dreapta pot fi acceptabile, si intre punctele cunoscute se poate face o interpolare
liniara.
O ilustrare a procedeului este data de programul urmator, in care se face
reprezentarea cu interpolare liniara a unei liste de valori { yi } determinate
experimental ale unei functii necunoscute pentru anumite valori ale argumentului
{ xi } .
Clear@"`∗"D;
x = 80, 1.5, 2, 3.5, 4, 5, 6, 7, 8, 9<; H∗ abscisele in care este cunoscuta functia ∗L
y = 8−2, 1, 3, 5, 2, 1, 2, −1, 3, −1<;H∗ valorile functiei in punctele cunoscute ∗L
g1 = ListPlot@lstD; H∗ reprezentarea prin puncte distincte ∗L
g2 = ListPlot@lst, PlotJoined → TrueD; H∗Reprezentarea prin segmente∗L
Show@g1, g2D; H∗ Graficul compus din puncte unite prin segmente ∗L

Rezulta urmatoarea reprezentare grafica:

5

2 4 6 8

-1

-2
 In schimb, daca functia are o variatie rapida iar numarul de puncte in care ea
este cunoscuta este mic erorile pot creste foarte mult. In figura urmatoare este
reprezentata o aproximare prin segmente de dreapta a unei functii sinusoidale
cunoscuta doar in 6 puncte. Se observa ca la anumite valori ale argumentului (de
exemplu la 0,5 sau 2,5) erorile care apar la o interpolare liniara sunt foarte mari, de
ordinul zecilor de procente.
1

0.5

1 2 3 4 5

-0.5

-1

Chiar si prin marirea numarului de puncte la 12 erorile sunt in cele mai multe
cazuri inacceptabile, dupa cum se poate observa in figura urmatoare:
1

0.5

2 4 6 8 10

-0.5

-1

Desigur, daca s-ar folosi pentru aproximarea intre doua puncte a unui polinom
de grad superior, rezultatele ar fi mult mai bune, dar procedeul ar putea fi prea
complicat din punct de vedere numeric.
In practica, cea mai avantajoasa metoda o constituie folosirea unor polinoame
de gradul 3, ceea ce conduce la un compromis optim intre precizie si complexitate.
Aceasta este metoda interpolarii spline cubice, denumire provenita de la forma
functiei pe fiecare portiune portiune a domeniului si de la gradul polinomului folosit.
Metoda presupune mai intai impartirea domeniului [ x0 , xn ] in mai multe
subintervale:
[ x0 , xn ] = [ x0 , x1 ] ∪ [ x1 , x2 ] ∪ ... ∪ [ xi , xi +1 ] ∪ ... ∪ [ xn −1 , xn ] (2.62)
Pe fiecare subinterval [ xi , xi +1 ] , functia este aproximata printr-un polinom de
gradul trei Si ( x) de forma
⎧⎪a ( x − xi )3 + bi ( x − xi ) 2 + ci ( x − xi ) + di , x ∈ [ xi , xi +1 ]
Si ( x ) = ⎨ i (2.63)
⎪⎩0, x ∉ [ xi , xi +1 ]
 i = 0,1,..., n − 1
Pe ansamblu functia este aproximata printr-o suma S ( x) de astfel de functii cu
domenii de definitie disjuncte
n −1
f ( x ) ≈ S ( x ) = ∑ Si ( x ) (2.64)
i =0
Pentru a se putea calcula cei 4n coeficienti din relatiile (2.63) se vor impune
70H

conditii de racordare in punctele intermediare: functiile si primele lor doua derivate

trebuie sa fie continue in aceste n − 1 puncte. Rezulta 3(n − 1) ecuatii, care impreuna
cu cele n + 1 ecuatii care dau valorile cunoscute ale functiei in cele n + 1 puncte
(incluzand si cele de la capetele intervalului) dau un total de 4n − 2 ecuatii.
Pentru ca sistemul sa fie determinat se aleg (subiectiv) doua valori pentru
derivata a doua la capetele intervalului, sistemul devenind astfel complet. Cea mai
uzuala alegere este ca derivata a doua sa fie 0 la capetele intervalului, caz in care
avem asa numita interpolare spline bicubica naturala.

Capitolul 3
Rezolvarea ecuatiilor neliniare
Problema gasirii acelor valori ale argumentului unei functii reale de argument
real care anuleaza functia, adica a gasirii punctelor de intersectie a functiei cu abscisa
(fig 3.1), este echivalenta cu rezolvarea ecuatiei:
F ( x) = 0 , F : → (3.1)
Valorile respective sunt numite radacinile ecuatiei (3.1). 71H

Problema admite solutii analitice pentru toate cazurile in care functia este un
polinom de grad cel mult 4, pentru anumite functii elementare sau speciale, si anumite
combinatii ale acestora. In foarte multe cazuri insa, nu este posibila gasirea unor
solutii analitice si este necesar sa se utilizeze metode numerice pentru obtinerea
radacinilor.
Exista o mare varietate de metode generale pentru rezolvarea numerica a
acestei probleme, caracterizate prin diverse grade de complexitate, eficienta si viteza.
Sunt de asemenea cunoscute o serie de metode speciale, recomandabile pentru cazuri
in care se pot exploata favorabil o serie de particularitati ale functiei in cauza.
Pentru optimizarea numarului de operatii necesare pentru obtinerea cu precizie
ridicata a tuturor radacinilor, se procedeaza de cele mai multe ori in doua etape, in
special atunci cand se presupune ca exista mai multe radacini distincte. Este
recomandabil de asemenea ca, daca mediul de programare folosit permite usor acest
lucru, sa se vizualizeze mai intai graficul functiei respective pe intervalul de interes,
ceea ce ofera o imagine asupra pozitiei radacinilor si numarului lor.
Intr-o prima etapa, se vor stabili subintervalele in care sunt plasate radacinile,
pentru a restrange procesul de cautare a valorilor exacte in domenii cat mai inguste, in
care gasirea valorilor cu o precizie mare sa se faca printr-un numar cat mai redus de
iteratii. De exemplu, in figura 3.1 se observa ca radacina x3 este in intervalul
[4h,5h].Vom numi in continuare acest proces etapa gasirii intervalelor (engl.
braketing).
 In a doua etapa, se va considera fiecare interval in parte, si prin metode
numerice cat mai eficiente se va gasi radacina din intervalul respectiv cu precizia
dorita. Vom numi in continuare acest proces etapa gasirii radacinilor (engl.
refining).
Dat fiind ca unele metode considerate foarte eficiente in majoritatea cazurilor
pot da rezultate eronate pentru anumite cazuri particulare, este obligatoriu ca dupa
gasirea radacinilor sa se faca verificarea acestora, introducand valorile respective in
ecuatia (3.1) si observand abaterea fata de 0.
72H

3.1 Gasirea intervalelor prin esantionare

Intervalele in care se afla cate o radacina pot fi determinate vizual pe baza
analizei graficului functiei, daca acesta poate fi reprezentat. Analitic, acest lucru este
posibil pe baza urmatoarei consecinte a teoremei valorii medii:
• Daca o functie este monotona si continua pe un interval [a,b] si indeplineste
conditia F (a ) F (b) < 0 , atunci in intervalul respectiv se afla o radacina.
Prin urmare, este posibila gasirea intervalelor in care se afla cate o radacina
daca se esantioneaza abscisa cu pasul h, se calculeaza valorile functiei in punctele
respective si produsul acestor valori in puncte consecutive. Daca un astfel de produs
este negativ, intre cele doua puncte de esantionare respective se afla cel putin o
radacina. Se observa ca acest procedeu este de fapt echivalent construirii sirului Rolle.
Radacina este unica daca functia este monotona (conditie suficienta, nu si
necesara), dar acest lucru nu este intotdeauna cunoscut. In cazul in care in intevral ar
exista mai multe radacini, majoritatea metodelor numerice vor gasi numai una,
celelalte pierzandu-se. Pentru a evita aceasta exista mai multe posibilitati: se poate
reprezenta grafic functia, se poate studia si semnul derivatei daca se dispune de
aceasta, sau se poate alege pasul de esantionare suficient de fin pentru a se surprinde
toate trecerile prin zero ale functiei. Uneori, numarul radacinilor poate fi dedus din
considerente teoretice (de exemplu in cazul functiilor polinomiale, desi nu
intotdeauna radacinile sunt reale si distincte). De subliniat insa ca nu exista o metoda
generala sigura pentru gasirea tuturor radacinilor ecuatiilor neliniare, dar o analiza
riguroasa a problemei da o mare probabilitate de rezolvare in majoritatea cazurilor.

f(x)

-h h 2h 3h 4h 5h 6h 7h 8h 9h
X1 0 X2 X3 X4
x

Figura 3.1 – Esantionarea abscisei pentru gasirea intervalelor care contin radacini
 Programul 3.1 pentru gasirea intervalelor cu radacini prin esantionare este
scris in Mathematica si are:
• intrari: functia, intervalul si pasul de esantionare
• iesiri: graficul functiei (figura 3.2) si lista intervalelor pe care se afla radacini
Este exemplificata gasirea intervalelor care contin radainile ecuatiei transcendente
3x − 5 x 2 + 2 = 0

Programul 3.1. Gasirea intervalelor cu radacini prin esantionare

Clear@"`∗"D;
f@x_D = 3x − 5 x2 + 2.; H∗ Definirea functiei ∗L
a = −2.;
b = 5; H∗ Stabilirea intervalului global∗L
Plot@f@xD, 8x, a, b<D; H∗ Vizualizarea graficului ∗L

h = 10−2; H∗ Pasul de esantionare ∗L

lst = 8<; H∗ Initializare lista intervalelor ∗L
For@x = a − h, x ≤ b + h, x = x + h,
yn = f@x + hD;
If@yn yv ≤ 0, AppendTo@lst, 8x, x + h<DD; H∗ Adauga intervalul in lista ∗L
yv = yn; H∗ Noua ordonata devine ordonata veche pentru urmatorul interval∗L
D;
Print@lstD;
30

20

10

-2 -1 1 2 3 4 5

-10

Figura 3.2. Graficul functiei F ( x) = 3x − 5 x 2 + 2

Intervalele pentru radacini rezultate din program sunt urmatoarele:

{{-0.71,-0.7},{0.99,1.},{3.93,3.94}}

Dupa cum se observa, conditia de existenta a radacinii in interval a fost pusa

sub forma F ( xv ) F ( xn ) ≤ 0 , deoarece exista posibilitatea ca o radacina sa fie chiar in
punctul de esantionare (ca in cazul radacinii x=1 in exemplul considerat). De
asemenea, ciclul a fost inceput cu un esantion inaintea capatului inferior al
intervalului pentru eventualitatea ca o radacina sa se afle chiar in acel capat.
 3.2 Metoda bisectiei
Una dintre cele mai simple si sigure metode pentru aflarea radacinii dintr-un
subinterval o costituie metoda bisectiei. Desi are o convergenta scazuta (cum se va
arata mai departe) ea este totusi preferata in multe aplicatii datorita usurintei de
implementare si datorita faptului ca viteza nu mai este o problema pentru
calculatoarele actuale intr-o clasa larga de probleme.
Principiul metodei este urmatorul: se imparte intervalul [a,b] in doua jumatati
si se studiaza semnul functiei in cele trei puncte obtinute. Se alege intervalul in care
semnele sunt opuse si apoi procesul se reia.
Astfel, noua valoare a variabilei este:
⎧ xk − a
⎪⎪ 2 , daca F (a) F ( xk ) < 0
xk +1 = ⎨ (3.2)
⎪ b − xk , daca F (a) F ( x ) > 0
⎪⎩ 2 k

iar noul intervat este [ xk +1 , a ] daca F (a) F ( xk +1 ) < 0 si [b, xk +1 ] in caz contrar.
Rezulta ca dupa n bisectii, intervalul se micsoreaza de 2n ori, deci si eroarea se
micsoreaza in aceeasi proportie:
1 1
ε = xn − x0 = n x1 − x0 ≤ n b − a (3.3)
2 2
unde xn este valoarea obtinuta pentru radacina dupa n bisectii x0 este valoarea reala iar
x1 este valoarea initiala aleasa pentru radacina (in intervalul respectiv, de regula a sau
b).
Numarul de pasi necesari obtinerii unei anumite erori se obtine prin
logaritmarea relatiei (3.3):
73H

b − a 10 b−a
n ≥ Log 2 ≈ Log10 , n∈ (3.4)
ε 3 ε
Convergenta gasirii radacinii poate fi evaluata din relatia intre erorile obtinute
in doi pasi consecutivi. Deoarece la fiecare pas intervalul se injumatateste, se poate
scrie:
1
ε n +1 = ε n = Const × ε n (3.5)
2
O astfel de relatie caracterizata printr-o proportionalitate intre erorile
consecutive defineste o convergenta liniara a metodei. Exista, dupa cum se va arata,
si metode care au o convergenta supraliniara, definita printr-o relatie de tipul:
ε n +1 = Const × ε nα , α > 1, α ∈ (3.6)
Evident, un procedeu este cu atat mai eficient cu cat α este mai mare, in
practica existand metode cu exponent cuprins intre 1 si 2.

Programul 3.2, scris in Mathematica, pentru rezolvarea unei ecuatii neliniare

prin metoda bisectiei foloseste ecuatia (3.2), parametrii de intrare fiind definitia
74H

ecuatiei, intervalul si precizia dorita.

 Programul 3.2. Rezolvarea unei ecuatii neliniare prin metoda bisectiei

Clear@"`∗"D;
f@x_D = Sin@xD − Log@xD; H∗ Ecuatia neliniara de rezolvat ∗L
a = 0; b = 3; H∗ Intervalul pe care se afla o radacina reala ∗L
h = 10−15; H∗ Precizia dorita ∗L
WhileAAbs@b − aD > h, H∗Ciclul repetat pana la atingerea preciziei∗L
; H∗ Noua valoare a solutiei ∗L
a+ b
c=
If@f@aD f@cD > 0, a = c, b = cDE; H∗Decizia de alegere a intervalului ∗L
2.

x = a;H∗Alegerea unei valori din intervalul final ∗L

Print@x, " ", f@xDD;

3.3 Metoda punctului fix

Metoda bisectiei este foarte sigura pentru determinarea radacinilor reale, dar
avand o rata de convergenta scazuta nu este recomandabila pentru rezolvarea de
ecuatii cu un grad de complexitate crescut. Se vor prezenta in paragrafele urmatoare
metode care au o convergenta mai rapida, dar, de cele mai multe ori, mai putin sigura.
Fie ecuatia neliniara:
F ( x) = 0 (3.7)
avand o solutie reala α in intervalul [ a, b ] stabilit anterior printr-o metoda de gasire a
intervalelor cu solutii.
Daca descompunem functia initiala F ( x) in forma:
F ( x) = x − f ( x) (3.8)
ecuatia (3.7) este echivalenta cu
75H

x − f ( x) = 0 (3.9)
O solutie α pentru ecuatia (3.9) se numeste punct fix, care este si solutie a
76H

ecuatiei (3.7) deoarece conform (3.9) are proprietatea:

77H 78H

α = f (α ) (3.10)
si conform (3.8) F (α ) = α − f (α ) = 0 .
79H

Se va demonstra in continuare ca un proces iterativ de forma

xk +1 = f ( xk ) (3.11)
poate asigura convergenta catre punctul fix α , solutie a ecuatiei (3.7) daca sunt
80H

indeplinite anumite conditii.

Intr-adevar, eroarea la iteratia k + 1 poate fi exprimata in functie de eroarea la
iteratia precedenta k , tinand cont de relatiile (3.10) si (3.11):
81H 82H

ε k +1 ≡ xk +1 − α = f ( xk ) − f (α ) ≤ f '(ξ ) xk − α = f '(ξ )ε k (3.12)

Inegalitatea din ecuatia (3.12) a fost scrisa in baza teoremei mediei, ξ fiind o
83H

valoare (necunoscuta) din intervalul [ xk , α ] ⊂ [ a, b ] .

 Aceasta ecuatie ne da si conditia de convergenta. Intr-adevar, eroarea scade de
la o iteratie la urmatoarea, daca derivata functiei in punctul dat de teorema mediei este
subunitara, adica
f ( xk ) − f (α )
f '(ξ ) = <1 (3.13)
xk − α
Rezulta, tinand cont de (3.11) si (3.8), urmatoarea regula de iteratie:
84H 85H

xk +1 = xk − F [ xk ] (3.14)
Aceasta poate fi interpretata in felul urmator: valoarea necunoscutei dupa
iteratie este egala cu valoarea de la iteratia precedenta, corectata cu valoarea
functiei chiar in punctul de la iteratia precedenta.
De subliniat insa ca nu intotdeauna convergenta acestei metode este asigurata,
fiind necesara indeplinirea conditiei (3.13), care insa nu se poate verifica apriori,
86H

nefiind cunoscuta valoarea ξ . Convergenta se poate verifica prin scaderea lui ε sub o
anumita limita impusa, si obligatoriu prin verificarea solutiei in ecuatia (3.7). Metoda
87H

are totusi o importanta practica redusa, in schimb conduce catre o generalizare

teoretica deosebit de importanta, care va fi prezentata in cele ce urmeaza.

3.4 Metoda Newton-Raphson

Metoda punctului fix poate fi generalizata prin modificarea relatiei finale

(3.14) sub forma:
88H

xk +1 = xk − β F [ xk ] (3.15)
Aceasta ar insemna ca prin inmultirea corectorului F [ xk ] cu un factor
corespunzator β se poate obtine in urma iteratiei (3.15) o valoare mult mai apropiata
89H

de radacina, deci procesul poate deveni mult mai rapid convergent.

Tinand seama in (3.15) de definita (3.8) obtinem:
90H 91H

xk +1 = xk + β [ f ( xk ) − xk ] (3.16)
Factorul β poate fi determinat punand conditia ca eroarea la iteratia k + 1 sa
fie minima. Rezulta succesiv:
f ( xk ) − xk
ε k +1 = xk +1 − α = xk − α + β [ f ( xk ) − xk ] = xk − α 1 + β
xk − α
Adunand si scazand f (α ) la numaratorul fractiei din relatia precedenta si
aplicand din nou teorema mediei obtinem:
ε k +1 = ε k 1 + β [ f '(ξ k ) − 1] (3.17)
Daca am cunoaste valoarea ξ am putea alege coeficientul β astfel ca eroarea
sa se anuleze:
1
β= (3.18)
1 − f '(ξ k )
In practica insa nu cunoastem pe ξ , deci trebuie sa ii atribuim una din valorile
cunoscute. Considerand ca cea mai apropiata valoare este xk , vom scrie
f (ξ k ) = f ( xk ) iar relatia (3.16) devine:
92H

1
xk +1 = xk + [ f ( xk ) − xk ] (3.19)
1 − f '( xk )
 Daca tinem seama de relatia (3.8) si de derivata acesteia, relatia (3.19) poate
93H 94H

fi scrisa sub forma:

F ( xk )
xk +1 = xk − (3.20)
F '( xk )
care defineste metoda de recurenta Newton-Raphson.
Se poate demonstra ca prin aceasta metoda se obtine o convergenta foarte
mare, ceea ce justifica larga sa utilizare. Intr-adevar, daca se dezvolta in serie Taylor
functia F ( x) si derivata sa in jurul radacinii α se obtin relatiile:
1
F ( xk ) = F (α ) + ε k F '(α ) + ε k2 F ''(α ) + ... (3.21)
2
F '( xk ) = F '(α ) + ε k F ''(α ) + ... (3.22)
In relatia (3.21), deoarece α este radacina, rezulta F (α ) = 0 ; din relatia
95H

(3.22), daca neglijam termenii care cuprind derivata de ordinul doi si cei superiori
96H

(lucru valabil daca suntem suficient de aproape de solutia corecta), obtinem

F '( xk ) = F '(α ) (3.23)
Daca scadem α din ambii membri ai relatiei (3.20) si tinem seama de (3.21)
97H 98H

si (3.23), rezulta:
99H

F ( xk ) 1 F ''(α )
xk +1 − α = xk − α − = ε k − ε k − ε k2 (3.24)
F '( xk ) 2 F '(α )
ε k +1 εk

care este echivalenta cu:

ε k +1 = Const.ε k2 (3.25)
Dupa cum se observa din relatia (3.25), daca alegerea initiala este suficient de
100H

apropiata de valoarea corecta a radacinii, convergenta este de ordinul 2, valoare care

indica o eficienta foarte mare a procedeului. In schimb daca alegerea initiala nu este
corecta, convergenta scade foarte mult, si exista chiar o probabilitate mult mai mare
decat la alte metode ca iteratia sa fie divergenta. De aceea, metoda Newton-Raphson
reclama gasirea in prealabil a unor intervale cat mai inguste cu radacini.
Daca alegerea initiala este un numar complex, (si desigur daca programul
foloseste numere complexe), metoda Newton-Raphson poate gasi si radacinile
complexe. Acest lucru precum si convergenta mare, o recomanda ca principala
metoda in rezolvarea ecuatiilor neliniare,
Programul 3.3, scris in Mathematica, pentru rezolvarea unei ecuatii neliniare
prin metoda Newton-Raphson foloseste ecuatia (3.20), parametrii de intrare fiind
101H

definitia ecuatiei, intervalul, precizia dorita si numarul maxim de iteratii. Acesta din
urma este necesar pentru eventualitatea ca metoda nu converge, caz in care programul
nu ramane in bucla ci iese dupa terminarea numarului maxim de iteratii. Alegerea
initiala este un capat al intervalului pe care se presupune ca exista o radacina, iar
celalalt capat este folosit ca valoare initiala de comparatie.
Programul 3.3. Rezolvarea unei ecuatii neliniare prin metoda Newton-
Raphson
 Clear@"`∗"D;
f@x_D = Sin@xD − Log@xD;
a = 0; b = 3.0;
ε = 10−12; imax = 100; H∗Precizia si numarul maxim de iteratii∗L
c = a; H∗Valoare initiala pentru radacina∗L
cv = b; H∗Valoare initiala pentru comparatie∗L
f@cvD
DoAc = cv − ; H∗ Calcul valoare noua ∗L
f'@cvD
If@Abs@c − cvD < ε, Break@D D; H∗Iese din ciclu daca s−a atins precizia∗L
cv = c, 8imax<E; H∗Ciclul Do se repeta de maxim imax ori ∗L
sol = c;
Print@sol, " ", f@solDD;
Null

Aici am ales un ciclu cu “Do” pentru familiarizarea cititorului cu tipuri

altenative de cicluri.

3.5 Metoda secantei

Metoda Newton-Raphson, desi rapid convergenta, este uneori dificil de utilizat
deoarece necesita cunoasterea derivatei functiei. De aceea este foarte raspandita o
varianta a acesteia, si anume metoda secantei.
Derivata functiei din ecuatia (3.20) poate fi inlocuita prin formula:
102H

ΔF ( xn ) F ( xn ) − F ( xn −1 )
F '( xn ) ≈ = (3.26)
Δxn xn − xn −1
daca iteratia n se face intre punctele xn −1 si xn .
In metoda secantei, fiecare punct nou gasit devine inceput de interval pentru
urmatoarea iteratie, indiferent daca derivata astfel calculata intersecteaza abscisa in
interiorul sau in exteriorul intervalului curent. Deoarece este posibil ca intersectia sa
apara in exteriorul intervalului, si acest lucru sa se mentina si la iteratiile urmatoare,
metoda poate deveni, ca si cea Newton-Raphson, divergenta.
Inlocuind (3.26) in ecuatia (3.20), si daca aproximatia radacinii la iteratia n
103H 104H

este xn , rezulta pozitia noului punct xn +1 :

F ( xn )( xn − xn −1 )
xn +1 = xn − (3.27)
F ( xn ) − F ( xn −1 )
sau
x F ( xn ) − xn F ( xn −1 )
xn +1 = n −1 (3.28)
F ( xn ) − F ( xn −1 )

Relatia (3.28) defineste metoda secantei, denumire datorata faptului ca

105H

derivata este inlocuita de dreapta care uneste punctele de la capetele intervalului de la

iteratia n.
 Deoarece derivata functiei este calculata cu aproximatie, metoda are o
convergenta mai scazuta decat Newton-Raphson, dar totusi supraliniara, ordinul de
1+ 5
convergenta fiind ≈ 1.618 , adica asa numitul “golden number”.
2
Programul 8.4, scris in Mathematica, pentru rezolvarea unei ecuatii neliniare
prin metoda secantei foloseste ecuatia (3.28), parametrii de intrare fiind:
106H

• Functia f [ x _] definita in prealabil ca o functie pura;

• intervalul pe care se cauta radacini;
• precizia dorita;
• numarul maxim de iteratii, Nmax.
Acesta din urma este necesar pentru eventualitatea ca metoda nu converge, caz
in care programul nu ramane in bucla ci iese dupa terminarea numarului maxim de
iteratii.

Programul 8.4. Rezolvarea unei ecuatii neliniare prin metoda secantei

Clear@"`∗"D;
f@x_D = Sin@xD − Log@xD;
a = 1; b = 3.0; H∗ Intervalul care contine radacinile ∗L
ε = 10−7; Nmax = 20; H∗ Conditii de oprire ∗L
cv = a; H∗Valoare initiala a radacinii, pentru comparatie∗L
ForAn = 1, n < Nmax, n++,
a f@bD − b f@aD
; H∗Valoarea noua a radacinii∗L
f@bD − f@aD
c=

If@f@aD f@cD > 0, a = c, b = cD

If@Abs@c − cvD < ε, Break@D D;
cv = cE; H∗Alegerea noului interval∗L
Print@c, " ", f@cDD;

Demonstarea ordinului de convergenta al metodei secantei este destul de dificila si poate fi

omisa la o prima lectura.
Pentru demonstratie, se pleaca de la formula (3.28), scazandu-se radacina α din ambii
107H

membri si tinand cont ca F (α ) = 0 . Rezulta succesiv:

( xn −1 − α ) [ F ( xn ) − F (α ) ] − ( xn − α ) [ F ( xn −1 ) − F (α ) ]
ε n +1 ≡ xn +1 − α = =
F ( xn ) − F ( xn −1 )
( xn −1 − α ) F '( xn )( xn − α ) − ( xn − α ) F '( xn −1 )( xn −1 − α )
=
F ( xn ) − F ( xn −1 )
F '( xn ) − F '( xn −1 ) xn − xn −1
= ( xn −1 − α ) ( xn − α )
xn − xn −1 F ( xn ) − F ( xn −1 )
ε n−1 εn
F ''(ξ ) 1/ F '(η )

F ''(ξ n )
ε n +1 = ε nε n −1 (3.29)
F '(η n )
S-a aplicat teorema mediei, ξn [
si η n fiind doua valori cuprinse in intervalul xn −1 , xn . ]
Rescriind ecuatia (3.29) pentru
108H ε n + 2 , si logaritmand-o obtinem:
 ⎛ F ''(ξ n ) ⎞
log(ε n + 2 ) = log(ε n +1 ) + log(ε n ) + log ⎜ ⎟
⎝ F '(η n ) ⎠
Notand
⎛ F ''(ξ n ) ⎞
log(ε n + 2 ) = ν n + 2 , log(ε n +1 ) = ν n +1 , log(ε n ) = ν n , log ⎜ ⎟ = cn (3.30)
⎝ F '(η n ) ⎠
se obtine:
ν n + 2 −ν n +1 −ν n = cn (3.31)
Definim operatorul de deplasare care face trecerea de la eroarea in iteratia n la eroarea din
iteratia n + 1 sub forma:
Eν n = ν n +1
care introdus in ecuatia (3.31) conduce la:
109H

( E 2 − E − 1)ν n = cn
1+ 5 1− 5
sau ( E − q )( E − p )ν n = cn unde am notat p = ,q =
2 2
Notand
( E − p)ν n = un (3.32)
rezulta ( E − q )un = cn si schimband indicele ( E − q )un −1 = cn −1 . Tinand cont de actiunea
operatorului E, aceasta relatie devine
un = cn −1 + qun −1 (3.33)
Rescriind (3.33) cu schimbarea indicelui n → n − 1 obtinem
110H

un −1 = cn − 2 + qun − 2 (3.34)
care introdusa in (3.33) duce la:
111H

un = cn −1 + qcn − 2 + q 2un − 2
Repetand procedeul prin inlocuirea lui un − 2 printr-o relatie analoaga cu (3.34), obtinem in 112H

final:
un = cn −1 + qcn − 2 + q 2 cn −3 + ... + q n −1c0 + q nu0 (3.35)
Tinand cont de actiunea operatorului E, si de definitia (3.32), membrul stang al ecuatiei
113H

(3.35) devine: un = ( E − p )ν n = ν n +1 − pν n .
114H

Membrul drept contine o serie alternata cn −1 + qcn − 2 + q 2 cn −3 + ... + q n −1c0

1− 5
(deoarece q = < 0 ), avand deci o valoare finita, precum si cantitatea finita
2
q nu0 = q n ( E − p)ν 0 = q n (ν 1 − pν 0 ) . Notand membrul drept al ecuatiei (3.35) cu L < ∞ , aceasta
115H

relatie devine:
ν n +1 − pν n = L (3.36)
si tinand cont de notatia (3.30) rezulta
116H log(ε n +1 ) − log(ε np ) = L , de unde se deduce relatia care da
ordinul de convergenta al metodei secantei:
1+ 5
ε n +1 = e Lε n 2
(3.37)
 3.6 Metoda Muller de interpolare cu parabola
Dupa cum s-a aratat, metoda secantei foloseste ultimele doua puncte
cunoscute, printre care se traseaza o dreapta a carei intersectie cu abscisa determina
noul punct.
Metoda Muller este o varianta a metodei secantei, folosind ultimele trei puncte
cunoscute, prin care se traseaza o parabola, a carei intersectie cu abscisa determina
noul punct. Este de asteptat sa se obtina o convergenta mai buna decat in metoda
secantei, dar formulele care iau in consideratie trei puncte pot fi mai complicate.
Fie xk , xk −1 si xk − 2 , k ≥ 2 cele trei valori consecutive pentru radacina, obtinute
pana la iteratia precedenta k-1, sau alese in intervalul stabilit pentru radacina inaintea
primei iteratii, ca in figura 3.2. O parabola p( x) care trece prin punctele
[ xk −2 , F ( xk −2 )] , [ xk −1 , F ( xk −1 )] si [ xk , F ( xk )] va avea ecuatia:
aν 2 + bν + c = p ( x) (3.38)
F(x)

h1
ν
x k-2 x k-1 xk x
0
h2

Figura 3.3 Metoda parabolei (Muller)

Notand ν = x − xk −1 (ne raportam la punctul aflat in mijloc, prin care va trece

deci si ordonata, conform figurii), si h1 = xk − xk −1 si h2 = xk −1 − xk − 2 , pentru ca
parabola p( x) sa treaca prin cele trei puncte prin care trece si functia F ( x) a carei
radacina o cautam, trebuie satisfacut sistemul de ecuaţii:
⎧a 02 + b0 + c = F ( xk −1 )
⎪ 2
⎨ah1 + bh1 + c = F ( xk ) (3.39)
⎪ 2
⎩ah2 − bh2 + c = F ( xk − 2 )
unde s-a tinut cont că în aceste puncte variabila ν a parabolei are valorile 0, h1
respectiv − h2 , conform figurii.
Pentru aflarea parametrilor a , b si c ai parabolei, rezolvam acest sistem de
ecuatii. Din prima ecuatie rezulta:
c = F ( xk −1 ) (3.40)
Inmultind a doua ecuatie cu h2 si a treia cu h1 si adunandu-le rezulta
 F ( xk ) F ( xk − 2 ) F ( xk −1 )
a= + − (3.41)
h1 (h1 + h2 ) h2 (h1 + h2 ) h1h2

Inmultind a doua ecuatie cu h22 si a treia cu h12 si scazandu-le rezulta

h F ( xk ) h F ( xk − 2 ) (h1 − h2 ) F ( xk −1 )
b= 2 − 1 + (3.42)
h1 (h1 + h2 ) h2 (h1 + h2 ) h1h2
Parabola intersecteaza abscisa in noul punct xk +1 care rezulta din ecuatia
p(ν ) = 0 cu parametrii dati de ecuatiile (3.40), (3.41) si (3.42).
117H 118H 119H

−b ± b 2 − 4ac 2c
ν 1,2 = =
2a −b ∓ b 2 − 4ac
Tinand cont de notatia ν = x − xk −1 , noua valoare pentru radacina va fi deci:

2c
xk +1 = xk −1 − (3.43)
b ± b 2 − 4ac
Semnul de la numitor se ia in asa fel incat sa se obtina valoarea maxima a
modulului acestuia, astfel incat sa se inainteze catre cea mai apropiata radacina fata de
punctele anterioare (deci plus daca b este pozitiv si minus daca b este negativ).
Se poate demonstra ca ordinul de convergenta al metodei Muller este 1,84,
foarte apropiat de metoda Newton, fara a fi necesara insa cunoasterea derivatei. De
asemenea, este de remarcat ca functia poate fi reala sau complexa si se pot obtine
atat radacinile reale cat si cele complexe, lucru care nu era posibil cu metodele
anterioare. Alegerea initiala necesita in schimb trei puncte pe abscisa, care trebuie sa
fie suficient de apropiate de radacina.
Programul 8.5, pentru calculul radacinilor unei ecuatii prin metoda Muller
foloseste relatiile (3.40) - (3.43), si esantioneaza aleator planul complex intr-un
120H 121H

numar N st de triplete de puncte localizate intr-un patrat de latura L centrat in origine,

pentru a gasi cele mai apropiate radacini de acestea. Daca se obtine o radacina (fapt
semnalat de o valoare suficient de mica a functiei in punctul respectiv) se analizeaza
partea sa imaginara. Daca aceasta parte imaginara este suficient de mica, se considera
ca radacina este reala si se afiseaza doar partea reala a radacinii respective.
Nu se poate cunoaste ordinul de multiplicitate al radacinilor, acesta trebuind sa
fie dedus prin alte metode.
Programul are intrarile:
• Functia f [ x _] definita in prealabil ca o functie pura;
• Eroarea maxima admisa in valoarea functiei, ε ;
• Numarul de starturi din puncte initiale, N st , suficient de mare pentru a se
asigura gasirea tuturor radacinilor;
• Latura patratului in planul complex in care se cauta radacinile, L ;
• Numarul maxim de iteratii N max daca nu se obtine convergenta;
Iesirea este o lista cu radacini gasite, dintre care unele se repeta fara a se sti
daca sunt multiple sau nu. Se afiseaja de asemenea si valorile functiei in radacaini,
pentru verificare.
Se observa ca punctele initiale se aleg aleator de N st ori, in planul complex in
patratul [ −0.5 L ⋅ i, 0.5 L ⋅ i ] × [ −0.5 L, 0.5 L ] .
 Programul 8.5 Calculul radacinilor unei ecuatii prin metoda Muller
f@x_D = 10 + 7 Sin@xD − 3 x2 − 2 x3 + x4 + 2 x5;
X = 8<; H∗ Vectorul radacinilor ∗L;
Nst = 40; H∗ Numarul de repetari ale startului∗L
L = 10; H∗Latura patratului in planul complex in care se cauta radacini∗L
ForAn = 0, n < Nst, n++,
x2 = L HRandom@D − 0.5L + L I HRandom@D − 0.5L; H∗ Puncte de start aleatoare ∗L
x1 = L HRandom@D − 0.5L + L I HRandom@D − 0.5L;
x0 = L HRandom@D − 0.5L + L I HRandom@D − 0.5L;
H∗ Se esantioneaza aleator spatiul complex ∗L
ε = 10−13; H∗ Eroarea admisibila in valoarea functiei ∗L
Nmax = 100; H∗ Numarul maxim de iteratii daca nu se obtine convergenta ∗L
k = 1;
DoAh1 = x0 − x1;
H∗ Ciclul pentru cautarea unei radacini in apropierea punctelor curente ∗L
h2 = x1 − x2;
c = f@x1D;
f@x0D f@x2D f@x1D
h1 Hh1 + h2L h2 Hh1 + h2L
a= + − ;
h2 f@x0D h1 f@x2D Hh1 − h2L f@x1D
h1 h2

h1 Hh1 + h2L h2 Hh1 + h2L

b= − + ;

x2 = x1; H∗ Se reactualizeaza punctele x2 si x1∗L

h1 h2

x1 = x0;
H∗ Pentru reactualizarea lui x0 se calculeaza pasul h= ∗L
b+"################# b+"#################
2c sau h= 2c

è!!!!!!!!!!!!!!!!!
b2−4 a c b2−4 a c

è!!!!!!!!!!!!!!!!!
d1 = b + b2 − 4 a c ;
d2 = b − b2 − 4 a c ;
IfAAbs@d1D > Abs@d2D, h = E;
2c 2c
,h=
x0 = x2 − h; H∗x2 a capatat mai sus valoarea x1, fata de care am scris pe h ∗L
d1 d2

IfAAbs@f@x0DD < ε ,
IfAAbs@Im@x0DD > 10−7, AppendTo@X, x0D,
AppendTo@X, Re@x0DDE; Break@D E;
H∗ Daca x0 da o valoare suficient de mica a functiei

se adauga doar partea reala la vectorul radacinilor ∗L

se opresc iteratiile. Daca partea imaginara a radacinii este prea mica

k++, 8Nmax<E;
E;
Y = Sort@XD;
For@n = 1, n ≤ Length@YD, n++, Print@Y@@nDD, " ", f@Y@@nDDDDD;
H∗ Afiseaza radacinile; unele se repeta, fara a fi a se sti daca sunt multiple ∗L
 3.7 Metoda Lobacevski-Graeffe de calculare a radacinilor reale
ale polinoamelor

Pentru aflarea radacinilor polinoamelor exista metode specifice avand in

general o convergenta si stabilitate mai buna decat metodele generale valabile pentru
orice ecuatii transcendente. Una dintre acestea este metoda Lobacevski-Graeffe care,
prin prelucrari algebrice ridica la o putere superioara toate radacinile permitand
separarea radacinii dominante si aplicand apoi relatiile Viete.
Fie ecuatia (3.1) in care membrul stang este un polinom de grad n care poate
122H

fi scris sub doua forme:

n n
Pn ( x) = ∑ an − j x j = a0 ∏ ( x − x j ), a0 ≠ 0 (3.44)
j =0 j =1

Intre radacinile si coeficientii ecuatiei polinomiale exista relatiile Viete:

a
x1 + x2 + ... + xn = − 1
a0
a2
x1 x2 + x1 x3 + ... + xn −1 xn =
a0
............................................
(3.45)
aj
x1 x2 ...x j −1 x j + x1 x2 ...x j −1 x j +1 + ... = (−1) j

a0
............................................
an
x1 x2 ...xn = (−1) n
a0
Daca una dintre radacini, s-o numim x1 ar fi dominanta, adica:
x1 x2 , x1 x3 ,...., x1 xn (3.46)
folosind prima relatie Viete in care neglijam toate celelalte radacini, putem sa obtinem
valoarea radacinii dominante:
a
x1 = − 1 (3.47)
a0
Evident, in general nu avem aceasta situatie foarte favorabila, dar ea se poate
obtine prin prelucrari algebrice.
Sa consideram radacina cu cea mai mare valoare absoluta si sa o notam cu x1 ,
ea avand proprietatea:
x1 > x2 , x1 > x3 ,...., x1 > xn (3.48)
Daca am reusi sa ridicam toate radacinile la o putere suficient de mare, relatia (3.48) 123H

ar deveni o relatie de tipul (3.46), ceea ce ar permite calculul simplu al radacinii

124H

dominante si apoi a celorlalte.

Pentru aceasta, sa consideram mai intai polinomul in variabila − x :
n n
Pn (− x) = ∑ an − j (− x) j = (−1) n a0 ∏ ( x + x j ) (3.49)
j =0 j =1

Facand produsul (−1) n Pn ( x) Pn (− x) , obtinem un polinom in variabila

x 2 , pe care sa-l notam cu Pn(1) ( x 2 ) :
 n n
Pn(1) ( x 2 ) = (−1) n Pn ( x) Pn (− x) = a02 ∏ ( x 2 − x 2j ) = ∑ A(1)
j x
2j
(3.50)
j =1 j =0

unde am notat cu A(1)

j coeficientii noului polinom Pn(1) ( x 2 ) .
Presupunand ca x1 > x2 , cu atat mai mult x12 > x22 , ceea ce incepe sa ne
apropie de situatia unei radacini dominante ca in (3.46). 125H

Repetand procedeul obtinem un polinom in x 4 :

n n
Pn(2) ( x 4 ) = (−1) n Pn(1) ( x 2 ) Pn(1) (− x 2 ) = a04 ∏ ( x 4 − x 4j ) = ∑ A(2)
j x
4j

j =1 j =0

Dupa aplicarea procedeului de s ori, se obtine:

n n

∏ ( x 2 − x2j ) = ∑ A(j s ) x2
s −1 s −1
Pn( s ) ( x 2 ) = (−1) n Pn( s −1) ( x 2 ) Pn( s −1) (− x 2 ) = a02
s s s s s
j
(3.51)
j =1 j =0

Se observa ca puterile radacinilor cresc foarte rapid cu s, ajungandu-se la o

relatie de ordine de tipul (3.46). 126H

Notand 2s = m , avem din prima relatie Viete:

As
x1m = − 1s (3.52)
A0
Aceasta permite calculul primei radacini dupa formula:
1
s m
A
x1 = ± 1
s
(3.53)
A 0

Semnul se stabileste prin calculul valorii functiei pentru un semn si celalalt,

alegand-l pe acela care da valoarea cea mai mica in modul.
Considerand ca dupa un numar suficient de ridicari la puteri radacinile sunt
ordonate dupa o relatie de forma: x1 x2 x3 ,...., xn , putem calcula si
celelalte radacini. Intr-adevar, a doua relatie Viete devine:
A2( s )
x1 x1 = ( s )
m m
(3.54)
A0
si impartind-o la relatia Viete precedenta (3.52) obtinem a doua radacina:
127H

1
A( s ) m
x2 = ± 2( s ) (3.55)
A1
Analog se obtin si celelalte radacini
1
(s) m
A
xj = ±
j
(s)
(3.56)
A j −1

semnul alegandu-se de asemenea prin testare in ecuatie.

s
Coeficientii Asj ai polinomului final Pn( s ) ( x 2 ) , necesari in relatiile de calcul
ale radacinilor, trebuie calculati din coeficientii polinomului initial Pn ( x ) . Pentru
aceasta se poate folosi o relatie de recurenta. Din relatia (3.51) rescrisa dentru 128H

s → s + 1 si tinand cont de expresia polinoamelor care se inmultesc, rezulta

n ⎛ n ( s ) 2s j ⎞ ⎛ n ⎞
∑ A( s +1) 2 s +1 j
n− j x = ( − 1) n (s)
Pn ( x 2s
) Pn
(s)
( − x 2s
) = ⎜ ∑ An − j x ⎟ ⎜ ∑ (−1) An − j x ⎟ (3.57)
j ( s ) 2s j

j =0 ⎝ j =0 ⎠ ⎝ j =0 ⎠
 Se poate deci deduce o relatie de recurenta intre coeficientul A(j s +1) al
polinomului din etapa s + 1 si coeficientii Ak( s ) polinoamelor din etapa s prin
identificarea termenilor cu aceeasi putere a lui x din primul si ultimul membru al
acestei relatii. Pentru aceasta se tine cont de urmatoarele:
• Termenul la puterea 2s j dintr-un polinom initial poate fi inmultit doar cu
termenul de aceleasi putere 2s j din cel de al doilea polinom initial pentru a
obtine termenul cu puterea 2 s +1 j din polinomul final;
• Termenul cu puterea 2 s +1 j din polinomul final mai poate fi obtinut si din
inmultirea termenului la puterea 2s ( j − k ) dintr-un polinom initial poate fi
inmultit doar cu termenul de aceleasi putere 2s ( j + k ) , k = 1,.., j − 1 din cel de
al doilea polinom initial (cazurile k = 0 si k = j au fost deja tratate la punctul
anterior) ;
• Datorita simetriei, termenii de la punctul 2 apar de doua ori, deci suma lor
trebuie dublata;
• Termenii de forma 2 s ( j − k ) din unul din polinoame au semnul (−1) .
j

Rezulta deci relatia de recurenta pentru calculul coeficientilor A(j s +1) sub
forma:
j −1
A(j s +1) = ⎡⎣ A(j s ) ⎤⎦ + 2∑ (−1) k A(j −s )k A(j +s )k
2
(3.58)
k =1
Alternativ, aceasta relatie se poate scrie si sub o alta forma, echivalenta. Intr-
adevar, se poate scrie:
Pn( s ) ( x m ) = (−1) n Pn( s −1) ( x m ) Pn( s −1) (− x m ) (3.59)
⎡ n ( s −1) m 2l ⎤ ⎡ n ( s −1) m k2 ⎤ n ⎡ n ( s −1) ( s −1) m⎛⎜⎝ 2 + 2 ⎞⎟⎠ ⎤
n l k

∑ A x = ⎢ ∑ An −l x ⎥ ⎢ ∑ An − k x ⎥ = ∑ ⎢∑ An −l An − k x
(s)
n− j
mj
⎥ (3.60)
j =0 ⎣ l =0 ⎦ ⎣ k =0 ⎦ l = 0 ⎢⎣ k = 0 ⎥⎦
Se observa ca termenul in x mj din primul membru trebuie identificat cu
⎛l k⎞
m⎜ + ⎟
⎝2 2⎠
termenul in x din ultimul membru ceea ce conduce la conditia:
⎛l k⎞
m ⎜ + ⎟ = mj (3.61)
⎝2 2⎠
deci
l = 2j−k (3.62)
iar relatia (3.60) se poate scrie:
129H

n n
⎡ n ⎤
∑A
j =0
(s)
n− j x mj =∑ x mj ⎢∑ An( s−−l 1) An( s−−k1) ⎥
l =0 ⎣ k =0 ⎦
(3.63)

ceea ce inseamna ca
n
An( −s )j = ∑ An( −s −l 1) An( −s −k1) , j = 1, 2,..., n − 1 (3.64)
k =0
Multimea de valori a lui j s-a stabilit tinand cont ca deja se cunoaste
A0( s ) = A02 iar termenii An( s ) sunt termini liberi si deci nu pot ridica gradul monomului
s

cu care se face produsul. Acelasi lucru se poate spune si despre indicii coeficientilor
din membrul drept al realtiei (3.64), deci: 130H

n − l ≤ n −1 (3.65)
 Tinand cont de (3.62) rezulta:
131H

n − (2 j − k ) ≤ n − 1 ⇒ k ≤ 2 j − 1 (3.66)
De asemena pentru indicile celui de al doilea coeficient se poate scrie:
n − k ≤ n −1 ⇒ k ≥ 1 (3.67)
Inlocuind pe l dat de (3.62) si tinand cont de limitele de variatie a lui k dat de
132H

relatiile (3.66) si (3.67), ecuatia (3.64) devine

133H 134H 135H

2 j −1
An( s−)j = ∑A
k =1
( s −1)
n − (2 j − k ) An( −s −k1) , j = 1, 2,..., n − 1 (3.68)

La utilizarea oricareia dintre aceste relatii de recurenta se va tine seama ca

A0(1) = a02
evidenta conform relatiei (3.50).136H

Pentru verificarea conditiei de oprire a iteratiilor, s-ar putea compara valoarea

functiei in x j sau in − x j , valoare care in unul dintre aceste cazuri (cel cu semnul care
va fi ales in final) trebuie sa difere de zero cu o cantitate suficient de mica ε , impusa.
Caracteristica cea mai importanta a metodei este ca nu sunt necesare puncte de
start, acesta fiind un caz foarte rar in metodele de rezolvare a ecuatiilor neliniare.
Totusi, deoarece prin ridicari succesive la patrat se pot obtine destul de repede valori
foarte mari pentru coeficientii polinoamelor obtinute, metoda trebuie folosita doar cu
un numar mic de iteratii, altfel este posibil ca ca apara depasire de registru (numerele
sa fie prea mari pentru mediul de programare folosit). De aceea, de cele mai multe ori
se impune ca intrare numarul maxim de iteratii si nu eroarea admisibila. Deoarece
eroarea la un numar de iteratii mic poate fi totusi destul de importanta, metoda este
recomandabila in special pentru gasirea intervalelor care contin radacini (braketing)
fiind apoi urmata de o alta metoda pentru rafinarea acestora.
Programul Mathematica permite valori destul de mari ale coeficientilor,
putand fi folosit pentru 20-30 de iteratii in majoritatea cazurilor, ceea ce asigura si
rafinarea radacinilor, dar metoda se recomanda totusi pentru polinoame de ordin
moderat.
Programul 8.6, pentru calculul radacinilor reale ale unui polinom prin metoda
Lobacevski-Graeffe foloseste relatiile Error! Reference source not found. si (3.56) 137H

Programul are intrarile:

• Polinomul f [ x _] definit ca o functie pura;
• Numarul maxim de iteratii N max =25;

Programul 8.6. Calculul radacinilor reale ale unui polinom prin metoda
Lobacevski-Graeffe
 f@x_D = −13824 + 27648 x + 11808 x2 − 44560 x3 − 3392 x4 + 30632 x5 + 1746 x6 −
10869 x7 − 1215 x8 + 1911 x9 + 297 x10 − 159 x11 − 29 x12 + 5 x13 + x14;
H∗ Polinomul ale carui radacini se cauta ∗L
a = N@CoefficientList@f@xD, xDD; H∗ Lista coeficientilor ∗L
n = Length@aD; H∗ Gradul polinomului +1 ∗L
A = Flatten@Append@a, Table@0, 8n<DDD;
H∗ Se completeaza lista coeficientilor cu n zerouri ∗L
B = A; H∗ Lista auxiliara ∗L
X = Table@0, 8n − 1<D; H∗ Lista radacinilor ∗L
ForAm = 1, m < 25, m ++, H∗ Iteratiile ∗L
ForAj = 1, j ≤ n, j++,

B@@jDD = H−1Lj ‚ H−1Lk+1 A@@2 j − kDD A@@kDD;E; H∗ Coeficientii auxiliari ∗L

2 j−1

A = B; H∗ devin coeficientii polinomului in urmatoarea iteratia ∗L

k=1

E;
ForAj = 1, j <= n − 1, j++,
B@@jDD y
r= i
j
jAbsA Ez
z 2m−1 ; H∗ calculul radacinilor ∗L
1

k B@@j + 1DD {
If@Abs@f@rDD <= Abs@f@−rDD, X@@jDD = r, X@@jDD = −rD; H∗ Alegerea semnului ∗L
E;
Print@Sort@XDD; H∗ Ordonarea radacinilor si afisarea ∗L

Este interesant ca se poate deduce prin aceasta metoda si existenta radacinilor

multiple si a radacinilor complexe (acestea din urma neputand insa fi calculate cu
exactitate).
In cazul radacinilor de multiplicitate M, prima relatie Viete (3.52) presupune
138H

adunarea a M radacini, deci se obtine:

As
Mx1m = 1s (3.69)
A0
iar radacina se calculeaza cu prima relatie Vieta,dupa formula:
1
s m
A
x1 = ± 1
(3.70)
MA0s
sau cu relatia Vieta M, dupa formula
1
As mM
x1 = ± Ms (3.71)
A0
In cazul radacinilor complexe, acestea fiind conjugate, se poate scrie:
x1 = ρ eiϕ , x2 = ρ e −iϕ
Astfel, dupa s ridicari la patrat, prima relatie Vieta devine:
A( s )
ρ m eimϕ + ρ m e− imϕ + x3m + ... + xnm = − 1( s )
A0
deci
 A1( s ) A1( s )
ρ m eimϕ + ρ m e −imϕ ≈ − , 2 ρ m
cos mϕ ≈ − (3.72)
A0( s ) A0( s )
2
⎡⎣ A(j s −1) ⎤⎦
Astfel, daca se calculeaza rapoartele de tipul se constata
A(j s )
urmatoarele:
• Pentru radacini simple aceste rapoarte tind spre 1. Intr-adevar, ridicand la
patrat ecuatia (3.52) scrisa pentru s − 1 se obtine:
139H

2
⎡ A(j s −1) ⎤
( x ) = ( x ) = ⎢− A( s−1) ⎥ s −1 2
j
m/2 2
j
⎣⎢ j ⎦⎥
(s)
A
Pe de alta parte aceasta este egala cu x sj = (js ) deci:
Aj −1
2
⎡ A(j s −1) ⎤ A(j s )
⎢ − ( s −1) ⎥ = s
⎣⎢ Aj −1 ⎦⎥ Aj −1
2 2
⎡⎣ A(j s −1) ⎤⎦ ⎡⎣ A(j s−−11) ⎤⎦
Deci: = care poate fi scrisa si pentru j − 1, j − 2,..., 0
A(j s ) Asj −1
2 2 2 2 2
⎡⎣ A(j s −1) ⎤⎦ ⎡⎣ A(j −s −11) ⎤⎦ ⎡⎣ A(j −s −21) ⎤⎦ ⎡⎣ A1( s −1) ⎤⎦ ⎡⎣ A0( s −1) ⎤⎦
= = = ... = = =1 (3.73)
A(j s ) Asj −1 A(j −s )2 A1s A0s

2
unde s-a tinut cont ca A0( s ) = ⎡⎣ A0( s −1) ⎤⎦
Rezulta :
2
⎡⎣ A(j s −1) ⎤⎦
rj = ≈1.
A(j s )
Deci, daca prin calculul acestui raport pentru radacina j se obtine o valoare
apropiata de 1, radacina este simpla.
• Pentru radacini de ordin de multiplicitate M aceste rapoarte tind spre M.
2
⎡⎣ A(j s −1) ⎤⎦
rj = ≈M
A(j s )
Demonstratia este identica, dar in locul relatiei (3.52), se pleaca de la relatia 140H

(3.69). Deci, daca prin calculul acestui raport pentru radacina j se obtine o
141H

valoare apropiata de M, radacina este de ordinul de multiplicitate M.

• Pentru radacini complex conjugate, acest raport este oscilant. Intr-adevar,
tinand cont de relatia (3.72), putem scrie: 142H

A( s )
x1m + x2m = 2 ρ m cos mϕ = 1( s ) ; ⇒ A1( s ) = 2 A0( s ) ρ m cos mϕ (3.74)
A0

(x )
2
m m 2
m ⎡ A( s −1) ⎤ 2 2 m
1
2
+ x2 2
= 4 ρ cos ϕ = ⎢ 1( s −1) ⎥ ; ⇒ ⎡⎣ A1( s −1) ⎤⎦ = 4 ⎡⎣ A0( s −1) ⎤⎦ ρ m cos 2 ϕ
m 2

2 ⎣ A0 ⎦ 2
 Deci, acest raport este:
m
⎡⎣ A(j s −1) ⎤⎦
2
2 cos 2 ϕ
rj = = 2 .
A(j s ) cos mϕ
Deci, daca prin calculul acestui raport pentru radacina j se obtine o valoare
oscilanta in functie de s, exista doua radacini compex conjugate.
In literatura sunt mentionate numeroase alte reguli de aplicare a metodei
Lobacevski-Graeffe pentru radacini multiple si complex conjugate. Aplicarea acestora
in algoritm poate imbunatati substantial performantele metodei.
Importanta metodei este data in principal de posibilitatea de determinare a
ordinului de multiplicitate a radacinilor. Cum s-a specificat anterior, dezavantajul
principal consta in posibilitatea aparitiei unor numere foarte mari prin ridicari
succesive la patrat, conducand la depasire de registru, ceea ce limiteaza precizia.

3.8 Metoda Bairstow de calculare a radacinilor complexe ale

polinoamelor
In cazul in care polinomul are si radacini complexe, metodele elementare
prezentate anterior (cu exceptia metodei Muller) nu pot calcula radacinile respective
datorita faptului ca graficul functiei nu intersecteaza abscisa la aceste valori ale
variabilei. Una dintre metodele cele mai eficiente de calculare a unor astfel de
radacini este metoda Bairstow, care consta in factorizarea polinomului cu un trinom
de gradul doi.
Fie un polinom de grad n, in variabila reala sau complexa z:
Pn ( z ) = an z n + an −1 z n −1 + an − 2 z n − 2 + .... + a2 z 2 + a2 z + a0 (3.75)
Se incearca impartirea cu un trinom de gradul doi de forma z 2 − pz − q , care
are radacinile:
p+ p 2 + 4q
z1 =
2 (3.76)
p−
p + 4q 2
z2 =
2
care pot fi reale sau complex conjugate dupa cum p 2 + 4q este pozitiv, respectiv
negativ.
Problema se reduce deci intr-o prima instanta la determinarea coeficientilor p
si q care sa realizeze o factorizare exacta a polinomului dat, dupa care catul obtinut
este din nou supus procedeului Bairstow s.a.m.d. pana la obtinerea tuturor radacinilor.
Dupa cum se va vedea, metoda este iterativa si are unele din caracteristicile
metodei Newton-Raphson, deci este necesara o aproximare initiala suficient de buna a
coeficientilor p0 si q0 , ceea ce asigura o convergenta rapida (de ordinul doi). In
schimb, daca alegerea initiala nu este suficient de apropiata de cea corecta,
convergenta este slaba sau poate sa nu existe deloc.
 Deoarece alegerea initiala este diferita intr-o oarecare masura de solutia
exacta, la impartire apare un rest sub forma unui polinom de gradul unu. O varianta de
scriere a factorizarii este urmatoarea 1: 0F

Pn ( z ) = ( z 2 − p0 z − q0 )Qn − 2 ( z ) + r1 ( p0 − z ) + r0 (3.77)
unde Qn − 2 ( z ) este un polinom de grad n − 2 reprezentand catul impartirii
Qn − 2 ( z ) = rn z n − 2 + rn −1 z n −3 + rn − 2 z n − 4 + ... + rk z k − 2 + ...r4 z 2 + r3 z + r2 (3.78)
iar r1 ( p0 − z ) + r0 este un polinom de gadul intai reprezentand restul impartirii. Daca
p0 si q0 ar avea valorile corecte, impartirea s-ar face exact, restul fiind nul. Prin
urmare, problema se reduce la minimizarea (anularea) coeficientilor restului, r1 si r0 .
Deoarece acesti coeficienti depind de valorile coeficientilor p si q
r0 = r0 ( p, q )
r1 = r1 ( p, q )
ei se pot dezvolta in serie Taylor in jurul unui punct ales initial ( p0 , q0 ) , urmarindu-se
ca in punctul final sa aiba valori cat mai mici, apropiate de zero. Daca alegerea initiala
este suficient de corecta, termenii de ordin superior din serie se pot neglija si se poate
scrie:
∂r ∂r
r0 ( p, q) = r0 ( p0 , q0 ) + 0 ( p − p0 ) + 0 (q − q0 ) = 0
∂p ∂q
(3.79)
∂r1 ∂r1
r1 ( p, q ) = r1 ( p0 , q0 ) + ( p − p0 ) + (q − q0 ) = 0
∂p ∂q
Prin urmare, se pot calcula valorile corecte p si q din sistemul de ecuatii
(3.79) in functie de valorile initiale p0 si q0 si de elementele matricii Jacobiene:
143H

⎡ ∂r0 ∂r0 ⎤
⎢ ∂p ∂q ⎥
J =⎢ ⎥ (3.80)
⎢ ∂r1 ∂r1 ⎥
⎢ ∂p ∂q ⎥
⎣ ⎦
calculate in punctele initiale p0 si q0 .
Intr-adevar, ecuatia (3.79) se poate scrie matricial sub forma:
144H

⎡ ∂r0 ∂r0 ⎤
⎢ ∂p ∂q ⎥ ⎡ p − p ⎤ ⎡ −r ( p , q ) ⎤
⎢ ⎥⋅⎢ 0
⎥ =⎢ 0 0 0 ⎥ (3.81)
⎢ ∂r1 ∂r1 ⎥ ⎣ − 0 ⎦ ⎣ −r1 ( p0 , q0 ) ⎦
q q
⎢ ∂p ∂q ⎥
⎣ ⎦
de unde rezulta noile valori ale lui p si q :
⎡ p ⎤ ⎡ p0 ⎤ −1 ⎡
−r0 ( p0 , q0 ) ⎤
⎢ q ⎥ = ⎢ q ⎥ + J ⋅ ⎢ −r ( p , q ) ⎥ (3.82)
⎣ ⎦ ⎣ 0⎦ ⎣ 1 0 0 ⎦
Deoarece in general nu se pot obtine dintr-un singur pas chiar valorile corecte
ale lui p si q , relatia (3.79) fiind doar aproximativa, procedeul se reia si se obtin prin
145H

1
Uneori se foloseste forma mai simpla si mai naturala r1 z + r0 pentru restul impartirii, dar aceasta
conduce la calcule ceva mai complicate si la relatii de recurenta mai complexe. In unele cazuri ea este
preferabila, dar pentru simplitatea programului am ales forma prezentata.
iteratii valori din ce in ce mai apropiate de cele corecte. Relatia (3.82) este deci 146H

calculata in mod repetat, sub forma:

⎡ p j +1 ⎤ ⎡ p j ⎤ −1 ⎡
− r0 ( p j , q j ) ⎤
⎢q ⎥ ⎢q ⎥ = + J j ⋅ ⎢ −r ( p , q ) ⎥ (3.83)
⎣ j +1 ⎦ ⎣ j ⎦ ⎣ 1 j j ⎦
astfel ca, daca procedeul converge, la iteratia j + 1 se obtin valori mai corecte decat
cele de la iteratia anterioara j .
Dupa cum se observa, metoda este o generalizare in doua dimensiuni a
metodei Newton-Raphson si deci prezinta avantajele si dezavantajele acesteia, avand
si o importanta similara.
Totusi, unul dintre dezavantajele metodei Newton Raphson, si anume
necesitatea cunoasterii derivatelor functiei poate fi inlaturat prin metoda Bairstow-
Linn. Se va arata in continuare ca, prin doua relatii de recurenta, se pot calcula aceste
derivate precum si coeficientii polinomului cat Q( z ) chiar din coeficientii
polinomului initial P( z ) .
Pentru determinarea polinomului cat, se introduce relatia (3.78) in (3.77).
147H 148H

Rezulta:
Pn ( z ) = z n rn + z n −1 (rn −1 − prn ) + z n − 2 (rn − 2 − prn −1 − qrn ) + ...
(3.85)
+ z k (rk − prk +1 − qrk + 2 ) + ... + z (r1 − pr2 − qr3 ) + (r0 − pr1 − qr2 )
Prin identificare relatiei (3.85) cu (3.75) se obtin relatiile de recurenta prin
149H 150H

care se pot calcula coeficientii polinomului cat:

rn = an
rn −1 = an −1 + prn (3.86)
rk = ak + prk +1 + qrk + 2 , k = 0...n − 2
Pentru determinarea derivatelor partiale din componenta Jacobianului, se vor
face urmatoarele notatii:
∂r
ck +1 ≡ k
∂p
(3.87)
∂rk
dk +2 ≡
∂q
Acum, prin derivarea in raport cu p a relatiilor (3.86) si tinand seama de
151H

definitiile (3.87), obtinem succesiv:

152H

∂rn ∂an
= =0
∂p ∂p
∂r ∂r
cn ≡ n −1 = p n + rn (3.88)
∂p ∂p
∂r ∂r ∂r
ck +1 ≡ k = rk +1 + p k +1 + q k + 2 , k = 0...n − 2 (3.89)
∂p ∂p ∂p
ck + 2 ck +3

Din relatia (3.88) si din (3.89), renotand k + 1 → k rezulta relatiile de

153H 154H

recurenta din care se pot determina elementele Jacobianului:

cn = rn
(3.90)
ck = rk + pck +1 + qck + 2 , k = 1...n − 1
 Se observa deplasarea cu o unitate a limitelor indicelui k , datorita definitiei
∂r ∂r
(3.87). Pentru Jacobian este nevoie de valorile c1 ≡ 0 si c2 ≡ 1 care sunt ultima si
∂p ∂p
155H

repectiv penultima derivata calculata prin recurenta (3.90). 156H

Pentru determinarea celorlalte doua elemente, nu mai sunt necesare calcule

suplimentare, deoarece se poate demonstra ca, tinand cont de definitiile (3.87), avem 157H

relatiile:
∂r
d 2 ≡ 0 = c2
∂q
(3.91)
∂r1
d3 ≡ = c3
∂q
Intr-adevar, daca derivam relatiile (3.86) in raport cu q si tinem seama de
158H

definitiile (3.87), obtinem succesiv:

159H

∂rn ∂an
= =0
∂q ∂q
∂rn −1 ∂an −1
= =0
∂q ∂q

∂rn − 2 ∂r ∂r
dn ≡ = p n −1 + q n + rn (3.92)
∂q ∂q ∂q
∂r ∂r ∂r
dk +2 ≡ k = p k +1 + rk + 2 + q k + 2 , k = 0...n − 2 (3.93)
∂p ∂q ∂q
d k +3 d k +4

Din (3.92) si (3.93), renotand k + 2 → k se obtin relatiile:

160H 161H

d n = an
(3.94)
d k = rk + pd k +1 + qd k + 2 , k = 2...n − 2
Se observa ca relatiile (3.94) si (3.90) sunt similare, astfel ca ck si d k avand
162H 163H

acelasi punct de plecare ( an ) si aceleasi relatii de recurenta, vor avea valori egale,
pentru k = 2...n , ceea ce demonstreaza relatiile (3.91). 164H

Jacobianul va avea deci forma:

⎡c c ⎤
J = ⎢ 1 2⎥ (3.95)
⎣c2 c3 ⎦
Din relatia (3.83), tinand seama de regula lui Cramer, solutia sistemului va
165H

putea fi scrisa sub forma:

− r0 ( p j , q j ) c2
−r1 ( p j , q j ) c3
p j +1 = p j + (3.96)
c1 c2
c2 c3
c1 − r0 ( p j , q j )
c2 − r1 ( p j , q j )
q j +1 = q j + (3.97)
c1 c2
c2 c3
 Iteratiile pot fi oprite atunci cand atat variatiile lui p cat si ale lui q in relatiile
(3.96) si (3.97) scad sub o limita de precizie impusa.
166H 167H

Cunoscand valorile corecte ale lui p si q , din relatia (3.76) se determina

168H

radacinile cele mai apropiate de cele determinate de valorile alese initial p0 si q0 , iar
din coeficientii rk se poate determina polinomul cat Qn − 2 ( z ) . Acestuia i se va aplica
din nou procedeul Bairstow daca are grad mai mare de 2, sau se vor calcula prin
formulele analitice radacinile daca gradul sau este 2 sau 1. De remarcat insa ca la
fiecare deflatie a polinomului, deoarece restul impartirii nu este totusi exact zero,
informatia despre radacini din polinomul initial este intr-o oarecare masura alterata in
polinomul cat Qn − 2 ( z ) , astfel ca ultimele radacini gasite pot fi destul de departate de
cele reale. Fenomenul este mai frecvent in cazul radacinilor cu ordin de multiplicitate
peste 3. Aceasta face ca deflatia sa nu fie recomandabila decat pentru un numar redus
de radacini si cu ordin de multiplicitate mic.
O varianta mai corecta de determinare a tuturor radacinilor ar fi deteminarea
intr-o prima faza a intervalelor pe care se afla radacini (prin metoda Lobacevski-
Graffe de exemplu) si apoi aplicarea cate o singura data pe fiecare interval a metodei
Bairstow. In felul acesta erorile aparute la determinare catului Qn − 2 ( z ) nu se propaga
prin folosirea repetata a metodei Bairstow, in schimb este necesara o procedura in doi
pasi, care complica programul.
Programul 8.6 pentru calculul radacinilor reale ale unui polinom prin metoda
Bairstow foloseste relatiile (3.86), (3.90), (3.95)- (3.97) si (3.76). Dupa calculul
169H 170H 171H 172H 173H

radacinilor complex conjugate cele mai apropiate de alegerea initiala a trinomului, se

reia procedeul folosindu-se polinomul cat si asa mai departe pana la obtinerea unui
polinom cat de grad cel mult 2. Se calculeaza apoi radacinile finale (doua sau una,
dupa cum gradul acestui polinom este 2 sau 1).
Programul are intrarile:
• Polinomul ale carui radacini se cauta;
• Eroarea admisa ε ;
• Numarul maxim de iteratii admise, max;
• Coeficientii alesi pentru trinomul cu care se face prima factorizare, p0 si q0.
Rezulta radacinile reale si complexe, valorile functiei in aceste puncte si
numarul de iteratii folosite.

Programul 8.7. Calculul radacinilor unui polinom prin metoda Bairstow

Clear@"`∗"D;
f@x_D = x22 + 12 x9 − 7 x7 − 2 x2 + 2 x + 3; H∗ Polinomul ale carui radacini le cautam∗L
Plot@f@xD, 8x, −1.5, 2<, AxesOrigin → 80, 0<, PlotRange → 8−10, 10<D; H∗ Trasarea graficului ∗L
a = N@CoefficientList@f@xD, xDD; H∗ Lista coeficientilor, cu evitarea calculelor cu intregi ∗L
Print@aD; n = Length@aD − 1;
p0 = 1. + 0.2 I; q0 = −.8 − 0.9 I; H∗ Alegerea initiala a sumei si produsului a doua radacini ∗L
εp = εq = 1;
c = r = Table@0, 8i, n + 1<D; H∗Initializarea matricilor coloana care vor fi calculate ∗L
max = 200; H∗ Numarul maxim admis de iteratii ∗L
ε = 10−8; H∗ Eroarea admisa ∗L
j = 0;
X = 8<; H∗ Initializarea matricii radacinilor ∗L
WhileAn > 3, H∗ Se va reduce gradul polinomului pana cand vor fi 3 sau 2 termeni ∗L
WhileA j < max,
n = Length@aD − 1;
rP n+1T = aP n+1T;
rP nT = aP nT + p0 rP n+1T;
For@ k = n − 1, 1 ≤ k, k−−,
rP kT = aP kT + p0 rP k+1T + q0 rP k+2T; D; H∗ Calculul coeficientilor polinomului cat ∗L
cP n+1T = rP n+1T;
cP nT = rP nT + p0 cP n+1T;
For@ k = n − 1, 2 ≤ k, k−−,
cP kT = rP kT + p0 cP k+1T + q0 cP k+2T; D; H∗ Calculul derivatelor care vor da matricea Jacobiana ∗L
d0 = DetAj i
j
cP2T cP3T z y
zE; d1 = DetAj i
j
−rP1T cP3T y
z i
zE; d2 = DetAj
j
cP2T −rP1T zy
zE;
k cP 3T cP 4T { k −rP2T cP 4T { k cP3T −rP 2T {
; H∗ Corectorii la valorile anterioare ale sumei si produsului radacinilor ∗L
d1 d2
εp = ; εq =
p1 = p0 + εp; q1 = q0 + εq; H∗ Noile valori ale sumei si produsului radacinilor ∗L
d0 d0

j = j + 1;
If@Abs@εpD < ε && Abs@εqD < ε, Break@DD; H∗ Oprirea daca ambii corectori sunt suficient de mici ∗L
p0 = p1; H∗ Noile valori ale sumei si produsului radacinilor devin valori anterioare∗L
q0 = q1; E;
è!!!!!!!!!!!!!!!!!! è!!!!!!!!!!!!!!!!!!
; H∗ Cele 2 radacini ∗L
p1 + p12 + 4 q1 p1 − p12 + 4 q1
x1 = ; x2 =

IfAAbs@Im@x1DD < 10 , x1 = Re@x1DE; IfAAbs@Im@x2DD < 10−5, x2 = Re@x2DE;

2 2
−5

H∗ S−au eliminat partile imaginare daca sunt prea mici∗L

AppendTo@X, 8x1, f@x1D<D; H∗ Daca s−a atins precizia, caluleaza radainile∗L
AppendTo@X, 8x2, f@x2D<D; H∗si adauga in lista impreuna cu valorile functiei∗L
p1 = r@@nDD; q1 = r@@n − 1DD;
n−−; H∗ Se scade ordinul polinomului ∗L
a = Table@0, 8i, n<D;
For@i = n + 1, i ≥ 2, i−−, a@@i − 1DD = r@@i + 1DD D; H∗ Se stabilesc noii coeficienti ai polinomului ∗L
E;
è!!!!!!!!!!!!!!!!!! è!!!!!!!!!!!!!!!!!!
IfAn > 2, p1 = a@@2DD; q1 = a@@1DD; x1 =
− p1 + p12 − 4 q1 −p1 − p12 − 4 q1
; x2 = ;

IfAAbs@Im@x1DD < 10 , x1 = Re@x1DE; IfAAbs@Im@x2DD < 10−5, x2 = Re@x2DE;

2 2
−5

H∗ Daca polinomul cat are gradul 2 se calculeaza cele doua radacini ∗L

AppendTo@X, 8x1, f@x1D<D; AppendTo@ X, 8x2, f@x2D<D
, AppendTo@X, 8− p1, f@−p1D<DE; H∗ Daca polinomul cat are grad 1, s−a calculat radacina ramasa ∗L
Print@TableForm@Sort@XDDD;
H∗ Se afiseaza radacinile si valorile functiei in punctele respective pentru verificare∗L
Print@jD; H∗ Se afisaza numarul maxim de iteratii ∗L
 3.9 Metoda Laguerre de calculare a radacinilor complexe ale
polinoamelor

Una dintre metodele cele mai raspandite si sigure pentru calculul radacinilor
reale sau complexe ale polinoamelor cu coeficienti reali sau complecsi este metoda
Laguerre, care foloseste un singur punct de plecare in planul complex si poate gasi
radacina cea mai apropiata de acesta . Prin alegerea unui numar sufficient de mare de
puncte de start se pot gasi toate radacinile reale sau complexe ale ecuatiei
polinomiale.
Pentru stabilirea relatiilor de aplicare a acestei metode sa pornim de la
expresia scrisa sub forma de produs a polinomului:
n
Pn ( x) = ∏ ( x − xk ) (3.98)
k =1
in care s-a impartit cu coeficientul nenul al termenului cu cea mai mare putere, ceea
ce nu modifica ecuatia Pn ( x) = 0 .
Logaritmul natural al polinomului este:
n
ln Pn ( x) = ∑ ln( x − xk ) (3.99)
k =1
Derivata intai a acestui logaritm o notam cu G si este, conform relatiei precedente
d 1 1 1
ln Pn ( x) = + + ... + ≡G (3.100)
dx x − x1 x − x2 x − xn
iar derivate a doua a logaritmului o notam cu H si este
d2 1 1 1
ln Pn ( x) = + + ... + ≡H (3.101)
( x − x1 ) ( x − x2 ) ( x − xn )
2 2 2 2
dx
Tinand cont de membrul stang al relatiilor (3.100) si (3.101) putem scrie pe
174H 175H

G si pe H sub formele:
P' ( x)
G= n (3.102)
Pn ( x)
P' ( x) Pn' ( x) − Pn ( x) Pn'' ( x)
H= n (3.103)
Pn2 ( x)
Se pot face acum unele aproximari, care chiar daca par fortate intr-o prima
instanta conduc spre un punct mai apropiat de solutie decat punctul ales initiala. Prin
urmare, daca se repeta iterativ procesul, se ajunge din ce in ce mai aproape de solutia
corecta.
Vom presupune ca punctul ales initial se afla la distanta a fata de una din
radacini, sa spunem x1 , si la distante egale b fata de toate celelalte radacini, deci
x − x1 = a, x − xi = b, i = 2,3,..., n (3.104)
Evident aceasta ultima aproximatie este foarte grosiera, ea ar fi adevarata doar
in cazuri exeptionale (celelalte radacini sunt multiple sau plasate pe un cerc in planul
complex), dar va conduce la o aproximare mai buna decat alegerea initiala care nu se
bazeaza pe nici o presupunere.
Introducand relatiile (3.104) in (3.100) si (3.101) rezulta sistemul de ecuatii
176H 177H 178H
 ⎧ 1 n −1
⎪⎪G = a + b
⎨ (3.105)
⎪H = + 1 n − 1
⎪⎩ a2 b2
Rezolvand acest sistem in raport cu a si b , rezulta corectorul
n
a= (3.106)
G ± (n − 1)(nH − G 2 )
Evident, acesta nu este exact, dar adunat cu vechiul punct ne da o valoare mai
apropiata de radacina decat acesta. Procedeul se repeta iterativ, deci daca la iteratia k ,
aveam valoarea xk aproximand radacini x1 , la iteratia urmatoare avem o valoare mai
apropiata, corectata cu valoarea ak . Considerand xk +1 ≈ x1 , conform relatiei (3.104)
179H

avem
xk +1 = xk − ak (3.107)
Semnul din ecuatia (3.106) se alege astfel incat numitorul sa aiba valoarea cea
180H

mai mare pentru micsorarea erorilor de rotunjire si pentru a nu se exagera corectorul,

ceea ce ar putea conduce in afara intervalului de convergenta.
Se constata experimental ca in majoritatea situatilor se obtine o convergenta
foarte rapida, fiind necesare putine iteratii, ceea ce recomanda metoda ca una dintre
cele mai bune pentru ecuatii polinomiale.
Programul 8.8 de calcul pentru radacinile unei ecuatii polinomiale prin metoda
Laguerre foloseste relatiile (3.102), (3.103), (3.106) si (3.107) si are urmatoarele
181H 182H 183H 184H

intrari:
• Polinomul ale carui radacini se cauta;
• Valoarea maxima a functiei intr-o radacina, ε ;
• Numarul maxim de puncte de start, Nr ;
• Numarul maxim de iteratii admise, Ni ;
Rezulta radacinile reale si complexe, valorea maxima a functiei in aceste
puncte, numarul de puncte de start si numarul total de iteratii folosite.
In acest program s-au exemplificat in plus unele metode care pot fi folosite si
la programele prezentate anterior pentru a asigura gasirea tuturor radacinilor:
• Punctul de plecare se alege aleator intr-un patrat in planul complex cu centrul
in originea axelor si de latura 10. S-a folosit instructiunea Random[ ] care
genereaza un numar aleator intre 0 si 1;
• Dupa o oprire a iteratiilor, se testeaza daca radacina respectiva a mai fost
gasita anterior, cu negatul functiei MemberQ[Xa,x2] care da valoarea logica
True in cazul in care x 2 ∈ X . Deoarece pentru diverse puncte de plecare se
pot obtine valori usor diferite pentru aceeasi radacina (de exemplu 1.1112223
si 1.1112224), se formeaza variabila auxiliara x 2 care trunchiaza radacina x
luandu-se numai trei cifre exacte, si cu ea se formeaza multimea radacinilor
trunchiate Xa . Daca se constata ca radacina trunchiata nu a mai fost gasita, se
introduce in X radacina netrunchiata, taindu-se eventual partea imaginara daca
aceasta este prea mica (insemnand ca radacina este de fapt reala).
• Inainte de afisare, radacinile au fost ordonate crescator folosind functia Sort.
• Intregul proces se reia cu un nou punct de plecare aleator, pana cand numarul
radacinilor gasite este egal cu gradul polinomului.
Programul 8.8. Calculul radacinilor unui polinom prin metoda Laguerre
f@x_D = 4 x10 − 5 x8 − 4 x7 − H6 − 2 IL x5 + 3 x2 + H2 + 3 IL x + 1;
H∗ Polinomul ale carui radacini se cauta ∗L
ε = 10−9; H∗ Eroarea maxima admisa in functie ∗L
Nr = 200; H∗ Nr maxim de puncte de start pentru iteratii ∗L
Ni = 200; H∗ Nr maxim de iteratii la fiecare punct ∗L
A = N@CoefficientList@f@xD, xDD; H∗ Coeficientii polinomului initial ∗L
X = 8<; H∗ Multimea radacinilor ∗L
Xa = 8<; H∗ Multimea radacinilor rotunjite cu 2 cifre ∗L
n = Length@AD − 1; H∗ Gradul polinomului ∗L
j = 0; H∗ Nr de puncte de start ∗L
i = 0; H∗ Nr total de iteratii ∗L
WhileALength@ XD < n, H∗Repeta pana cand nr radacinilor este egal cu gradul polinomului∗L
j++;
If@j > Nr, Print@"Peste ", Nr, " de incercari"D; Break@DD;
H∗Daca sunt prea multe incercari iese∗L
x1 = 10 HH Random@D − 0.5L + I H Random@D − 0.5LL; H∗ Puncte de start complexe aleatoare ∗L
ForAk = 1, k ≤ Ni, k++, H∗ Iteratiile pentru fiecare punct de start ∗L
i++;
If@Abs@f@x1DD < ε, Break@DD; H∗ Se intrerup daca functia e suficient de mica ∗L

fd@x_D = ‚ A@@k + 1DD k xk−1; H∗ Derivata intai a polinomului ∗L

n

k=1

fdd@x_D = ‚ A@@k + 2DD k Hk + 1L xk−1; H∗Derivata a doua a polinomului ∗L

n−1

fd@x1D
k=1

f@x1D
G= ;

i fd@x1D z
ij y2 fdd@x1D z y
H= jjj z − z;
kk f@x1 D { f@x1D {

è!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
n
G − Hn − 1L Hn H − G2L
a1 = ;

è!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
; H∗ cele 2 posibilitati de alegere a corectorului ∗L
n
G + Hn − 1L Hn H − G2L
a2 =

If@Abs@a1D <= Abs@a2D, a = a1, a = a2D; H∗ Se alege cel care are numitorul mai mare ∗L
x1 = x1 − a; H∗ Noua valoare a radacinii dupa corectie ∗L
E;
x2 = AbsA NA10−2 RoundA102 x1EEE; H∗ Radacina rotunjita cu doua cifre ∗L
If@! MemberQ@ Xa, x2D, H∗ Testare daca radacina trunchiata a mai fost gasita ∗L
AppendTo@ Xa, x2D; If@Abs@Im@x1DD > ε, AppendTo@X, x1D, AppendTo@ X, Re@x1DDD;D;E;
H∗ Daca radacina rotunjita nu este in lista Xa,

in lista X, taind partea imaginara daca aceasta este neglijabila ∗L

atunci se introduce radacina nerotunjita

Print@"Nr de incercari=", jD;

Print@"Nr total de iteratii=", iD;
Print@"Radacini:"D;
Print@TableForm@Sort@ XDDD; H∗ Aranjeaza radacinile crescator si le afiseaza ca tabel ∗L
Print@"Valoarea maxima a functiei in radacini=", Max@Abs@f@ XDDDD;
NSolve@f@xD 0, xD

O varianta de rezultat obtinut prin acest progaram este urmatoarea:

Nr de incercari=41
Nr total de iteratii=204
Radacini:
−1.05057 − 0.551074
−1.04202 + 0.657252
−0.6718 − 0.430802
−0.376545 + 0.732874
−0.10401 + 0.22578
−0.0140766 − 0.859236
0.435058 + 0.918046
0.438331 − 0.779951
0.87771 + 0.134909
1.50792 − 0.0477977
Valoarea maxima a functiei in radacini=1.27737 × 10−12

3.10 Comparatie a metodelor de determinare a radacinilor ecuatiilor

neliniare

Nici una dintre metodele de determinare a radacinilor ecuatiilor neliniare nu

acopera toate situatiile posibile. De aceea, cel mai corect ar fi sa fie folosite doua
metode in tandem, care sa se completeze una pe cealalta. Pentru alegerea lor se poate
folosi tabelul 8.1 in care sunt prezentate caracteristici ale metodelor expuse.

Tabelul 8.1 Caracteristicile si utilizarea metodelor de determinare a radacinilor

ecuatiilor neliniare

Metoda Bisectie Secanta Newton Muller Lobacevski Bairstow Laguerre

Raphson Graeffe
Functie Oarecare Oarecare Oarecare Oarecare Polinom Polinom Polinom
Radacini Nu Nu Da Da Nu Da Da
complexe
Multiplicitate Nu Nu Nu Nu Da Da Nu
radacini
Puncte start Interval Interval 1 3 Nu 2 1
Convergenta Minima Medie Foarte Mare Medie Mare Mare
mare
Dezavantaj Convergenta Radacini Necesita Start Precizie Erori la Necesita
principal mica reale derivata complicat redusa deflatie derivata
Utilizare Braketing Refining Refining Refining Braketing Refining Refining
recomandata

In special atunci cand nu se cunoaste nimic despre pozitia radacinilor

polinoamelor, este recomandabil sa se procedeze in doua etape: mai intai localizarea
radacinilor (braketing) printr-o metoda care nu necesita puncte de start, cum ar fi
metoda Lobacevski-Graeffe (care permite si determinarea ordinului de multiplicitate),
urmata de etapa de rafinare (refining) printr-o metoda care permite atat determinarea
radacinilor reale cat si a celor complexe.
 Capitolul 4
Sisteme de ecuatii neliniare
In capitolul precedent s-a exemplificat deja rezolvarea unui sistem de doua
ecuatii neliniare, intalnit la metoda Bairstow. Intr-un spatiu tridimensional, problema
sistemului de doua ecatii neliniare consta in gasirea unei perechi de valori ( x1 , x2 ) din
planul x1 0 x2 pentru care doua functii de aceste variabile F1 ( x1 , x2 ) respectiv
F2 ( x1 , x2 ) sa se anuleze simultan.
Cele doua functii reprezinta suprafete intr-un spatiu tridimensional, ale caror
intersectii cu planul x1 0 x2 sunt curbe formate de punctele de anulare ale functiilor.
Solutiile sistemului de doua ecuatii sunt punctele de intersectie ale acestor curbe in
planul x1 0 x2 . Evident, este posibil sa existe mai multe solutii, sau ca multimea
solutiilor sa fie vida, daca aceste curbe nu se intersecteaza in nici un punct. Problema
este prin urmare mai dificila decat in cazul unidimensional, si gradul de dificultate
creste puternic cu dimensiunea sistemului. Aceasta in primul rand deoarece, spre
deosebire de cazul unidimensional, nu exista metode simple de localizare a
radacinilor. De asemenea, daca exista mai multe solutii, alegerea lor trebuie facuta pe
baza analizarii semnificatiei acestora in cadrul modelului (de exemplu fizic) care a
generat problema. Nu in ultimul rand, trebuie remarcat ca volumul de calcule poate sa
creasca exponential cu numarul de ecuatii, ceea ce face ca problema convergentei sa
fie deseori critica. De aceea, in cele ce urmeaza se vor analiza doar acele metode,
relativ putine, care asigura o convergenta acceptabila si se vor omite metodele care
prezinta doar un interes academic.
Exemplul precedent se poate generaliza pentru un sistem de n ecuatii neliniare,
problema formulandu-se astfel: trebuie gasit setul de valori ( x1 , x2 ,..., xn ) pentru care
functiile F1 ( x1 , x2 ,..., xn ) , F2 ( x1 , x2 ,..., xn ) ,..., Fn ( x1 , x2 ,..., xn ) sa se anuleze simultan:
⎧ F1 ( x1 , x2 ,..., xn ) = 0
⎪ F ( x , x ,..., x ) = 0
⎪ 2 1 2
⎨
n
(4.1)
⎪ ..............................
⎪⎩ Fn ( x1 , x2 ,..., xn ) = 0
Acest sistem poate fi scris matriceal sub forma:
F ( x) = 0 (4.2)
unde F = ( F1 , F2 ,..., Fn ) este vectorul coloana al valorilor functiei iar x = ( x1 , x2 ,..., xn )
este vectorul coloana al variabilelor.

4.1 Metoda Newton-Raphson pentru sisteme de ecuatii

neliniare

Aceasta metoda este preferata in majoritatea cazurilor, poate in primul rand

datorita popularitatii metodei analoage pentru cazul unidimensional. Totusi, s-au
impus si sunt recomandabile si metode alternative, datorita unor dezavantaje ale
metodei Newton-Raphsn care vor reiesi pe parcurs.
Se poate face o asemanare formala cu ecuatiile intalnite in cazul
unidimensional, pornindu-se de la ecuatia (4.2) si facandu-se dezvotarea in serie
185H

Taylor a functiei din membrul stang in jurul unui punct intr-un spatiu n-dimensional.
 n
∂Fi (x)
Fi (x + δ x) = Fi (x) + ∑ δ x + 0 (δ x) 2 , i = 1, 2,..., n (4.3)
j =1 ∂x j

Deoarece derivatele partiale ale tuturor componentelor functiei in raport cu variabilele

formeaza matricea Jacobiana
∂F (x)
J ij = i (4.4)
∂x j
ecuatia (4.3) se poate scrie sub forma matriceala:
186H

F(x + δ x) = F(x) + J ⋅ δ x + 0 (δ x) 2 (4.5)

Deoarece se urmareste ca dupa corectie F(x + δ x) sa fie cat mai aproape de 0,
din (4.5) rezulta corectorul:
187H

δ x = −J −1 ⋅ F(x) (4.6)
Metoda se aplica iterativ, astfel ca la iteratia k+1 se obtine noul vector
radacina:
x k +1 = x k + δ x k (4.7)
Convergenta este foarte rapida, ordinul de convergenta fiind 2, dar numai daca
alegerea initiala este suficient de apropiata de solutia corecta, in caz contrar metoda
fiind sau slab convergenta sau chiar divergenta. Deoarece alegerea initiala presupune
atatea valori cate ecuatii sunt in sistem, ea este uneori mult mai dificila decat in cazul
unidimensional. Aceasta limiteaza utilizarea metodei Newton–Raphson in cazul
sistemelor de ecuatii neliniare, existand o serie de metode alternative care pot fi
aplicate cu succes in diverse situatii.

Capitolul 5
Derivarea numerica
5.1 Derivarea directa

5.2 Derivarea prin interpoalre

Exista numeroase metode alternative pentru derivarea numerica, dintre care
sunt mai importante cele care folosesc decompunerea polinomului de interpolare dupa
un set de functii ortogonale, cu aplicatii directe in metodele spectrale de rezolvare a
ecuatiilor diferentiale.
Sa consideram un set de valori ale functiei (necunoscuta analitic) ale carei
derivate dorim sa le calculam. O posibilitate evidenta de obtinere a derivatelor o
constituie gasirea unei functii de interpolare, in forma analitica, aproximand functia
necunoscuta, printr-o metoda de tipul celor prezentate anterior. Noua functie poate fi
deci derivata analitic si apoi se pot determina valorile derivatei in punctele dorite, care
pot fi si altele decat cele initiale.
Interpolarea se poate face folosind diverse seturi de functii, fiind preferate cele
cu proprietati de ortogonalitate, care permit calculul simplu al coeficientilor
dezvoltarii. De exemplu, se pot folosi functiile trigonometrice sau polinoame
ortogonale, cele mai raspandite fiind polinoamele Chebyshev datorita calculului
simplu si proprietatilor de minimizare a erorilor sugerate de teoreme mini-max.
 Pe langa polinoamele Chebyshev de speta 1, prezentate deja in capitolul II si
definite prin formula:
Tn ( x) = cos[n arccos( x)] (5.1)
se mai definesc si polinoamele Chebyshev de speta 2, prin formula:
sin[( n + 1) arccos( x)]
U n ( x) = (5.2)
sin x
S-a aratat ca una dintre metodele de interpolare a functiilor se bazeaza pe
descompunerea acestora intr-o serie de N polinoame Chebyshev:
N −1
f ( x) = ∑ cnTn ( x) (5.3)
n =0

unde coeficientii ck se calculeaza tinand seama de relatiile de ortogonalitate:

⎧
1
⎪0, m ≠ n
Tn ( x)Tm ( x) ⎪
∫−1 1 − x 2 dx = ⎨π , m = n = 0 (5.4)
⎪π
⎪ ,m = n ≠ 0
⎩2
Inmultind relatia (5.3) cu Tm ( x) integrand si tinand cont de ralatia (5.4),
188H 189H

rezulta coeficientii dezvoltarii ca produse scalare ale functiei cu vectorii ortogonali

Tn ( x) :
1
1 f ( x)Tk ( x)
ck =
αk ∫ −1 1 − x 2
dx

⎧π , k = 0 (5.5)
⎪
α k = ⎨π
⎪⎩ 2 , k ≠ 0
Aceasta relatie nu este insa practica deoarece necesita o integrare numerica,
introducand erori suplimentare si consumand timp de calcul. Coeficientii dezvoltarii
(5.3) se pot calcula insa si fara a fi necesara integrarea, tinand seama de relatia de
190H

ortogonalitate discreta a polinaomelor Chebyshev:

⎧0, n ≠ m
⎪N
N
⎪
∑ Tn ( xk )T m ( xk ) = ⎨ ,n = m ≠ 0
2
(5.6)
k =1 ⎪
⎪⎩ N , n = m = 0
unde xk sunt radacinile polinomului Chebyshev de ordinul N , obtinute din relatia de
definitie:
TN ( xk ) = cos( N arccos xk ) = 0
Rezulta
π
arccos xk = (2k − 1)
2N
(5.7)
⎡⎛ 1⎞π ⎤
xk = cos ⎢⎜ k − ⎟ ⎥
⎣⎝ 2⎠ N ⎦
Intr-adevar, membrul stang al relatiei (5.6) se scrie succesiv:
191H
 N N
⎡ ⎛ 1⎞π ⎤ ⎡ ⎛ 1⎞π ⎤
S = ∑ Tn ( xk )Tm ( xk ) = ∑ cos ⎢ n ⎜ k − ⎟ ⎥ cos ⎢ m ⎜ k − ⎟ ⎥ =
k =1 k =1 ⎣ ⎝ 2⎠ N⎦ ⎣ ⎝ 2⎠ N ⎦
1 N ⎡ ⎛ 1⎞π ⎤ 1 N ⎡ ⎛ 1⎞ π ⎤
∑
2 k =1
cos ⎢
⎣
( n + m ) ⎜
⎝
k − ⎟ ⎥ + ∑ cos ⎢(n − m) ⎜ k − ⎟ ⎥ =
2 ⎠ N ⎦ 2 k =1 ⎣ ⎝ 2⎠ N ⎦
(5.8)

1 N ⎡ i ( n + m )⎛⎜ k − 12 ⎞⎟ πN ⎤ 1 N ⎡ i ( n − m )⎛⎜ k − 12 ⎞⎟ πN ⎤
∑ Re ⎢e
2 k =1 ⎢⎣
⎝ ⎠
⎥ + ∑ Re ⎢ e
2
⎝ ⎠
⎥
⎥⎦ k =1 ⎢⎣ ⎥⎦
S1 S2

Cei doi termeni pot fi calculati ca sume de progresii geometrice. De exemplu,

π π
i ( n+m) i ( n+m)
termenul S1 are primul termen a1 = e 2N
si ratia r = e N
, de unde rezulta:

1 ⎡⎢ i ( n + m ) 2πN 1 − ei ( n + m )π ⎤
⎥
S1 = Re e π
2 ⎢ i ( n+m) ⎥
⎣ 1− e N
⎦
π
i(n+m)
Impartind la e 2N
se poate obtine la numitor un sinus, si tinand cont ca
i ( n + m )π
e = (−1) n + m , rezulta:
⎡ ⎤
1 ⎢ i 1 − (−1) n + m ⎥
S1 = Re ⎢ ⎥ (5.9)
2 ⎢ 2 sin(n + m) π ⎥
⎣ 2N ⎦
Similar se obtine si celalalt termen:
⎡ ⎤
1 ⎢ i 1 − (−1) n − m ⎥
S 2 = Re ⎢ ⎥ (5.10)
2 ⎢ 2 sin(n − m) π ⎥
⎣ 2N ⎦

Se observa ca avem trei cazuri:

a) Daca n = m = 0 , ambele sume din ultimul membru al relatiei (5.8) au 192H

exponentul nul si deci suma totala este N ;

b) Daca n = m ≠ 0 , prima suma va fi nula conform relatiei (5.9) in care in 193H

paranteza dreapta avem un numar pur imaginar iar a doua suma se poate
calcula din termenul final al relatiei (5.8) in care exponentul e nul, suma totala
194H

N
fiind deci ;
2
c) Daca n ≠ m , si prima si a doua suma sunt nule conform relatiilor (5.9) si 195H

(5.10), in care in parantezele dreapte avem numere pur imaginare, deci suma
196H

totala este nula.

Relatia (5.6) este deci demonstrata si poate fi folosita pentru calculul
197H

coeficientilor dezvoltarii (5.3). Pentru inceput se scrie aproximarea functiei initiale

198H

prin polinoame Chebyshev conform relatiei (5.3), pentru punctele xk , 199H

N −1
f ( xk ) = ∑ cnTn ( xk ) (5.11)
n =0

Inmultind aceasta relatie cu Tm ( xk ) , sumand dupa k si tinand cont de relatia

de ortogonalitate discreta (5.6), obtinem:
200H
 N N N −1

∑
k =1
f ( xk )Tm ( xk ) = ∑∑ cnTn ( xk )Tm ( xk ) =
k =1 n = 0

N −1 ⎧ Nc0 , m = 0
N
⎪
∑ cn ∑ Tn ( xk )Tm ( xk ) = ⎨ N
n=0 k =1 ⎪⎩ 2 cm , m ≠ 0
Coeficientii dezvoltarii vor fi deci:

1 N
c0 =
N
∑ f ( x )T ( x )
k =1
k 0 k (5.12)

2 N
cm = ∑ f ( xk )Tm ( xk ), m = 1, 2,..., N − 1
N k =1
(5.13)

sau, tinand cont de definitia (5.7) a radacinilor polinoamelor Chebyshev

201H

2 N ⎧ ⎡⎛ 1 ⎞ π ⎤⎫ ⎡ ⎛ 1⎞π ⎤
cm = ∑ f ⎨cos ⎢⎜ k − ⎟ ⎥ ⎬ cos ⎢ m ⎜ k − ⎟ ⎥, m = 1, 2,..., N − 1 (5.14)
N k =1 ⎩ ⎣⎝ 2 ⎠ N ⎦⎭ ⎣ ⎝ 2⎠ N⎦
Aceasta relatie defineste asa-numita transformare cosinus discreta si pentru
calculul sau exista si metode rapide bazate pe transformata Fourier rapida (FFT).
Pentru a nu fi necesare doua formule (adica (5.12) si (5.13)) pentru calculul 202H 203H

coeficientilor dezvoltarii, in practica se prefera o forma usor modificata pentru relatia

(5.11):
204H

N −1
1
f ( x) ≈ ∑ cnTn ( x) + c0 (5.15)
n =1 2
in care toti coeficientii, inclusiv c0 , se calculeaza cu aceeasi relatie (5.13), dar se tine 205H

seama ca din aceasta relatie c0 rezulta de doua ori mai mare decat din (5.12); de 206H

aceea se separa termenul respectiv in relatia (5.15) si se incepe indexarea in suma de

207H

la 1. S-a tinut cont si de faptul ca T0 ( x) = 1 .

Derivata functiei f ( x) se poate calcula din relatia (5.3), in care functiile 208H

dependente de x din membrul drept sunt polinoamele Chebyshev:

N −1
f '( x) = ∑ cnTn' ( x) (5.16)
n =0
cu coeficientii calculati dupa (5.12) si (5.13), sau rescrisa cu coeficienti calculati
209H 210H

numai dupa (5.13), sub o forma similara cu (5.15):

211H 212H

N −1
1
f '( x) = ∑ cnTn' ( x) + c0 (5.17)
n =1 2
Aceste forme au dezavantajul ca trebuiesc cunoscute si derivatele
polinoamelor Chebyshev, ceea ce consuma timp suplimentar. Exista insa posibilitatea
ca aceste derivate ale polinoamelor sa se exprime prin relatii de recurenta in functie de
polinoamele Chebashev nederivate, de care se presupune ca se dispune deja pentru
calculul coeficientilor cn .
Se pot demonstra urmatoarele relatii:
d
Tn ( x) = nU n −1 ( x) (5.18)
dx
1
Tn ( x) = [U n ( x) − U n − 2 ( x) ] (5.19)
2
Prima relatie se demonstreaza prin calcul direct:
 d n sin( n arccos x)
cos(n arccos x ) = =
dx 1 − x2
sin(n arccos x) sin( n arccos x)
=n =n = nU n −1 ( x)
1 − [cos(arccos x)]2 sin(arccos x )
A doua relatie se demonstreaza pe baza identitatii trigonometrice:
1
cos a sin b = [sin(a + b) − sin(a − b) ]
2
Notand θ = arccos x si alegand a = nθ , b = (n − 1)θ se obtine succesiv:
1
cos nθ sin θ = [sin(n + 1)θ − sin(n − 1)θ ]
2
1 ⎧ sin [ (n + 1)θ ] sin [ (n − 1)θ ] ⎫
cos( n arccos x ) = ⎨ − ⎬=
2⎩ sin θ θ ⎭
1 ⎧ sin [ ( n + 1) arccos x ] sin [ (n − 1) arccos x ] ⎫
= ⎨ − ⎬
2 ⎩ sin(arccos x) sin(arccos x) ⎭
obtinandu-se relatia (5.19). 213H

Rescriind relatia (5.18) pentru n + 1 si n − 1 se obtine

214H

d
Tn +1 ( x) = (n + 1)U n ( x)
dx
d
Tn −1 ( x) = (n − 1)U n − 2 ( x)
dx

Impartind prima relatie cu n + 1 si a doua cu n − 1 , adunandu-le si tinand

seama de relatia (5.19), rezulta:
215H

Tn'+1 ( x) Tn'−1 ( x)
2Tn ( x) = − (5.20)
n +1 n −1
Aceasta relatie este fundamentala pentru metodele spectrale, permitand
calculul prin recurenta al coeficientilor seriilor derivate.
Astfel, se pot exprima derivatele polinoamelor Chebyshev din relatia (5.16) in 216H

functie de polinoamele nederivate, obtinandu-se o noua serie cu alti coeficienti bn

care pot fi calculati prin recurenta din coeficientii deja cunoscuti cn . Pentru gasirea
relatiei intre cele doua seturi de coeficienti, se scrie (5.16) in cele doua forme si se
217H

egaleaza coeficientii polinoamelor Chebyshev de acelasi ordin.

N −1 N −1
f '( x) = ∑ cnTn' ( x) = ∑ bnTn ( x) (5.21)
n =0 n=0
Din relatia (5.20) rezulta derivata polinomului Chebyshev de ordin n + 1 in
218H

functie de polinomul nederivat de ordin n :

n +1 '
Tn'+1 ( x) = Tn −1 ( x) + 2(n + 1)Tn ( x)
n −1
care rescrisa trecand de la n + 1 la n este:
n
Tn' ( x) = Tn'− 2 ( x) + 2nTn −1 ( x) (5.22)
n−2

Tot din relatia (5.20) rezulta derivata polinomului Chebyshev de ordin n − 1

219H

in functie de polinomul nederivat de ordin n :

 n −1 '
Tn'−1 ( x) =
Tn +1 ( x) − 2(n − 1)Tn ( x)
n +1
care rescrisa trecand de la n − 1 la n este:
n
Tn' ( x) = Tn'+ 2 ( x) − 2nTn +1 ( x) (5.23)
n+2

Din (5.22)si (5.23) se observa ca la termenul in Tn' ( x) din suma din stanga
220H 221H

relatiei (5.21) contribuie termenii in Tn −1 ( x) si Tn +1 ( x) ai sumei din dreapta acestei

222H

relatii, avand coeficientii bn −1 si respectiv bn +1 :

cnTn' ( x) = bn −1Tn −1 ( x) + bn −1Tn −1 ( x) (5.24)
Tot din relatiile (5.22) si (5.23) rezulta ca factorii de proportionalitate intre
223H 224H

1 1
Tn −1 ( x) si Tn' ( x) , respectiv Tn +1 ( x) si Tn' ( x) sunt , respectiv − , deci inlocuind in
2n 2n
(5.24) obtinem:
225H

b b
cnTn' ( x) = n −1 Tn' ( x) − n −1 Tn' ( x)
2n 2n
Rezulta relatia de recurenta pentru calculul coeficientilor dezvoltarii derivatei:
bn −1 = bn +1 + 2ncn (5.25)
Ei se calculeaza cu conditiile de start cn +1 = 0 si cn = 0 , iar deoarece b0 se
calculeaza din c1 , el rezulta cu o valoare dubla fata de cea normala si trebuie
injumatatit. Se obtine deci formula de interpolare pentru derivata functiei:
N −1
1
f '( x) ≈ ∑ bnTn ( x) + b0 (5.26)
n =1 2
coeficientii fiind calculati prin recurenta cu relatia (5.25). Egalitatea este exacta
226H

numai in radacinile polinomului Chebyshev de ordinul N − 1 , in rest fiind cu atat mai

corecta cu cat N este mai mare.
Aceste relatii sunt valabile, ca si cele de interpolare, numai pe intervalul de
ortogonalitate al polinoamelor Chebyshev, adica [-1,1]. In cazul in care functia a carei
derivata o calculam are un domeniu de definitie [ xmin , xmax ] mai mare, trebuie facuta o
schimbare de variabila dupa formula:
x −x x + xmin
x = y max min + max (5.27)
2 2
aplicata pentru punctele de esantionare ale functiei, astfel incat relatia (5.13) devine:
227H

2 N x −x x + xmin
cm = ∑ f ( xk max min + max )Tm ( xk ), m = 0, 2,..., N − 1
N k =1 2 2
La calculul functiei de interpolare trebuie facuta trecerea inversa in variabila
polinoamelor Chebyshev, astfel incat formula (5.15) devine: 228H

N −1
x − xmin 1
f ( x) ≈ ∑ cnTn (2 − 1) + c0
n =1 xmax − xmin 2
De asemenea, la calculul functiei de interpolare a derivatei trebuie tinut cont
de factorul care apare la derivarea relatiei de schimbare de variabila (5.27), astfel ca
229H

relatia (5.26) devine:

230H

2 N −1
x − xmin 1
f '( x) ≈
xmax − xmin n =1
∑ bnTn (2
xmax − xmin
− 1) + b0
2
In[1]:= Needs@"Graphics`Colors`"D;
f@x_D = + CosAx2E;H∗Cuprinde functia Runge − interpolarea echidistanta nu converge! ∗L
1

xm = −5; xM = 5; H∗ Domeniul de definitie, extins in afara lui @−1,1D ∗L

1 + x2

g1 = Plot@f@xD, 8x, xm, xM<, PlotStyle → BlueD;

Na = 60; H∗ Numarul de termeni din dezvoltarea in polinoame Chebyshev∗L

H∗ Se calculeaza coeficientii de interpolare Chebyshev pentru f HxL ∗L

ForAk = 0, k <= Na, k++,
i H2 j − 1L y xM + xm i k H2 j − 1L
j z xM − xm j 2 y
z
c@kD = „ fAj
j
jCosA Ez
z Ej
j
jCosA Ez
z
zE;
Na π π

z
2. 2

k { k {
+
Na Na 2 2 Na
j=1

H∗ S−a facut schimbarea de variabila x=y xM−xm + xM+xm pentru extinderea intervalului ∗L

H∗p@x_D=‚ c@nD ChebyshevTAn,2 x−xm −1E+ 1 c@0D;

2 2
Na

Genereaza p HxL in forma polinomiala dar cu erori mai mari pentru Na mare∗L
n=1 xM−xm 2

p@x_D = ‚ c@nD CosAn ArcCosA2 − 1EE + c@0D;

Na x − xm 1
xM − xm 2
H∗ Forma preferabila in Mathematica 5.2, erori mai mici decat cu functia ChebyshevT@n,xD ∗L
n=1

H∗ S−a facut schimbarea inversa de variabila y=2 x−xm −1 ∗L

xM−xm

g2 = Plot@p@xD, 8x, xm, xM<, PlotStyle → GreenD;

Show@g1, g2D; H∗ Afiseaza functia initiala si aproximanta ∗L

H∗ Aici incepe derivarea numerica ∗L

b@Na + 1D = 0;
b@NaD = 0;
For@k = Na, k > 0, k−−,
b@k − 1D = b@k + 1D + 2 k c@kDD;
i
j y
dp@x_D = j
j
j ‚ b@nD CosAn ArcCosA2 − 1EE + b@0Dz
z
z
z;
2 Na x − xm 1
xM − xm k n=1 xM − xm 2 {
H∗ S−a inmultit cu dy din schimbarea de variabila ∗L
dx

df@x_D = ∂x f@xD; H∗ Derivata analitica pentru comparatie ∗L

g3 = Plot@df@xD, 8x, xm, xM<, PlotStyle → RedD;
g4 = Plot@dp@xD, 8x, xm, xM<, PlotStyle → GreenD;
Show@g3, g4D; H∗ Afiseaza derivata analitica si derivata numerica ∗L

Nm = 500; H∗ Nr de puncte in care se calculeaza abaterea ∗L

ForAk = 0, k < Nm, k++, r@kD = xm + k E;
xM − xm
Nm

‚ Abs@df@r@kDD − dp@r@kDDD;
1 Nm−1
e=
Nm k=0
Print@"Eroare=", eD;

Programul prezentat calculeaza functia de interpolare a unei functii date si a

derivatei acesteia, avand urmatoarele intrari:
• Functia initiala f ( x) ;
• Domeniul de definitie [ xm, xM ] ;
• Numarul de termeni ai polinomului de interpolare.
Iesirile sunt:
 • Polinomul de aproximare Chebyshev al functiei, p( x) ;
• Polinomul de aproximare Chebyshev al derivatei functiei, dp( x) ;
• Graficele functiei initiale, polinomului de aproximare (figura 5.1) si al
derivatelor acestora (derivata functiei initiala, considerata analitica in
exemplul de fata se calculeaza analitic pentru comparatie) – figura 5.2;
• Eroarea definita ca media modulurilor abaterilor intr-un numar de puncte de
comparare, Nm .

Pentru o aproximare cu 60 de termeni de exemplu, eroarea obtinuta pentru

aceasta functie destul de dificila, este de aproximativ 3.2 × 10−5 .

3.5

2.5

1.5

-4 -2 2 4
4

3.5

2.5

1.5

-4 -2 2 4
Figura 5.1 Graficul functiei f ( x) si al polinomului de aproximare p( x)
 10

-4 -2 2 4

-5

-10
10

-4 -2 2 4

-5

-10
Figura 5.2 Graficul derivatei analitice a functiei, f '( x) , si al polinomului de
aproximare numerica a derivatei, dp ( x) .
 Capitolul 6
Integrarea numerica
b

∫ f ( x)dx
a
(6.1)

6.1 Integrarea Gauss

Fie integrala definita pe un domeniu finit a unei functii integrabile f ( x) :

b

∫ f ( x)dx
a
(6.2)

Metodele elementare prezentate anterior folosesc pentru calculul valorilor

integrandului o divizare echidistanta a intervalului de integrare. Printr-o alegere
avantajoasa a punctelor de evaluare a integrantului (puncte numite si noduri), se pot
obtine metode mult mai eficiente de integrare, printre care si metoda de cuadratura
Gauss.
Functia f ( x) poate fi aproximata prin una dintre metodele cunoscute, de
exemplu prin dezvoltare Lagrange cu n termeni:
n
f ( x) ∑ f ( x ) L ( x)
i =1
i i (6.3)

unde polinomul Lagrange are forma cunoscuta:

n
x − xi
Li ( x) = ∏ (6.4)
j =1 x j − xi
j ≠i

si deoarece in (6.4) exista conditia ca j ≠ i , acest polinom este de grad n − 1 .

231H

Introducand (6.4) in (6.3) si apoi in (6.2), rezulta:

232H 233H 234H

x − xi
b n b n

∫a f ( x ) dx ∑
i =1
f ( xi ∫∏
)
x j − xi
dx (6.5)
a j =1
j ≠i

Rezulta ca prin aceasta formula se va aproxima integrala functiei print

integrala unui polinom de grad n − 1 , care trece prin n puncte xi . In general, cu cat n
este mai mare cu atat aproximatia functiei si deci si a integralei va fi mai buna. Intr-
adevar, conform formule pentru restul aproximarii unei functii printr-un polinom de
grad n care trece prin n puncte xi care a fost dedusa la studiul interpolarii polinomiale:
f ( n ) (ξ ) n
Rn −1 ( x) =
n ! i =1
∏ ( x − xi ) (6.6)

Se observa o scadere a erorii cu n ! , ceea ce in unele cazuri este suficient.

Dupa cum se va arata in continuare, exista insa posibilitatea de a aproxima
integrala printr-un polinom de grad 2n − 1 folosind o dexvoltare cu numai n termeni,
daca acest polinom face parte dintr-o baza de polinoame ortogonale.
Sa consideram o aproximare a functiei f ( x ) printr-o functie polinomiala de
grad 2n − 1 , f 2 n −1 ( x) . Aceasta functie o factorizam cu un polinom de grad n , Pn ( x) ,
obtinandu-se un cat q ( x) de grad (2n − 1) − n = n − 1 si un rest r ( x) de asemenea de
grad n − 1 :
f 2 n −1 ( x) = Pn ( x)qn −1 ( x) + rn −1 ( x) (6.7)
 Daca impunem ca polinomul Pn ( x ) sa faca parte dintr-o baza de polinoame
{Pn ( x)} ortogonale pe intervalul [ a, b] , el satisface o relatie de ortogonalitate de tipul:
b

∫ P ( x) P ( x)w( x)dx = c δ
2
n m n nm (6.8)
a

unde w( x) este functia pondere ale bazei de polinoame ortogonale considerate.

Inmultind relatia (6.7) cu w( x) si integrand pe intervalul [ a, b ] obtinem:
235H

b b b

∫a
f 2 n −1 ( x) w( x)dx = ∫ Pn ( x)qn −1 ( x) w( x)dx + ∫ rn −1 ( x) w( x)dx
a a
(6.9)

Exprimand catul ca o dezvoltare dupa polinoamele din baza ortogonala

{Pn ( x)}
n −1
qn −1 ( x) = ∑ ak Pk ( x) (6.10)
k =1
Se observa ca deoarece catul are gradul n − 1 , si polinoamele folosite din baza
au gradul cel mult n − 1 . Ele sunt deci toate ortogonale pe polinomul Pn ( x) , si deci
conform relatiei (6.8), prima integrala din membrul drept al relatiei (6.9) se anuleaza:
236H 237H

⎡ n −1 ⎤
b b

∫a n n−1
P ( x ) q ( x ) w( x ) dx = ∫a n ⎢⎣∑
P ( x )
k =1
ak Pk ( x) ⎥ w( x)dx =
⎦
n −1 b
= ∑ ak ∫ Pn ( x)Pk ( x)w( x)dx = cn2δ nk = 0
k =1 a

deoarece k ≤ n − 1 , deci k ≠ n .
Relatia (6.9) devine:
238H

b b

∫
a
f 2 n −1 ( x) w( x)dx = ∫ rn −1 ( x) w( x)dx
a
(6.11)

Sa consideram acum radacinile polinomului ortogonal Pn ( x ) , adica acele

puncte xi pentru care avem relatia:
Pn ( xi ) = 0, i = 1, 2,..., n (6.12)
Scriind relatia (6.7) in aceste puncte rezulta:
239H

f 2 n −1 ( xi ) = rn −1 ( xi ) (6.13)
Restul rn −1 ( x ) poate fi interpolat Lagrange printr-o relatie de tipul (6.3) si 240H

tinand cont de relatia (6.13) se poate scrie:

241H

n n
rn −1 ( x) ∑ rn−1 ( xi ) Li ( x) = ∑ f 2n−1 ( xi ) Li ( x)
i =1 i =1
(6.14)

Introducand relatia (6.14) in (6.11) rexulta 242H 243H

b n b

∫ f 2 n −1 ( x) w( x)dx = ∑ f 2 n −1 ( xi ) ∫ Li ( x) w( x)dx (6.15)

a i =1 a

Tinand cont ca f 2 n −1 ( x) aproximeaza pe f ( x) , relatia (6.15) se poate scrie: 244H

b n

∫ f ( x) w( x)dx ∑ f (x ) W
i =1
i i (6.16)
a

unde marimile
 b
Wi = ∫ Li ( x) w( x)dx , i = 1, 2,..., n (6.17)
a

se numesc ponderi de integrare Gauss. Ele sunt independente de functia de integrat,

dupa cum arata relatia (6.17) si pot fi calculate cu aceasta relatie din functia pondere
245H

w( x) a bazei de polinoame ortogonale considerate si definitia (6.4) a polinomului 246H

Lagrange Li ( x) pentru nodurile xi radacini ale polinomului ortogonal de ordin maxim

folosit, sau se pot obtine din tabele.
Aproximatia data de aceasta metoda corespunde unei erori de interpolare a
functiei cu un polimom de grad 2n − 1 , desi functia se calculeaza numai in n puncte.
Restul aproximarii va fi analog relatiei (6.6), dar pentru un grad 2n − 1 :
247H

f ( n ) (ξ ) n
Rn −1 ( x) = ∏ ( x − xi )
(2n − 1)! i =1
(6.18)

find deci de (2n − 1)!/ n ! = (n + 1)(n + 2)...(2n − 1) ori mai mic. De exemplu, pentru
n = 10 , restul poate fi de 1.4 ×1012 ori mai mic, ceea ce indica o metoda
hipereficienta.
Relatia (6.16) defineste metoda de quadratura Gauss in mod general. Se
248H

remarca faptul ca integrandul contine functia pondere a polinomului ortogonal, de

aceea este necesar ca functia de integrat, de exemplu g ( x) sa fie in prealabil
factorizata in aceasta forma
g ( x) = f ( x) w( x) (6.19)
De aici rezulta

g ( x)
f ( x) = (6.20)
w( x)
iar relatia (6.16) devine
249H

b
g ( xi ) n

a
∫ g ( x)dx = ∑
i =1 w ( xi )
Wi

Trecand din nou la functia f ( x) , obtinem in final

b n
f ( xi )
∫ f ( x)dx = ∑ Wi (6.21)
a i =1 w( xi )

6.2 Integrarea de tip Gauss-Legendre

Daca se folosesc polinoamele Legendre ca baza, avem avantajul ca functia lor

pondere este 1, simplificand relatia de cuadratura. Rezulta formula de cuadratura
Gauss-Legendre:
b n

∫ f ( x)dx = ∑ f ( xi )Wi (6.22)

a i =1

unde ponderile sunt:

x − xi
b b n
Wi = ∫ Li ( x)dx = ∫ dx∏ (6.23)
a a j =1 x j − xi
j ≠i
nodurile xi fiind radacinile polinomului Legendre de ordin n . Pentru calculul
nodurilor xi trebuie rezolvata ecuatia polinomiala (6.12) pentru polinoame Legendre,
250H

sau se pot lua din tabele. De asemenea, exista tabele si pentru aceste functii pondere.
Se mai poate folosi de asemenea si formula:
P ( x), Pn −1 ( x) An
Wi = n −1 (6.24)
Pn −1 ( xi ) Pn'−1 ( xi ) An −1
unde An si An −1 sunt coeficientii termenilor de grad maxim din polinoamele Pn ( x) ,
respectiv Pn −1 ( x) .
O problema care apare in practica este datorata faptului ca domeniul de
ortogonalitate al polinoamelor Legendre este [ −1,1] . Daca integrala este definita pe
un domeniu oarecare, dar finit, [ a, b ] , trebuie facuta schimbarea de variabila:
b−a b+a
x+ z= (6.25)
2 2
astfel ca formula de cuadratura Gauss Legendre devine:

b−a n ⎛b−a b+a⎞

b

∫ f ( x)dx = ∑
2 i =1
f⎜
⎝ 2
xi + ⎟Wi
2 ⎠
(6.26)
a

iar ponderile raman neschimbate.

6.3 Cuadratura Clenshaw Curtis

Tinand cont de argumentul dat de relatia (6.18), se pare ca metoda de

251H

cuadratura Gauss este cea mai eficienta din cate se cunose. Totusi, in anumite cazuri
exista si argumente in favoarea unor alte metode, cum ar fi Clenshaw-Curtis. Astfel,
se pot mentiona urmatoarele:
• Cuadratura Clenshaw-Curtis are formule mult mai simple de calcul a nodurilor
si ponderilor decat Gauss-Legendre si deci erorile introduse de acestea sunt
mai mici la un numar de puncte de esantionare mare;
• In varianta de esantionare Gauss-Lobatto metodele adaptive pot reutiliza
jumatate din punctele anterior calculate la trecerea de la n puncte la 2n puncte
de esantionare.
Aceste argumente sunt importante in special atunci cand trebuiesc efectuate
numeroase cuadraturi pentru rezolvarea unei probleme, cum ar fi in cazul ecuatiilor
integrale si integro-diferentiale, asa cum se va arata mai tarziu.
Pentru prezentarea algoritmului Clenshaw-Curtis se porneste de la exprimarea
functiei de integrat ca o dezvoltare in serie dupa polinoame Chebyshev Tr ( x) de ordin
cel mult n , conform relatiei (5.11) prezentata in capitolul anterior, unde semnul prim
252H

al sumei arata ca primul termen se ia injumatatit

n
f ( x) = ' a T ( x)
∑ r r (6.27)
r =0

Coeficientii ar ai dezvoltarii rezulta din (5.13) inclusiv pentru r = 0

253H

2 n
ar = ∑ f ( xk )Tr ( xk ), r = 0,1,..., n
n k =1
(6.28)
 Vom considera in continuare ca esantionarea se face in punctele xk date de
radacinile polinomului Chebyshev de ordinul n , deci
Tn ( xk ) = cos n arccos xk = 0, n arccos xk = (2k − 1)π / 2, k = 1, 2,..n
(2k − 1)π / 2
xk = cos . (6.29)
n
In anumite cazuri se poate considera si esantionarea de tip Lobatto care se face
kπ
in extremele polinomului Chebyshev de ordinul n adica xk = cos . Aceasta
n
metoda are avantajul ca, daca se dubleaza numarul de puncte de esantionare (in cadrul
asa-numitelor metode adaptive de cuadratura), jumatate din punctele obtinute coincid
cu cele vechi si deci valorile respective pot fi reutilizate. Totusi se constata in calcule
ca erorile introduse in schema de cuadratura cu esantionarea in extreme sunt mai mari
decat cele introduse de esantionarea in zerouri astfel ca daca nu se folosesc algoritmi
adaptivi varianta data de relatia (6.29) este mai buna.
254H

Pe baza dezvoltarii functiei in serie de polinoame Chebyshev (6.27) se poate 255H

calcula integrala definita (6.1) pe intervalul [ −1,1] in forma

256H

1 1 1 1 1
' a T ( x)dx = ' a T ( x)dx = a T ( x)dx + 1 a T ( x)dx (6.30)
n n n

∫−1 f ( x)dx = −∫1 ∑

r =0
r r ∑
r =0
r ∫ r
−1
∑
r =1
r ∫ r
−1
2 −∫1
0 0

deci problema se reduce la calculul coeficientilor dezvoltarii prin relatia (6.28) si al 257H

integralelor polinoamelor Chebyshev pana la ordinul n .

Ca si in cazul cuadraturii Gauss, aceste integrale nu depind de functia initiala,
fiind universal valabile si deci pot fi gasite in tabele. Chiar si intr-un caz mai general,
cu limita superioara de integrare variabila (ca in cazul ecuatiilor integrale Volterra
care vor fi prezentate mai tarziu), integrala initiala se poate calcula si fara a se
cunoaste aceste noi integrale, folosindu-se doar coeficientii dezvoltarii functiei in
serie de polinoame Chebyshev.
Intr-adevar, integralele de forma
t

∫ T ( x)dx
−1
r cu r = 0,1,..., n se pot calcula analitic cu ajutorul relatiei (5.20). Rezulta
258H

t t
t t
⎡ Tr'+1 ( x) Tr'−1 ( x) ⎤ Tr +1 ( x) Tr −1 ( x) Tr +1 (t ) Tr −1 (t ) (−1)r +1
∫−1Tr (x)dx = −∫1 ⎢⎣ 2(r +1) − 2(r −1) ⎥⎦ dx = 2(r +1) − 2(r −1) = 2(r +1) − 2(r −1) + r2 −1 (6.31)
−1 −1

unde s-a tinut cont ca Tr −1 (−1) = (−1) = (−1) (−1) = (−1) r +1 = Tr +1 (−1) .
r −1 2 r −1

Relatia (6.31) este utilizabila pentru r = 2,3,..., n , astfel ca integralele pentru

259H

n = 0,1 trebuiesc calculate separat:

t
n = 0 ⇒ ∫ dx = t − (−1) = T1 (t ) + 1
−1
(6.32)
t 2 (−1) 2 1
t
1
n = 1 ⇒ ∫ xdx = − = T2 (t ) −
−1
2 2 4 4
unde s-a tinut cont ca primele polinoame Chebyshev
sunt T0 (t ) = 1, T1 (t ) = t , T2 (t ) = 2t − 1 .
2

Introducand relatiile (6.31) si (6.32) in ultimul membru al relatiei (6.30) si

260H 261H 262H

ordonand dupa ordinul polinoamelor Chebyshev care apar, rezulta

 (−1) r +1 ⎛a a ⎞ ⎛a a ⎞
t
a0 a1 n
∫ f ( x)dx = − + ∑ ar 2 + T1 (t ) ⎜ 0 − 2 ⎟ + T2 (t ) ⎜ 1 − 3 ⎟
−1
2 4 r =2 r −1 ⎝ 2 2⎠ ⎝4 4⎠
⎛a a ⎞ ⎛a a ⎞
+T3 (t ) ⎜ 2 − 4 ⎟ + ... + Tn (t ) ⎜ n −1 − n +1 ⎟
⎝ 6 6⎠ ⎝ 2n 2n ⎠
Daca se ia intervalul de integrare [ −1, t ] , tinand cont ca Tr (1) = 1 , se obtine
formula de cuadratura Clenshaw-Curtis
a0 a1 n −1 (−1) r +1 n a −a
t

∫−1 f ( x ) dx = − + ∑ ar 2 + ∑ Tr (t ) r −1 r +1 (6.33)
2 4 r =2 r − 1 r =1 2r
in care coeficientii dezvoltarii se cacluleaza cu relatia (6.28) care poate fi explicitata
263H

ca o transformare cosinus discreta

2 n ⎛ (2k − 1)π / 2 ⎞ r (2k − 1)π / 2
ar = ∑ f ⎜ cos ⎟ cos , r = 0,1,..., n
n k =1 ⎝ n ⎠ n
iar coeficientul ar +1 = 0 , deoarece seria (6.27) se considera trunchiata la n .
264H

Ca si in cazul integrarii Gauss Legendre, daca integrala este definita pe un

domeniu oarecare, dar finit, [ a, b ] , trebuie facuta schimbarea de variabila:
b−a b+a
z=
x+ (6.34)
2 2
Prin urmare, daca limita superioara de integrare este variabila, pentru un
interval [ a, t ] formula de cuadratura Clenshaw-Curtis devine:
t − a ⎛ a0 a1 n −1 (−1) r +1 n a −a ⎞
t

∫a ⎜ − + ∑ ar 2
f ( x)dx = + ∑ Tr (t ) r −1 r +1 ⎟ (6.35)
2 ⎝ 2 4 r =2 r − 1 r =1 2r ⎠
iar relatia de calcul pentru coeficienti va fi
2 n ⎛ (2k − 1)π / 2 t − a t + a ⎞ r (2k − 1)π / 2
ar = ∑ f ⎜ cos + ⎟ cos , r = 0,1,..., n (6.36)
n k =1 ⎝ n 2 2 ⎠ n

Daca limita superioara de integrare nu este neaparat variabila este

avantajos sa se faca schimbarea de variabila (6.34) catre intervalul de integrare
265H

[ −1,1] . In acest caz, tinand cont ca Tr (1) = 1 pentru orice r , pornind direct de la
relatia (6.31) si facand observatia ca termenii cu r impar devin nuli, se obtine
266H

urmatoarea formula, mult mai practica, de cuadratura Clenshaw-Curtis

⎛ n
−1 ⎞
b−a⎜
b 2
2 ⎟
∫a f ( x)dx = 2 ⎜⎜ a0 + ∑
r =1
ar
1 − (2 r ) 2 ⎟
⎟
(6.37)
⎝ ⎠
unde coeficientii ar sunt dati de relatia de transformare cosinus discreta:
2 n ⎛ (2k − 1)π / 2 b − a b + a ⎞ r (2k − 1)π / 2 n
ar = ∑ f ⎜ cos
n k =1 ⎝ n 2
+
2 ⎠
⎟ cos
n
, r = 0,1,..., − 1 (6.38)
2
Se obtine o dublare a vitezei de calcul datorita injumatatirii numarului de
coeficienti care trebuiesc calculati.
Programul 6.5 caluleaza integrala definita a unei functii (in exemplu functia
Runge) pe un interval finit oarecare, pe baza relatiilor (6.37) si (6.38).
267H 268H
 Programul 6.5. Calculul unei integrale definite prin cuadratura Clenshaw-
Curtis
f@x_D = ;H∗ Functia de integrat ∗L
1
1 + x2
a = −3.; b = 3.;H∗ Intervalul ∗L
Na = 21;
ForAr = 0, r ≤ Na, r += 2, H∗ se calculeaza doar coeficientii pari ∗L
H2 j − 1L r H2 j − 1L
A@rD = „ fACosA E E CosA EE;
Na π π
2. 2 b− a b + a 2
+
Na Na 2 2 Na
j=1

i Na
j y
z
j
j z
z
j
j ‚ A@2 rD + A@0Dz
z
−1

j z
b− a 2

j 1 − H2 rL z
2
j r=1 z
Icc = ;
k {
2 2

Deoarece nodurile sunt mai simplu de calculat decat in cazul cuadraturii

Gauss-Legendre, in care ele trebuiesc calculate ca radacini ale unei ecuatii
polinomiale si deci au deja anumite erori inerente, practic numarul de termeni care
trebuiesc considerati pentru obtinerea unei precizii date prin cuadratura Clenshaw-
Curtis este comparabil cu cel pentru cuadratura Gauss-Legendre. Erorile de calcul
pentru nodurile cuadraturii Gauss-Legendre crescand cu ordinul polinomului folosit,
exista un ordin peste care cuadratura Clenshaw-Curtis devine mai precisa, in special
in cazul unor functii dificil de integrat.
Printre avantajele cuadraturii Clenshaw-Curtis fata de Gauss-Legendre se
citeaza si posibilitatea de a se obtine transformata cosinus discreta (partea cea mai
cronofaga a algoritmului) printr-o transformare Fourier rapida (Fast Fourier
Transform- FFT) care reduce numarul de calcule de la n 2 la n log 2 n daca n este o
putere a lui 2. Astfel, de exemplu pentru n = 32 , volumul de calcule scade de peste 6
ori, iar pentru n = 64 de peste 10 ori.

Capitolul 7
Ecuatii diferentiale ordinare
 7.1 Metoda Euler de ordinul I

Sa consideram forma explicita a unei ecuatii diferentiale ordinare de ordinul I

y ' ( x) = f [ x, y ( x)] (7.1)
Daca pornim dintr-un punct initial x0 , in care valoarea lui y ( x) este cunoscuta
ca o conditie initiala y0 , se poate dezvolta y ( x) in serie Taylor la distanta δ x de
acest punct. Limitand-ne la primii doi termeni, avem
y ( x0 + δ x) = y ( x0 ) + δ x y ' ( x0 ) + 0 (δ x 2 ) (7.2)
Notand δ x cu h , (pasul de esantionare) si pe y ' ( x0 ) cu k1 , neglijand
termenii de ordin superior putem afla valoarea lui y ( x) in punctul x0 + h
y ( x0 + h) = y ( x0 ) + h k1 (7.3)
in care, conform relatiei (7.1)
269H

k1 = f ( xn , yn ) (7.4)

Figura 7.1 Metoda Euler de ordinul I pe un subinterval

Putem repeta procedeul plecand acum din punctul x0 + h si determinam

valoarea lui y ( x) intr-un punct x0 + 2h , si asa mai departe.
Notand in general x0 + nh = xn si y ( x0 + nh) = yn rezulta ca , pornind de la
conditia initiala cunoscuta y ( x0 ) = y0 , putem calcula succesiv valorile lui y ( x) intr-
un numar oarecare de puncte echidistante dintr-un interval [a, b] prin formule similare
cu (7.3) si (7.4):
270H 271H

yn +1 = yn + h k1 (7.5)
k1 = f ( x0 , y0 ) (7.6)
Pasul este dat de
 b−a
h= (7.7)
N
unde N este numarul de diviziuni ale intervalului, iar relatiile sunt cu atat mai exacte
cu cat acest numar este mai mare. Intr-adevar, conform relatiei (7.2), eroarea la
272H

⎛ (b − a ) 2 ⎞
fiecare interval este 0 (h 2 ) = 0 ⎜ 2 ⎟ si deoarece erorile se cumuleaza la fiecare
⎝ N ⎠
interval si avem N intervale, eroarea totala este
⎛ (b − a ) 2 ⎞ 1
ε t = N 0 (h 2 ) = 0 ⎜ ⎟ ∼ 0 (h) ∼ (7.8)
⎝ N ⎠ N

Figura 7.2 Metoda Euler de ordinul I pe intregul interval

Datorita faptului ca eroarea este proportionala cu puterea intai a pasului de

esantionare, deci destul de mare, metoda nu prezinta importanta practica, dar poate
conduce la generalizari cu caracteristici superioare, care vor fi prezentate in cele ce
urmeaza.

7.2 Metode Euler de ordinul II

Din figura 7.1 se observa ca eroarea de estimare a valorii lui y ( x) in punctul

urmator celui curent poate fi destul de mare daca derivata intai variaza rapid. De aceea
ar fi avantajos sa se ia in consideratie si derivata a doua a lui y ( x) , obtinand-se astfel
formule Euler de ordinul II. Deoarece derivata a doua nu se cunoaste si nici nu se
poate calcula direct deoarece chiar functia de derivat este necunoscuta, se va proceda
pe o cale indirecta.
Astfel, se observa ca daca am lua in considerare un punct intermediar xn + α h ,
cu α ∈ (0,1) , si daca s-ar folosi in relatia (7.5) derivata k2 in acest punct, s-ar putea
273H

ca pentru un α corespunzator sa se obtina o eroare de semn invers fata de cea data de

derivata in punctul xn , ca in figura 7.3. O medie ponderata a acestor doua derivate ar
putea da o panta
k = w1k1 + w2 k2 (7.9)
(cu ponderile w1 si w2 cuprinse intre 0 si 1) care sa conduca chiar in punctul corect,
y ( xn +1 ) ,anuland eroarea. Din pacate, nici panta k2 in punctul xn + α h , data conform
ecuatiei explicite (7.1) de f [ xn + α h, y ( xn + α h) ] , nu se cunoaste deoarece nu se
274H

cunoaste valoarea functiei y ( xn + α h) in acest punct. Am putea insa aproxima pe

y ( xn + α h) dupa formula de ordinul I
y ( xn + α h) ≈ yn + α hk1
Totusi, nici aceasta nu da o valoare exacta pentru y ( xn + α h) ; cea corecta s-ar
obtine daca am lua o alta abscisa, xn + β h , in care apare o noua necunoscuta, β , dar
asigura egalitatea y ( xn + α h) = yn + β hk1 . Rezumand, vom folosi urmatoarul sistem:
yn +1 = yn + h( w1k1 + w2 k2 ) (7.10)
k1 = f [ xn , yn ] (7.11)
k2 = f [ xn + α h, yn + β hk1 ] (7.12)
 Figura7.3 Pantele necesare pentru metoda Euler de ordinul II

Exista patru necunoscute, w1 , w2 , α si β , pentru aflarea carora sa punem

urmatoarea conditie: dezvoltarea in serie Taylor a lui y ( xn + h) in punctul xn
trebuie sa coincida cu relatia (7.10), in care pe k2 il dezvoltam de asemenea in
275H

serie in punctul xn .
Sa dezvoltam mai intai pe y :
h 2 ''
yn +1 = yn + hy ( xn ) +
'
y ( xn ) + 0 (h3 ) (7.13)
2!
Conform relatiilor (7.1) si (7.4), avem y ' ( xn ) = f ( xn , yn ) = k1 , iar derivata a
276H 277H

doua se poate scrie pornind de la aceleasi relatii:

∂f ( x, y ) ∂f ( x, y )
y '' ( xn ) = f ' ( xn , yn ) = + ⋅ y ' ( xn )
∂x x = xn ∂y y = yn
si notand derivatele partiale in raport cu x si y cu f xn respectiv f yn obtinem
y '' ( xn ) = f xn + f yn k1 (7.14)
Rezulta ca dezvoltarea in serie Taylor are forma
h2 h2
yn +1 = yn + hk1 + f xn + f yn k1 (7.15)
2 2
 Sa consideram acum relatia (7.10) in care de asemenea dezvoltam in serie
278H

Taylor pe k2 = f [ xn + α h, yn + β hk1 ] , conform relatiei (7.12), ca functie de doua

279H

variabile, in jurul punctului ( xn , yn )

∂f ( x, y ) ∂f ( x, y )
k2 = f [ xn , yn ] + α h + β hk1 + ... (7.16)
∂x x = xn ∂y y = yn
Limitandu-ne la acesti primi termeni, si tinand seama de definitiile precedenta
obtinem
k2 = k1 + α hf xn + β hk1 f yn
care introdusa in relatia (7.10) conduce la
280H

yn +1 = yn + h ( w1k1 + w2 k1 ) + h 2 w2α f xn + h 2 w2 β k1 f yn (7.17)

Identificand asa cum am propus relatiile (7.13) si (7.17) dupa puterile lui h
281H 282H

obtinem urmatoarele conditii:

⎧
⎪ w1 + w2 = 1
⎪
⎪ 1
⎨ w2α = (7.18)
⎪ 2
⎪ 1
⎪⎩ w2 β = 2
Avand numai trei ecuatii si patru necunoscute, sistemul este nedetermiant si va
trebui sa alegem o valoare pentru una dintre necunoscute, ceea ce va face insa
imposibila anularea erorii, realizand totusi o micsorare puternica a acesteia. Sunt mai
raspandite doua alegeri posibile, care vor fi prezentate in continuare.
1
Alegand ponderi egale, rezulta w1 = w2 = , de unde α = β = 1 , deci
2
algoritmul devine:
k +k
yn +1 = yn + h 1 2 (7.19)
2
k1 = f [ xn , yn ] (7.20)
k2 = f [ xn + h, yn + hk1 ] (7.21)
definind o metoda Euler de ordinul II pe care o putem numi metoda pantei de mijoc.
1 1
Alegand pe α = , rezulta w2 = 1 , w1 = 0 si β = , deci algoritmul devine:
2 2
yn +1 = yn + hk2 (7.22)
k1 = f [ xn , yn ] (7.23)
⎡ h h ⎤
k2 = f ⎢ xn + , yn + k1 ⎥ (7.24)
⎣ 2 2 ⎦

definind o metoda Euler de ordinul II pe care o putem numi metoda punctului de

mijoc.
Conform relatiei (7.13) termenii neglijati sunt de ordinul h3 , acestia dand
283H

b−a
eroarea pe fiecare interval, astfel ca eroarea totala pe N = intervale va fi:
h
 1
ε t ∼ h2 ∼ (7.25)
N2
deci mult mai mica (de N ori) decat la metoda de ordinul I.

Clear@"`∗"D;
H∗ Problema Cauchy: y'@xD=f@x,yD ∗L
f@x_, y_D = −y Hx + 1L − Cos@xD;
yi = 1; H∗ Conditia initiala ∗L
xi = 0.2; xf = 4.; H∗ Intervalul ∗L
Np = 100; H∗ Numar de puncte de evaluare ∗L
;H∗ Pasul∗L
xf − xi
h=
Np
y@0D = yi;
Y = 88xi, yi<<; H∗ Lista de interpolare a solutiei ∗L
H∗ membrul drept al ecuatiei ∗L
ForAm = 0, m ≤ Np − 1, m ++,
xm = xi + m h; H∗ Punctul curent pe abscisa ∗L
K1 = f@xm, y@mDD;
K2 = f@xm + h, y@mD + h K1D;
y@m + 1D = y@mD + HK1 + K2L; H∗ Formula Euler cu panta medie ∗L
h

AppendTo@Y, 8xm + h, y@ m + 1D<D;

2

E;
g1 = ListPlot@Y, PlotJoined → True, PlotStyle → BlueD;
Print@YD;

H∗ Verificarea, prin comparatie cu solutia analitica, daca aceasta exista ∗L

sol@x_D = DSolve@8y'@xD −y@xD Hx + 1L − Cos@xD, y@xiD yi<, y@xD, xD@@1, 1, 2DD
H∗ Solutia analitica ∗L
g2 = Plot@sol@xD, 8x, xi, xf<, PlotStyle → BlackD;
Show@g1, g2D;
erabs = 8<;
For@n = 1, n ≤ Np, n++, AppendTo@erabs, 8Y@@n, 2DD − N@sol@xi + Hn − 1L hDD<D;D;
Print@Max@Abs@erabsDDD;
 H∗ Problema: y'@xD=−y@xD Hx+1L−Cos@xD ∗L
f@x_, y_D = −y Hx + 1L − Cos@xD;
yi = 1; H∗ Conditia initiala ∗L
xi = 0.2; xf = 4.; H∗ Intervalul ∗L
Np = 100; H∗ Numar de puncte de evaluare ∗L
; H∗ Pasul∗L
xf − xi
h=
Np
y@0D = yi;
Y = 88xi, yi<<; H∗ Lista de interpolare a solutiei ∗L
ForAm = 0, m ≤ Np − 1, m ++,
xm = xi + m h; H∗ Punctul curent pe abscisa ∗L
K1 = f@xm, y@mDD;
K2 = fAxm + , y@mD + K1E;
h h

y@m + 1D = y@mD + h H K2L; H∗ Formula Euler cu punct mediu ∗L

2 2

AppendTo@Y, 8xm + h, y@ m + 1D<D;

E;
g3 = ListPlot@Y, PlotJoined → True, PlotStyle → BlueD;
Show@g2, g3D;
erabs = 8<;
For@n = 1, n ≤ Np, n++, AppendTo@erabs, 8Y@@n, 2DD − N@sol@xi + Hn − 1L hDD<D;D;
Print@Max@Abs@erabsDDD;

7.3 Metode Runge-Kutta

Metoda prezentata anterior si pe care am denumit-o Euler de ordinul 2 mai

este denumita si metoda Runge-Kutta de ordinul 2. De fapt ea este o generalizare a
metodei Euler si in acelasi timp un caz particular al domeniului mult mai vast al
metodelor de tip Runge Kutta.
Ideea fundamentala a metodei Euler de ordinul 2, dupa cum s-a vazut, are
doua componente:
a) Impartirea unui interval in subintervale, aproximarea pantei cautate printr-o
medie ponderata a mai multor pante. Astfel daca se folosesc numai doua
pante, se obtine o metoda numita de ordinul 2, cu panta data de relatia (7.9) 284H

k = w1k1 + w2 k2 .
Daca s-ar folosi mai multe pante, s-ar obtine metode de ordin superior, m ,
printr-o generalizare de tipul:
m
k = ∑ wi ki (7.26)
i =1
b) Fiecare panta se determina din cea precedenta, fiind o aproximare din ce in ce
mai buna a pantei corecte. Astfel, la metoda de ordinul 2, prima panta se
determina direct conform metodei Euler de ordinul 1, prin relatia (7.4) 285H

k1 = f ( xn , yn ) , iar a doua panta se determina pe baza acesteia, prin relatia

(7.12) k2 = f [ xn + α h, yn + β hk1 ] .
286H
 Generalizand, o panta de ordinul m se poate obtine din panta de ordin
m − 1 printr-o relatie de tipul:

⎡ i −1 ⎤
ki = f ⎢ xn + α i h, yn + ∑ β ij hki −1 ⎥ (7.27)
⎣ j =1 ⎦
unde i = 1, m, j = 1, i − 1, α1 = β10 = 0 .
In functie de valoarea lui m se obtin metode de diverse ordine, denumite
metode de tip Runge-Kutta. Desigur, pot exista multe metode de acest tip, dar in
practica s-au impus pe scara larga metodele de ordinul 4, abreviate de obicei ca
RK4.
Trebuie observat ca pe masura ce ordinul creste, dezvoltarile in serie Taylor
(7.13) si (7.17) folosite in paragraful precedent pentru gasirea
287H 288H

coeficientilor α , β , w si k trebuie sa cuprinda termeni pana la ordinul m , astfel ca

expresiile finale de tipul (7.17) devin din ce in ce mai complexe. De exemplu,
289H

pentru deducerea metodei RK4, tinand seama de expresiile derivatelor functiilor

compuse pana la ordinul 4, expresiile finale au aproximativ dimensiunea unei
pagini. Acest calcul poate fi facut prin algebra computerizata folosind programul
Mathematica, dar chiar si prelucrarea si interpretarea rezultatului este o operatiune
laborioasa; efortul pe care ar trebui sa il depuna cursantul pentru parcurgerea
acestei demonstaratii este prea mare in comparatie cu beneficiul pedagogic, astfel
ca ea nu va fi prezentata amanuntit in acest material.
Pe de alta parte, numarul de necunoscute creste din ce in ce mai mult in raport
cu numarul de ecuatii pe masura ce ordinul creste, astfel ca sunt necesare alegeri
arbitrare pentru multe dintre acestea, generand pentru fiecare ordin in parte o clasa
de metode RK.
Deoarece precizia obtinuta nu mai creste liniar peste ordinul 4, asa cum am
mentionat acesta este ordinul cel mai des utilizat, iar alegerea parametrilor se face
utilizand in general valori de tip α i = 1/ 2 . Rezulta printr-un procedeu analog cu
cel prezentat in paragraful precedent, urmatoarea schema larg raspandita de tip
RK4, bazata pe astfel de alegeri:
h
yn +1 = yn + (k1 + 2k2 + 2k3 + k4 )
6
k1 = f ( xn , yn )
⎛ h h ⎞
k2 = f ⎜ xn + , yn + k1 ⎟ (7.28)
⎝ 2 2 ⎠
⎛ h h ⎞
k3 = f ⎜ xn + , yn + k2 ⎟
⎝ 2 2 ⎠
k4 = f ( xn + h, yn + hk3 )
Clear@"`∗"D;
H∗ y'@xD=−y@xD Hx+1L−Cos@xD; ODE sub forma explicita ∗L
f = −y Hx + 1L − Cos@xD; H∗ Forma explicita fara argumente ∗L
y@0D = 1; H∗ Conditia initiala ∗L
xM = 4.; x0 = 0; H∗ Intervalul ∗L

Np = 40; H∗ Numarul de puncte ales pe interval ∗L

; H∗ Pasul ∗L
xM − x0
h=
Np
Y = 880, y@0D<<; H∗ Punctul initial al graficului solutiei ∗L

ForAm = 0, m ≤ Np − 1, m ++, H∗ algoritmul RK4 in fiecare punct ∗L

xm = m h;
K1 = f ê. 8x → xm, y → y@ mD<;
K2 = f ê. 9x → xm + , y → y@ mD + K1=;
h h
2 2
K3 = f ê. 9x → xm + , y → y@ mD + K2=;
h h

K4 = f ê. 8x → xm + h, y → y@ mD + h K3<;
2 2

y@ m + 1D = y@ mD + H K1 + 2 K2 + 2 K3 + K4L;
h

AppendTo@ Y, 8xm + h, y@ m + 1D<D; E;

6

H∗ Y contine lista cu punctele solutiei ∗L

g1 = ListPlot@ Y, PlotJoined → FalseD; H∗ reprezentarea grafica a solutiei ∗L

H∗ Se interpoleaza lista construind polinomul de interpolare Lagrange ∗L

ii−1
j x − Y@@j, 1DD yj
z i Np+1 x − Y@@j, 1DD y
z
p@x_D := „ Y@@i, 2DD j
j
j‰
z
z j
j‰ z
z
Np+1

kj=1
z
Y@@i, 1DD − Y@@j, 1DD { j
kj=i+1 Y@@i, 1 DD − Y@@j, 1DD z;
{
i=1

g2 = Plot@ Evaluate@ p@xDD, 8x, x0, xM<, PlotStyle → RedD;

Show@g1, g2D;
Print@"PHxL=", N@ Expand@p@xDDDD;

H∗ In continuare se calculeaza eroarea maxima ∗L

Clear@x, yD;
H∗ Aceasta ecuatie are solutie analitica exacta ∗L
sol = DSolve@8y'@xD −y@xD Hx + 1L − Cos@xD, y@0D 1<, y@xD, xD@@1, 1, 2DD;
nl = 10 Np; H∗ Se esantioneaza intr−un numar suficient de puncte ∗L
; H∗ Pas mai fin de esantionare ∗L
xM − x0
h1 =
ert = 8<; H∗ Lista erorilor ∗L
nl

z@x_D = sol;
For@n = 1, n ≤ nl, n++, AppendTo@ert, 8 p@x0 + h1 nD − z@x0 + h1 nD<D;D;
Print@Max@Abs@ertDDD;
H∗ Se afiseaza eroarea absoluta maxima ∗L
 Metoda RK4 este inca dominanta in rezolvarea ecuatiilor diferentiale cu
conditii initiale, fiind uneori combinata cu alte metode in cazul problemelor mai
dificile sau care solicita o precizie foarte mare.

7.4 Algoritmi adaptivi

Pentru imbunatatirea performatelor metodelor directe de rezolvare a ecuatiilor

diferentiale prezentate anterior se recomanda utilizarea acestora in conjunctie cu
algoritmi adaptivi, cu pas variabil. Aceste procedee tin seama de compromisul care
trebuie in general efectuat in metodele de aproximatie numerica la esantionarea unui
interval: daca pasul de esantionare este prea mare, viteza este de asemena mare, dar
precizia este mica si invers.
Principiul metodei consta in calculul initial al unui punct cu un pas relativ
mare, care ar asigura o viteza ridicata, apoi recalcularea punctului respectiv cu un pas
injumatatit. Dupa compararea celor doua rezultate se va lua o decizie: daca diferenta
este sub o anumita limita (de ordinul de marime al erorii admise), se trece la
urmatorul punct, iar daca diferenta este mai mare decat cea admisa, se injumatateste
din nou pasul si se reia calculul aceluiasi punct. In acest fel, pe portiunile cu variatii
lente, pasul va putea sa aiba valori suficient de mari pentru a nu se pierde timp in mod
inutil, iar pe portiunile critice pasul va fi ajustat automat pana la obtinerea preciziei
dorite. De obicei se limiteaza numarul de injumatatiri admise, deoarece este posibil pe
de o parte ca in unele cazuri sa nu se obtina convergenta, iar in altele calculul sa
dureze prea mult. De remarcat ca dupa o prima injumatatire, pentru a se ajunge in
acelasi punct xn +1 este necesar calculul in 2 puncte, dupa inca o injumatatire in 4
puncte, apoi in 8, s.a.m.d. ceea ce poate accentua problemele de viteza in cazul unor
functii neconvenabile.

7.5 Extrapolarea Richardson prin metoda Burlich-Stoer

Un alt procedeu puternic de imbunatatire a performantelor metodelor de

rezolvare a ecuatiilor diferentiale este dat de metoda Burlich-Stoer de extrapolare
Richardson. Ideea de baza a extrapolarii Richardson in acest caz este de a calcula
b−a
valoarea functiei in punctul urmator yn +1 ( xn + h1 ) folosind mai intai un pas h1 =
Np
(unde a, b si N p sunt respectiv capetele intervalului si numarul de puncte din
interval), si apoi repetarea calculului folosind valori hk din ce in ce mai mici pentru
pas. Rezultatele ynhk+1 obtinute intr-un anumit punct xn +1 pot forma perechi cu valorile
respective ale lui hk sub forma: { ynh1+1 , h1},{ ynh2+1 , h2 },...,{ ynhm+1 , hm }, iar aceste puncte pot
fi folosite pentru o interpolare, de exemplu polinomiala, care permite calculul valorii
functiei pentru h = 0 , (care ar corespunde valorii corecte). Valorile pasilor trebuie
alese in asa fel incat de fiecare data, pornind din xn sa se ajunga la acelasi xn +1 la care
se ajunsese cu pasul initial. Propunerea Burlich-Stoer pentru marirea numarului de
h h h h h h
subintervale este: hk = 1 , 1 , 1 , 1 , 1 , 1 ,... deci nu o simpla injumatatire care ar
2 4 6 8 12 16
creste prea rapid numarul de subintervale. De asemenea, Burlich si Stoer au aratat ca
in general este mai avantajoasa interpolarea cu functii rationale decat cu polinoame.
 Desi acest procedeu se poate folosi cu metode de orice ordin, de cele mai
multe ori el este asociat cu metode de ordinul 2, precizia obtinuta crescand rapid cu
numarul de subdivizari ale intervalelor.

12.4 Metode spectrale de rezolvare a ecuatiilor diferentiale ordinare