Metode munerice

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:

11 1 12 2 1 1

21 1 22 2 2 2

1 1 2 2

......

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

n n

n n

n n nn n n

a x a x a x ba x a x a x b

a x a x a x b

+ + + =⎧⎪ + + + =⎪⎨⎪⎪ + + + =⎩

(1.1)

sau concis:

1

, 1, 2,..n

ij j ij

a x b i n=

= =∑ (1.2)

Considerand trei matrici: A matricea coeficientilor variabilelor (patrata), X matricea necunoscutelor (coloana) si B matricea termenilor liberi (coloana), avand urmatoarele forme:

1 1111 12

2 221 22 2

1 2

.... .

; ;. .. .

n

n

n n nn

n n

x baa a

x ba a a

a a a

x b

⎛ ⎞ ⎛ ⎞⎛ ⎞ ⎜ ⎟ ⎜ ⎟⎜ ⎟ ⎜ ⎟ ⎜ ⎟⎜ ⎟ ⎜ ⎟ ⎜ ⎟⎜ ⎟= = =⎜ ⎟ ⎜ ⎟⎜ ⎟ ⎜ ⎟ ⎜ ⎟⎜ ⎟ ⎜ ⎟ ⎜ ⎟⎜ ⎟ ⎜ ⎟ ⎜ ⎟⎝ ⎠ ⎜ ⎟ ⎜ ⎟

⎝ ⎠ ⎝ ⎠

A X B

…

(1.3)

se observa ca in ecuatiile 0H(1.2) membrul stang este echivalent cu inmultirea 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 1−A , de unde rezulta direct solutiile: 1.−=X A 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:

0, pentru

; , 1, 2,..0, pentru

ij

ij

a i ji j n

a i j

= >⎧⎪ =⎨ ≠ ≤⎪⎩ (1.6)

si urmatoarea forma

111 12

22 20 ...

0 0

n

n

nn

aa aa a

a

⎛ ⎞⎜ ⎟⎜ ⎟⎜ ⎟=⎜ ⎟⎜ ⎟⎜ ⎟⎝ ⎠

A

…

(1.7)

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

1n − variabile sunt nuli, rezulta:

nn

nn

bxa

= (1.8)

In penultima ecuatie, numai ultimele doua necunoscute au coeficienti nenuli. Dar deoarece nx este deja cunoscut, se poate determina 1nx − :

1 1,1

1, 1

n n n nn

n n

b a xx

a− −

−− −

−= (1.9)

Continuind in acest mod, pentru ecuatia i , coeficientii primelor 1i − necunoscute sunt nuli iar necunoscutele 1 2, ,...,i i nx x x+ + au fost deja determinate. Ecuatia se scrie deci , 1 1 ...ii i i i i in n ia x a x a x b+ ++ + + = (1.10) de unde se obtine necunoscuta ix :

1

n

i ij jj i

iii

b a xx

a= +

−=

∑ (1.11)

Procedeul se va repeta pana la obtinerea tuturor necunoscutelor folosind formula generala 1H(1.11). 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 2H(1.1) si care in mod evident pastreaza validitatea 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 2m si sa o scadem din a doua, apoi sa inmultim prima linie cu o alta constanta 3m 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: 1 11 0i ia a m− = (1.12) de unde rezulta constanta cu care trebuie inmultita prima linie pentru anularea elementului din prima coloana a liniei i :

1

11

ii

ama

= (1.13)

Elementul 11a 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: '

1 ; 2,3..., ; 1, 2...,ij ij j ia a a m i n j 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.

111 12' '22 2' '32 3

' '2

0 ......0

0

n

n

n

n nn

aa a

a a

a a

a a

⎛ ⎞⎜ ⎟⎜ ⎟⎜ ⎟⎜ ⎟⎜ ⎟=⎜ ⎟⎜ ⎟⎜ ⎟⎜ ⎟⎜ ⎟⎝ ⎠

A

…

(1.14)

In continuare, pentru a forma a doua coloana, trebuiesc anulate elementele ' ' '32 42 2, ,..., na a a . Aceasta se va face in mod similar, dar folosindu-se ca pivot elementul

de pe diagonala celei de a doua linii, '22a . Conditia de anulare a unui element din linia

i a celei de a doua coloane este: ' '

2 22 0i ia a m− = (1.15) Constanta de multiplicare a celei de a doua linii este deci:

'2'22

ii

ama

= (1.16)

iar operatiile asupra celorlalte linii: '' ' '

2 ; 3, 4..., ; 2,3,...,ij ij j ia a a m i n j n= − = = (1.17) Este evident ca elementele nule din prima coloana raman nule deoarece din ele

se scade constanta inmultita cu '1 ja 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.

111 12' '22 2

'3

'

0 ......0 0

0 0

n

n

n

nn

aa a

a a

a

a

⎛ ⎞⎜ ⎟⎜ ⎟⎜ ⎟⎜ ⎟⎜ ⎟=⎜ ⎟⎜ ⎟⎜ ⎟⎜ ⎟⎜ ⎟⎝ ⎠

A

…

(1.18)

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 1n + coloane, iar operatiile descrise sa se faca pentru 1,2,..., 1j n= + . 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:

'

' ; 1, 2, 1iki

kk

am k na

= = − (1.19)

'' ' ' ; 1, 2,..., 1; 1, 2..., ; , 1..., 1ij ij kj ia a a m k n i k k n j k k n= − = − = + + = + + (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 ∗Lb= 82, −1, 1, 3<; H∗ Matricea termenilor liberi ∗Ln= Length@bD; H∗Dimensiunea matricii ∗La= Table@0, 8n<, 8n+ 1<D; H∗ Formarea matricii extinse ∗Lx= 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 ∗LPrint@MatrixForm@aDD;ForAk= 1, k≤ n−1, k++,

ForAi= k+ 1, i ≤ n, i++,

m =a@@i, kDDa@@k, kDD ; H∗ Constanta pentru linia k ∗L

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

E;

Print@MatrixForm@aDD;H∗ Matricea extinsa, superior triunghiulara∗LH∗ Rezolvarea prin substitutie inversa ∗L

x@@nDD =a@@n, n+ 1DD

a@@n, nDDForAj= n− 1, j ≥ 1, j−−,

x@@jDD =a@@j, n+ 1DD −⁄k=j+1

n a@@j, kDD x@@kDDa@@j, jDD ;

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

ai.xH∗ Verificarea solutiilor ∗Lb

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 3H(1.19) va da rezultat infinit, 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 ∗Lb= 82, −1, 1, 3<; H∗ Matricea termenilor liberi ∗Ln= Length@bD; H∗Dimensiunea matricii ∗La= Table@0, 8n<, 8n+ 1<D;x= Table@0, 8n<D;H∗ Formarea directa a matricii extinse ∗LFor@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 ∗LFor@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++,

m =a@@i, kDDa@@k, kDD ; H∗ Constanta pentru linia k ∗L

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

Print@MatrixForm@aDD; H∗ Matricea extinsa, superior triunghiulara∗LH∗ Rezolvarea prin substitutie inversa ∗L

x@@nDD =a@@n, n+ 1DD

a@@n, nDDForAj= n− 1, j ≥ 1, j−−,

x@@jDD =a@@j, n+ 1DD −⁄k=j+1

n a@@j, kDD x@@kDDa@@j, jDD ;

Print@x@@jDD D;E H∗ Afisarea solutiilor ∗Lai.xH∗ Verificarea solutiilor ∗Lb

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 4H(1.19) se va face cu erorile de rotunjire cele mai mici posibile.

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, 1−A , care inmultita cu A sa conduca la o matrice unitara I: -1A .A = I (1.22) unde

1 0 ... 00 1... 0.............0 0 ...1

⎛ ⎞⎜ ⎟⎜ ⎟=⎜ ⎟⎜ ⎟⎝ ⎠

I (1.23)

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 5H(1.4) se formeaza mai intai o matrice extinsa de dimensiuni (2 1)n n× + , prin reunirea coloanelor matricii coeficientilor, a matricii termenilor liberi si a coloanelor unei matrici unitare, e =A 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 1.− =A A I , 1.− =A B X si 1 1.− −=A I A , obtinem: ( )1 1 1. .e

− − −= =A A A A B I I X A∪ ∪ ∪ ∪ (1.25) Se observa ca se obtine in primele n coloane o matrice unitara si in plus, ca

“produse secundare” se obtin: in coloana 1n + 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 ∗LAE1= Table@0, 8n<, 8n+ 1<D;For@i= 1, i≤ n, i++, AE1@@iDD = Flatten@Join@A@@iDD, 8b@@iDD<DDD;

H∗Formarea matricei extinse totale ∗La= Transpose@Join@Transpose@AE1D, IdentityMatrix@nDDD;Print@MatrixForm@aDD;

H∗Permutarea liniilor matricei extinse partiale ∗Lp= 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, cDD = tD DD;DPrint@pD;H∗ numarul de permutari de linii ∗L;Print@MatrixForm@aDD;

d = H−1Lp; H∗Valoare initiala a determinantului, daca este necesar∗LForAi= 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=", dDx= 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;

A.Aiêê MatrixForm H∗ Verificarea matricii inverse obtinute ∗L

A.x H∗ Verificarea solutiilor ∗Lb

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:

'

' ; 1, 2, ;iki

kk

am k n k ia

= = ≠ (1.26)

'' ' ' ; 1, 2,..., ; 1, 2..., ; ; 1, 2..., 2 1ij ij kj ia a a m i n k n k i j n= − = = ≠ = + (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)

111 12

22 2 12

1 2

1 0 00 ... 1 ... 0

;0 0 1

n

n

nn n n

uu uu u l

u l l

⎛ ⎞ ⎛ ⎞⎜ ⎟ ⎜ ⎟⎜ ⎟ ⎜ ⎟⎜ ⎟ ⎜ ⎟= =⎜ ⎟ ⎜ ⎟⎜ ⎟ ⎜ ⎟

⎜ ⎟⎜ ⎟ ⎝ ⎠⎝ ⎠

U L

… …

(1.29)

Elementele acestora au deci proprietatile:

0, pentru0, pentru

; 1 ; , 1, 2,..0, pentru

0, pentru

ijij

iiij

ij

l i ju i j

l i j nu i j

l i j

= >⎧= >⎧ ⎪⎪ = =⎨ ⎨≠ ≤⎪ ⎪⎩ ≠ <⎩

(1.30)

Explicitand relatia 6H(1.28) pe elemente sub forma:

1

; , 1, 2,...,n

ij ik kjk

a l u i j n=

= =∑ (1.31)

si tinand cont de proprietatile 7H(1.30), rezulta relatiile:

1

1 1; , 1, 2,...,

j j

ij ik kj ik kj ij jjk k

a l u l u l u i j n−

= =

= = + =∑ ∑ (1.32)

1

1 1; , 1, 2,...,

i i

ij ik kj ik kj ijk k

a l u l u u i j n−

= =

= = + =∑ ∑ (1.33)

In prima relatie s-a tinut cont ca primul indice al lui kju , 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 ikl , care este chiar indice de sumare, poate fi doar mai mic sau egal cu primul si in plus 1iil = .

Rezulta urmatoarele relatii care permit determinarea succesiva a tuturor elementelor matricilor L si U:

1

1 ; , 1, 2,...,

j

ij ik kjk

ijjj

a l ul i j n

u

−

=

−= =

∑ (1.34)

1

1; , 1, 2,...,

i

ij ij ik kjk

u a l u i j n−

=

= − =∑ (1.35)

Astfel, pentru 1i = se obtin:

11 11 11

12 12

1 1

1;

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

n n

l u au a

u a

= ==

=

(1.36)

Pentru 2i = obtinem:

2121 22 22 21 12

11

22 23 22 21 12

2 2 21 1

;

1;.............

n n n

al u a l uu

l u a l u

u a l u

= = −

= = −

= −

(1.37)

Pentru 3i = obtinem:

3131 33 33 31 13 32 23

11

32 11 1222 34 34 31 13 32 24

22

33

3 3 31 1 32 2

;

;

1;.............

n n n n

al u a l u l uua l ul u a l u l u

ul

u a l u l u

= = − −

−= = − −

=

= − −

(1.38)

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∗LB= 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 ∗LForAk= 1, k≤ n, k++,

LPk,kT = 1; H∗ Elementele de pe diagonala lui L sunt 1 ∗LForAj= k, j≤ n, j++,

UPk,jT = APk,jT − ‚m=1

k−1LPk,mT UPm,jT; H∗ Pentru linia k se calculeaza elemntele lui U ∗L

H∗ pentru k=1, rezulta U@@1,jDD=A@@1,jDD ∗LH∗ U@@m,jDD se cunosc, deoarece suma este cu m<k deci s−au calculat anterior∗LH∗ L@@k,mDD se calculasera anterior din bucla urmatoare ∗L

ForAi= k+ 1, i ≤ n, i++,

LPi,kT =i

kjjjjAPi,kT − ‚

m=1

k−1LPi,mT UPm,kT

y

{zzzz ì UPk,kT;H∗ Pentru linia i, se calculeaza elementele lui L ∗L

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

Print@MatrixForm@LDD;Print@MatrixForm@UDD; MatrixForm@A− L.UD H∗ Verificari factorizare ∗LMax@Abs@A− L.UDD Y= Table@0, 8n<D;H∗ Forward Substitution pentru a afla Y ∗L

YPnT =BPnT

LPn,nT;

ForAi= 1, i≤ n, i++,

YPiT =i

kjjjjBPiT −‚

j=1

i−1LPi,jT YPjT

y

{zzzz ì LPi,iT;E;

H∗ Back substitution pentru a afla X ∗LX= Table@0, 8n<D;

XPnT =YPnT

UPn,nT;

ForAi= n− 1, i ≥ 1, i−−,

XPiT =i

kjjjjYPiT − ‚

j=i+1

nUPi,jT XPjT

y

{zzzz ì UPi,iT;E;

Print@MatrixForm@XDD;A.XH∗ Verificari solutie ∗LBMax@[email protected]− BDD

Pentru rezolvarea unui sistem de ecuatii liniare de tipul 8H(1.4) se scrie: ⇒ ⇒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 inversa.Matricile 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 1112

2 12 22

1 2

01 0... ... 00 1

;0 0 1

n

n

n n nn

u luu l l

l l l

⎛ ⎞ ⎛ ⎞⎜ ⎟ ⎜ ⎟⎜ ⎟ ⎜ ⎟⎜ ⎟ ⎜ ⎟= =⎜ ⎟ ⎜ ⎟⎜ ⎟ ⎜ ⎟

⎜ ⎟⎜ ⎟ ⎝ ⎠⎝ ⎠

U L

… …

(1.41)

Elementele celor doua matrici au deci proprietatile

,

1;0,

kjkk kj

u k ju u

k j

<⎧= = ⎨

>⎩ (1.42)

,

0,ik

ik

l k il

k i<⎧

= ⎨ >⎩ (1.43)

Conform relatiilor 9H(1.40)-10H(1.43), separand termenul diagonal 11 1u = , se poate scrie

1

1 1 1

i i

ii ij ji ij ji ii iij j

a l u l u l u−

= =

= = +∑ ∑ (1.44)

de unde se poate obtine termenul diagonal al matricei U

1

1

i

ii ii ij jij

l a l u−

=

= −∑ (1.45)

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

1

1 1

i i

ij ik kj ik kj ii ijk k

a l u l u l u−

= =

= = +∑ ∑ (1.46)

de unde rezulta elementul omolog al matricii U

1

1

i

ij ik kjk

ijii

a l uu

l

−

=

−=

∑ (1.47)

Pentru acelasi element al matricii A, tinand cont de relatia 12H(1.42)

1

1 1

j j

ij ik kj ik kj ij jjk k

a l u l u l u−

= =

= = +∑ ∑ (1.48)

Deoarece 1jju = , putem afla si elementul corespunzator al matricii L prin relatia

1

1

j

ij ij ik kjk

l a l u−

=

= −∑ (1.49)

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++,

LPk,kT = APk,kT − ‚m=1

k−1LPk,mT UPm,kT;

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

UPk,jT =i

kjjjjAPk,jT − ‚

m=1

k−1LPk,mT UPm,jT

y

{zzzz ì LPk,kT;

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++,

LPi,kT =i

kjjjjAPi,kT − ‚

m=1

k−1LPi,mT UPm,kT

y

{zzzz;

E;E;E

Print@MatrixForm@LDD;Print@MatrixForm@UDD; MatrixForm@A− L.UDMax@Abs@A− L.UDD Y= Table@0, 8n<D;

YPnT =BPnT

LPn,nT;


YPiT =i

kjjjjBPiT −‚

j=1

i−1LPi,jT YPjT

y

{zzzz ì LPi,iT;E;

X= Table@0, 8n<D;

XPnT =YPnT

UPn,nT;

ForAi= n− 1, i ≥ 1, i−−,

XPiT =i

kjjjjYPiT − ‚

j=i+1

nUPi,jT XPjT

y

{zzzz ì UPi,iT;E

Print@MatrixForm@XDD;Max@[email protected]− 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) ij jiu l= (1.51) Din relatiile 13H(1.50) si 14H(1.51), pentru un element diagonal, prin separarea ultimului termen se deduce

1

2 2 2

1 1 1

n k k

kk km mk km km kkm m m

a l u l l l−

= = =

= = = +∑ ∑ ∑ (1.52)

de unde rezulta elementele de pe diagonale

1

2

1

k

kk kk kmm

l a l−

=

= −∑ (1.53)

De asemenea pentru un element nediagonal se poate scrie

1

1 1 1

k k k

ik im mk im km im km ik kkm m m

a l u l l l l l l−

= = =

= = ⋅ = +∑ ∑ ∑ (1.54)

de unde rezulta elementele nediagonale ale matricilor triunghiulare

1

1

k

ik im kmm

ikkk

a l ll

l

−

=

−=

∑ (1.55)

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 ∗LH∗ Se construieste matricea ca matrice simetrica ∗LFor@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 ∗LForAk= 1, k≤ n, k++,

LPk,kT =$%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%APk,kT − ⁄m=1

k−1HLPk,mTL2;

ForAi= k+ 1, i ≤ n, i++,

LPi,kT =i

kjjjjAPi,kT − ‚

m=1

k−1LPi,mT LPk,mT

y

{zzzz ì LPk,kT;

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

YPnT =BPnT

LPn,nT;


YPiT =i

kjjjjBPiT −‚

j=1

i−1LPi,jT YPjT

y

{zzzz ì LPi,iT;E;

H∗ Back substitution ∗LX= Table@0, 8n<D;

XPnT =YPnT

UPn,nT;

ForAi= n− 1, i ≥ 1, i−−,

XPiT =i

kjjjjYPiT − ‚

j=i+1

nUPi,jT XPjT

y

{zzzz ì UPi,iT;E

Print@MatrixForm@XDD;Max@[email protected]− BDDH∗ 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

,

0,

,

0,

ij ij ij

ijij

ijij

d a

a i jl

i j

a i ju

i j

δ=

>⎧= ⎨

<⎩<⎧

= ⎨>⎩

(1.58)

Tinand seama de relatiile 15H(1.56) si 16H(1.57) se poate scrie urmatoarea ecuatie 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 1k + in functie de matricea necunoscutelor la iteratia anterioara k 1 -k k+( ) -1 ( ) -1X = D .(L+U).X +D .B (1.60) Relatia 17H(1.59) scrisa pe componente este

1

( 1) ( ) ( )

1 1 1

n i nk k k

ij j i ij j ij jj j j i

d x b a x a x−

+

= = = +

= − −∑ ∑ ∑ (1.61)

Tinand cont de definitia elementelor matricii diagonale ij ij ijd a δ= (1.62) rezulta ecuatia metodei iterative Jacobi 18H(1.60) scrisa pe componente

1

( 1) ( ) ( )

1 1

1 i nk k k

i i ij j ij jj j iii

x b a x a xa

−+

= = +

⎛ ⎞= − −⎜ ⎟

⎝ ⎠∑ ∑ (1.63)

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 19H(1.56) si 20H(1.57) se poate scrie urmatoarea ecuatie 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 1k + in functie de matricea necunoscutelor la iteratia anterioara k ( 1)k k+ = +-1 ( ) -1X -(D+L) .U.X (D+L) .B (1.65) Relatia 21H(1.64) scrisa pe componente este

( 1) ( )

1 1

( ) ( )i n

k kij j i ij j

j j i

x b x+

= = +

= −∑ ∑D+L U (1.66)

unde s-a tinut cont ca ( ) 0ij =D+L pentru j i> si ( )ijU =0 pentru j i≤ . Aceasta se poate scrie si sub forma

1

( 1) ( 1) ( )

1 1

i nk k k

ii i ij j i ij jj j i

a x a x b a x−

+ +

= = +

+ = −∑ ∑ (1.67)

unde s-a tinut cont de definitiile22H(1.58) ale elementelor matricii diagonale si ale matricilor triunghilare.

Rezulta ecuatia metodei iterative Gauss-Siedel scrisa pe componente

1

( 1) ( 1) ( )

1 1

1 i nk k k

i i ij j ij jj j iii

x b a x a xa

−+ +

= = +

⎛ ⎞= − −⎜ ⎟

⎝ ⎠∑ ∑ (1.68)

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 ( 1)k

jx + si nu pe cele de la iteratia anterioara ( )kjx .

Acest lucru este posibil, deoarece daca in fiecare iteratie valorile necunoscutelor se calculeaza pe rand cu indice crescator, la calculul necunoscutei ( 1)k

ix + sunt cunoscute deja valorile ( 1)

0kx + , ( 1)

1kx + ,…, ( 1)

1k

ix +− . 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 ∗LX= Xo;H∗ Noile valori produse dupa iteratie ∗LNk= 2500; H∗ Numar maxim de iteratii ∗Lk= 0;WhileAk≤ Nk,


XPiT =i

kjjjjBPiT −‚

j=1

i−1APi,jT XPjT − ‚

j=i+1

nAPi,jT XoPjT

y

{zzzz ì APi,iT;

E; H∗ Gasirea noilor aproximatii ∗LXo= X; H∗ Pentru iteratia urmatoare ele devin valori vechi ∗Lk++E;


O imbunatatire in continuare a vitezei de convergenta se poate obtine folosind metoda suprarelaxarii. Daca in relatia 23H(1.68) se aduna si se scade ( )k

ix , incluzind termenul scazut in ultima suma se obtine relatia echivalenta ( )k

ix

1

( 1) ( ) ( 1) ( ) ( ) ( )

1

1 i nk k k k k k

i i i ij j ij j i ij j iii

x x b a x a x x Ra

−+ +

= =

⎛ ⎞= + − − = +⎜ ⎟

⎝ ⎠∑ ∑ (1.69)

Semnificatia acestei relatii este ca la valoarea din iteratia precedenta se aduna un corector dat de restul membrului drept, ( )k

iR . 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 ( 1) ( ) ( ) ; [0, 2]k k k

i i ix x wR w+ = + ∈ (1.70) sau desfasurat

1

( 1) ( ) ( 1) ( )

1

i nk k k k

i i i ij j ij jj j iii

wx x b a x a xa

−+ +

= =

⎛ ⎞= + − −⎜ ⎟

⎝ ⎠∑ ∑ (1.71)

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

1.633w = 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 ∗LX= Xo; H∗ Noile valori produse dupa iteratie ∗LNk= 2000; H∗ Numar maxim de iteratii ∗Lk= 0;WhileAk≤ Nk,


XPiT = XoPiT + 1.633 i

kjjjjBPiT −‚

j=1

i−1APi,jT XPjT −‚

j=i

nAPi,jT XoPjT

y

{zzzz ì APi,iT;

E; H∗ Gasirea noilor aproximatii ∗LXo= X; H∗ Pentru iteratia urmatoare ele devin valori vechi ∗Lk++E;


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

1 1 1 1

1 2 2 2 2

1

0............0............0

...............................0 0 0... n n n n

d c x ba d c x b

a d x b−

⎛ ⎞ ⎛ ⎞ ⎛ ⎞⎜ ⎟ ⎜ ⎟ ⎜ ⎟⎜ ⎟ ⎜ ⎟ ⎜ ⎟⋅ =⎜ ⎟ ⎜ ⎟ ⎜ ⎟⎜ ⎟ ⎜ ⎟ ⎜ ⎟⎝ ⎠ ⎝ ⎠ ⎝ ⎠

(1.72)

Se observa ca matricea coeficientilor sistemului este

1, 1

1 11 , 1

............................................0...... ................00............... .....0.............................................

ii iii i i

i i ii i i

i i i

dd aa d c

A c aa d c

a a

−+

+ +− −

⎛ ⎞ ⎧ =⎜ ⎟ ⎪⎜ ⎟= =⎨⎜ ⎟ ⎪⎜ ⎟ =⎩⎝ ⎠

1 1, 1

1,

1 1, 2

i i

i i i

i i i

aa ac a

+ + +

+

+ + +

⎧ =⎪

=⎨⎪ =⎩

(1.73)

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

'

10i i i ia a m d+= = − (1.74) rezulta multiplicator pentru linia 1i +

1i

ii

amd+ = (1.75)

unde id este elementul pivot. Tinand cont ca operatia trebuie facuta pe intraga linie 1i + '

1 1 1i i i id d m c+ + += − ⋅ (1.76) rezulta ca celelalte elemente care se modifica pe linia 1i + sunt (1.77)

'1 1

ii i i

i

ad d cd+ += − (1.78)

'1 1

ii i i

i

ab b bd+ += − (1.79)

In aceste relatii avem 1, 2,..., 1i n= − . Matricea devine si mai rara, si in plus superior triunghiulara

'1'2

''

''1 1

'

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

............0 ..........' ;

............0 0

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

i i

ii i

n

b

bd c

A Bbd c

b

+ +

⎛ ⎞⎜ ⎟

⎛ ⎞ ⎜ ⎟⎜ ⎟ ⎜ ⎟⎜ ⎟ ⎜ ⎟= =⎜ ⎟ ⎜ ⎟⎜ ⎟ ⎜ ⎟⎜ ⎟⎝ ⎠ ⎜ ⎟

⎜ ⎟⎝ ⎠

(1.80)

Prin urmare necunoscutele se pot determina din relatii cu numai trei termeni

' '1i i i i id x c x b+⋅ + ⋅ = (1.81)

Deci si substitutia inversa se realizeaza mult mai rapid deoarece in relatiile finale nu mai apar sume ca in cazul matricilor dense

'

1'

i i ii

i

b c xxd

+−= (1.82)

unde 1, 2,...,1i n n= − − . Evident, ultima necunoscuta este data de

'

'n

nn

bxd

= (1.83)

Programul 1.9 rezolva un sistem de ecuatii liniare cu matrice tridiagonala prin metoda Gauss folosind relatiile 24H(1.78), 25H(1.79), 26H(1.82) si 27H(1.83), unde s-au notat cu

[ ]d i elementele diagonale ale matricii superior triunghiulare obtinute. Programul 1.9 Rezolvarea unui sistem de ecuatii liniare cu matrice tridiagonala prin metoda Gauss

A=

i

k

jjjjjjjjjjjjj

2 5 0 0 03 4 2 0 00 1 5 2 00 0 3 3 10 0 0 3 4

y

{

zzzzzzzzzzzzz; H∗ Matricea tridiagonala ∗L

b= 83, 2, 1, 4, 5<; H∗ Matricea termenilor liberi ∗LB= b; H∗ Copie pentru verificare ∗Ln= 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

d@i+ 1D = A@@i+ 1, i+ 1DD−A@@i+1, iDD

d@iD A@@i, i+ 1DD;

b@@i+ 1DD = B@@i+ 1DD −A@@i+ 1, iDD

d@iD b@@iDD;

E;

x@@nDD =b@@nDDd@nD ;

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

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

d@iD ;

E;

x H∗ Afisarea solutiilor ∗LA.x− BH∗ 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

11 12

21 22 23

1 1

0..............01 0............0 01............0 0 0 ............0

;

0 0............ 1 0 0 0..............n n nn

u ul u u

l u− −

⎛ ⎞⎛ ⎞⎜ ⎟⎜ ⎟⎜ ⎟⎜ ⎟= =⎜ ⎟⎜ ⎟⎜ ⎟⎜ ⎟

⎝ ⎠ ⎝ ⎠

L U (1.86)

si deci elementele lor sunt

/ . , 1 / . , 1

;0 / . , 1 0 / . , 1ij ij

ij ij

l i j i j u i j i jl u

i j i j i j i j

= = + = = −⎧ ⎧= =⎨ ⎨

< > + > < −⎩ ⎩ (1.87)

Din relatia 28H(1.85) rezulta pentru elementele de deasupra diagonalei

, 1 , 1 , 11 1

,n

k k km m k kk k km

a l u l u+ + +=

= =∑ (1.88)

unde s-a tinut cont de 29H(1.87) din care rezulta ca { 1, } { 1, } { }m k k k k k∈ − ∩ + = . Prin urmare, elementele nediagonale ale matricii U sunt chiar elementele nediagonale ale matricii initiale A. , 1 , 1k k k ku a+ += (1.89) Tot din relatia 30H(1.85) rezulta pentru elementele de dedesubtul diagonalei

, 1 , 1 , 1 1, 11

n

k k km m k k k k km

a l u l u− − − − −=

= =∑ (1.90)

unde s-a tinut cont de 31H(1.87) din care rezulta ca { 1, } { 2, 1} { 1}m k k k k k∈ − ∩ − − = − . Prin urmare, elementele nediagonale ale matricii L sunt

, 1, 1

1, 1

k kk k

k k

al

u−

−− −

= (1.91)

Din relatia 32H(1.85) se observa imediat ca primul element diagonal al matricii U este egal cu primul element diagonal al matricii A. 11 11a u= ⇒ 11 11u a= (1.92) Din relatia pentru elementele diagonale

, 1 1,1 1

n

kk km mk k k k k kk kkm

a l u l u l u− −=

= = +∑ (1.93)

unde s-a tinut cont de 33H(1.87) din care rezulta { 1, } { 1, } { 1, }m k k k k k k∈ − − = −∩ , rezulta celelalte elemente diagonale

1,

, 1 1,

k k

kk kk k k k k

a

u a l u−

− −= − ⋅ (1.94)

si deci , 1 1,kk kk k k k ku a l a− −= − ⋅ (1.95) Relatiile de factorizare 34H(1.89),35H(1.91),36H(1.92) si 37H(1.95) implica mult mai putine 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 ( ) ( ) ( ){ }1 2, ,..., nY f x f x f x= , pentru anumite valori

(alese sau impuse) ale argumentului { }1 2, ,..., nX x x x= , ca in figura 2.1. Daca dorim sa determinam prin calcul valoarea functiei intr-un punct ix care nu apartine multimii de valori experimentale, ix 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.

x1

x

f (x )a

f(x )

f(x)

i

i

i

Functie aproximataFunctia exacta

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

-0.75 -0.5 -0.25 0.25 0.5 0.75

0.9

0.95

1.05

1.1

1.15

1.2

1.25

Figura 2.2 Aproximarea functiei hypergeometrice 2 1(1/ 3,1/ 3,2 / 3; )F 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ϕ

1

( ) ( )n

i ii

f x a xϕ=

=∑ (2.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 38H(2.1) o constituie functiile de tip monom ix , astfel ca se 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 :

0 1 1 2 1 2

1 1 2 1 1

( ) ( ) ( )( )... ( )( )...( ) ( ) ( ) ( )n n n n n

f x a a x x a x x x xa x x x x x x R x P x R x− − −

= + − + − − ++ + − − − + = +

(2.2)

unde 1( )nP x− este polinomul aproximant, de grad 1n− , iar ( )nR x este restul aproximarii, presupus foarte mic in modul si care este un polinom de grad n .

Coeficientii ia ai dezvoltarii sunt dati de urmatoarele formule care definesc asa numitele diferente divizate , 1 1( , ,..., )k kDD x x x− :

0 1

2 11 2 1

2 1

3 2 2 12 3 2 1

3 1

1 2 1 2 11

1

( );( ) ( ) ( , );

( , ) ( , ) ( , , );

...................................................................( , ,... ) ( , ,... ) (n n n n

n nn

a f xf x f xa DD x x

x xDD x x DD x xa DD x x x

x x

DD x x x DD x x xa DD xx x

− − −−

=−

= =−

−= =

−

−= =

− , 1,..., 1);nx x−

(2.3)

Prin urmare, formula de interpolare Newton poate fi scrisa sub forma: 1

1 , 1 1 1 12 1 1

( ) ( ) ( , ,..., ) ( ) ( , ,..., , ) ( );k nn

k k i n n ik i i

f x f x DD x x x x x DD x x x x x x−

− −= = =

= + − + −∑ ∏ ∏ (2.4)

Demonstratia se poate face prin inductie completa. Astfel pentru 1n = formula este adevarata deoarece ea se poate scrie:

11

1

( ) ( )( ) ( ) f x f xf x f xx x

− += +

− + 1( )x x− ( );f x= (2.5)

Presupunand acum ca ea este valabila pentru n , putem sa o scriem pentru 1n n= + :

1 11

1 1 1 1 12 1 1

( ) ( ) ( , ,...., ) ( ) ( , ,..., , ) ( )k nn

k k i n n ik i i

f x f x DD x x x x x DD x x x x x x− ++

− += = =

= + − + −∑ ∏ ∏ (2.6)

Separand factorul cu indice 1i n= + din produs si explicitand coeficientul acestuia conform definitiilor diferentelor divizate 39H(2.3), termenul final din 40H(2.6) devine:

1

1 11

1 1 1 11

11

( , ,...., , ) ( )

( , ,...., ) ( , ,...., , ) ( ) ( )

n

n n ii

nn n n n

n iin

DD x x x x x x

DD x x x DD x x x x x x x xx x

+

+=

+ −+

=+

−

−= − −

−

∏

∏ (2.7)

Separand si ultimul termenul cu indice 1k n= + din suma din relatia 41H(2.6) si schimband semnul si la numitorul si la numaratorul fractiei din 42H(2.7) obtinem

1

1 1 1 1 11

( ) ( ) ( , ,..., ) ( ) ( , ,..., ) ( )n

k k i n n ii

f x f x DD x x x x x DD x x x x x+

− +=

= + − + −1−1

2 1

1 1 1 1

1

( , ,...., ) ( , ,...., , )

kn

k i

n n n n

n

DD x x x DD x x x xx x

−

= =

+ −

+

− ++

− +

∑ ∏ ∏

1( )nx x+− + 11

( )n

i

x x=

−∏(2.8)

Se observa ca aceasta relatie se reduce la relatia 43H(2.4) pentru n , pe care am 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:

1

1 , 1 12 1

( ) ( ) ( , ,..., ) ( )kn

k k ik i

f x f x DD x x x x x−

−= =

≈ + −∑ ∏ (2.9)

Restul interpolarii polinomiale este prin urmare ultimul termen din dezvoltarea 44H(2.4)

1 11

( ) ( , ,..., , ) ( );n

n n n ii

R x DD x x x x x x−=

= −∏ (2.10)

Acesta se poate estima folosind functia ajutatoare:

1 1 11

( ) ( ) ( ) ( , ,..., , ) ( )n

n n n ii

Q t f t P t DD x x x x t x− −=

= − − −∏ (2.11)

Aceasta functie are pe intervalul [ ],a b un numar de 1n+ radacini deoarece pentru t x= ea este nula conform relatiei 45H(2.4), iar pentru cele n valori it x= produsul se anuleaza si functia ( )iQ x se anuleaza de asemenea conform relatiei 46H(2.4)

; ( ) 0;; ( ) 0; 1, 2,...,i i

t x Q tt x Q x i n= == = =

(2.12)

Daca ( )Q t are 1n+ radacini, inseamna ca derivata de ordinul intai '( )Q t are n radacini, derivata de ordinul doi ''( )Q t are 1n− radacini, si asa mai departe, iar derivata de ordinul n are o radacina. Fie [ ],a bξ ∈ aceasta radacina

( ) ( ) 0nQ ξ = (2.13) Luad t ξ= in relatia 47H(2.11) si tiand cont ca derivata de ordinul n in raport cu t a produsului din 48H(2.11) este !n , obtinem ( ) ( )

1( ) ( )n nnf Pξ ξ−− 1 1( , ,..., , ) ! 0n nDD x x x x n−− = (2.14)

unde al doilea termen este nul deoarece este derivata de ordin n a unui polinom de grad 1n− . Rezulta

( )

1 1( )( , ,..., , )!

n

n nfDD x x x x

nξ

− = (2.15)

si inlocuind in relatia 49H(2.10) se obtine formula

( )

1

( )( ) ( )!

n n

n ii

fR x x xnξ

=

= −∏ (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 SinAx2

− 3E;H∗ functia de aproximat ∗L

g1= Plot@f@xD, 8x, 0, 2 π<, PlotStyle → BlueD; H∗ graficul functiei ∗Ln= 16; H∗ numarul de termeni ai polinomului ∗L

ForAk= 1, k<= n, r@kD =2. π

n k; k++E; H∗ punctele de esantionare∗L

ForAi= 1, i≤ n,

DD@i+ 1, iD =f@r@i+ 1DD− f@r@iDD

r@i+ 1D− r@iD ; H∗ diferentele finite directe ∗L

i++E;

ForAj= 3, j<= n,

ForAi= 1, i< j−1,

DD@j, j−i− 1D =DD@j, j−iD − DD@j− 1, j−i− 1D

r@jD − r@j− i− 1D ;H∗diferente finite recursive∗L

i++E;

j++E;

p@x_D = f@r@1DD+ „k=2

n

DD@k, 1D ‰l=1

k−1Hx− r@lDL; H∗ Polinomul de interpolare ∗L

g2= Plot@p@xD, 8x, 0, 2 π<, PlotStyle → RedD;Show@g1, g2D; H ∗graficeleambelorfunctii ∗L

Print@"PHxL=", Expand@p@xDDD;

Cele doua functii afisate in final au grafice perfect suprapuse:

1 2 3 4 5 6

-1

1

2

3

4

5

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 ∗Ly= 82, −1, 0, 3, 5, 7, 8, 7, 5<; H∗ valorile corespunzatoare ale functiei ∗Ln= 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,

DD@i+ 1, iD =y@@i+ 1DD− y@@iDDr@@i+ 1DD− r@@iDD ;

i++E;

ForAj= 3, j<= n,

ForAi= 1, i< j−1,


r@@jDD − r@@j− i− 1DD ;

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

j++E;

pN@x_D = y@@1DD + „k=2

n

DD@k, 1D ‰l=1

k−1Hx− r@@lDDL;

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

2

4

6

8

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 1n− 1( ) ( )nf x P x−≈ (2.17) folosind formula

11 1

( ) ( ) ( )nn

jn i

i j i jj i

x xf x P x f x

x x−= =

≠

−≈ =

−∑ ∏ (2.18)

Aceasta relatie se mai poate scrie si sub forma

1

( ) ( ) ( )n

i ii

f x f x L x=

≈ ∑ (2.19)

unde functiile ( )iL x sunt polinoamele Lagrange avand expresia

1

( )n

ji

j i jj i

x xL x

x x=≠

−=

−∏ (2.20)

Ca si in cazul interpolarii Newton, trebuie sa existe egalitate stricta in relatia 50H(2.19) in punctele in care functia este cunoscuta, ix , iar in afara lor sa existe o diferenta cat mai mica, reprezentata de acelasi rest dat de relatia 51H(2.16).

1

( ) ( ) ( )n

i i i ii

f x f x L x=

= ∑ (2.21)

1

( ) ( ) ( ) ( )n

i i ni

f x f x L x R x=

= +∑ (2.22)

Pentru demonstarea formulei, scriem mai intai forma desfasurata a polinoamelor Lagrange

1 2 1 1

1 2 1 1

( )( )...( )...( )( )...( )( )( )( )...( )...( )( )...( )

k i i ni

i i i k i i i i n

x x x x x x x x x x x xL xx x x x x x x x x x x x

− +

− +

− − − − − −=

− − − − − − (2.23)

Se observa ca pentru ix x= rezulta ( ) 1i iL x = , deoarece fiecare factor de la numarator se simplifica cu unul de la numitor. Pentru kx x= unde k i≠ , exista la numarator un factor ( ) 0k kx x− = , astfel ca

( ) 0i kL x = . Rezulta urmatoarea proprietate a polinoamelor Lagrange: ( )i k ikL x δ= (2.24) Calculand cu relatia 52H(2.18) valoarea functiei polinomiale intr-un punct jx , rezulta:

11

( ) ( ) ( )n

n j i ij ji

P x f x f xδ−=

= =∑ (2.25)

astfel ca este demonstrata relatia 53H(2.21). 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 + 1E;

a= 3;b= 9;g1= Plot@f@xD, 8x, a, b<, PlotStyle → BlueD;n= 21;r= Table@0, 8n<D;

AbsoluteTimingAForAk= 1, k<= n, k++, r@@kDD = a+b−a

n kE;

p@x_D = „i=1

n

f@r@@iDDD i

kjjjj‰

j=1

i−1 x− r@@jDDr@@iDD − r@@jDD

y

{zzzz i

kjjjj ‰

j=i+1

n x− r@@jDDr@@iDD − r@@jDD

y

{zzzz;E

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 1P polinomul de gradul zero (o constanta) care trece prin

punctul {{ }1 1( ),f x x , deci 1 1( )P f x= , cu 2P polinomul de gradul zero (o constanta)

care trece prin punctul {{ }2 2( ),f x x , deci 2 2( )P f x= , si asa mai departe pana la ( )n nP f x= . Din aceste polinoame se pot forma polinoame de grad unu (o dreapta) notate

astfel: cu 12P polinomul care trece prin punctele { }1 1( ),f x x si { }2 2( ),f x x , cu 23P

polinomul care trece prin punctele { }2 2( ),f x x si { }3 3( ),f x x , si asa mai departe pana

la 1,n nP − polinomul care trece prin punctele { }1 1( ),n nf x x− − si { }( ),n nf x x . Se continua acest proces prin formarea de polinoame din grad din ce in ce mai

mare, pana se ajunge la 1,2,..., 1,n nP − polinomul care trece prin punctele

{ }1 1( ),f x x ,{ }2 2( ),f x x ,…, { }1 1( ),n nf x x− − si { }( ),n nf x x , adica tocmai polinomul cautat.

Se poate alcatui deci o schema arborescenta in care fiecare coloana se calculeaza pe rand pe baza coloanei precedente.

1 1 1

12

2 2 2 123

23 1234

3 3 3 234

34

4 4 4

: ( )

: ( )

: ( )

: ( )

x f x PP

x f x P PP P

x f x P PP

x f x P

=

=

=

=

(2.26)

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

grad 1k − dupa urmatoarea formula obtinuta pe baza diferentelor divizate:

( 1)...( 1) ( 1)( 2)...( )( 1)...( )

( ) ( )i k i i i k i i i i ki i i k

i i k

x x P x x PP

x x+ + + − + + +

+ ++

− + −=

− (2.27)

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, 8M<D;AbsoluteTimingAForAk= 1, k≤ M, k++,

r@@kDD = a+b−a

M k; d@k, kD = f@r@@kDDDE;

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

d@i, i+ mD = ExpandA Hx− r@@i+ mDDL d@i, i+ m − 1D−Hx− r@@iDDL d@i+ 1, i+ mDr@@iDD − r@@i+ mDD EE;EE

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 =1

x2 + 1;H∗ functia Runge ∗L

a= −3.0;b= 3.0;g1= Plot@f@xD, 8x, a, b<, PlotStyle → BlueD; H∗ graficul functiei ∗Ln= 10; H∗ numarul de termeni ai polinomului ∗L

ForAk= 1, k<= n, r@kD = a+b− a

n k; k++E; H∗ punctele de esantionare∗L

ForAi= 1, i≤ n,



i++E;

ForAj= 3, j<= n,

ForAi= 1, i< j−1,



i++E;

j++E;

p@x_D = f@r@1DD+ „k=2

n

DD@k, 1D ‰l=1


g2= Plot@p@xD, 8x, a, b<, PlotStyle → RedD;Show@g1, g2D; H∗ graficele ambelor functii ∗LPrint@"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:

-3 -2 -1 1 2 3

0.5

1

1.5

2

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

in puncte echidistante a functiei Runge este urmatorul:

-3 -2 -1 1 2 3

-1

-0.5

0.5

1

1.5

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 ( )nR x dat de relatia

54H(2.16)

( )

1

( )( ) ( )!

n n

n ii

fR x x xnξ

=

= −∏

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 ix astfel incat produsul

1

( ) ( )n

n ii

p x x x=

= −∏ (2.28)

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 ix sunt chiar radacinile unui polinom monic Chebyshev de grad n , produsul 55H(2.28) este un astfel de polinom si deci restul 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:

( ) cos( arccos( ))nT z n z= (2.29) Notand [ ]arccos ; 1,1z zθ = ∈ − (2.30) rezulta forma echivalenta: ( ) cosnT z nθ= (2.31) unde cosz θ= (2.32) Pornind de la aceasta definitie rezulta ca primele doua polinoame Chebyshev sunt:

0

1

( ) cos0 1( ) cos

T zT z z

θθ

= == =

(2.33)

Restul polinoamelor Chebyshev se pot calcula si prin relatia de recurenta: 1 2( ) 2 ( ) ( ); 2n n nT z zT z T z n− −= − ≥ (2.34) Relatia 56H(2.34) se poate demonstra pe baza urmatoarei identitati trigonometrice:

1cos( ) cos( ) [cos( ) cos( )]2

α β α β α β= + + − (2.35)

Luand in aceasta relatie α θ= si ( 1)nβ θ= − se obtine:

1 2

( )( ) ( )

1cos cos( 1) [cos cos( 2) ]2

nn nz T zT z T z

n n nθ θ θ θ− −

− = + − (2.36)

de unde se deduce relatia 57H(2.34). Se poate demonstra ca functiile ( )nT z sunt polinoame de grad n , avand urmatoarea forma generala: 1

1( ) 2 ( )n nn nT z z R z−

−= + (2.37) unde 1( )nR z− este un polinom de grad 1n − . Vom folosi notatia { }1 1( ) ( )n nR z z− −∈ ,

unde { }1( )n z− este multimea polinoamelor de grad cel mult 1n − . Demonstratia se poate face prin inductie completa. Intr-adevar pentru 1n = relatia 58H(2.37) este evidenta. Presupunand acum ca ea este valabila pentru n , in cazul

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

1 1 1 2

( )

( ) 2 ( ) ( ) 2 2 2 ( ) 2 ( ) 2 ( )n

n n n n n nn n n n n n

R z

T z zT z T z z z zR z z R z z R z− − − ++ − − −= − = + − − = + (2.38)

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

( ) 0; cos 0; (2 1) ; 1, 2,...,2n k k k

nradacini

T z n n k k nπθ θ= ⇒ = = − = (2.39)

1 ;2k k

nπθ ⎛ ⎞= −⎜ ⎟

⎝ ⎠ (2.40)

1cos2kz k

nπ⎛ ⎞= −⎜ ⎟

⎝ ⎠ (2.41)

Extremele polinomului au valoarea 1 sau -1: ( ) cos [ 1,1]nT z nθ= ∈ − (2.42) Ele se pot obtine simplu punand conditia de extrem in definitia 60H(2.31):

| cos | 1 ;k k kn n k knπθ θ π θ= ⇒ = = (2.43)

1 extreme

cos ; 0,1,2,...,k

n

z k k nnπ

+

⎛ ⎞= =⎜ ⎟⎝ ⎠

(2.44)

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

n kT z = − (2.45) Conform celor precizate in alineatul precedent, rezulta ca la extrema dreapta a domeniului de definitie avem extremul 0 1z = , pentru care valoarea polinomului este intotdeauna 0( 1) 1− = , in timp ce la extrema stanga avem extremul 1nz = − , 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

1n + extreme. Programul 2.6 Calculul si reprezentarea polinoamelor Chebyshev

T@0D = 1;T@1D = x;T@n_D:= 2 xT@n− 1D − T@n−2D;gr= 8<;ForAn= 0, n≤ 6, n++, Print@Expand@T@nDDD

AppendToAgr, PlotAT@nD, 8x, −1, 1<,

PlotStyle→ RGBColorA n10

,n10

,10− n

10EEEE;

Show@grD; Reprezentarea pe un singur grafic a primelor 7 polinoame Chebyshev, generata de acest program este urmatoarea:

-1 -0.5 0.5 1

-1

-0.5

0.5

1

De asemenea programul genereaza urmatoarele reprezentari pe grafice individuale pentru aceste polinoame:

1

-1 -0.5 0.5 1

0.5

1

1.5

2

x

-1 -0.5 0.5 1

-1

-0.5

0.5

1

−1 + 2 x2

-1 -0.5 0.5 1

-1

-0.5

0.5

1

−3 x + 4 x3

-1 -0.5 0.5 1

-1

-0.5

0.5

1

1− 8 x2 + 8 x4

-1 -0.5 0.5 1

-1

-0.5

0.5

1

5 x− 20 x3 + 16 x5

-1 -0.5 0.5 1

-1

-0.5

0.5

1

−1 + 18 x2 − 48 x4 + 32 x6

-1 -0.5 0.5 1

-1

-0.5

0.5

1

Numim polinoame monice acele polinoame care au coeficientul termenului de gradul cel mai inalt egal cu 1. Din relatia 62H(2.37) se observa ca se poate obtine simplu un polinom monic Chebyshev de gradul n , pe care il vom nota in continuare ( )nT z , prin impartirea unui polinom Chebyshev cu 12n−

1 111

( )( ) 2 ( ) 2 ( )2

n n nnn n nn

T zT z T z z R z− −−−= = = + (2.46)

Conform definitiei 63H(2.31) se poate scrie ( ) 1nT z ≤ (2.47) deci polinomul monic Chebyshev are proprietatea

11

1( ) 22

nn nT z −

−≤ = (2.48)

Demonstratia teoremei mini-max. Aceasta demonstratie se va face prin reducere la absurd. Sa presupunem ca ar exista un polinom monic ( )nf z cu proprietatea

1( ) ( ) 2 nn nf z T z −< ≤ (2.49)

Daca ( )nf z este un polinom monic se poate scrie 1( ) ( )n

n nf z z F z−= + (2.50) cu { }1 1( ) ( )n nF z z− −∈ Sa formam functia ajutatoare ( ) ( ) ( )n nQ z T z f z= − (2.51) Tinand seama de relatiile 64H(2.46) si 65H(2.50) se observa ca avem 1

1 1( ) 2 ( ) ( )nn nQ z R z F z−− −= −

{ }1( ) ( )nQ z z−∈ (2.52) deci functia ( )Q z este un polinom de grad cel mult 1n − . Sa studiem comportarea acestei functii intre doua extreme consecutive ale polinomului monic Chebyshev. Evisent, acest polinom monic are tot 1n + 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 2 1k ± . In extremul cu numar par, conform relatiei 66H(2.45), rezulta ca polinomul monic este pozitiv

2

2 1

( 1)( ) 02

k

n k nT ξ −

−= > (2.53)

deci 2 2( ) ( )n k n kT Tξ ξ= (2.54)

iar functia ajutatoare este pozitiva 2 2 2 2 2( ) ( ) ( ) ( ) ( ) 0k n k n k n k n kQ T f T fξ ξ ξ ξ ξ= − = >∓ (2.55)

deoarece conform ipotezei modulul lui ( )nf z este mai mic decat al polinomului monic Chebyshev. Pe de alta parte, in extremul cu numar impar, conform relatiei 67H(2.45), rezulta ca polinomul monic este negativ

2 1

2 1 1

( 1)( ) 02

k

n k nT ξ±

± −

−= < (2.56)

deci 2 1 2 1( ) ( )n k n kT Tξ ξ± ±= − (2.57)

iar functia ajutatoare este negativa 2 1 2 1 2 1 2 1 2 1( ) ( ) ( ) ( ) ( ) 0k n k n k n k n kQ T f T fξ ξ ξ ξ ξ± ± ± ± ±= − = − <∓ (2.58)

deoarece conform ipotezei modulul lui ( )nf 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 1n + 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 1n − , ceea ce arata ca ipoteza este falsa si deci este adevarata contrara ei: 1( ) ( ) 2 n

n nf z T z −≥ ≤ (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 68H(2.41). Daca avem un interval oarecare [ ],a b , se poate utiliza o schimbare de variabila de tipul

2 2

b a b az x− += + (2.60)

astfel ca nodurile de interpolare sunt date de relatia

1cos , 1,2,...,2 2 2

b a b az k k nnπ+ − ⎡ ⎤⎛ ⎞= + − =⎜ ⎟⎢ ⎥⎝ ⎠⎣ ⎦

(2.61)

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 69H(2.61), conform teoremei mini-max.

Programul 2.7- Interpolare polinomiala globala a functiei Runge cu esantionare in noduri Chebashev

f@x_D =1

x2 + 1;H∗ functia Runge ∗L

a= −3.0;b= 3.0;g1= Plot@f@xD, 8x, a, b<, PlotStyle → BlueD;H∗ graficul functiei ∗Ln= 10; H∗ numarul de termeni ai polinomului ∗L

ForAk= 1, k<= n, r@kD =a+ b

2+

b− a2

CosAJk−12N

π

nE; k++E; H∗ punctele de esantionare∗L

ForAi= 1, i≤ n,



i++E;

ForAj= 3, j<= n,

ForAi= 1, i< j−1,



i++E;

j++E;

p@x_D = f@r@1DD+ „k=2

n

DD@k, 1D ‰l=1


g2= Plot@p@xD, 8x, a, b<, PlotStyle → RedD;Show@g1, g2D; H∗ graficele ambelor functii ∗LPrint@"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:

-3 -2 -1 1 2 3

0.2

0.4

0.6

0.8

1

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:

-3 -2 -1 1 2 3

0.2

0.4

0.6

0.8

1

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 { }iy determinate experimental ale unei functii necunoscute pentru anumite valori ale argumentului { }ix .

Clear@"̀ ∗"D;x= 80, 1.5, 2, 3.5, 4, 5, 6, 7, 8, 9<;H∗ abscisele in care este cunoscuta functia ∗Ly= 8−2, 1, 3, 5, 2, 1, 2, −1, 3, −1<;H∗ valorile functiei in punctele cunoscute ∗Lg1= ListPlot@lstD;H∗ reprezentarea prin puncte distincte ∗Lg2= ListPlot@lst, PlotJoined→ TrueD;H∗Reprezentarea prin segmente∗LShow@g1, g2D; H∗ Graficul compus din puncte unite prin segmente ∗L

Rezulta urmatoarea reprezentare grafica:

2 4 6 8

-2

-1

1

2

3

4

5

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 2 3 4 5

-1

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

2 4 6 8 10

-1

-0.5

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 0[ , ]nx x in mai multe subintervale: 0 0 1 1 2 1 1[ , ] [ , ] [ , ] ... [ , ] ... [ , ]n i i n nx x x x x x x x x x+ −= ∪ ∪ ∪ ∪ ∪ (2.62) Pe fiecare subinterval 1[ , ]i ix x + , functia este aproximata printr-un polinom de gradul trei ( )iS x de forma

3 2

1

1

( ) ( ) ( ) , [ , ]( )

0, [ , ]i i i i i i i i i

ii i

a x x b x x c x x d x x xS x

x x x+

+

⎧ − + − + − + ∈⎪= ⎨∉⎪⎩

(2.63)

0,1,..., 1i n= − Pe ansamblu functia este aproximata printr-o suma ( )S x de astfel de functii cu domenii de definitie disjuncte

1

0

( ) ( ) ( )n

ii

f x S x S x−

=

≈ =∑ (2.64)

Pentru a se putea calcula cei 4n coeficienti din relatiile 70H(2.63) se vor impune conditii de racordare in punctele intermediare: functiile si primele lor doua derivate trebuie sa fie continue in aceste 1n − puncte. Rezulta 3( 1)n − ecuatii, care impreuna cu cele 1n + ecuatii care dau valorile cunoscute ale functiei in cele 1n + puncte (incluzand si cele de la capetele intervalului) dau un total de 4 2n − 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: ( ) 0F x = , :F → (3.1) Valorile respective sunt numite radacinile ecuatiei 71H(3.1). 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 72H(3.1) si observand abaterea fata de 0. 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 ( ) ( ) 0F a F b < , 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.

X1 X2X3 X4

f(x)

x0

h 2h 3h 4h 5h 6h 7h 8h 9h-h

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

23 5 2 0x x− + =

Programul 3.1. Gasirea intervalelor cu radacini prin esantionare

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

h= 10−2; H∗ Pasul de esantionare ∗Llst= 8<; H∗ Initializare lista intervalelor ∗LFor@x= a− h, x ≤ b+ h, x = x+ h,yn= f@x+ hD;If@ynyv≤ 0, AppendTo@lst, 8x, x+ h<DD; H∗ Adauga intervalul in lista ∗Lyv= yn;H∗ Noua ordonata devine ordonata veche pentru urmatorul interval∗LD;Print@lstD;

-2 -1 1 2 3 4 5

-10

10

20

30

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

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 ( ) ( ) 0v nF x F x ≤ , deoarece exista posibilitatea ca o radacina sa fie chiar in punctul de esantionare (ca in cazul radacinii 1x= 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:

1

, ( ) ( ) 02

, ( ) ( ) 02

kk

kk

k

x a daca F a F xx

b x daca F a F x+

−⎧ <⎪⎪= ⎨ −⎪ >⎪⎩

(3.2)

iar noul intervat este [ ]1,kx a+ daca 1( ) ( ) 0kF a F x + < si [ ]1, kb x + in caz contrar.

Rezulta ca dupa n bisectii, intervalul se micsoreaza de 2n ori, deci si eroarea se micsoreaza in aceeasi proportie:

0 1 01 12 2n n nx x x x b aε = − = − ≤ − (3.3)

unde nx este valoarea obtinuta pentru radacina dupa n bisectii 0x este valoarea reala iar

1x 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 73H(3.3):

2 1010 ,3

b a b an Log Log nε ε− −

≥ ≈ ∈ (3.4)

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:

112n n nConstε ε ε+ = = × (3.5)

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: 1 , 1,n nConst αε ε α α+ = × > ∈ (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 74H(3.2), parametrii de intrare fiind definitia 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 ∗La= 0; b = 3;H∗ Intervalul pe care se afla o radacina reala ∗Lh= 10−15;H∗ Precizia dorita ∗LWhileAAbs@b− aD > h, H∗Ciclul repetat pana la atingerea preciziei∗L

c=a+ b2.

; H∗ Noua valoare a solutiei ∗L

If@f@aDf@cD > 0, a= c, b= cDE;H∗Decizia de alegere a intervalului ∗Lx= a;H∗Alegerea unei valori din intervalul final ∗LPrint@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: ( ) 0F x = (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 75H(3.7) este echivalenta cu ( ) 0x f x− = (3.9)

O solutie α pentru ecuatia 76H(3.9) se numeste punct fix, care este si solutie a ecuatiei 77H(3.7) deoarece conform 78H(3.9) are proprietatea: ( )fα α= (3.10) si conform 79H(3.8) ( ) ( ) 0F fα α α= − = . Se va demonstra in continuare ca un proces iterativ de forma 1 ( )k kx f x+ = (3.11) poate asigura convergenta catre punctul fix α , solutie a ecuatiei 80H(3.7) daca sunt indeplinite anumite conditii.

Intr-adevar, eroarea la iteratia 1k + poate fi exprimata in functie de eroarea la iteratia precedenta k , tinand cont de relatiile 81H(3.10) si 82H(3.11): 1 1 ( ) ( ) '( ) '( )k k k k kx f x f f x fε α α ξ α ξ ε+ +≡ − = − ≤ − = (3.12) Inegalitatea din ecuatia 83H(3.12) a fost scrisa in baza teoremei mediei, ξ fiind o valoare (necunoscuta) din intervalul [ ] [ ], ,kx 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

( ) ( )

'( ) 1k

k

f x ff

xα

ξα

−= <

− (3.13)

Rezulta, tinand cont de 84H(3.11) si 85H(3.8), urmatoarea regula de iteratie: 1 [ ]k k kx x F x+ = − (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 86H(3.13), care insa nu se poate verifica apriori, nefiind cunoscuta valoarea ξ . Convergenta se poate verifica prin scaderea lui ε sub o anumita limita impusa, si obligatoriu prin verificarea solutiei in ecuatia 87H(3.7). Metoda 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 88H(3.14) sub forma: 1 [ ]k k kx x F xβ+ = − (3.15) Aceasta ar insemna ca prin inmultirea corectorului [ ]kF x cu un factor corespunzator β se poate obtine in urma iteratiei 89H(3.15) o valoare mult mai apropiata de radacina, deci procesul poate deveni mult mai rapid convergent. Tinand seama in 90H(3.15) de definita 91H(3.8) obtinem: [ ]1 ( )k k k kx x f x xβ+ = + − (3.16) Factorul β poate fi determinat punand conditia ca eroarea la iteratia 1k + sa fie minima. Rezulta succesiv:

[ ]1 1( )( ) 1 k k

k k k k k kk

f x xx x f x x xx

ε α α β α βα+ +

−= − = − + − = − +

−

Adunand si scazand ( )f α la numaratorul fractiei din relatia precedenta si aplicand din nou teorema mediei obtinem: [ ]1 1 '( ) 1k k kfε ε β ξ+ = + − (3.17)

Daca am cunoaste valoarea ξ am putea alege coeficientul β astfel ca eroarea sa se anuleze:

11 '( )kf

βξ

=−

(3.18)

In practica insa nu cunoastem pe ξ , deci trebuie sa ii atribuim una din valorile cunoscute. Considerand ca cea mai apropiata valoare este kx , vom scrie

( ) ( )k kf f xξ = iar relatia 92H(3.16) devine:

[ ]11( )

1 '( )k k k kk

x x f x xf x+ = + −

− (3.19)

Daca tinem seama de relatia 93H(3.8) si de derivata acesteia, relatia 94H(3.19) poate fi scrisa sub forma:

1( )'( )

kk k

k

F xx xF x+ = − (3.20)

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:

21( ) ( ) '( ) ''( ) ...2k k kF x F F Fα ε α ε α= + + + (3.21)

'( ) '( ) ''( ) ...k kF x F Fα ε α= + + (3.22) In relatia 95H(3.21), deoarece α este radacina, rezulta ( ) 0F α = ; din relatia 96H(3.22), daca neglijam termenii care cuprind derivata de ordinul doi si cei superiori (lucru valabil daca suntem suficient de aproape de solutia corecta), obtinem '( ) '( )kF x F α= (3.23) Daca scadem α din ambii membri ai relatiei 97H(3.20) si tinem seama de 98H(3.21) si 99H(3.23), rezulta:

1

21

( ) 1 ''( )'( ) 2 '( )

k k

kk k k k k

k

F x Fx xF x F

ε ε

αα α ε ε εα

+

+ − = − − = − − (3.24)

care este echivalenta cu: 2

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

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 101H(3.20), parametrii de intrare fiind 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∗Lc= a;H∗Valoare initiala pentru radacina∗Lcv= b;H∗Valoare initiala pentru comparatie∗L

DoAc= cv−f@cvD

f'@cvD ;H∗ Calcul valoare noua ∗L

If@Abs@c− cvD < ε, Break@D D; H∗Iese din ciclu daca s−a atins precizia∗Lcv= 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 102H(3.20) poate fi inlocuita prin formula:

1

1

( ) ( ) ( )'( ) n n nn

n n n

F x F x F xF xx x x

−

−

Δ −≈ =

Δ − (3.26)

daca iteratia n se face intre punctele 1nx − si nx . 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 103H(3.26) in ecuatia 104H(3.20), si daca aproximatia radacinii la iteratia n este nx , rezulta pozitia noului punct 1nx + :

11

1

( )( )( ) ( )

n n nn n

n n

F x x xx xF x F x

−+

−

−= −

− (3.27)

sau

1 11

1

( ) ( )( ) ( )

n n n nn

n n

x F x x F xxF x F x

− −+

−

−=

− (3.28)

Relatia 105H(3.28) defineste metoda secantei, denumire datorata faptului ca 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

convergenta fiind 1 5 1.6182+

≈ , adica asa numitul “golden number”.

Programul 8.4, scris in Mathematica, pentru rezolvarea unei ecuatii neliniare prin metoda secantei foloseste ecuatia 106H(3.28), parametrii de intrare fiind:

• 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 ∗Lcv= a;H∗Valoare initiala a radacinii, pentru comparatie∗LForAn= 1, n< Nmax, n++,

c=af@bD − bf@aD

f@bD− f@aD ; H∗Valoarea noua a radacinii∗L

If@f@aDf@cD > 0, a= c, b= cDIf@Abs@c− cvD < ε, Break@D D;

cv= cE;H∗Alegerea noului interval∗LPrint@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 107H(3.28), scazandu-se radacina α din ambii membri si tinand cont ca ( ) 0F α = . Rezulta succesiv:

[ ] [ ]

1

1 11 1

1

1 1 1

1

1 11

1

''( )

( ) ( ) ( ) ( ) ( ) ( )( ) ( )

( ) '( )( ) ( ) '( )( )( ) ( )

'( ) '( )( ) ( )( )

n n

n n n nn n

n n

n n n n n n

n n

n n n nn n

n n n

F

x F x F x F x Fx

F x F xx F x x x F x x

F x F xF x F x x xx x

x x F x Fε ε

ξ

α α α αε α

α α α α

α α−

− −+ +

−

− − −

−

− −−

−

− − − − −≡ − = =

−− − − − −

=−

− −= − −

− − 1

1/ '( )

( )n

F

xη

−

1 1''( )'( )

nn n n

n

FF

ξε ε εη+ −= (3.29)

S-a aplicat teorema mediei, nξ si nη fiind doua valori cuprinse in intervalul [ ]1,n nx x− .

Rescriind ecuatia 108H(3.29) pentru 2nε + , si logaritmand-o obtinem:

2 1''( )log( ) log( ) log( ) log'( )

nn n n

n

FF

ξε ε ε

η+ +

⎛ ⎞= + + ⎜ ⎟

⎝ ⎠

Notand

2 2 1 1''( )log( ) , log( ) , log( ) , log'( )

nn n n n n n n

n

F cF

ξε ν ε ν ε ν

η+ + + +

⎛ ⎞= = = =⎜ ⎟

⎝ ⎠ (3.30)

se obtine: 2 1n n n ncν ν ν+ +− − = (3.31) Definim operatorul de deplasare care face trecerea de la eroarea in iteratia n la eroarea din iteratia 1n + sub forma: 1n nEν ν += care introdus in ecuatia 109H(3.31) conduce la: 2( 1) n nE E cν− − =

sau ( )( ) n nE q E p cν− − = unde am notat 1 5 1 5,

2 2p q+ −= =

Notand ( ) n nE p uν− = (3.32)

rezulta ( ) n nE q u c− = si schimband indicele 1 1( ) n nE q u c− −− = . Tinand cont de actiunea operatorului E, aceasta relatie devine 1 1n n nu c qu− −= + (3.33)

Rescriind 110H(3.33) cu schimbarea indicelui 1n n→ − obtinem 1 2 2n n nu c qu− − −= + (3.34) care introdusa in 111H(3.33) duce la: 2

1 2 2n n n nu c qc q u− − −= + + Repetand procedeul prin inlocuirea lui 2nu − printr-o relatie analoaga cu 112H(3.34), obtinem in final: 2 1

1 2 3 0 0... n nn n n nu c qc q c q c q u−

− − −= + + + + + (3.35) Tinand cont de actiunea operatorului E, si de definitia 113H(3.32), membrul stang al ecuatiei 114H(3.35) devine: 1( )n n n nu E p pν ν ν+= − = − .

Membrul drept contine o serie alternata 2 11 2 3 0... n

n n nc qc q c q c−− − −+ + + +

(deoarece1 5 0

2q −= < ), avand deci o valoare finita, precum si cantitatea finita

0 0 1 0( ) ( )n n nq u q E p q pν ν ν= − = − . Notand membrul drept al ecuatiei 115H(3.35) cu L < ∞ , aceasta relatie devine: 1n np Lν ν+ − = (3.36)

si tinand cont de notatia 116H(3.30) rezulta 1log( ) log( )pn n Lε ε+ − = , de unde se deduce relatia care da

ordinul de convergenta al metodei secantei:

1 5

21

Ln neε ε

+

+ = (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 kx , 1kx − si 2kx − , 2k ≥ 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 [ ]2 2, ( )k kx F x− − , [ ]1 1, ( )k kx F x− − si [ ], ( )k kx F x va avea ecuatia:

2 ( )a b c p xν ν+ + = (3.38)

x

F(x)

xkxk-1xk-2

0

h1

h2

ν

Figura 3.3 Metoda parabolei (Muller) Notand 1kx xν −= − (ne raportam la punctul aflat in mijloc, prin care va trece deci si ordonata, conform figurii), si 1 1k kh x x −= − si 2 1 2k kh x x− −= − , 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:

21

21 1

22 2 2

0 0 ( )

( )

( )

k

k

k

a b c F x

ah bh c F x

ah bh c F x

−

−

⎧ + + =⎪

+ + =⎨⎪ − + =⎩

(3.39)

unde s-a tinut cont că în aceste puncte variabila ν a parabolei are valorile 0, 1h respectiv 2h− , conform figurii.

Pentru aflarea parametrilor a , b si c ai parabolei, rezolvam acest sistem de ecuatii. Din prima ecuatie rezulta: 1( )kc F x −= (3.40) Inmultind a doua ecuatie cu 2h si a treia cu 1h si adunandu-le rezulta

2 1

1 1 2 2 1 2 1 2

( ) ( ) ( )( ) ( )

k k kF x F x F xah h h h h h h h

− −= + −+ +

(3.41)

Inmultind a doua ecuatie cu 2

2h si a treia cu 21h si scazandu-le rezulta

2 1 2 1 2 1

1 1 2 2 1 2 1 2

( ) ( ) ( ) ( )( ) ( )

k k kh F x h F x h h F xbh h h h h h h h

− −−= − +

+ + (3.42)

Parabola intersecteaza abscisa in noul punct 1kx + care rezulta din ecuatia ( ) 0p ν = cu parametrii dati de ecuatiile 117H(3.40), 118H(3.41) si 119H(3.42).

2

1,2 2

4 22 4

b b ac ca b b ac

ν − ± −= =

− −∓

Tinand cont de notatia 1kx xν −= − , noua valoare pentru radacina va fi deci:

1 1 2

24

k kcx x

b b ac+ −= −

± − (3.43)

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 120H(3.40) - 121H(3.43), si esantioneaza aleator planul complex intr-un numar stN 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, stN , 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 maxN 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 stN ori, in planul complex in patratul [ ] [ ]0.5 ,0.5 0.5 ,0.5L i L i L L− ⋅ ⋅ × − .

Programul 8.5 Calculul radacinilor unei ecuatii prin metoda Muller f@x_D = 10+ 7 Sin@xD − 3x2 − 2x3 + x4 + 2 x5;X= 8<;H∗ Vectorul radacinilor ∗L;

Nst= 40; H∗ Numarul de repetari ale startului∗LL= 10; H∗Latura patratului in planul complex in care se cauta radacini∗LForAn= 0, n< Nst, n++,

x2= L HRandom@D− 0.5L + LIHRandom@D −0.5L; H∗ Puncte de start aleatoare ∗Lx1= L HRandom@D− 0.5L + LIHRandom@D −0.5L;x0= L HRandom@D− 0.5L + LIHRandom@D −0.5L;H∗ Se esantioneaza aleator spatiul complex ∗Lε = 10−13; H∗ Eroarea admisibila in valoarea functiei ∗LNmax= 100; H∗ Numarul maxim de iteratii daca nu se obtine convergenta ∗Lk= 1;DoAh1= x0− x1;

H∗ Ciclul pentru cautarea unei radacini in apropierea punctelor curente ∗Lh2= x1− x2;c= f@x1D;

a=f@x0D

h1 Hh1+ h2L +f@x2D

h2 Hh1+ h2L −f@x1Dh1h2

;

b=h2f@x0D

h1 Hh1+ h2L −h1f@x2D

h2 Hh1+ h2L +Hh1− h2L f@x1D

h1h2;

x2= x1; H∗ Se reactualizeaza punctele x2 si x1∗Lx1= x0;

H∗ Pentru reactualizarea lui x0 se calculeaza pasul h= 2 c

b+"#### ## #### ## #####b2−4 a csau h= 2 c

b+"#### ## ###### #####b2−4 a c∗L

d1= b+è!!!!!!!!!!!!!!!!!

b2 − 4 ac;

d2= b−è!!!!!!!!!!!!!!!!!

b2 − 4 ac;

IfAAbs@d1D > Abs@d2D, h=2 cd1

, h =2 cd2

E;

x0= x2− h; H∗x2 a capatat mai sus valoarea x1, fata de care am scris pe h ∗LIfAAbs@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 opresc iteratiile. Daca partea imaginara a radacinii este prea mica

se adauga doar partea reala la vectorul radacinilor ∗Lk++, 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 122H(3.1) in care membrul stang este un polinom de grad n care poate fi scris sub doua forme:

0 00 1

( ) ( ), 0nn

jn n j j

j j

P x a x a x x a−= =

= = − ≠∑ ∏ (3.44)

Intre radacinile si coeficientii ecuatiei polinomiale exista relatiile Viete:

11 2

0

21 2 1 3 1

0

1 2 1 1 2 1 10

1 20

...

...

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

... ... ... ( 1)

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

... ( 1)

n

n n

jjj j j j

n nn

ax x xa

ax x x x x xa

ax x x x x x x x

a

ax x xa

−

− − +

+ + + = −

+ + + =

+ + = −

= −

(3.45)

Daca una dintre radacini, s-o numim 1x ar fi dominanta, adica: 1 2 1 3 1, ,...., nx x x x x x (3.46) folosind prima relatie Viete in care neglijam toate celelalte radacini, putem sa obtinem valoarea radacinii dominante:

11

0

axa

= − (3.47)

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 1x , ea avand proprietatea: 1 2 1 3 1, ,...., nx x x x x x> > > (3.48) Daca am reusi sa ridicam toate radacinile la o putere suficient de mare, relatia 123H(3.48) ar deveni o relatie de tipul 124H(3.46), ceea ce ar permite calculul simplu al radacinii dominante si apoi a celorlalte. Pentru aceasta, sa consideram mai intai polinomul in variabila x− :

00 1

( ) ( ) ( 1) ( )nn

j nn n j j

j j

P x a x a x x−= =

− = − = − +∑ ∏ (3.49)

Facand produsul ( 1) ( ) ( )nn nP x P x− − , obtinem un polinom in variabila

2x , pe care sa-l notam cu (1) 2( )nP x :

(1) 2 2 2 2 (1) 20

01

( ) ( 1) ( ) ( ) ( )n n

n jn n n j j

jj

P x P x P x a x x A x==

= − − = − =∑∏ (3.50)

unde am notat cu (1)jA coeficientii noului polinom (1) 2( )nP x .

Presupunand ca 1 2x x> , cu atat mai mult 2 21 2x x> , ceea ce incepe sa ne

apropie de situatia unei radacini dominante ca in 125H(3.46). Repetand procedeul obtinem un polinom in 4x :

(2) 4 (1) 2 (1) 2 4 4 4 (2) 40

01

( ) ( 1) ( ) ( ) ( )n n

n jn n n j j

jj

P x P x P x a x x A x==

= − − = − =∑∏

Dupa aplicarea procedeului de s ori, se obtine:

1 1( ) 2 ( 1) 2 ( 1) 2 2 2 2 ( ) 2

001

( ) ( 1) ( ) ( ) ( )s s s s s s s

n ns n s s s j

n n n j jjj

P x P x P x a x x A x− −− −

==

= − − = − =∑∏ (3.51)

Se observa ca puterile radacinilor cresc foarte rapid cu s, ajungandu-se la o relatie de ordine de tipul 126H(3.46). Notand 2s m= , avem din prima relatie Viete:

11

0

sm

s

AxA

= − (3.52)

Aceasta permite calculul primei radacini dupa formula:

1

11

0

s m

s

AxA

= ± (3.53)

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: 1 2 3 ,...., nx x x x , putem calcula si celelalte radacini. Intr-adevar, a doua relatie Viete devine:

( )2

1 1 ( )0

sm m

s

Ax xA

= (3.54)

si impartind-o la relatia Viete precedenta 127H(3.52) obtinem a doua radacina:

1( )2

2 ( )1

s m

s

AxA

= ± (3.55)

Analog se obtin si celelalte radacini

1( )

( )1

s mj

j sj

Ax

A −

= ± (3.56)

semnul alegandu-se de asemenea prin testare in ecuatie. Coeficientii s

jA ai polinomului final ( ) 2( )ss

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

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

1( 1) 2 ( ) 2 ( ) 2 ( ) 2 ( ) 2

0 0 0( 1) ( ) ( ) ( 1)

s s s s sn n n

s j n s s s j j s jn j n n n j n j

j j jA x P x P x A x A x

++− − −

= = =

⎛ ⎞⎛ ⎞= − − = −⎜ ⎟⎜ ⎟

⎝ ⎠⎝ ⎠∑ ∑ ∑ (3.57)

Se poate deci deduce o relatie de recurenta intre coeficientul ( 1)sjA + al

polinomului din etapa 1s + si coeficientii ( )skA 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 12s j+ din polinomul final;

• Termenul cu puterea 12s j+ din polinomul final mai poate fi obtinut si din inmultirea termenului la puterea 2 ( )s j k− dintr-un polinom initial poate fi inmultit doar cu termenul de aceleasi putere 2 ( )s j k+ , 1,.., 1k j= − din cel de al doilea polinom initial (cazurile 0k = 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 ( 1)s

jA + sub forma:

12( 1) ( ) ( ) ( )

12 ( 1)

js s k s s

j j j k j kk

A A A A−

+− +

=

⎡ ⎤= + −⎣ ⎦ ∑ (3.58)

Alternativ, aceasta relatie se poate scrie si sub o alta forma, echivalenta. Intr-adevar, se poate scrie: ( ) ( 1) ( 1)( ) ( 1) ( ) ( )s m n s m s m

n n nP x P x P x− −= − − (3.59)

( ) ( 1) ( 1) ( 1) ( 1) 2 22 2

0 0 0 0 0

l kl kn n n n n mm ms mj s s s sn j n l n k n l n k

j l k l k

A x A x A x A A x⎛ ⎞+⎜ ⎟− − − − ⎝ ⎠

− − − − −= = = = =

⎡ ⎤⎡ ⎤ ⎡ ⎤= = ⎢ ⎥⎢ ⎥ ⎢ ⎥

⎢ ⎥⎣ ⎦ ⎣ ⎦ ⎣ ⎦∑ ∑ ∑ ∑ ∑ (3.60)

Se observa ca termenul in mjx din primul membru trebuie identificat cu

termenul in 2 2l km

x⎛ ⎞+⎜ ⎟⎝ ⎠ din ultimul membru ceea ce conduce la conditia:

2 2l km mj⎛ ⎞+ =⎜ ⎟

⎝ ⎠ (3.61)

deci 2l j k= − (3.62) iar relatia 129H(3.60) se poate scrie:

( ) ( 1) ( 1)

0 0 0

n n ns mj mj s s

n j n l n kj l k

A x x A A− −− − −

= = =

⎡ ⎤= ⎢ ⎥⎣ ⎦∑ ∑ ∑ (3.63)

ceea ce inseamna ca

( ) ( 1) ( 1)

0

, 1, 2,..., 1n

s s sn j n l n k

k

A A A j n− −− − −

=

= = −∑ (3.64)

Multimea de valori a lui j s-a stabilit tinand cont ca deja se cunoaste ( ) 20 0

ssA A= iar termenii ( )snA sunt termini liberi si deci nu pot ridica gradul monomului

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

Tinand cont de 131H(3.62) rezulta: (2 ) 1 2 1n j k n k j− − ≤ − ⇒ ≤ − (3.66)

De asemena pentru indicile celui de al doilea coeficient se poate scrie: 1 1n k n k− ≤ − ⇒ ≥ (3.67)

Inlocuind pe l dat de 132H(3.62) si tinand cont de limitele de variatie a lui k dat de relatiile 133H(3.66) si 134H(3.67), ecuatia 135H(3.64) devine

2 1

( ) ( 1) ( 1)(2 )

1

, 1, 2,..., 1j

s s sn j n j k n k

k

A A A j n−

− −− − − −

=

= = −∑ (3.68)

La utilizarea oricareia dintre aceste relatii de recurenta se va tine seama ca (1) 2

0 0A a= evidenta conform relatiei 136H(3.50).

Pentru verificarea conditiei de oprire a iteratiilor, s-ar putea compara valoarea functiei in jx sau in jx− , 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 137H(3.56)

Programul are intrarile: • Polinomul [ _]f x definit ca o functie pura; • Numarul maxim de iteratii maxN =25;

Programul 8.6. Calculul radacinilor reale ale unui polinom prin metoda Lobacevski-Graeffe

f@x_D = −13824+ 27648x+11808x2 − 44560x3 − 3392x4 + 30632x5 + 1746x6 −

10869x7 − 1215x8 + 1911x9 + 297x10 − 159x11− 29x12 + 5x13 +x14;H∗ Polinomul ale carui radacini se cauta ∗La= N@CoefficientList@f@xD, xDD; H∗ Lista coeficientilor ∗Ln= Length@aD; H∗ Gradul polinomului +1 ∗LA= Flatten@Append@a, Table@0, 8n<DDD;H∗ Se completeaza lista coeficientilor cu n zerouri ∗LB= A; H∗ Lista auxiliara ∗LX= Table@0, 8n−1<D; H∗ Lista radacinilor ∗LForAm = 1, m < 25, m++, H∗ Iteratiile ∗LForAj= 1, j≤ n, j++,

B@@jDD = H−1Lj ‚k=1

2 j−1H−1Lk+1 A@@2 j− kDD A@@kDD;E; H∗ Coeficientii auxiliari ∗L

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

ForAj= 1, j<= n−1, j++,

r=ikjjAbsA B@@jDD

B@@j+ 1DD Ey{zz

12m−1 ; H∗ calculul radacinilor ∗L

If@Abs@f@rDD <= Abs@f@−rDD, X@@jDD = r, X@@jDD = −rD;H∗ Alegerea semnului ∗LE;

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 138H(3.52) presupune adunarea a M radacini, deci se obtine:

11

0

sm

s

AMxA

= (3.69)

iar radacina se calculeaza cu prima relatie Vieta,dupa formula:

1

11

0

s m

s

AxMA

= ± (3.70)

sau cu relatia Vieta M, dupa formula

1

10

s mMMs

AxA

= ± (3.71)

In cazul radacinilor complexe, acestea fiind conjugate, se poate scrie: 1 2,i ix e x eϕ ϕρ ρ −= = Astfel, dupa s ridicari la patrat, prima relatie Vieta devine:

( )1

3 ( )0

...s

m im m im m mn s

Ae e x xA

ϕ ϕρ ρ −+ + + + = −

deci

( ) ( )

1 1( ) ( )0 0

, 2 coss s

m im m im ms s

A Ae e mA A

ϕ ϕρ ρ ρ ϕ−+ ≈ − ≈ − (3.72)

Astfel, daca se calculeaza rapoartele de tipul 2( 1)

( )

sj

sj

AA

−⎡ ⎤⎣ ⎦ se constata

urmatoarele: • Pentru radacini simple aceste rapoarte tind spre 1. Intr-adevar, ridicand la

patrat ecuatia 139H(3.52) scrisa pentru 1s − se obtine:

( ) ( )2( 1)

2 21 / 2( 1)

sjs m

j j sj

Ax x

A

−−

−

⎡ ⎤= = −⎢ ⎥

⎢ ⎥⎣ ⎦

Pe de alta parte aceasta este egala cu ( )

( )1

sjs

j sj

Ax

A −

= deci:

2( 1) ( )

( 1)1 1

s sj js s

j j

A AA A

−

−− −

⎡ ⎤− =⎢ ⎥⎢ ⎥⎣ ⎦

Deci:2 2( 1) ( 1)

1( )

1

s sj j

s sj j

A AA A

− −−

−

⎡ ⎤ ⎡ ⎤⎣ ⎦ ⎣ ⎦= care poate fi scrisa si pentru 1, 2,...,0j j− −

2 2 2 2 2( 1) ( 1) ( 1) ( 1) ( 1)

1 2 1 0( ) ( )

1 2 1 0

... 1s s s s s

j j js s s s s

j j j

A A A A AA A A A A

− − − − −− −

− −

⎡ ⎤ ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ ⎡ ⎤⎣ ⎦ ⎣ ⎦ ⎣ ⎦ ⎣ ⎦ ⎣ ⎦= = = = = = (3.73)

unde s-a tinut cont ca

2( ) ( 1)0 0

s sA A −⎡ ⎤= ⎣ ⎦ Rezulta :

2( 1)

( ) 1s

jj s

j

Ar

A

−⎡ ⎤⎣ ⎦= ≈ .

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( 1)

( )

sj

j sj

Ar M

A

−⎡ ⎤⎣ ⎦= ≈

Demonstratia este identica, dar in locul relatiei 140H(3.52), se pleaca de la relatia 141H(3.69). Deci, daca prin calculul acestui raport pentru radacina j se obtine o 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 142H(3.72), putem scrie:

( )

( ) ( )11 2 1 0( )

0

2 cos ; 2 coss

m m m s s ms

Ax x m A A mA

ρ ϕ ρ ϕ+ = = ⇒ = (3.74)

( )2( 1)2 2 22 ( 1) ( 1) 212 2

1 2 1 0( 1)0

4 cos ; 4 cos2 2

sm m m s s ms

Am mx x A AA

ρ ϕ ρ ϕ−

− −−

⎡ ⎤⎡ ⎤ ⎡ ⎤+ = = ⇒ =⎢ ⎥ ⎣ ⎦ ⎣ ⎦

⎣ ⎦

Deci, acest raport este:

22( 1)

( )

2cos2

cos

sj

j sj

mA

rA 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: 1 2 2

1 2 2 2 0( ) ....n n nn n n nP z a z a z a z a z a z a− −

− −= + + + + + + (3.75) Se incearca impartirea cu un trinom de gradul doi de forma 2z pz q− − , care are radacinile:

2

1

2

2

42

42

p p qz

p p qz

+ +=

− +=

(3.76)

care pot fi reale sau complex conjugate dupa cum 2 4p q+ 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 0p si 0q , 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 urmatoarea0F

1: 2

0 0 2 1 0 0( ) ( ) ( ) ( )n nP z z p z q Q z r p z r−= − − + − + (3.77) unde 2 ( )nQ z− este un polinom de grad 2n − reprezentand catul impartirii 2 3 4 2 2

2 1 2 4 3 2( ) ... ...n n n kn n n n kQ z r z r z r z r z r z r z r− − − −− − −= + + + + + + + (3.78)

iar 1 0 0( )r p z r− + este un polinom de gadul intai reprezentand restul impartirii. Daca

0p si 0q ar avea valorile corecte, impartirea s-ar face exact, restul fiind nul. Prin urmare, problema se reduce la minimizarea (anularea) coeficientilor restului, 1r si 0r . Deoarece acesti coeficienti depind de valorile coeficientilor p si q

0 0

1 1

( , )( , )

r r p qr r p q==

ei se pot dezvolta in serie Taylor in jurul unui punct ales initial 0 0( , )p q , 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:

0 00 0 0 0 0 0

1 11 1 0 0 0 0

( , ) ( , ) ( ) ( ) 0

( , ) ( , ) ( ) ( ) 0

r rr p q r p q p p q q

p qr r

r p q r p q p p q qp q

∂ ∂= + − + − =

∂ ∂∂ ∂

= + − + − =∂ ∂

(3.79)

Prin urmare, se pot calcula valorile corecte p si q din sistemul de ecuatii 143H(3.79) in functie de valorile initiale 0p si 0q si de elementele matricii Jacobiene:

0 0

1 1

r rp q

Jr rp q

∂ ∂⎡ ⎤⎢ ⎥∂ ∂⎢ ⎥=∂ ∂⎢ ⎥⎢ ⎥∂ ∂⎣ ⎦

(3.80)

calculate in punctele initiale 0p si 0q . Intr-adevar, ecuatia 144H(3.79) se poate scrie matricial sub forma:

0 0

0 0 0 0

0 1 0 01 1

( , )( , )

r rp p r p qp qq q r p qr r

p q

∂ ∂⎡ ⎤⎢ ⎥ − −∂ ∂ ⎡ ⎤ ⎡ ⎤⎢ ⎥ ⋅ =⎢ ⎥ ⎢ ⎥− −∂ ∂⎢ ⎥ ⎣ ⎦ ⎣ ⎦⎢ ⎥∂ ∂⎣ ⎦

(3.81)

de unde rezulta noile valori ale lui p si q :

0 0 0 01

0 1 0 0

( , )( , )

p r p qpJ

q r p qq− −⎡ ⎤ ⎡ ⎤⎡ ⎤

= + ⋅⎢ ⎥ ⎢ ⎥⎢ ⎥ −⎣ ⎦ ⎣ ⎦ ⎣ ⎦ (3.82)

Deoarece in general nu se pot obtine dintr-un singur pas chiar valorile corecte ale lui p si q , relatia 145H(3.79) fiind doar aproximativa, procedeul se reia si se obtin prin

1 Uneori se foloseste forma mai simpla si mai naturala 1 0r z r+ 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 146H(3.82) este deci calculata in mod repetat, sub forma:

1 01

1 1

( , )( , )

j j j jj

j j j j

p p r p qJ

q q r p q+ −

+

−⎡ ⎤ ⎡ ⎤ ⎡ ⎤= + ⋅⎢ ⎥ ⎢ ⎥ ⎢ ⎥−⎣ ⎦ ⎣ ⎦ ⎣ ⎦

(3.83)

astfel ca, daca procedeul converge, la iteratia 1j + 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 147H(3.78) in 148H(3.77). Rezulta:

1 2

1 2 1

1 2 1 2 3 0 1 2

( ) ( ) ( ) ...

( ) ... ( ) ( )

n n nn n n n n n n

kk k k

P z z r z r pr z r pr qr

z r pr qr z r pr qr r pr qr

− −− − −

+ +

= + − + − − +

+ − − + + − − + − − (3.85)

Prin identificare relatiei 149H(3.85) cu 150H(3.75) se obtin relatiile de recurenta prin care se pot calcula coeficientii polinomului cat:

1 1

1 2 , 0... 2

n n

n n n

k k k k

r ar a prr a pr qr k n− −

+ +

== +

= + + = − (3.86)

Pentru determinarea derivatelor partiale din componenta Jacobianului, se vor face urmatoarele notatii:

1

2

kk

kk

rc

pr

dq

+

+

∂≡∂∂

≡∂

(3.87)

Acum, prin derivarea in raport cu p a relatiilor 151H(3.86) si tinand seama de definitiile 152H(3.87), obtinem succesiv:

0n nr ap p

∂ ∂= =

∂ ∂

1n nn n

r rc p r

p p−∂ ∂

≡ = +∂ ∂

(3.88)

2 3

1 21 1 , 0... 2

k k

k k kk k

c c

r r rc r p q k np p p

+ +

+ ++ +

∂ ∂ ∂≡ = + + = −∂ ∂ ∂

(3.89)

Din relatia 153H(3.88) si din 154H(3.89), renotand 1k k+ → rezulta relatiile de recurenta din care se pot determina elementele Jacobianului:

1 2 , 1... 1

n n

k k k k

c rc r pc qc k n+ +

== + + = −

(3.90)

Se observa deplasarea cu o unitate a limitelor indicelui k , datorita definitiei

155H(3.87). Pentru Jacobian este nevoie de valorile 01

rcp

∂≡∂

si 12

rcp∂

≡∂

care sunt ultima si

repectiv penultima derivata calculata prin recurenta 156H(3.90). Pentru determinarea celorlalte doua elemente, nu mai sunt necesare calcule suplimentare, deoarece se poate demonstra ca, tinand cont de definitiile 157H(3.87), avem relatiile:

02 2

13 3

rd cqrd cq

∂≡ =∂∂

≡ =∂

(3.91)

Intr-adevar, daca derivam relatiile 158H(3.86) in raport cu q si tinem seama de definitiile 159H(3.87), obtinem succesiv:

0n nr aq q

∂ ∂= =

∂ ∂

1 1 0n nr aq q− −∂ ∂= =

∂ ∂

2 1n n nn n

r r rd p q rq q q− −∂ ∂ ∂

≡ = + +∂ ∂ ∂

(3.92)

3 4

1 22 2 , 0... 2

k k

k k kk k

d d

r r rd p r q k np q q

+ +

+ ++ +

∂ ∂ ∂≡ = + + = −∂ ∂ ∂

(3.93)

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

1 2 , 2... 2

n n

k k k k

d ad r pd qd k n+ +

== + + = −

(3.94)

Se observa ca relatiile 162H(3.94) si 163H(3.90) sunt similare, astfel ca kc si kd avand acelasi punct de plecare ( na ) si aceleasi relatii de recurenta, vor avea valori egale, pentru 2...k n= , ceea ce demonstreaza relatiile 164H(3.91). Jacobianul va avea deci forma:

1 2

2 3

c cJ

c c⎡ ⎤

= ⎢ ⎥⎣ ⎦

(3.95)

Din relatia 165H(3.83), tinand seama de regula lui Cramer, solutia sistemului va putea fi scrisa sub forma:

0 2

1 31

1 2

2 3

( , )( , )

j j

j jj j

r p q cr p q c

p pc cc c

+

−−

= + (3.96)

1 0

2 11

1 2

2 3

( , )( , )

j j

j jj j

c r p qc r p q

q qc cc c

+

−−

= + (3.97)

Iteratiile pot fi oprite atunci cand atat variatiile lui p cat si ale lui q in relatiile 166H(3.96) si 167H(3.97) scad sub o limita de precizie impusa. Cunoscand valorile corecte ale lui p si q , din relatia 168H(3.76) se determina radacinile cele mai apropiate de cele determinate de valorile alese initial 0p si 0q , iar din coeficientii kr se poate determina polinomul cat 2 ( )nQ 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 2 ( )nQ 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 2 ( )nQ 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 169H(3.86), 170H(3.90), 171H(3.95)-172H(3.97) si 173H(3.76). Dupa calculul 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∗LPlot@f@xD, 8x, −1.5, 2<, AxesOrigin→ 80, 0<, PlotRange → 8−10, 10<D; H∗ Trasarea graficului ∗La= N@CoefficientList@f@xD, xDD; H∗ Lista coeficientilor, cu evitarea calculelor cu intregi ∗LPrint@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 ∗Lmax= 200; H∗ Numarul maxim admis de iteratii ∗Lε = 10−8; H∗ Eroarea admisa ∗Lj= 0;X= 8<; H∗ Initializarea matricii radacinilor ∗LWhileAn> 3, H∗ Se va reduce gradul polinomului pana cand vor fi 3 sau 2 termeni ∗LWhileAj< max,

n= Length@aD− 1;rPn+1T = aPn+1T;rPnT = aPnT + p0rPn+1T;For@k= n− 1, 1 ≤ k, k−−,rPkT = aPkT + p0rPk+1T + q0rPk+2T;D; H∗ Calculul coeficientilor polinomului cat ∗L

cPn+1T = rPn+1T;cPnT = rPnT + p0cPn+1T;For@k= n− 1, 2 ≤ k, k−−,cPkT = rPkT + p0cPk+1T + q0cPk+2T;D; H∗ Calculul derivatelor care vor da matricea Jacobiana ∗L

d0 = DetAikjjcP2T cP3T

cP3T cP4Ty{zzE; d1 = DetAi

kjj −rP1T cP3T

−rP2T cP4Ty{zzE; d2 = DetAi

kjjcP2T −rP1T

cP3T −rP2Ty{zzE;

εp=d1d0

;εq =d2d0

; H∗ Corectorii la valorile anterioare ale sumei si produsului radacinilor ∗L

p1 = p0 + εp; q1 = q0 + εq; H∗ Noile valori ale sumei si produsului radacinilor ∗Lj = j+ 1;If@Abs@εpD < ε && Abs@εqD < ε, Break@DD; H∗ Oprirea daca ambii corectori sunt suficient de mici ∗Lp0 = p1; H∗ Noile valori ale sumei si produsului radacinilor devin valori anterioare∗Lq0 = q1; E;

x1=p1+

è!!!!!!!!!!!!!!!!!!p12 + 4 q12

;x2 =p1−

è!!!!!!!!!!!!!!!!!!p12 +4 q12

; H∗ Cele 2 radacini ∗L

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

H∗ S−au eliminat partile imaginare daca sunt prea mici∗LAppendTo@X, 8x1, f@x1D<D; H∗ Daca s−a atins precizia, caluleaza radainile∗LAppendTo@X, 8x2, f@x2D<D;H∗si adauga in lista impreuna cu valorile functiei∗Lp1= r@@nDD;q1 = r@@n− 1DD;n−−; H∗ Se scade ordinul polinomului ∗La= 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 q12

; x2=−p1−

è!!!!!!!!!!!!!!!!!!p12 − 4 q12

;

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

H∗ Daca polinomul cat are gradul 2 se calculeaza cele doua radacini ∗LAppendTo@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∗LPrint@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:

1

( ) ( )n

n kk

P x x x=

= −∏ (3.98)

in care s-a impartit cu coeficientul nenul al termenului cu cea mai mare putere, ceea ce nu modifica ecuatia ( ) 0nP x = . Logaritmul natural al polinomului este:

1

ln ( ) ln( )n

n kk

P x x x=

= −∑ (3.99)

Derivata intai a acestui logaritm o notam cu G si este, conform relatiei precedente

1 2

1 1 1ln ( ) ...nn

d P x Gdx x x x x x x

= + + + ≡− − −

(3.100)

iar derivate a doua a logaritmului o notam cu H si este

( ) ( ) ( )

2

2 2 221 2

1 1 1ln ( ) ...nn

d P x Hdx x x x x x x

= + + + ≡− − −

(3.101)

Tinand cont de membrul stang al relatiilor 174H(3.100) si 175H(3.101) putem scrie pe G si pe H sub formele:

'( )( )

n

n

P xGP x

= (3.102)

2

' ' ''( ) ( ) ( ) ( )( )

n n n n

n

P x P x P x P xHP x−

= (3.103)

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 1x , si la distante egale b fata de toate celelalte radacini, deci 1 , , 2,3,...,ix x a x x b i 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 176H(3.104) in 177H(3.100) si 178H(3.101) rezulta sistemul de ecuatii

2 2

1 1

1 1

nGa b

nHa b

−⎧ = +⎪⎪⎨ −⎪ = +⎪⎩

(3.105)

Rezolvand acest sistem in raport cu a si b , rezulta corectorul

2( 1)( )

naG n nH G

=± − −

(3.106)

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 kx aproximand radacini 1x , la iteratia urmatoare avem o valoare mai apropiata, corectata cu valoarea ka . Considerand 1 1kx x+ ≈ , conform relatiei 179H(3.104) avem 1k k kx x a+ = − (3.107) Semnul din ecuatia 180H(3.106) se alege astfel incat numitorul sa aiba valoarea cea 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 181H(3.102), 182H(3.103), 183H(3.106) si 184H(3.107) si are urmatoarele 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 2x 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 2x 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 ∗LNr= 200; H∗ Nr maxim de puncte de start pentru iteratii ∗LNi= 200; H∗ Nr maxim de iteratii la fiecare punct ∗LA= N@CoefficientList@f@xD, xDD; H∗ Coeficientii polinomului initial ∗LX= 8<; H∗ Multimea radacinilor ∗LXa= 8<; H∗ Multimea radacinilor rotunjite cu 2 cifre ∗Ln= Length@AD− 1; H∗ Gradul polinomului ∗Lj= 0; H∗ Nr de puncte de start ∗Li= 0; H∗ Nr total de iteratii ∗LWhileALength@XD < n, H∗Repeta pana cand nr radacinilor este egal cu gradul polinomului∗Lj++;If@j> Nr, Print@"Peste ", Nr, " de incercari"D; Break@DD;H∗Daca sunt prea multe incercari iese∗Lx1= 10 HHRandom@D− 0.5L + IHRandom@D −0.5LL; H∗ Puncte de start complexe aleatoare ∗LForAk= 1, k≤ Ni, k++, H∗ Iteratiile pentru fiecare punct de start ∗Li++;If@Abs@f@x1DD < ε, Break@DD; H∗ Se intrerup daca functia e suficient de mica ∗L

fd@x_D = ‚k=1

nA@@k+ 1DD kxk−1; H∗ Derivata intai a polinomului ∗L

fdd@x_D = ‚k=1

n−1A@@k+ 2DD k Hk+ 1Lxk−1; H∗Derivata a doua a polinomului ∗L

G=fd@x1Df@x1D ;

H=ikjjikjjfd@x1D

f@x1Dy{zz

2−

fdd@x1Df@x1D

y{zz;

a1=n

G−è!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!Hn−1L Hn H− G2L

;

a2=n

G+è!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!Hn−1L Hn H− G2L

; H∗ cele 2 posibilitati de alegere a corectorului ∗L

If@Abs@a1D <= Abs@a2D, a= a1, a= a2D; H∗ Se alege cel care are numitorul mai mare ∗Lx1= x1− a;H∗ Noua valoare a radacinii dupa corectie ∗LE;

x2= AbsANA10−2 RoundA102 x1EEE; H∗ Radacina rotunjita cu doua cifre ∗LIf@! MemberQ@Xa, x2D, H∗ Testare daca radacina trunchiata a mai fost gasita ∗LAppendTo@Xa, x2D;If@Abs@Im@x1DD > ε, AppendTo@X, x1D, AppendTo@X, Re@x1DDD;D;E;

H∗ Daca radacina rotunjita nu este in lista Xa,

atunci se introduce radacina nerotunjita

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

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 ∗LPrint@"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.8592360.435058+ 0.9180460.438331− 0.7799510.87771+ 0.1349091.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

Raphson Muller Lobacevski

Graeffe Bairstow Laguerre

Functie Oarecare Oarecare Oarecare Oarecare Polinom Polinom Polinom Radacini complexe

Nu Nu Da Da Nu Da Da

Multiplicitate radacini

Nu Nu Nu Nu Da Da Nu

Puncte start Interval Interval 1 3 Nu 2 1 Convergenta Minima Medie Foarte

mare Mare Medie Mare Mare

Dezavantaj principal

Convergenta mica

Radacini reale

Necesita derivata

Start complicat

Precizie redusa

Erori la deflatie

Necesita derivata

Utilizare recomandata

Braketing Refining Refining Refining Braketing Refining Refining

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 1 2( , )x x din planul 1 20x x pentru care doua functii de aceste variabile 1 1 2( , )F x x respectiv

2 1 2( , )F x x sa se anuleze simultan. Cele doua functii reprezinta suprafete intr-un spatiu tridimensional, ale caror intersectii cu planul 1 20x x sunt curbe formate de punctele de anulare ale functiilor. Solutiile sistemului de doua ecuatii sunt punctele de intersectie ale acestor curbe in planul 1 20x x . 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 1 2( , ,..., )nx x x pentru care functiile 1 1 2( , ,..., )nF x x x , 2 1 2( , ,..., )nF x x x ,..., 1 2( , ,..., )n nF x x x sa se anuleze simultan:

1 1 2

2 1 2

1 2

( , ,..., ) 0( , ,..., ) 0

..............................( , ,..., ) 0

n

n

n n

F x x xF x x x

F x x x

=⎧⎪ =⎪⎨⎪⎪ =⎩

(4.1)

Acest sistem poate fi scris matriceal sub forma: ( ) 0=F x (4.2) unde 1 2( , ,..., )nF F F=F este vectorul coloana al valorilor functiei iar 1 2( , ,..., )nx x x=x 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 185H(4.2) si facandu-se dezvotarea in serie Taylor a functiei din membrul stang in jurul unui punct intr-un spatiu n-dimensional.

2

1

( )( ) ( ) ( ) , 1, 2,...,n

ii i

j j

FF F i nx

δ δ δ=

∂+ = + + =

∂∑ xx x x x x (4.3)

Deoarece derivatele partiale ale tuturor componentelor functiei in raport cu variabilele formeaza matricea Jacobiana

( )iij

j

FJx

∂=

∂x (4.4)

ecuatia 186H(4.3) se poate scrie sub forma matriceala: 2( ) ( ) ( )δ δ δ+ = + ⋅ +F x x F x J x x (4.5) Deoarece se urmareste ca dupa corectie ( )δ+F x x sa fie cat mai aproape de 0, din 187H(4.5) rezulta corectorul: 1 ( )δ −= − ⋅x J F x (4.6) Metoda se aplica iterativ, astfel ca la iteratia k+1 se obtine noul vector radacina: 1k k kδ+ = +x x x (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: ( ) cos[ arccos( )]nT x n x= (5.1) se mai definesc si polinoamele Chebyshev de speta 2, prin formula:

sin[( 1)arccos( )]( )sinn

n xU xx

+= (5.2)

S-a aratat ca una dintre metodele de interpolare a functiilor se bazeaza pe descompunerea acestora intr-o serie de N polinoame Chebyshev:

1

0( ) ( )

N

n nn

f x c T x−

=

= ∑ (5.3)

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

1

21

0,( ) ( ) , 01

, 02

n m

m nT x T x dx m n

xm n

ππ−

⎧⎪ ≠⎪

= = =⎨− ⎪

⎪ = ≠⎩

∫ (5.4)

Inmultind relatia 188H(5.3) cu ( )mT x integrand si tinand cont de ralatia 189H(5.4), rezulta coeficientii dezvoltarii ca produse scalare ale functiei cu vectorii ortogonali

( )nT x :

1

21

( ) ( )11

, 0

, 02

kk

k

k

f x T xc dxx

k

k

α

πα π

−

=−

=⎧⎪= ⎨

≠⎪⎩

∫ (5.5)

Aceasta relatie nu este insa practica deoarece necesita o integrare numerica, introducand erori suplimentare si consumand timp de calcul. Coeficientii dezvoltarii 190H(5.3) se pot calcula insa si fara a fi necesara integrarea, tinand seama de relatia de ortogonalitate discreta a polinaomelor Chebyshev:

1

0,

( ) ( ) , 02

, 0

N

n k m kk

n mNT x T x n m

N n m=

≠⎧⎪⎪= = ≠⎨⎪

= =⎪⎩

∑ (5.6)

unde kx sunt radacinile polinomului Chebyshev de ordinul N , obtinute din relatia de definitie: ( ) cos( arccos ) 0N k kT x N x= = Rezulta

arccos (2 1)

21cos2

k

k

x kN

x kN

π

π

= −

⎡ ⎤⎛ ⎞= −⎜ ⎟⎢ ⎥⎝ ⎠⎣ ⎦

(5.7)

Intr-adevar, membrul stang al relatiei 191H(5.6) se scrie succesiv:

1

1 1

1 1

1 1( ) ( )2

1

1 1( ) ( ) cos cos2 2

1 1 1 1cos ( ) cos ( )2 2 2 2

1 1Re Re2 2

N N

n k m kk k

N N

k k

N i n m k i n m kN

k

S

S T x T x n k m kN N

n m k n m kN N

e eπ

π π

π π= =

= =

⎛ ⎞+ − − −⎜ ⎟⎝ ⎠

=

⎡ ⎤ ⎡ ⎤⎛ ⎞ ⎛ ⎞= = − − =⎜ ⎟ ⎜ ⎟⎢ ⎥ ⎢ ⎥⎝ ⎠ ⎝ ⎠⎣ ⎦ ⎣ ⎦⎡ ⎤ ⎡ ⎤⎛ ⎞ ⎛ ⎞+ − + − − =⎜ ⎟ ⎜ ⎟⎢ ⎥ ⎢ ⎥⎝ ⎠ ⎝ ⎠⎣ ⎦ ⎣ ⎦⎡ ⎤

+⎢ ⎥⎢ ⎥⎣ ⎦

∑ ∑

∑ ∑

∑2

2

1

NN

k

S

π⎛ ⎞⎜ ⎟⎝ ⎠

=

⎡ ⎤⎢ ⎥⎢ ⎥⎣ ⎦

∑

(5.8)

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

termenul 1S are primul termen ( )

21

i n mNa eπ

+= si ratia

( )i n mNr eπ

+= , de unde rezulta:

( )( )

21

( )

1 1Re2

1

i n mi n mN

i n mN

eS ee

π π

π

++

+

⎡ ⎤−⎢ ⎥=⎢ ⎥

−⎣ ⎦

Impartind la ( )

2i n m

Neπ

+ se poate obtine la numitor un sinus, si tinand cont ca

( ) ( 1)i n m n me π+ += − , rezulta:

11 1 ( 1)Re2 2 sin( )

2

n miSn m

Nπ

+⎡ ⎤⎢ ⎥− −

= ⎢ ⎥⎢ ⎥+⎣ ⎦

(5.9)

Similar se obtine si celalalt termen:

21 1 ( 1)Re2 2 sin( )

2

n miSn m

Nπ

−⎡ ⎤⎢ ⎥− −

= ⎢ ⎥⎢ ⎥−⎣ ⎦

(5.10)

Se observa ca avem trei cazuri:

a) Daca 0n m= = , ambele sume din ultimul membru al relatiei 192H(5.8) au exponentul nul si deci suma totala este N ;

b) Daca 0n m= ≠ , prima suma va fi nula conform relatiei 193H(5.9) in care in paranteza dreapta avem un numar pur imaginar iar a doua suma se poate calcula din termenul final al relatiei 194H(5.8) in care exponentul e nul, suma totala

fiind deci 2N ;

c) Daca n m≠ , si prima si a doua suma sunt nule conform relatiilor 195H(5.9) si 196H(5.10), in care in parantezele dreapte avem numere pur imaginare, deci suma totala este nula. Relatia 197H(5.6) este deci demonstrata si poate fi folosita pentru calculul

coeficientilor dezvoltarii 198H(5.3). Pentru inceput se scrie aproximarea functiei initiale prin polinoame Chebyshev conform relatiei 199H(5.3), pentru punctele kx ,

1

0

( ) ( )N

k n n kn

f x c T x−

=

= ∑ (5.11)

Inmultind aceasta relatie cu ( )m kT x , sumand dupa k si tinand cont de relatia de ortogonalitate discreta 200H(5.6), obtinem:

1

1 1 0

01

0 1

( ) ( ) ( ) ( )

, 0( ) ( )

, 02

N N N

k m k n n k m kk k n

N N

n n k m kn k m

f x T x c T x T x

Nc mc T x T x N c m

−

= = =

−

= =

= =

=⎧⎪= ⎨

≠⎪⎩

∑ ∑∑

∑ ∑

Coeficientii dezvoltarii vor fi deci:

0 01

1 ( ) ( )N

k kk

c f x T xN =

= ∑ (5.12)

1

2 ( ) ( ), 1, 2,..., 1N

m k m kk

c f x T x m NN =

= = −∑ (5.13)

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

1

2 1 1cos cos , 1, 2,..., 12 2

N

mk

c f k m k m NN N N

π π=

⎧ ⎫⎡ ⎤ ⎡ ⎤⎛ ⎞ ⎛ ⎞= − − = −⎨ ⎬⎜ ⎟ ⎜ ⎟⎢ ⎥ ⎢ ⎥⎝ ⎠ ⎝ ⎠⎣ ⎦ ⎣ ⎦⎩ ⎭∑ (5.14)

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 202H(5.12) si 203H(5.13)) pentru calculul coeficientilor dezvoltarii, in practica se prefera o forma usor modificata pentru relatia 204H(5.11):

1

01

1( ) ( )2

N

n nn

f x c T x c−

=

≈ +∑ (5.15)

in care toti coeficientii, inclusiv 0c , se calculeaza cu aceeasi relatie 205H(5.13), dar se tine seama ca din aceasta relatie 0c rezulta de doua ori mai mare decat din 206H(5.12); de aceea se separa termenul respectiv in relatia 207H(5.15) si se incepe indexarea in suma de la 1. S-a tinut cont si de faptul ca 0 ( ) 1T x = . Derivata functiei ( )f x se poate calcula din relatia 208H(5.3), in care functiile dependente de x din membrul drept sunt polinoamele Chebyshev:

1

'

0'( ) ( )

N

n nn

f x c T x−

=

= ∑ (5.16)

cu coeficientii calculati dupa 209H(5.12) si 210H(5.13), sau rescrisa cu coeficienti calculati numai dupa 211H(5.13), sub o forma similara cu 212H(5.15):

1

'0

1

1'( ) ( )2

N

n nn

f x c T x c−

=

= +∑ (5.17)

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

Se pot demonstra urmatoarele relatii:

1( ) ( )n nd T x nU xdx −= (5.18)

[ ]21( ) ( ) ( )2n n nT x U x U x−= − (5.19)

Prima relatie se demonstreaza prin calcul direct:

2

12

sin( arccos )cos( arccos )1

sin( arccos ) sin( arccos ) ( )sin(arccos )1 [cos(arccos )]

n

d n n xn xdx x

n x n xn n nU xxx

−

= =−

= = =−

A doua relatie se demonstreaza pe baza identitatii trigonometrice:

[ ]1cos sin sin( ) sin( )2

a b a b a b= + − −

Notand arccos xθ = si alegand , ( 1)a n b nθ θ= = − se obtine succesiv:

[ ]1cos sin sin( 1) sin( 1)2

n n nθ θ θ θ= + − −

[ ] [ ]

[ ] [ ]

sin ( 1) sin ( 1)1cos( arccos )2 sin

sin ( 1)arccos sin ( 1)arccos12 sin(arccos ) sin(arccos )

n nn x

n x n xx x

θ θθ θ

⎧ ⎫+ −= − =⎨ ⎬

⎩ ⎭⎧ ⎫+ −

= −⎨ ⎬⎩ ⎭

obtinandu-se relatia 213H(5.19). Rescriind relatia 214H(5.18) pentru 1n + si 1n − se obtine

1

1 2

( ) ( 1) ( )

( ) ( 1) ( )

n n

n n

d T x n U xdxd T x n U xdx

+

− −

= +

= −

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

seama de relatia 215H(5.19), rezulta:

' '

1 1( ) ( )2 ( )1 1

n nn

T x T xT xn n+ −= −+ −

(5.20)

Aceasta relatie este fundamentala pentru metodele spectrale, permitand calculul prin recurenta al coeficientilor seriilor derivate. Astfel, se pot exprima derivatele polinoamelor Chebyshev din relatia 216H(5.16) in functie de polinoamele nederivate, obtinandu-se o noua serie cu alti coeficienti nb care pot fi calculati prin recurenta din coeficientii deja cunoscuti nc . Pentru gasirea relatiei intre cele doua seturi de coeficienti, se scrie 217H(5.16) in cele doua forme si se egaleaza coeficientii polinoamelor Chebyshev de acelasi ordin.

1 1

'

0 0'( ) ( ) ( )

N N

n n n nn n

f x c T x b T x− −

= =

= =∑ ∑ (5.21)

Din relatia 218H(5.20) rezulta derivata polinomului Chebyshev de ordin 1n + in functie de polinomul nederivat de ordin n :

' '1 1

1( ) ( ) 2( 1) ( )1n n n

nT x T x n T xn+ −+

= + +−

care rescrisa trecand de la 1n + la n este:

' '2 1( ) ( ) 2 ( )

2n n nnT x T x nT x

n − −= +−

(5.22)

Tot din relatia 219H(5.20) rezulta derivata polinomului Chebyshev de ordin 1n − in functie de polinomul nederivat de ordin n :

' '1 1

1( ) ( ) 2( 1) ( )1n n n

nT x T x n T xn− +−

= − −+

care rescrisa trecand de la 1n − la n este:

' '2 1( ) ( ) 2 ( )

2n n nnT x T x nT x

n + += −+

(5.23)

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

relatiei 222H(5.21) contribuie termenii in 1( )nT x− si 1( )nT x+ ai sumei din dreapta acestei relatii, avand coeficientii 1nb − si respectiv 1nb + : '

1 1 1 1( ) ( ) ( )n n n n n nc T x b T x b T x− − − −= + (5.24) Tot din relatiile 223H(5.22) si 224H(5.23) rezulta ca factorii de proportionalitate intre

1( )nT x− si ' ( )nT x , respectiv 1( )nT x+ si ' ( )nT x sunt 12n

, respectiv 12n

− , deci inlocuind in

225H(5.24) obtinem:

' ' '1 1( ) ( ) ( )2 2n n

n n n nb bc T x T x T x

n n− −= −

Rezulta relatia de recurenta pentru calculul coeficientilor dezvoltarii derivatei: 1 1 2n n nb b nc− += + (5.25) Ei se calculeaza cu conditiile de start 1 0nc + = si 0nc = , iar deoarece 0b se calculeaza din 1c , el rezulta cu o valoare dubla fata de cea normala si trebuie injumatatit. Se obtine deci formula de interpolare pentru derivata functiei:

1

01

1'( ) ( )2

N

n nn

f x b T x b−

=

≈ +∑ (5.26)

coeficientii fiind calculati prin recurenta cu relatia 226H(5.25). Egalitatea este exacta numai in radacinile polinomului Chebyshev de ordinul 1N − , 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 min max[ , ]x x mai mare, trebuie facuta o schimbare de variabila dupa formula:

max min max min

2 2x x x xx y − +

= + (5.27)

aplicata pentru punctele de esantionare ale functiei, astfel incat relatia 227H(5.13) devine:

max min max min

1

2 ( ) ( ), 0, 2,..., 12 2

N

m k m kk

x x x xc f x T x m NN =

− += + = −∑

La calculul functiei de interpolare trebuie facuta trecerea inversa in variabila polinoamelor Chebyshev, astfel incat formula 228H(5.15) devine:

1

min0

1 max min

1( ) (2 1)2

N

n nn

x xf x c T cx x

−

=

−≈ − +

−∑

De asemenea, la calculul functiei de interpolare a derivatei trebuie tinut cont de factorul care apare la derivarea relatiei de schimbare de variabila 229H(5.27), astfel ca relatia 230H(5.26) devine:

1

min0

1max min max min

2 1'( ) (2 1)2

N

n nn

x xf x b T bx x x x

−

=

−≈ − +

− −∑

In[1]:= Needs@"Graphics̀ Colors̀ "D;

f@x_D =1

1+ x2+ CosAx2E;H∗Cuprinde functia Runge − interpolarea echidistanta nu converge! ∗L

xm = −5; xM= 5; H∗ Domeniul de definitie, extins in afara lui @−1,1D ∗Lg1= 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 ∗LForAk= 0, k<= Na, k++,

c@kD =2.Na

„j=1

Na

fAi

kjjjjCosA

H2 j− 1L π

2

NaEy

{zzzz

xM− xm2

+xM+ xm

2E i

kjjjjCosA

kH2 j− 1L π

2

NaEy

{zzzzE;

H∗ S−a facut schimbarea de variabila x=y xM−xm2

+ xM+xm2

pentru extinderea intervalului ∗L

H∗p@x_D=‚n=1Na c@nD ChebyshevTAn,2 x−xm

xM−xm−1E+ 1

2 c@0D;

Genereaza p HxL in forma polinomiala dar cu erori mai mari pentru Na mare∗L

p@x_D = ‚n=1

Nac@nD CosAnArcCosA2

x−xmxM− xm

− 1EE +12

c@0D;

H∗ Forma preferabila in Mathematica 5.2, erori mai mici decat cu functia ChebyshevT@n,xD ∗LH∗ S−a facut schimbarea inversa de variabila y=2 x−xm

xM−xm−1 ∗L

g2= Plot@p@xD, 8x, xm, xM<, PlotStyle → GreenD;Show@g1, g2D; H∗ Afiseaza functia initiala si aproximanta ∗L

H∗ Aici incepe derivarea numerica ∗Lb@Na+ 1D = 0;b@NaD = 0;For@k= Na, k> 0, k−−,

b@k− 1D = b@k+ 1D +2 kc@kDD;

dp@x_D =2

xM− xm i

kjjjj‚

n=1

Nab@nD CosAnArcCosA2

x− xmxM− xm

− 1EE +12

b@0Dy

{zzzz;

H∗ S−a inmultit cu dydx

din schimbarea de variabila ∗L

df@x_D = ∂x f@xD; H∗ Derivata analitica pentru comparatie ∗Lg3= 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 xM−xm

NmE;

e=1

Nm ‚

k=0

Nm−1Abs@df@r@kDD − dp@r@kDDD;

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 53.2 10−× .

-4 -2 2 4

1.5

2

2.5

3

3.5

4

-4 -2 2 4

1.5

2

2.5

3

3.5

4

Figura 5.1 Graficul functiei ( )f x si al polinomului de aproximare ( )p x

-4 -2 2 4

-10

-5

5

10

-4 -2 2 4

-10

-5

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

a

f x dx∫ (6.1)

6.1 Integrarea Gauss

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

( )b

a

f x dx∫ (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:

1

( ) ( ) ( )n

i ii

f x f x L x=∑ (6.3)

unde polinomul Lagrange are forma cunoscuta:

1

( )n

ii

j j ij i

x xL xx x=

≠

−=

−∏ (6.4)

si deoarece in 231H(6.4) exista conditia ca j i≠ , acest polinom este de grad 1n − . Introducand 232H(6.4) in 233H(6.3) si apoi in 234H(6.2), rezulta:

1 1

( ) ( )b b nn

ii

i j j ia aj i

x xf x dx f x dxx x= =

≠

−−∑ ∏∫ ∫ (6.5)

Rezulta ca prin aceasta formula se va aproxima integrala functiei print integrala unui polinom de grad 1n − , care trece prin n puncte ix . 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 ix care a fost dedusa la studiul interpolarii polinomiale:

( )

11

( )( ) ( )!

n n

n ii

fR x x xnξ

−=

= −∏ (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 2 1n − 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 2 1n − , 2 1( )nf x− . Aceasta functie o factorizam cu un polinom de grad n , ( )nP x , obtinandu-se un cat ( )q x de grad (2 1) 1n n n− − = − si un rest ( )r x de asemenea de grad 1n − : 2 1 1 1( ) ( ) ( ) ( )n n n nf x P x q x r x− − −= + (6.7)

Daca impunem ca polinomul ( )nP x sa faca parte dintr-o baza de polinoame

{ }( )nP x ortogonale pe intervalul [ ],a b , el satisface o relatie de ortogonalitate de tipul:

2( ) ( ) ( )b

n m n nma

P x P x w x dx c δ=∫ (6.8)

unde ( )w x este functia pondere ale bazei de polinoame ortogonale considerate. Inmultind relatia 235H(6.7) cu ( )w x si integrand pe intervalul [ ],a b obtinem:

2 1 1 1( ) ( ) ( ) ( ) ( ) ( ) ( )b b b

n n n na a a

f x w x dx P x q x w x dx r x w x dx− − −= +∫ ∫ ∫ (6.9)

Exprimand catul ca o dezvoltare dupa polinoamele din baza ortogonala { }( )nP x

1

11

( ) ( )n

n k kk

q x a P x−

−=

= ∑ (6.10)

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

1

11

12

1

( ) ( ) ( ) ( ) ( ) ( )

( ) ( ) ( ) 0

b b n

n n n k kka a

bn

k n k n nkk a

P x q x w x dx P x a P x w x dx

a P x P x w x dx c δ

−

−=

−

=

⎡ ⎤= =⎢ ⎥

⎣ ⎦

= = =

∑∫ ∫

∑ ∫

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

2 1 1( ) ( ) ( ) ( )b b

n na a

f x w x dx r x w x dx− −=∫ ∫ (6.11)

Sa consideram acum radacinile polinomului ortogonal ( )nP x , adica acele puncte ix pentru care avem relatia: ( ) 0, 1, 2,...,n iP x i n= = (6.12) Scriind relatia 239H(6.7) in aceste puncte rezulta: 2 1 1( ) ( )n i n if x r x− −= (6.13) Restul 1( )nr x− poate fi interpolat Lagrange printr-o relatie de tipul 240H(6.3) si tinand cont de relatia 241H(6.13) se poate scrie:

1 1 2 11 1

( ) ( ) ( ) ( ) ( )n n

n n i i n i ii i

r x r x L x f x L x− − −= =

=∑ ∑ (6.14)

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

2 1 2 11

( ) ( ) ( ) ( ) ( )b bn

n n i iia a

f x w x dx f x L x w x dx− −=

= ∑∫ ∫ (6.15)

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

1

( ) ( ) ( )b n

i iia

f x w x dx f x W=∑∫ (6.16)

unde marimile

( ) ( )b

i ia

W L x w x dx= ∫ , 1, 2,...,i n= (6.17)

se numesc ponderi de integrare Gauss. Ele sunt independente de functia de integrat, dupa cum arata relatia 245H(6.17) si pot fi calculate cu aceasta relatie din functia pondere

( )w x a bazei de polinoame ortogonale considerate si definitia 246H(6.4) a polinomului Lagrange ( )iL x pentru nodurile ix 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 2 1n − , desi functia se calculeaza numai in n puncte. Restul aproximarii va fi analog relatiei 247H(6.6), dar pentru un grad 2 1n − :

( )

11

( )( ) ( )(2 1)!

n n

n ii

fR x x xn

ξ−

=

= −− ∏ (6.18)

find deci de (2 1)!/ ! ( 1)( 2)...(2 1)n n n n n− = + + − ori mai mic. De exemplu, pentru 10n = , restul poate fi de 121.4 10× ori mai mic, ceea ce indica o metoda

hipereficienta. Relatia 248H(6.16) defineste metoda de quadratura Gauss in mod general. Se 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 xf xw x

= (6.20)

iar relatia 249H(6.16) devine

1

( )( )( )

b ni

ii ia

g xg x dx Ww x=

= ∑∫

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

1

( )( )( )

b ni

ii ia

f xf x dx Ww x=

=∑∫ (6.21)

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:

1

( ) ( )b n

i iia

f x dx f x W=

= ∑∫ (6.22)

unde ponderile sunt:

1

( )b b n

ii i

j j ia aj i

x xW L x dx dxx x=

≠

−= =

−∏∫ ∫ (6.23)

nodurile ix fiind radacinile polinomului Legendre de ordin n . Pentru calculul nodurilor ix trebuie rezolvata ecuatia polinomiala 250H(6.12) pentru polinoame Legendre, sau se pot lua din tabele. De asemenea, exista tabele si pentru aceste functii pondere.

Se mai poate folosi de asemenea si formula:

1 1'

1 1 1

( ), ( )( ) ( )

n n ni

n i n i n

P x P x AWP x P x A

− −

− − −

= (6.24)

unde nA si 1nA − sunt coeficientii termenilor de grad maxim din polinoamele ( )nP x , respectiv 1( )nP 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:

2 2

b a b az x− += + (6.25)

astfel ca formula de cuadratura Gauss Legendre devine:

1

( )2 2 2

b n

i iia

b a b a b af x dx f x W=

− − +⎛ ⎞= +⎜ ⎟⎝ ⎠

∑∫ (6.26)

iar ponderile raman neschimbate. 6.3 Cuadratura Clenshaw Curtis Tinand cont de argumentul dat de relatia 251H(6.18), se pare ca metoda de 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 ( )rT x de ordin cel mult n , conform relatiei 252H(5.11) prezentata in capitolul anterior, unde semnul prim al sumei arata ca primul termen se ia injumatatit

0

'( ) ( )n

r rr

f x a T x=

= ∑ (6.27)

Coeficientii ra ai dezvoltarii rezulta din 253H(5.13) inclusiv pentru 0r =

1

2 ( ) ( ), 0,1,...,n

r k r kk

a f x T x r nn =

= =∑ (6.28)

Vom considera in continuare ca esantionarea se face in punctele kx date de radacinile polinomului Chebyshev de ordinul n , deci ( ) cos arccos 0, arccos (2 1) / 2, 1,2,..n k k kT x n x n x k k nπ= = = − =

(2 1) / 2coskkx

nπ−

= . (6.29)

In anumite cazuri se poate considera si esantionarea de tip Lobatto care se face

in extremele polinomului Chebyshev de ordinul n adica coskkxnπ

= . Aceasta

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 254H(6.29) este mai buna. Pe baza dezvoltarii functiei in serie de polinoame Chebyshev 255H(6.27) se poate calcula integrala definita 256H(6.1) pe intervalul [ ]1,1− in forma

1 1 1 1 1

0 00 0 11 1 1 1 1

1' '( ) ( ) ( ) ( ) ( )2

n n n

r r r r r rr r r

f x dx a T x dx a T x dx a T x dx a T x dx= = =− − − − −

= = = +∑ ∑ ∑∫ ∫ ∫ ∫ ∫ (6.30)

deci problema se reduce la calculul coeficientilor dezvoltarii prin relatia 257H(6.28) si al 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

1

( )t

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

' ' 11 1 1 1 1 1

21 1 1 1

( ) ( ) ( ) ( ) ( ) ( ) ( 1)( )2( 1) 2( 1) 2( 1) 2( 1) 2( 1) 2( 1) 1

t tt t rr r r r r r

rT x T x T x T x T t T tT x dx dx

r r r r r r r

++ − + − + −

− − − −

⎡ ⎤ −= − = − = − +⎢ ⎥+ − + − + − −⎣ ⎦

∫ ∫ (6.31)

unde s-a tinut cont ca 1 2 1 11 1( 1) ( 1) ( 1) ( 1) ( 1) ( 1)r r r

r rT T− − +− +− = − = − − = − = − .

Relatia 259H(6.31) este utilizabila pentru 2,3,...,r n= , astfel ca integralele pentru 0,1n = trebuiesc calculate separat:

1

1

2 2

21

0 ( 1) ( ) 1

( 1) 1 11 ( )2 2 4 4

t

t

n dx t T t

tn xdx T t

−

−

= ⇒ = − − = +

−= ⇒ = − = −

∫

∫ (6.32)

unde s-a tinut cont ca primele polinoame Chebyshev sunt 2

0 1 2( ) 1, ( ) , ( ) 2 1T t T t t T t t= = = − . Introducand relatiile 260H(6.31) si 261H(6.32) in ultimul membru al relatiei 262H(6.30) si ordonand dupa ordinul polinoamelor Chebyshev care apar, rezulta

10 0 31 2 1

1 2221

1 12 43

( 1)( ) ( ) ( )2 4 1 2 2 4 4

( ) ... ( )6 6 2 2

t rn

rr

n nn

a a aa a af x dx a T t T tr

a aa aT t T tn n

+

=−

− +

− ⎛ ⎞ ⎛ ⎞= − + + − + −⎜ ⎟ ⎜ ⎟− ⎝ ⎠ ⎝ ⎠

⎛ ⎞⎛ ⎞+ − + + −⎜ ⎟ ⎜ ⎟⎝ ⎠ ⎝ ⎠

∑∫

Daca se ia intervalul de integrare [ ]1, t− , tinand cont ca (1) 1rT = , se obtine formula de cuadratura Clenshaw-Curtis

11

0 1 1 12

2 11

( 1)( ) ( )2 4 1 2

t rn nr r

r rr r

a a a af x dx a T tr r

+−− +

= =−

−−= − + +

−∑ ∑∫ (6.33)

in care coeficientii dezvoltarii se cacluleaza cu relatia 263H(6.28) care poate fi explicitata ca o transformare cosinus discreta

1

2 (2 1) / 2 (2 1) / 2cos cos , 0,1,...,n

rk

k r ka f r nn n n

π π=

− −⎛ ⎞= =⎜ ⎟⎝ ⎠

∑

iar coeficientul 1 0ra + = , deoarece seria 264H(6.27) se considera trunchiata la n . Ca si in cazul integrarii Gauss Legendre, daca integrala este definita pe un

domeniu oarecare, dar finit, [ ],a b , trebuie facuta schimbarea de variabila:

2 2

b a b az x− += + (6.34)

Prin urmare, daca limita superioara de integrare este variabila, pentru un interval [ ],a t formula de cuadratura Clenshaw-Curtis devine:

11

0 1 1 12

2 1

( 1)( ) ( )2 2 4 1 2

t rn nr r

r rr ra

a a a at af x dx a T tr r

+−− +

= =

⎛ ⎞−− −= − + +⎜ ⎟−⎝ ⎠

∑ ∑∫ (6.35)

iar relatia de calcul pentru coeficienti va fi

1

2 (2 1) / 2 (2 1) / 2cos cos , 0,1,...,2 2

n

rk

k t a t a r ka f r nn n n

π π=

− − + −⎛ ⎞= + =⎜ ⎟⎝ ⎠

∑ (6.36)

Daca limita superioara de integrare nu este neaparat variabila este avantajos sa se faca schimbarea de variabila 265H(6.34) catre intervalul de integrare [ ]1,1− . In acest caz, tinand cont ca (1) 1rT = pentru orice r , pornind direct de la relatia 266H(6.31) si facand observatia ca termenii cu r impar devin nuli, se obtine urmatoarea formula, mult mai practica, de cuadratura Clenshaw-Curtis

1

2

0 21

2( )2 1 (2 )

nb

rra

b af x dx a ar

−

=

⎛ ⎞− ⎜ ⎟= +⎜ ⎟−⎜ ⎟

⎝ ⎠∑∫ (6.37)

unde coeficientii ra sunt dati de relatia de transformare cosinus discreta:

1

2 (2 1) / 2 (2 1) / 2cos cos , 0,1,..., 12 2 2

n

rk

k b a b a r k na f rn n n

π π=

− − + −⎛ ⎞= + = −⎜ ⎟⎝ ⎠

∑ (6.38)

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 267H(6.37) si 268H(6.38).

Programul 6.5. Calculul unei integrale definite prin cuadratura Clenshaw- Curtis

f@x_D =1

1+ x2;H∗ Functia de integrat ∗L

a= −3.; b= 3.;H∗ Intervalul ∗LNa= 21;ForAr= 0, r≤ Na, r+= 2, H∗ se calculeaza doar coeficientii pari ∗L

A@rD =2.Na

„j=1

Na

fACosAH2 j− 1L π

2

NaE

b− a2

+b+ a

2E CosA

r H2 j− 1L π

2

NaEE;

Icc=b− a

2

i

k

jjjjjjjj‚r=1

Na2 −1

A@2 rD 2

1− H2 rL2+ A@0D

y

{

zzzzzzzz;

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 2n la 2logn n daca n este o putere a lui 2. Astfel, de exemplu pentru 32n = , volumul de calcule scade de peste 6 ori, iar pentru 64n = 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 0x , in care valoarea lui ( )y x este cunoscuta ca o conditie initiala 0y , se poate dezvolta ( )y x in serie Taylor la distanta xδ de acest punct. Limitand-ne la primii doi termeni, avem ' 2

0 0 0( ) ( ) ( ) ( )y x x y x x y x xδ δ δ+ = + + (7.2) Notand xδ cu h , (pasul de esantionare) si pe '

0( )y x cu 1k , neglijand termenii de ordin superior putem afla valoarea lui ( )y x in punctul 0x h+ 0 0 1( ) ( )y x h y x h k+ = + (7.3) in care, conform relatiei 269H(7.1) 1 ( , )n nk f x y= (7.4)

Figura 7.1 Metoda Euler de ordinul I pe un subinterval

Putem repeta procedeul plecand acum din punctul 0x h+ si determinam valoarea lui ( )y x intr-un punct 0 2x h+ , si asa mai departe.

Notand in general 0 nx nh x+ = si 0( ) ny x nh y+ = rezulta ca , pornind de la conditia initiala cunoscuta 0 0( )y x y= , putem calcula succesiv valorile lui ( )y x intr-un numar oarecare de puncte echidistante dintr-un interval [ , ]a b prin formule similare cu 270H(7.3) si 271H(7.4): 1 1n ny y h k+ = + (7.5) 1 0 0( , )k f x y= (7.6)

Pasul este dat de

b ahN−

= (7.7)

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 272H(7.2), eroarea la

fiecare interval este 2

22

( )( ) b ahN

⎛ ⎞−= ⎜ ⎟

⎝ ⎠ si deoarece erorile se cumuleaza la fiecare

interval si avem N intervale, eroarea totala este

2

2 ( ) 1( ) ( )tb aN h h

N Nε

⎛ ⎞−= = ⎜ ⎟

⎝ ⎠∼ ∼ (7.8)

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 nx hα+ , cu (0,1)α ∈ , si daca s-ar folosi in relatia 273H(7.5) derivata 2k in acest punct, s-ar putea ca pentru un α corespunzator sa se obtina o eroare de semn invers fata de cea data de derivata in punctul nx , ca in figura 7.3. O medie ponderata a acestor doua derivate ar putea da o panta 1 1 2 2k w k w k= + (7.9) (cu ponderile 1w si 2w cuprinse intre 0 si 1) care sa conduca chiar in punctul corect,

1( )ny x + ,anuland eroarea. Din pacate, nici panta 2k in punctul nx hα+ , data conform ecuatiei explicite 274H(7.1) de [ ], ( )n nf x h y x hα α+ + , nu se cunoaste deoarece nu se cunoaste valoarea functiei ( )ny x hα+ in acest punct. Am putea insa aproxima pe

( )ny x hα+ dupa formula de ordinul I 1( )n ny x h y hkα α+ ≈ + Totusi, nici aceasta nu da o valoare exacta pentru ( )ny x hα+ ; cea corecta s-ar obtine daca am lua o alta abscisa, nx hβ+ , in care apare o noua necunoscuta, β , dar asigura egalitatea 1( )n ny x h y hkα β+ = + . Rezumand, vom folosi urmatoarul sistem: 1 1 1 2 2( )n ny y h w k w k+ = + + (7.10) [ ]1 ,n nk f x y= (7.11)

[ ]2 1,n nk f x h y hkα β= + + (7.12)

Figura7.3 Pantele necesare pentru metoda Euler de ordinul II

Exista patru necunoscute, 1w , 2w ,α si β , pentru aflarea carora sa punem urmatoarea conditie: dezvoltarea in serie Taylor a lui ( )ny x h+ in punctul nx trebuie sa coincida cu relatia 275H(7.10), in care pe 2k il dezvoltam de asemenea in serie in punctul nx . Sa dezvoltam mai intai pe y :

2

' '' 31 ( ) ( ) ( )

2!n n n nhy y hy x y x h+ = + + + (7.13)

Conform relatiilor 276H(7.1) si 277H(7.4), avem '1( ) ( , )n n ny x f x y k= = , iar derivata a

doua se poate scrie pornind de la aceleasi relatii:

'' ' '( , ) ( , )( ) ( , ) ( )n n n nn n

f x y f x yy x f x y y xx x y yx y

∂ ∂= = + ⋅

= =∂ ∂

si notand derivatele partiale in raport cu x si y cu xnf respectiv ynf obtinem

''1( )n xn yny x f f k= + (7.14)

Rezulta ca dezvoltarea in serie Taylor are forma

2 2

1 1 12 2n n xn ynh hy y hk f f k+ = + + + (7.15)

Sa consideram acum relatia 278H(7.10) in care de asemenea dezvoltam in serie Taylor pe [ ]2 1,n nk f x h y hkα β= + + , conform relatiei 279H(7.12), ca functie de doua

variabile, in jurul punctului ( ),n nx y

[ ]2 1( , ) ( , ), ...n n

n n

f x y f x yk f x y h hkx x y yx y

α β∂ ∂= + + +

= =∂ ∂ (7.16)

Limitandu-ne la acesti primi termeni, si tinand seama de definitiile precedenta obtinem 2 1 1xn ynk k hf hk fα β= + + care introdusa in relatia 280H(7.10) conduce la ( ) 2 2

1 1 1 2 1 2 2 1n n xn yny y h w k w k h w f h w k fα β+ = + + + + (7.17) Identificand asa cum am propus relatiile 281H(7.13) si 282H(7.17) dupa puterile lui h obtinem urmatoarele conditii:

1 2

2

2

11212

w w

w

w

α

β

⎧⎪ + =⎪⎪ =⎨⎪⎪ =⎪⎩

(7.18)

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.

Alegand ponderi egale, rezulta 1 212

w w= = , de unde 1α β= = , deci

algoritmul devine:

1 21 2n n

k ky y h+

+= + (7.19)

[ ]1 ,n nk f x y= (7.20)

[ ]2 1,n nk f x h y hk= + + (7.21) definind o metoda Euler de ordinul II pe care o putem numi metoda pantei de mijoc.

Alegand pe 12

α = , rezulta 2 1w = , 1 0w = si 12

β = , deci algoritmul devine:

1 2n ny y hk+ = + (7.22) [ ]1 ,n nk f x y= (7.23)

2 1,2 2n nh hk f x y k⎡ ⎤= + +⎢ ⎥⎣ ⎦

(7.24)

definind o metoda Euler de ordinul II pe care o putem numi metoda punctului de mijoc. Conform relatiei 283H(7.13) termenii neglijati sunt de ordinul 3h , acestia dand

eroarea pe fiecare interval, astfel ca eroarea totala pe b aNh−

= intervale va fi:

22

1t h

Nε ∼ ∼ (7.25)

deci mult mai mica (de N ori) decat la metoda de ordinul I.

Clear@"̀ ∗"D;H∗ Problema Cauchy: y'@xD=f@x,yD ∗Lf@x_, y_D = −yHx+ 1L − Cos@xD;yi= 1;H∗ Conditia initiala ∗Lxi= 0.2;xf = 4.; H∗ Intervalul ∗LNp= 100;H∗ Numar de puncte de evaluare ∗L

h=xf− xi

Np;H∗ Pasul∗L

y@0D = yi;Y= 88xi, yi<<; H∗ Lista de interpolare a solutiei ∗LH∗ membrul drept al ecuatiei ∗LForAm = 0, m ≤ Np−1, m++,

xm = xi+ m h; H∗ Punctul curent pe abscisa ∗LK1= f@xm, y@mDD;K2= f@xm+ h, y@mD + h K1D;

y@m+ 1D = y@mD +h2

HK1+ K2L; H∗ Formula Euler cu panta medie ∗L

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

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

H∗ Verificarea, prin comparatie cu solutia analitica, daca aceasta exista ∗Lsol@x_D = DSolve@8y'@xD −y@xD Hx+ 1L − Cos@xD, y@xiD yi<, y@xD, xD@@1, 1, 2DDH∗ Solutia analitica ∗Lg2= 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− 1LhDD<D;D;Print@Max@Abs@erabsDDD;

H∗ Problema: y'@xD=−y@xD Hx+1L−Cos@xD ∗Lf@x_, y_D = −yHx+ 1L − Cos@xD;yi= 1;H∗ Conditia initiala ∗Lxi= 0.2;xf = 4.; H∗ Intervalul ∗LNp= 100; H∗ Numar de puncte de evaluare ∗L

h=xf− xi

Np; H∗ Pasul∗L

y@0D = yi;Y= 88xi, yi<<; H∗ Lista de interpolare a solutiei ∗LForAm = 0, m ≤ Np−1, m++,

xm = xi+ m h; H∗ Punctul curent pe abscisa ∗LK1= f@xm, y@mDD;

K2= fAxm+h2

, y@mD+h2

K1E;

y@m+ 1D = y@mD + h H K2L; H∗ Formula Euler cu punct mediu ∗LAppendTo@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− 1LhDD<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 284H(7.9)

1 1 2 2k w k w k= + . Daca s-ar folosi mai multe pante, s-ar obtine metode de ordin superior, m ,

printr-o generalizare de tipul:

1

m

i ii

k w k=

=∑ (7.26)

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 285H(7.4)

1 ( , )n nk f x y= , iar a doua panta se determina pe baza acesteia, prin relatia 286H(7.12) [ ]2 1,n nk f x h y hkα β= + + .

Generalizand, o panta de ordinul m se poate obtine din panta de ordin 1m − printr-o relatie de tipul:

1

11

,i

i n i n ij ij

k f x h y hkα β−

−=

⎡ ⎤= + +⎢ ⎥

⎣ ⎦∑ (7.27)

unde 1 101, , 1, 1, 0i m j i α β= = − = = .

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 287H(7.13) si 288H(7.17) folosite in paragraful precedent pentru gasirea coeficientilor , , wα β si k trebuie sa cuprinda termeni pana la ordinul m , astfel ca expresiile finale de tipul 289H(7.17) devin din ce in ce mai complexe. De exemplu, 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 1/ 2iα = . Rezulta printr-un procedeu analog cu cel prezentat in paragraful precedent, urmatoarea schema larg raspandita de tip RK4, bazata pe astfel de alegeri:

( )

( )

1 1 2 3 4

1

2 1

3 2

4 3

( 2 2 )6

,

,2 2

,2 2,

n n

n n

n n

n n

n n

hy y k k k k

k f x y

h hk f x y k

h hk f x y k

k f x h y hk

+ = + + + +

=

⎛ ⎞= + +⎜ ⎟⎝ ⎠⎛ ⎞= + +⎜ ⎟⎝ ⎠

= + +

(7.28)

Clear@"̀ ∗"D;H∗ y'@xD=−y@xD Hx+1L−Cos@xD; ODE sub forma explicita ∗Lf= −yHx+ 1L −Cos@xD; H∗ Forma explicita fara argumente ∗Ly@0D = 1; H∗ Conditia initiala ∗LxM= 4.;x0 = 0; H∗ Intervalul ∗L

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

h=xM− x0

Np; H∗ Pasul ∗L

Y= 880, y@0D<<; H∗ Punctul initial al graficului solutiei ∗L

ForAm = 0, m ≤ Np−1, m++, H∗ algoritmul RK4 in fiecare punct ∗Lxm = m h;K1= fê. 8x→ xm, y→ y@mD<;

K2= fê. 9x→ xm+h2

, y→ y@mD+h2

K1=;

K3= fê. 9x→ xm+h2

, y→ y@mD+h2

K2=;

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

y@m+ 1D = y@mD +h6

HK1+ 2 K2+ 2 K3+ K4L;

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

H∗ Y contine lista cu punctele solutiei ∗L

g1= ListPlot@Y, PlotJoined→ FalseD; H∗ reprezentarea grafica a solutiei ∗LH∗ Se interpoleaza lista construind polinomul de interpolare Lagrange ∗L

p@x_D:= „i=1

Np+1

Y@@i, 2DD i

kjjjj‰

j=1

i−1 x− Y@@j, 1DDY@@i, 1DD − Y@@j, 1DD

y

{zzzz i

kjjjj ‰

j=i+1

Np+1 x− Y@@j, 1DDY@@i, 1DD − Y@@j, 1DD

y

{zzzz;

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 ∗LClear@x, yD;H∗ Aceasta ecuatie are solutie analitica exacta ∗Lsol= 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

h1=xM− x0

nl; H∗ Pas mai fin de esantionare ∗L

ert= 8<; H∗ Lista erorilor ∗Lz@x_D = sol;For@n= 1, n≤ nl, n++, AppendTo@ert, 8p@x0+ h1nD − z@x0+ h1nD<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 1nx + 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

valoarea functiei in punctul urmator 1 1( )n ny x h+ + folosind mai intai un pas 1p

b ahN−

=

(unde ,a b si pN sunt respectiv capetele intervalului si numarul de puncte din interval), si apoi repetarea calculului folosind valori kh din ce in ce mai mici pentru pas. Rezultatele 1

khny + obtinute intr-un anumit punct 1nx + pot forma perechi cu valorile

respective ale lui kh sub forma: 1 21 1 1 2 1{ , },{ , },...,{ , },mhh h

n n n my h y h y h+ + + iar aceste puncte pot fi folosite pentru o interpolare, de exemplu polinomiala, care permite calculul valorii functiei pentru 0h = , (care ar corespunde valorii corecte). Valorile pasilor trebuie alese in asa fel incat de fiecare data, pornind din nx sa se ajunga la acelasi 1nx + la care se ajunsese cu pasul initial. Propunerea Burlich-Stoer pentru marirea numarului de

subintervale este: 1 1 1 1 1 1, , , , , ,...2 4 6 8 12 16kh h h h h hh = deci nu o simpla injumatatire care ar

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

Metode munerice

Documents

Transcript of Metode munerice