Metode Iterative de Rezolvare a Sistemelor Algebrice Liniare

Capitolul I. Introducere

Dezvoltarea metodelor de calcul numeric a fost influențată de progresele obținute în domeniul calculatoarelor electronice. Îmbunătățirile aduse în domeniul hardware, gestiunii resurselor și tehnicilor moderne de prelucrare a informației au condus la apariția unor metode de calcul, algoritmi, criterii de alegere între diverse metode ce se pot aplica aceleiași probleme date.

Dezvoltarea în domeniul hardware a permis construirea unor algoritmi de calcul numeric ce necesită un grad mare de paralelism. În acest sens se poate exemplifica rezolvarea sistemelor algebrice liniare prin metode iterative :

Metoda Gauss – Seidel care implică modificarea unei necunoscute la fiecare moment de timp este preferabilă metodei Jacobi care implică modificarea tuturor necunoscutelor în același timp, asta c=nd se dispune de un calculator serie obișnuit.

Metoda Jacobi este preferabilă metodei Gauss – Seidel c=nd se dispune de un calculator paralel, iar timpul necesar pentru o iterație va fi mult mai mic.

Trebuie menționat faptul că la rezolvarea numerică a unei probleme cu ajutorul calculatorului terbuie avută în vedere aritmetica calculatorului utilizat.

Este necesar a se menționa două aspecte din acest domeniu :

Analiza erorilor pentru operațiile executate în virgulă mobilă. În cazul unui calculator cu registre de lungime dublă, rezultatul operației de adunare sau scădere a două numere (a ^ b) cu a și b reprezentate în virgulă mobilă va fi de forma (a ^ b)(1 ^ ), unde 2-n iar n este numărul de cifre binare din mantisa numerelor reprezentate în virgulă mobilă.

Un alt aspect este intervalul aritmetic, orice număr x se poate reprezenta cu ajutorul unui interval x . C=nd se execută operații aritmetice are loc o recalculare a intervalelor corespunzătoare și în fiecare moment se dispune de o limită inferioară și superioară pentru toate numerele, totuși în foarte multe cazuri care implică șiruri mari de calcule, dimensiunea intervalelor implicate devine destul de mare.

Fiecare problemă se rezolvă separat iar combinarea rezultatelor dă soluția problemei originale.

În cadrul metodelor de calcul numeric mai trebuie considerată și problema complexității calculelor. Se pune problema determinării numărului de operații necesare pentru a rezolva o problemă dată.

În metoda Gauss de eliminare sunt necesare n3/3 operații de înmulțire și n3/6 operații de adunare, iar dacă facem o organizare adecvată a algoritmului aceeași metodă a lui Gauss de eliminare necesită numai n3/6 operații de înmulțire, dar 2n3/3 adunări. În metodele directe de rezolvare a sistemelor de ecuații liniare, în general, se folosesc n3/3 operații elementare, iar în metodele iterative se folosesc n2 operații aritmetice.

Din aceste exemple prezentate se vede că înainte de prelucrarea numerică cu ajutorul calculatorului a unei probleme date se poate alege metoda care necesită cel mai mic număr de operații, obțin=ndu-se astfel reducerea timpului de execuție pe calculator.

Din cele prezentate se poate spune că metodele de calcul numeric vor cunoaște o dezvoltare serioasă, dezvoltare ce este în legătură cu facilitățile oferite de calculatoarele eletronice foarte rapide precum și cu progresele privind tehnicile de prelucrare electronică a informației, conduc=nd astfel la noi probleme de cercetare în domeniu și în afara domeniului.

Algoritmi de calcul

Prin algoritm se înțelege un sistem de reguli, care, aplicat la o anumită clasă de probleme de același tip, conduce de la condițiile inițiale ale clasei de probleme la soluție, cu ajutorul unor operații succesive unic determinate.

Dezvoltarea noțiunii de algoritm a fost impusă de necesități practice și a luat amploare datorită apariției mașinilor de calcul. Scopul teoriei algoritmilor este de a găsi o metodă generală de pregătire a problemelor din punct de vedere matematic și logic și astfel apare necesitatea unei formulări corespunzătoare a problemelor. Această formalizare a problemelor se poate numi algoritmizarea aparatului matematic.

Un algoritm constituie un sistem de reguli care transformă o informație inițială într-o informație finală cu ajutorul unui șir de informații intermediare.

Caracteristicile unui algoritm :

a) Generalitatea.

Algoritmul descrie numai cum trebuie efectuate calculele : indiferent de cazul particular care ne interesează, calculele se vor desfășura întotdeauna după aceleași reguli, aplicate însă datelor caracteristice cazului particular considerat.

b) Finitudinea.

Prin aceasta se înțelege că numărul de transformări intermediare aplicate informației inițiale pentru a obține informația finală este finit.

c) Unicitatea.

Toate transformările intermediare aplicate informației inițiale pentru a obține informația finală sunt unic determinate de regulile algoritmului (după fiecare etapă a transformărilor, regulile algoritmului precizează care anume transformare va urma și de asemenea precizează c=nd se obține informația finală după care activitatea algoritmului se întrerupe).

În rezolvarea unei probleme cu ajutorul calculatorului pot fi puse în evidență următoarele etape :

Enunțarea problemei și formularea matematică. În această etapă se precizează (enunță) problema practică. Pentru rezolvarea acesteia apare necesitatea formulării matematice (modelării). Se exprimă matematic relații și restricții între parametrii problemei, se pun în evidență condițiile inițiale ale problemei, restricțiile referitoare la soluție.

Alegerea metodei numerice. Rezolvarea unei probleme presupune existența unui algoritm de calcul.. Deoarece problema se rezolvă cu ajutorul calculatorului, pentru alegerea metodei pot fi luate în considerare următoarele criterii :

precizie;

viteza de calcul;

cantitatea de date pentru a stabili memoria necesară;

simplitatea formulelor de calcul (necesară pentru întocmirea cu ușurință a programului de calcul).

Întocmirea programului de calcul. După ce algoritmul metodei numerice a fost pus în evidență, pentru a putea utiliza calculatorul, acesta trebuie transcris în limbaj de programare (se întocmește un program de calcul cu ajutorul instrucțiunilor).

Testarea rezultatelor. Sunt necesare diferite procedee de control care să permită detectarea eventualelor erori apărute în calcule, în funcție de acestea oprirea sau reluarea calculelor.

Interpretarea rezultatelor din punctul de vedere al problemei practice propuse.

Capitolul al IIlea. Noțiuni pregătitoare

1). Vectori și matrici. Clase speciale de matrici.

Fie X o mulțime nevidă.

Definiție : Aplicația d : XxX R^ se numește distanță (metrică) dacă îndeplinește următoarele condiții :

M1) d(x,y) 0, x,yX

d(x,y) # 0 x # y

M2) d(x,y) # d(y,x), x,yX

M3) d(x,z) d(x,y) ^ d(y,z), x,y,zX

Definiție : Un spațiu liniar real este o mulțime X înzestrată cu o operație internă de adunare și o operație externă de înmulțire cu numere reale (scalari), ce satisface axiomele :

x ^ y # y ^ x, x,yX

x ^ (y ^ z) # (x ^ y) ^ z, x,y,zX

0 X, x ^ 0 # x, xX

x X, (-x) X așa înc=t x ^ (-x) # 0

(x ^ y) # x ^ y, x,yX, R

( ^ ) x # x ^ x, xX, ,X

(x) # ()x, xX, ,X

1x # x, xX

Definiție : Fie M X. Mulțimea M se numește subspațiu liniar al spațiului liniar X dacă îndeplinește condițiile :

1). x ^ y M, x,yM

2). x M, xM, R

Fie Rn un spațiu liniar.

Elementele acestui spațiu le vom nota : x # (1, 2,…,n), cu i R, i # 1,n

În acest spațiu operațiile de adunare internă și înmulțire cu scalari reali se defnesc în mod natural prin operații pe componente.

Definiție : Elementele spațiului liniar Rn se numesc vectori (x # (1, 2,…,n)T)

Vom nota cu ei # (0,0,…,0,1,0,…,0)T, în care 1 este pe poziția i.

Acești vectori vor forma baza naturală a spațiului Rn, B # { ei }i = 1,n

În aceste condiții, vectorii x se vor scrie :

Fie vectorii : x # (1, 2,…,n)T și y # (1,2,…,n)

Definiție : Definim produsul scalar al vectorilor x și y ca fiind :

Structura de spațiu liniar euclidian a spațiului Rn este dată de produsul scalar.

Proprietățile produsului scalar :

P1). (x,y) # (y,x), x,yRn

P2). (x^y, z) # (x,z) ^ (y,z), x,y,zRn, , R

(x,y^z) # (x,y) ^ (x,z), x,y,zRn, , R

P3).

Vom nota cu ca fiind modulul vectorului x Rn (acesta este un număr strict pozitiv, pentru x 0, după P3).

Definiție : Dacă (x,y) # 0 atunci x și y se numesc vectori ortogonali.

Observație : Vectorul nul 0 # (0,0,…,0)T este ortogonal pe orice alt vector din spațiu.

Definiție : 1). Dacă vectorii unei baze a spațiului euclidian sunt ortogonali doi c=te doi atunci baza se numește bază ortogonală.

B = {yi}iN ; (yi , yj ) = 0, cu B se numește ortogonală.

2). Dacă o bază ortogonală are vectorii de modul egal cu 1 atunci ea se numește bază ortonormală.

B = {yi}iN ; yi = 0 , cu B se numește ortogonală.

Observație : Baza naturală {e1, e2, …, en} este bază ortonormală.

Matrici. Clase speciale de matrici.

Fie X, Y – două spații liniare.

Definiție : 1). O aplicație A : X Y se numește surjectivă dacă : yY, x X a.î.

y # A(x).

2). O aplicație A : X Y se numește injectivă dacă :

x,y X, A(x), A(y) Y a.î. A(x) # A(y) x # y.

3). O aplicație A : X Y se numește bijectivă dacă este injectivă și surjectivă.

4). O aplicație A : X Y cu proprietatea : x # y A(x) # A(y) se numește funcție sau operator.

5). Mulțimea D(A) X se numește domeniul de definiție a funcției A.

6). Mulțimea R(A) Y se numește codomeniul funcției A.

7). Aplicația A-1 : R(A) X, definită prin

y R(A), x # A-1(y) y # A(x) , se numește funcția inversă a funcției A.

8). O funcție A : X Y se numește operator liniar dacă :

A(x ^ y) # A(x) ^ A(y), x,y X, , R.

În particular : A(0) # 0.

Mulțimea operatorilor liniari este un spațiu liniar, motat cu L(X,Y), în care operațiile se definesc în mod natural cu ajutorul operațiilor lui Y.

Fie X # Rn și Y # Rm.

Fie B # {e1,e2,…,en} și B* # {f1,f2,…,fm} – bazele naturale ale celor două spații.

Un operator liniar A : X Y se va defini în modul următor :

Fie vectorul , atunci

Așadar, fiecărui operator liniar A : Rn Rm îi putem atașa matricea (aij) # A. Coloanele matricei A sunt vectorii .

Reciproc, o matrice A de dimensiune mxn definește în mod unic un operator liniar care duce orice vector x # (1, 2,…,n) Rn în vectorul

.

Adunarea, înmulțirea cu numere reale și compunerea a doi operatori liniari conduc la adunarea, înmulțirea cu numere reale și multiplicarea matricelor.

În acest mod putem să identificăm operatorii liniari în spații finit dimensionale cu matricele.

Cu aceste notații vectorii pot fi considerați și ei ca matrici cu o singură coloană, iar produsul scalar ca un produs matricial :

Definiție : 1). Mm,n – spațiul matricelor cu m – linii și n – coloane .

2). Pentru o matrice A #(aij), i #1,m și j # 1,n vom nota cu AT # (aji) Mn,m – matricea transpusă.

– matricea adjunctă (transpusă și conjugată)

În cazul real, A Mm,n, avem AH # AT.

Proprietăți :

P1). (AH)H # A, A Mm,n

P2). (AB)H # BH AH, A Mm,n și B Mn,p

P3). (AH)-1 # (A-1)H, A Mm,n – matricea A este pătratică și nesingulară (det A 0)

P4). (Ax,y) # (x, AHy), A Mm,n, x Cn, y Cm

Definiții : 1). Dacă A este pătratică și A# AH atunci A se numește matrice hermitiană (autoadjunctă).

2). Dacă A este pătratică și A# AT atunci A se numește matrice simetrică (elemenele simetrice față de diagonala principală a matricei sunt egale).

În particular (cazul R) : cele două noțiuni se confund.

Observație : Fie A – matrice complexă hermitiană. Atunci produratori liniari conduc la adunarea, înmulțirea cu numere reale și multiplicarea matricelor.

În acest mod putem să identificăm operatorii liniari în spații finit dimensionale cu matricele.

Cu aceste notații vectorii pot fi considerați și ei ca matrici cu o singură coloană, iar produsul scalar ca un produs matricial :

Definiție : 1). Mm,n – spațiul matricelor cu m – linii și n – coloane .

2). Pentru o matrice A #(aij), i #1,m și j # 1,n vom nota cu AT # (aji) Mn,m – matricea transpusă.

– matricea adjunctă (transpusă și conjugată)

În cazul real, A Mm,n, avem AH # AT.

Proprietăți :

P1). (AH)H # A, A Mm,n

P2). (AB)H # BH AH, A Mm,n și B Mn,p

P3). (AH)-1 # (A-1)H, A Mm,n – matricea A este pătratică și nesingulară (det A 0)

P4). (Ax,y) # (x, AHy), A Mm,n, x Cn, y Cm

Definiții : 1). Dacă A este pătratică și A# AH atunci A se numește matrice hermitiană (autoadjunctă).

2). Dacă A este pătratică și A# AT atunci A se numește matrice simetrică (elemenele simetrice față de diagonala principală a matricei sunt egale).

În particular (cazul R) : cele două noțiuni se confund.

Observație : Fie A – matrice complexă hermitiană. Atunci produsul scalar (Ax,x) este real.

3). Dacă, în plus, produsul scalar :

– (Ax,x) 0, atunci A se numește matrice pozitiv semidefinită.

– (Ax,x) > 0, x0, atunci A se numește matrice pozitiv definită.

4).O matrice care are coloanele e1,e2,…,en (baza naturală a spațiului euclidian) se numește matricea unitate (notată cu I ).

5). O matrice A se numește matrice unitară dacă AHA# I (pentru cazul real, dacă ATA#I atunci A se numește matrice ortogonală).

6). Dacă AHA#AAH atunci A se numește matrice normală.

7). Fie A # (aij) Mn,n(R).

Dacă aij#0, n i > j 1, atunci A se numește matrice superior triunghiulară.

În acest caz A se va scrie sub forma :

Dacă aij # 0, 1 i < j n, atunci A se numește matrice inferior triunghiulară.

În acest caz A se va scrie sub forma :

8). Dacă A este în același timp și inferior triunghiulară și superior triunghiulară atunci A se numește matrice diagonală.

În acest caz, A se va scrie sub forma : A # diag (a11,a22,….,ann).

Observație : O matrice normală și triunghiulară este diagonală.

9). Dacă matricea A are elementele diferite de 0 pe diagonala principală și pe cele două diagonale alăturate atunci ea se numește matrice triagonală.

În acest caz, A se va scrie sub forma :

10). Matricea A se numește matrice nesingulară dacă det A 0.

11). Matricea A se numește matrice singulară dacă det A # 0.

12). Pentru AMn(R) și H – nesingulară atunci matricele A și HAH-1 se numesc matrici asemenea.

13). Matrice A se numește matrice regulată dacă ea este asemenea cu o matrice B diagonală.

14). Matricea pentru care dimensiunea spațiului generat de vectorii ei proprii este mai mică dec=t ordinul matricii se numește matrice defectivă (de exemplu : atunci c=nd o matrice are valorile proprii multiple).

15). O matrice cu multe elemente egale cu zero se numește matrice rară.

2). Norme vectoriale. Norme matriceale.

Definiție 1 : Fie X un spațiu liniar. Aplicația : X R se numește normă dacă îndepli-nește condițiile :

N1) și (vectorul nul)

N2)

N3)

Spațiul liniar X, înzestrat cu norma , se numește spațiu liniar normat.

Consecință : Are loc inegalitatea :

x – y # x – y , x,yX

Definiție 2 : Spunem că șirul {xn} X este șir convergent la xX, dacă

Oservație : Din consecința de mai sus ne rezultă ușor continuitatea normei ca funcție de elementele spațiului X.

Definiție 3: Spunem că , 1 sunt norme echivalente dacă există constantele m,M >0 așa înc=t să avem :

mx 1 x M x 1 , xX

Observație : În acest caz, orice șir convergent în sensul uneia dintre norme este convergent și în sensul celeilalte norme.

Teorema 1 : Pe Rn toate normele vectoriale sunt echivalente și convergența, în sensul oricăror norme, este convergența pe componente.

Spațiul Rn este un spațiu liniar normat.

Pentru a dovedi acest lucru, vom demonstra că modulul vectorial : Rn R este o normă pe Rn.

: Rn R

N1) x 0, xRn

Acest lucru este adevărat deoarece (din proprietatea P3 a produsului scalar)

x # 0 (x,x) i # 0, x # 0.

x # (1, 2,…, n)T

N2) , xRn, R

N3) Fie x # (1, 2,…, n)T, y # (1, 2,…, n)T, x,yRn

Fie R

Conform proprietăților produsului scalar, vom avea :

Cum

Deoarece (care este inegalitatea lui Schwarz)

Am arătat că este o normă pe Rn, care va fi numită norma euclidiană a spațiului liniar normat Rn.

Această normă o vom mai nota cu E . Deci pentru xRn, vom avea :

x E #

Generalizat, vom avea : x p # , 1 p < .

sau x # (această normă vectorială se obține pentru p )

Lema 2 : Pentru orice normă vectorială și P o matrice pătratică nesingulară de ordin n, putem defini o normă vectorială prin :

x P # Px , xRn

Demonstrație :

Pentru a demonstra că x P # Px , xRn este o normă vectorială trebuie să arătăm că au loc axiomele din definiția normei vectoriale.

(N1) x P 0, xRn și x P # 0 x # 0.

x P # Px 0 (deoarece Px este normă). Deci x P 0, xRn.

Din x P # 0, xRn, avem că Px # 0 și, deoarece P este o matrice nesingulară (det P 0), avem că : x # 0.

(N2) x P # x P , R, x Rn.

Avem : x P # Px # # Px # # x P.

(N3) x^ y P x P^ y P , x,y Rn

Avem : x^ y P # P(x^ y) # Px^ Py Px ^ Py

Deci : x^ y P x P^ y P , x,y Rn

În particular : pentru P # diag (p1, p2, …, pn), pi > 0 putem obține normele ponderate :

x 1,P #

x 2,P #

x , P #

Fie A # (aij)nxm – o matrice pătratică

Definiție : Aplicația : Mn,m R se numește normă matriceală dacă satisface condițiile :

NM1) și (matricea nulă)

NM2)

NM3)

NM4)

Observație : Toate normele matriceale sunt echivalente și determină aceeași convergență a șirurilor de matrici (este convergentă după elementele matricei).

Exemplu de normă matriceală :

1). Norma euclidiană (a lui Fröbenius) :

{tim că matricile sunt operatori liniari.

Să considerăm norme matriceale care să satisfacă o proprietate de compatibilitate cu norme vectoriale.

Fie două norme : vectorială și una matriceală, notate cu .

A Mn,n, x Rn, vom avea : (cele două norme sunt compatibile)

În particular : E este compatibilă cu E : (aceasta se poate verifica cu inegalitatea lui Schwartz)

Observație : Putem defini o normă matriceală cu ajutorul unei norme vectoriale, dacă ele sunt compatibile.

Fie – o normă vectorială.

Ea este o funcție continuă pe mulțimea S # {y Cn; y = 1}

Fie A # – aceasta se numește norma matriceală indusă de norma vectorială sau normă naturală.

Teorema 3 : Norma matriceală indusă de norma vectorială :1 este definită pentru orice matrice A # (aij)nxn prin :

Demonstrație : fie x # (1, 2, …, n)T Rn, cu

Atunci vom avea :

()

Pentru a termina demonstrația este suficient să găsim un vetor y cu y# 1 așa înc=t să avem Ay1 # A1 .

Dacă , putem alege vectorul y # ek, pentru care inegalitățile de mai sus () devin egalități.

Deci : Aek 1 # A1 .

Am arătat că A1 este definită prin : și o vom nota cu AC .

Teorema 4 : Norma matriceală indusă de norma vectorială este definită pentru orice matrice A # (aij)nxn prin :

Demonstrație : fie x # (1, 2, …, n)T Rn, cu .

Atunci vom avea :

.

Fie și fie vectorul y # (1 ,2 , …, n)T , vector ale cărui componentele sunt date de :

Deoarece avem : urmează că avem : Ay # A, ceea ce demonstrează teorema.

Teorema 5 : Fie norma vectorială și matriceală notată prin , P – o matrice nesingulară. Atunci norma matriceală indusă de norma vectorială P este definită pentru orice matrice A Mn,n prin :

În particular : P # diag (p1, p2, …, pn), pi > 0.

Atunci normele matriceale ponderate, induse de normele vectoriale ponderate, sunt :

a). 1,P :

b). ,P :

c). norma euclidiană ponderată :

Teorema 6 : Pentru orice matrice pătratică complexă A, există o matrice unitară așa înc=t B # PHAP să fie superior triunghiulară.

Norma spectrală :

Fie – valoare proprie a matricei A Mn,n (Ax # x)

Definiții : 1). Atunci se numește raza spectrală a matricei A.

2). Norma indusă de modulul vectorial prin se numește norma spectrală.

Teorema razei spectrale : Raza spectrală nu depășește nici o normă a matricii, care este compatibilă cu o normă vectorială. Mai mult, pentru fiecare matrice A și pentru un > 0 dat, există o normă naturală , așa înc=t raza spectrală (A) satisface relația :

(A) A (A) ^

3). {iruri și serii de puteri matriceale

Propoziția 1 : Raza spectrală nu depășește nici o normă a matricii, care este compatibilă cu o normă vectorială. Mai mult, pentru fiecare matrice A și pentru un > 0 dat, există o normă naturală , așa înc=t raza spectrală (A) satisface relația :

(1) (A) A (A) ^

Deomonstrație :

Fie dată o normă vectorială, notată cu , și o normă matriceală compatibilă cu ea (notată cu ).

Fie o matrice A Mm,n (R).

Atunci, pentru orice valoare proprie și pentru orice x vector propriu a matricii A, avem :

x# x# Ax A x , x – vector propriu al lui A.

Deci, avem : A , adică (A) A , ceea ce demonstrează prima parte a teoremei.

Pentru a demonstra a doua parte a teoremei, să observăm că prima inegalitate din relația (1) este satisfăcută pentru orice normă naturală (acestea sunt compatibile cu norme vetoriale).

Aplic=nd teorema 6 (din paragraful precedent), teoremă ce asigură existența unei matrici unitare P astfel ca B # PHAP să fie superior triunghiulară, vom putea scrie :

B # ^ U, unde # diag (1 , 2 , …, n) ; i – valori proprii.

U # (uij ), cu uij # 0 pentru i j.

Fie, > 0 oarecare și fie D # diag (1, , 2, …, n-1).

Să notăm cu : C # D-1BD # D-1D ^ D-1UD # ^ E, unde E # (eij) cu elementele

.

Observăm că, dacă 0, atunci elementele matricii E tind și ele la zero ( eij 0).

Să notăm cu () # E e , atunci vom avea : .

Deoarece C # D-1BD # D-1 (PHAP)D, dacă notăm cu R # D-1PH și folosim lema 2 din paragraful precedent (pentru orice normă definită peste R și pentru P o matrice pătratică, de ordinul n, nesingulară, ave definită o normă vectorială prin xP # Px, xRn) vom avea în cazul nostru :

xR # R , xRn, atunci matricea indusă este, după teorema 5 (paragraful precedent) :

Fie z # (1 ,2 ,…,n)T oarecare cu .

Atunci :

Doarece , putem alege pe suficient de mic așa înc=t să avem C (A) ^ . Cu aceasta am demonstrat că a doua inegalitate din relația (1) are loc pentru A R .

Observație : Această teoremă ne arată că, pentru fiecare matrice, raza sa spectrală este marginea inferioară a tuturor normelor naturrale.

Pe baza acestei teoreme vom demonstra următoarea teoremă, teoremă care este rezultatul principal privind convergența la matricea nulă a șirurilor de puteri a unei matrici.

Teorema 2 : Pentru ca este suficient ca A < 1, măcar pentru o normă a matricei A și este necesar și suficient ca (A) < 1.

Demonstrație :

Dacă A < 1, după a patra axiomă a normei matriceale, avem : .

Să presupunem că , deoarece avem : dacă este valoare prorpie a matricii A, atunci n este autovaloare pentru An, deci urmează că , care tinde la zero (pentru n tinz=nd la ).

Deci am demonstrat necesitatea teoremei : (A) < 1.

Reciproc, dacă (A) < 1, fie (A) < < 1 și # – (A) > 0.

Conform proprietății 1 avem că există o nomră naturală cu proprietatea :

A (A) ^ # < 1, deci avem : .

Deci am demonstrat teorema.

Definiție : Fie A Mn,n, se definește (1) ca fiind serie de puteri matriceală.

Convergența seriei (1) se definește prin convergența șirului sumelor parțiale

pentru n .

Definiție : 1). Dacă șirul sumelor parțiale Sn este convergent atunci seria (1) este convergentă.

2). Dacă șirul sumelor parțiale Sn este divergent atunci seria (1) este divergentă.

Teorema 3 : Pentru convergența seriei (1) este suficient ca A < 1, măcar pentru o normă a matricei A și este necesar și suficient ca (A) < 1.

În caz de convergență, avem : det (I – A) 0 și .

Deomonstrație :

Pentru n, p N , avem :

Dacă A < 1, avem că : , p N

Deci șirul Sn este fundamental, adică Sn este convergent.

Să presupunem că seria este convergentă deci și conform teoremei 2 avem că (A) < 1.

Reciproc, dacă (A) < 1, toate valorile proprii i ale matricii A au modulul subunitar ( i < 1, i N), deci 1 nu este valoare proprie pentru matricea A. Ceea ce înseamnă că det (A – I) 0.

Deci, avem că I – A este nesingulară și deci există (I – A)-1.

Să demonstrăm egalitatea : (I – A)(I ^ A ^ A2 ^…^ An) # I – An^1, pentru n N.

Avem :

(I – A)(I ^ A ^ A2 ^…^ An) # I ^ A ^ A2 ^ … ^ An – A – A2 – …- An – An^1 # # I – An^1.

Deci are loc egalitatea de mai sus și o vom înmulții la st=nga cu (I – A)-1 și vom obține :

I ^ A ^ A2 ^…^ An # (I – A)-1 ( I – An^1).

Avem : Sn # (I – A)-1 – (I – A)-1 An^1.

Conform teoremei 2 avem că : . Deci avem .

Teorema 4 : Dacă A < 1 pentru o normă a matricei A, atunci I ^ A și I – A sunt matrici nesingulare și au loc inegalitățile :

Demonstrație :

Deoarece (A) A < 1, atunci valorile 1 și –1 nu pot fi valori proprii ale matricii A. În acest caz, avem : det (A I) 0. Deci matricile (I – A) și (I ^ A) sunt matrici nesingulare, ceea ce ne arată că avem existența matricilor (I – A)-1 și (I ^ A)-1.

Cum () I # (A I) (A I)-1 vom avea :

(1)

Tot din relația (), avem :

(2)

Din relațiile (1) și (2) obținem inegalitățile din enunțul teoremei :

Capitolul al IIIlea. Rezolvarea sistemelor de ecuații algebrice liniare

1). Noțiuni generale privind rezolvarea sistemelor liniare.

Prezentare generală

Rezolvarea sistemelor de ecuații liniare este una din cele mai frecvente aplicații care apar în calculul numeric. Există o mare clasă de fenomene care au ca model matematic sistemele de ecuații sau generează sisteme de ecuații. Metodele de rezolvare a sistemelor de ecuații liniare pot fi grupate în trei categorii :

cu ajutorul determinanților

cu ajutorul metodei eliminării

cu ajutorul tehnicilor iterative.

Din punct de vedere al calculelor efectuate cu ajutorul calculatorului prima metodă de rezolvare cu ajutorul determinanților este greoaie, în mod curent fiind utilizate metodele de eliminare și cele iterative.

În contextul de mai sus cuv=ntul “liniar” este de obicei luat în sensul că variabilele apar la puterea înt=ia în fiecare termen al funcției. Astfel

f(x,y) # a1x ^ a2y și g(x,y) # b1x^ b2y

Proprietățile pe care le oferă liniaritatea unei funcții sau ecuații sunt utilizate ca proprietăți de bază pentru metodele de eliminare utilizate la rezolvarea numerică a sistemelor de ecuații liniare. Ecuațiile sistemului inițial pot fi amplificate prin constante și combinate în așa fel, cu ajutorul operațiilor elementare, obțin=ndu-se în final o formă mult mai simplă a sistemului de ecuații considerat.

Există numeroase situații dependente de natura coeficienților unui sistem pentru care soluția nu există sau o soluție cu un anume grad de precizie nu poate fi găsită. În cazul unui sistem de două ecuații cu două necunoscute, în care reprezentarea celor două ecuații este materializată prin două drepte paralele în planul (x1 , x2) al necunoscutelor, conduce la existența unui șir infinit de soluții.

Există alte situații în care o soluție exactă a sistemului e practic imposibil de a fi determinată; astfel de sisteme se găsesc sub denumirea de sisteme “incorect condiționate” sau sisteme “singulare”. Calitatea unui sistem de a fi “incorect condiționat” sau “aproape singular” este pusă în evidență de natura coeficienților matricei asociate sistemului considerat.

Anticip=nd, se poate spune că un sistem de ecuații liniare permite obținerea unei soluții unice și cu un grad înalt de precizie dacă elementele diagonalei principale din matricea asociată sistemului sunt elemente dominante.

Pentru rezolvarea numerică a sistemelor de ecuații algebrice liniare, există o serie de metode, dar unele din aceste metode sunt impracticabile din punct de vedere al calculatorului, spre exemplu metoda matricei inverse care în cazul sistemelor mari de ecuații algebrice liniare, implică un volum mare de operații și o capacitate mare de memorie a sistemului de calcul. În plus, datorită secvenței mari de calcule ce trebuiesc executate este posibilă scăderea preciziei rezultatelor finale.

La alegerea unei metode de calcul pentru o anumită aplicație și un sistem de calcul dat trebuie avute în vedere o serie de criterii cum ar fi :

Care este numărul de operații aritmetice necesare pentru aplicația respectivă și metoda de calcul numeric aleasă.

Care va fi precizia rezultatelor finale pentru metoda aleasă.

Cum poate fi testată precizia calculelor prin verificări intermediare.

În funcție de aceste criterii, pentru o anume aplicație dată se alege metoda de calcul numeric care poate satisface în măsura cea mai mare toate cele trei deziderate.

Pentru rezolvarea sistemelor de ecuații algebrice liniare și neomogene se pot utiliza două tipuri de metode :

Metode iterative.

Metode directe.

2). Metode exacte de rezolvare a sistemelor de ecuații algebrice liniare. Prezentare generală

Am văzut în paragraful anterior că, metodele de rezolvare a sistemelor algebrice liniare pot fi grupate în două clase :

Metode iterative.

Metode directe.

Metodele directe de rezolvare a sistemelor liniare ne furnizează soluția într-un număr finit de pași. Ele se rezumă la rezolvarea sistemelor liniare algebrice scrise sub forma Ax # b, sisteme ce pot fi rezolvate direct x # A-1b sau cu metodele Cramer, Gauss, a radicalului, etc.

Metodele exacte permit obținerea soluției exacte a sistemului considerat, făc=nd abstracție de erorile de rotunjire și trunchiere, folosind un nunmăr finit de operații elementare.

Pentru sistemele scrise sub forma Ax # b cu soluția este x # A-1b, avem următoarele metode directe :

1). Metoda Gauss a eliminării;

2). Descompunerea LU (unde matricea A #LU) – matricea A-1 # U-1L-1 (metodă ce nu se poate aplica pentru lii # 0, i # 1,n-1)

3). Descompunerea QR (unde matricea A # QR) – matricea A-1 # R-1Q-1 (unde Q – matrice unitară, R – matrice superior triunghiulară)

4). Metoda lui Gauss–Jordan , directă sau prin pivotare.

5). Metoda lui Gauss – Newton.

3). Metode iterative de rezolvare a sistemelor

de ecuații algebrice liniare. Prezentare generală.

Deoarece metodele directe de rezolvare a sistemelor liniare nu ne permit să obținem precizia pe care o dorim, vom folosi așa numitele metode iterative care generează un șir de aproximații, soluția sistemului liniar algebric fiind limita acestui șir.

Metodele iterative construiesc șiruri de vectori ce converg la soluția exactă, care se obțin fără a modifica forma matricii. Astfel, dacă se calculează suficient de mulți vectori din șir, ne apropiem oric=t de mult de soluția exactă, obțin=nd aproximări ale acesteia în orice precizie.

Utilizarea calculatoarelor electronice ne oferă posibilitatea folosirii acestor metode iterative la rezolvarea anumitor tipuri de probleme. În cadrul unei metode iterative se alege o aproximație inițială și succesiv se îmbunătățește această aproximație a soluției (prin iterare) în așa fel ca șirul de soluții îmbunătățite să conveargă către soluția problemei considerate.

Astfel de metode sunt destul de simple și foarte atractive pentru a fi utilizate pe un calculator electronic, cu toate că implică destul de multe operații aritmetice, lucru care nu devine un impediment av=nd în vedere viteza de calcul pentru calculatoarele actuale.

O trăsătură importantă a metodelor iterative este că acestea sunt “auto-corectoare”. În timp ce metodele directe, cel puțin teoretic, converg într-un număr finit de etape, metodele iterative necesită un număr infinit de etape pentru convergență.

La realizarea metodelor iterative trebuie avute în vedere următoarele elemente :

metoda să fie convergentă,

să se poată determina viteza de convergență,

să stabilească criterii pentru stoparea procesului iterativ în momentul obținerii unei aproximări acceptabile pentru soluția problemei considerate.

Exemple de metode iterative :

Metoda Jacobi.

Metoda Gauss – Seidel.

Metoda gradientului (metoda celei mai rapide descreșteri).

Metoda gradientului conjugat.

Metoda Monte Carlo.

Capitolul al IV lea . Metode iterative de rezolvare a

sistemelor algebrice de ecuații algebrice liniare

Partea A.

1). Formule de iterare

O categorie de metode de rezolvare a sistemelor de ecuații liniare, frecvent utilizate, o constituie metodele de determinare a soluțiilor prin aproximări succesive sau metode iterative.

Fie aplicația f : A B, A Rn, B Rn

Fie sistemul (sau ecuația) : f (x) # 0 (1)

Sistemul (1) are, într-o vecinătate V , o rădăcină unică x # .

Fie x(0) – aproximarea inițială a rădăcinii sistemului (1)

Definiție : Se numește formulă de iterare (șir de iterare) pentru determinarea rădăcinii un șir de forma generală :

x(k) # k (x(0), x(1), …, x(k-1)), satisfăc=nd condiția :

Observație : Funcția k depinde de f și, uneori, de rangul termenilor din șirul de iterare.

Pentru determinarea rădăcinii x # a sistemului f(x) # 0 prin metode iterative, în mod frecvent se determină o ecuație de forma x # (x), echivalentă cu sistemul inițial cel puțin într-o vecinătate V , adică :

Convergența metodelor de iterare, existența și unicitatea soluției x # pentru o ecuație x # (x) se bazează pe teorema contracției.

2). Convergența șirurilor de iterare

Fie A : E E; A – operator, E – spațiu metric complet.

Definiție : Operatorul A realizează o transformare contractantă, dacă :

Ax E, x E, , 0 < <1 așa înc=t x,y E (Ax, Ay) (x,y).

Teorema contracției : Fie E spațiu metric complet, A : E E o transformare contractantă a spațiului E pe E; atunci

x E – unic, așa înc=t Ax # x

și

, unde x(k) # Ax(k-1); .

Observație : Teorema afirmă că transformarea contractantă A : E E are un singur punt fix x E.

Demonstrație :

Vom arăta că șirul {x(k)}kN este șir Cauchy.

Fie k > m, x(0) E.

Din A – transformare contractantă, avem :

Aplicăm succesiv și obținem :

(1)

Din relația : rezultă :

aplicăm acesteia rela-ția (1) și obținem :

(pentru 0 < < 1)

Deci am obținut : și aplic=nd-o pentru relația (1) vom obține :

.

Deoarece 0 < < 1 avem , rezultă :

(dar aceasta este tocmai definiția unui șir Cauchy). Deci șirul {x(k)}kN este șir Cauchy și cum E este spațiu metric complet rezultă că șirul este convergent : .

Din convergența șirului {x(k)}kN avem :

Dar : , atunci avem :

Trec=nd la limită în relația : x(k) # Ax(k-1) obținem : x# Ax.

Unicitatea :

Fie x,y E cu x y și y # Ay(adică (x,y) 0 )

Avem :

(contradicție cu ipoteza 0 < < 1). Deci x este unica soluție pentru ecuația x # Ax, x E.

Teorema 2 : Fie – spațiu metric complet, E Rn;

x’, x” E se satisface : () ((x’), (x”) (x’,x”), cu 0 < < 1,

iar () (x0,(x0) ) (1- ) r .

Atunci ecuația x # (x) are pe E soluție unică x # , , unde x(k) # (x(k-1)) , iar viteza de convergență a șirului (x(k))kN către este caracterizată de relația :

(x(k) , ) k (x(0) , ).

Observație : Relația (x(k) , ) k (x(0) , ) reprezintă evaluarea erorii.

Demonstrație :

Avem :

(s-a folosit relațiile () și () din ipoteză)

Deci : pentru x E, avem (x) E .

Aplicația : E E este o transformare contractantă deoarece :

(din relația () din ipoteză), x,y E

Din : – transformare contractantă și E – spațiu metric complet obținem :

{irul de iterare pentru rădăcina x # este dat de : x(k) # (x(k-1)), k N, iar convergența sa rezultă din teorema contracției :

cu (0 < < 1) și .

Trec=nd la limită în x(k) # (x(k-1)) rezultă # ().

Evaluarea erorii este dată de : (x(k) , ) k (x(0) , ) .

Ordinul unei formule de iterare.

Fie : E E, E R, E – spațiu metric complet.

Fie x # soluție a ecuației x # (x) dată de șirul : x(k) # (x(k-1)) (1) , k N.

Definiție : Iterarea definită prin șirul (1) are ordin m dacă

(2)

Observație : Viteza de convergență crește odată cu ordinul formulei.

Teoremă : Dacă (2) este iterare de ordinul m, atunci viteza de convergență a șirului irului {x(k)}kN către este caracterizată de relația :

Demonstrație :

Scriem dezvoltarea în serie Taylor a aplicației (x) în punctul x # :

{tiind că avem o iterare de ordinul m, atunci conform definiției acesteia avem :

Vom ține seama de # () și făc=nd înlocuirile x # x(k-1), (x(k-1) # (x), din ultima relație vom obține :

Vom presupune , pentru x E, atunci rezultă :

Transcriind sucesiv avem :

Am obținut :

Să vedem dacă < 1.

Luăm aproximarea inițială x(0) așa înc=t aceasta să fie suficient de apropiată de soluția x # , atunci din modul de alegere a lui obținem < 1. Cu acesta am demonstrat teorema.

3). Propagarea erorilor

În paragrafele precedente am văzut cum se definește un șir de iterare precum și în ce condiții acesta converge către o soluție unică x # , cu .

Am mai văzut că șirul de iterare (x(k)) dat de : (2) x(k) # k (x(0), x(1), …, x(k-1)) are viteza de convergență caracterizată prin relația :

(3) (x(k) , ) k (x(0) , ), cu 0 < < 1.

În aceste condiții, relația (3) reprezintă evaluarea erorii.

Av=nd în vedere că toate aceste calcule se execută cu ajutorul calculatorului electronic, nu este posibil ca funcția (x) să fie evaluată exact.

Pentru orice a se poate reprezenta aproximarea lui (x) prin (x) # (x) ^ (x), unde (x) este eroarea comisă în evaluarea funcției.

De obicei se cunoaște o margine pentru (x), adică .

În acest caz schema iterativă care se utilizează poate fi reprezentată astfel :

(4) x(k^1) # (x(k)) ^ (k), , unde x(k) sunt valorile obținute din calcul și (k) satisface relația : , .

Din figura de mai jos, se poate observa că, pentru un caz particular c=nd (x) # ^ (x – ) (cu 0 < < 1), eroarea în rădăcina este mărginită prin .

Se observă că, dacă 1, problema este slab condiționată.

În cazul în care schema de iterație este convergentă, prezența erorii în calculul funcției (x), de mărime mărginită prin , face ca schema iterativă să estimeze rădăcina cu o imprecizie mărginită prin .

Fie x(0) orice valoare astfel ca , unde .

În acest caz iterația x(k), calculată prin relația (4) cu o eroare mărginită prin , se află în intervalul și are loc :

(5) , unde .

Din afirmația precedentă se vede că eroarea de calcul care apare în evaluarea lui (x) este cel mult .

Numărul de iterații necesar este : (6) .

Desigur dacă eroarea acceptabilă este mai mare de , numărul iterațiilor dat de relațiile (6) este estimat adecvat.

4). Metode iterative staționare și nestaționare

În paragraful 1) am văzut cum definim o formulă de iterare (șir de iterare) :

(2) x(k) # k (x(0), x(1), …, x(k-1)), satisfăc=nd condiția :

Aplicația k depinde de funcția f și de k (rangul termenilor din șirul de iterare)

Definiție : 1). O astfel de metodă, în care s-a definit un șir de iterare (2), se numește metodă iterativă nestaționară.

2). Dacă pentru fiecare k, funcția k depinde cel mult de una din variabilele x(k-s^1), x(k-s^2),…., x(k-1), metoda iterativă se numește metodă iterativă cu un pas.

3). Dacă funcția k nu depinde de k, metoda iterativă se numește metodă iterativă staționară.

Exemple : 1). Metode iterative nestaționare : – metoda bisecției

– metoda poziției false

2). Metode iterative staționare : – metoda secantei (cu doi pași)

O clasă importantă de metode iterative este clasa metodelor staționare cu un pas.

În acest caz se alege o singură valoare de start x(0) (aproximație inițială) și o funcție (x), iar șirul x(1), x(2), x(3),… este calculat cu ajutorul relației :

(3) x(k^1) # (x(k)),

Funcția (x) se numește funcție de iterare.

Consider=ndu-se o funcție iterativă (x), există o serie de întrebări care trebuie puse în legătură cu metoda iterativă respectivă :

{irul construit prin (3) converge către o limită unică.

Dacă este soluția exactă și valoarea de start x(0) este suficient de aproape de , șirul construit prin (3) converge la .

Toate acestea le vom avea în vedere atunci c=nd vom prezenta exemple de metode iterative (metoda Jacobi, metoda Gauss-Seidel, ș.a.)

5). Schema generală a unei metode iterative

de rezolvare a sistemelor liniare algebrice

Am văzut, în paragrafele anterioare, ce înseamnă un șir de iterare pentru determinarea rădăcinii unui sistem, în ce condiții acesta converge la soluția unică a sistemului, precum și posibilitatea găsirii unei soluții aproximative.

{tim că pentru a folosi metode iterative trebuie să construim un șir de vectori care converg la soluția exactă, fără, însă, să modificăm forma matricei date.

Fie sistemul de ecuații liniare :

(1) , sau scris sub formă matriceală : (1’) Ax # b, unde matricea A este matrice pătratică nesingulară (det A 0).

În continuare vom folosi forma (1’) pentru a ușura calculele.

O clasă foarte largă de metode iterative putem obține dacă scriem matricea sub forma : (2) A # B – C , unde B este matrice nesingulară (de regulă diagonală sau triunghiulară).

Din relațiile (1’) și (2) avem : (B – C)x # b, pentru orice matrice B scrisă sub forma : (3) Bx # Cx ^ b.

Forma (3) ne va permite definirea unei metode iterative în modul următor :

Pentru orice vector x(0) Rn, se determină termenii succesivi ai unui șir {x(k)}k N prin recurența :

(4) Bx(k^1) # Cx(k) ^ b,

Deoarece matricea B este nesingulară și ușor de inversat, pentru fiecare x(k) deja determinat, se obține următorul vector x(k^1) rezolv=nd sistemul (4).

Dacă vom impune condiții care asigură convergența șirului {x(k)}k N pentru k , limita acestui șir va satisface (3) Bx # Cx ^ b, adică echivalent cu sistemul dat Ax # b, deci va fi soluția acestui sistem.

În acest mod sistemul Ax # b se va reduce la rezolvarea unui șir de sisteme cu matricea B.

În continuare vom vedea în ce condiții șirul {x(k)}k N este convergent. Vom prezenta următoarea teoremă :

Teoremă : Condiția necesară și suficientă ca șirul {x(k)}k N , obținut cu iterația Bx(k^1) # Cx(k) ^ b, , pentru x(0) Rn, să fie convergent, este ca raza spectrală a matricei M # B –1C să fie subunitară : (7) (M) < 1.

Dacă măcar o normă a matricei M satisface : (8) M < 1 atunci șirul {x(k)}k N , este convergent, x(0) Rn.

Demonstrație : Fie x – soluția exactă a sistemului Ax # b.

Atunci ecuațiile : Bx # Cx ^ b se vor scrie sub forma : x # Mx ^c

Bx(k^1) # Cx(k) ^ b x(k^1) # Mx(k) ^c

B-1 / Bx # Cx ^ b x # B-1Cx ^ B-1b.

Analog : B-1 / Bx(k^1) # Cx(k) ^ b x(k^1) # B-1Cx(k) ^ B-1b.

În aceste ultime două relații vom nota M # B-1C ; c # B-1b.

Scăz=nd cele două relații : : x # Mx ^c vom obține : x9k^1) – x # M (x(k) – x)

x(k^1) # Mx(k) ^c

Dăm valori lui k :

k # k – 1 x(k) – x # M (x(k-1) – x)

k # k – 2 x(k-1) – x # M (x(k-2) – x) x(k) – x # M2 (x(k-2) – x)

………..

k # k – k x(k) – x # M (x(0) – x)

Vom obține : (9) x(k) – x # M (x(k-1) – x) # M2 (x(k-2) – x) # … # Mk (x(0) – x)

Ultima relație obținută arată că pentru a avea convergența lui{x(k)}k N , , este necesar și suficient ca (ceea ce este adevărat pe baza proprietății 2 de la paragraful 3, capitolul II).

Definiție : Vom numi matricea M – matricea iterației (datorită modului de definire M # B-1C).

Convergența șirului iterațiilor {x(k)}k N , , nu e suficient pentru studiul unei metode iterative.

Pentru acest lucru, vom fi interesați în continuare și de rapiditatea convergeței (viteza de convergență) a șirurilor {x(k)}k N la x pentru k .

În acest scop vom evalua numărul de iterații necesare pentru a obține o anumită valoare relativă.

Vom alege o normă matriceală naturală M suficient de apropiată de (M).

Conform proprietății 1, de la paragraful 3, capitolul II, pentru un oarecare , există o astfel de normă așa înc=t să avem :

M (M) ^

Folosind o normă vectorială compatibilă cu această normă matriceală, vom obține relația : (10) x(k) – x # Mk (x(0) – x) (11)

Vom presupune că > 0 este destul de mic așa înc=t , dacă (12) ș(M)ț < 10 –m să avem și (13) ș(M) + țk 10-m de la un k încolo.

Deci, dacă relația (12) este satisfăcută, vom avea, din (11) : (13)

Sau , scrisă altfel : (14) . Ceea ce ne arată că o eroare relativă nu va depăși 10 –m .

Acest lucru ne va asigura că, pentru o iterație inițială exactă numai în ce privește prima zecimală, vom obține după k – iterații o aproximare cu m – cifre exacte.

Observație : 1). Condiția (12) ș(M)ț < 10 –m ne dă, pentru (Mț < 1(adică în condiții de convergență) :

(15) , unde se numește rata spectrală de convergență.

2). Cu c=t R este mai mare cu at=t sunt suficiente mai puține iterații pentru a obține precizia 10 – m .

În practică este, în general, mai greu de calculat raza spectrală (M) implicit și rata spectrală de convergență R.

O evaluare a erorii mai practică se obține astfel :

Din definiția iterațiilor

x(k) – x # M (x(k-1) – x) # M2 (x(k-2) – x) # … # Mk (x(0) – x) vom avea : x(k^j^1) – x(k^j) # M (x(k^j) – x(k^j-1) ) # M2 (x(k^j-1) – x(k^j-2)) # … # Mj (x(K^1) – x(k))

Deci :

În ipoteza că (M) < 1, pentru p , se va obține :

În ipoteza că M < 1vom avea evaluarea (pe baza proprietății 4, de la paragra-ful 3, capitolul II) :

, relație ce folosește numai cantități cunoscute : iterațiile x(k-1), x(k) și matricea M.

Partea B. Exemple de metode iterative

1). Metoda lui Jacobi

a). Algoritmul Jacobi

Fie un sistem algebric de ecuații liniare de forma : (1) Ax # B.

Notăm : – elementele matricei A.

– termenii liberi.

Pentru sistemul (1) avem o descompunere de forma (2) Bx # Cx ^ b.

În această descompunere alegem matricea B # diag (a11 , a22 , …, ann)

{tim că, după schema generală a metodei iterative, avem definită o metodă iterativă în modul următor :

(3) Bx(k^1) # Cx(k) ^ b,

Vom scrie relația (3) sub formă matriceală :

(3’)

În relațiile (3’), xi(k) și xi(k^1) sunt componentele vectorilor x(k) și x(k^1).

Dacă în relațiile (3’) facem calculele vom obține :

Analog, vom obține pentru valoarea i # 2 :

Prin inducție, vom obține pentru i # n, relațiile :

Deci, metoda iterativă corespunzătoare este de forma :

(4)

Vom pune o condiție foarte importantă :

(5) aii 0, – ceea ce ne asigură că matricea B este nesingulară (det B 0).

În aceste condiții, recurența (4) împreună cu condiția (5), vom obține relațiile :

(6) , .

Metoda iterativă dată de relațiile (6) se numește metoda lui Jacobi.

b). Condiții de convergență pentru metoda lui Jacobi

Teoremă : Condiția necesară și suficientă de convergență a metodei iterative (6) este ca modulele valorilor proprii ale matricei iterației să fie subunitare.

Demonstrație :

Conform teoremei de la paragraful 4, partea A, avem că :condiția necesară și suficientă ca șirul {x(k)}k N să fie convergent este ca : (M) < 1 (unde M # B-1 C).

Cum, ecuația caracteristică a matricei M este : det ( I – B-1 C) # 0, sau echivalent cu det (b-1 (B – C) # 0 și det B 0, atunci condiția necesară și suficientă de convergență a șirului {x(k)}k N este ca toate rădăcinile ecuației det (B – C) # 0 să se găsească în cercul de rază 1 din planul complex și cu centrul în origine ( z # 1).

{tiind că B # diag (a11 , a22 , …,ann), vom avea :

(7)

Deci, condiția necesară și suficientă de convergență a metodei iterative dată de relațiile (6) va fi : ca toate rădăcinile ecuației (7) să se găsească în cercul de rază 1 și cu centrul în O(0,0) din planul complex.

Aceste condiții au, însă, numai o valoare teoretică.

În practică vom folosi numai condiții de suficiență, care le vom obține din teorema de la paragraful 4, partea A :

Dacă măcar o normă a matricei M satisface : (8) M < 1 atunci șirul {x(k)}k N , este convergent, x(0) Rn – o aproximație inițială.

{tiind că ; Bx # Cx ^ b x # B-1 Cx ^ B-1b, avem o metodă iterativă de forma : x(k^1) # Mx(k) ^ B-1 b.

Această ultimă relație împreună cu relațiile (6) ne va da elementele matricei M :

Deci, elementele matricei M sunt date de relațiile : (7)

Cu ajutorul definițiilor normelor putem obține diferite condiții de convergență :

Din M 1 < 1 , se va obține condițiile de convergență : (9) , .

Iar pentru norma duală, condiția M *1 < 1 , vom obține condițiile de convergență :

(10) , – condiție care poartă numele de dominanța diagonalei pe linii.

Să vedem acest lucru : dominanța diagonalei pe linii se obține din condiția : M < 1 ( adică : M *1 < 1 ) :

Pentru norma euclidiană M E < 1 , vom obține condițiile de convergență :

(11)

Putem obține și alte criterii de convergență, folosind normele ponderate ( 1,P , 2,P , ,P ) :

Fie constantele p1, p2 , …, pn pozitive astfel ca :

(12) ,

sau

(13)

sau, încă

(14) ,

Condiția (14) se mai poate rescrie, renot=nd , sub forma :

(15) ,

Lu=nd, în particular, în criteriul (12), condiția ca , , vom obține relațiile :

(16) , – care se numește dominanța diagonalei pe coloane.

Se pot obține și alte criterii de convergență folosind diferite norme matriceale.

Observație : Să mai observăm că la metoda lui Jacobi, deoarece M # B-1 C, unde B # diag (a11 , a22 , …, ann ), putem scrie matricea M și sub forma : M # – (L ^ U), unde L și U sunt matrici inferioară, respectiv superioară triunghiulară.

1). Metoda Gauss – Seidel

a). Algoritmul Gauss – Seidel

Am văzut că metoda Jacobi (dată de relațiile (6) de la paragraful 1 al acestui capitol) folosește pentru calculul fiecărei componente xi(k^1) pe cale ale lui x(k), deși x1(k^1), x2(k^1), …, xi-1(k^1) sunt deja calculate.

Acest lucru va necesita stocarea în memoria calculatorului at=t a vectorului x(k), c=t și a lui x(k^1), pe măsură ce se calculează componentele lui.

Metoda lui Gauss – Seidel este o modificare a metodei Jacobi, în care componentele deja calculate vor fi folosite în această iterație.

Fie un sistem algebric de ecuații liniare de forma : (1) Ax # B.

Notăm : – elementele matricei A, punem condiția ca aii 0,

– termenii liberi.

Pentru sistemul (1) avem o descompunere de forma (2) Bx # Cx ^ b.

În această descompunere vom lua matricea B – matrice inferior triunghiulară :

și C # B – A – matrice superior triunghiulară cu zerouri pe diagonală.

Cu aceste alegeri, sistemul (2) se va scrie sub forma :

(2’)

În relațiile (2’), xi sunt componentele vectorului x.

Vom face calculele 1n relațiile (2’) și impărțim rezultatul la aii , vom obține :

Vom împărți rezultatul la a11 și vom avea :

Proced=nd analog, vom obține pentru valoarea i # 2 :

Să încercăm și pentru valoarea i # 3 :

Prin inducție, vom obține pentru i # n, relațiile :

Împărțind fiecare relație obținută la aii , vom obține formulele :

(3) , cu aii 0,

}in=nd cont că avem definită o metodă iterativă de forma Bx(k^1) # Cx(k) ^ b, , și relațiile (3), vom avea în acest caz o iterație dată de :

(4) , cu aii 0,

Metoda iterativă dată de relațiile (4) se numește metoda Gauss – Seidel.

Observație : Metoda lui Gauss – Seidel generează șiruri convergente pentru orice iterație inițială, dacă și numai dacă rădăcinile ecuației det(I – B-1 C) # 0 se găsesc în interiorul cercului unitate din planul complex.

Conform teoremei de la paragraful 4, partea A, avem că :condiția necesară și suficientă ca șirul {x(k)}k N să fie convergent este ca : (M) < 1 (unde M # B-1 C).

Cum, ecuația caracteristică a matricei M este : det ( I – B-1 C) # 0, sau echivalent cu det (b-1 (B – C) # 0 și det B 0, atunci condiția necesară și suficientă de convergență a șirului {x(k)}k N este ca toate rădăcinile ecuației det (B – C) # 0 să se găsească în cercul de rază 1 din planul complex și cu centrul în origine ( z # 1).

{tiind că B # este o matrice inferior triunghiulară vom avea :

Deci, condiția necesară și suficientă de convergență a metodei iterative dată de relațiile (6) va fi : ca toate rădăcinile ecuației (7) să se găsească în cercul de rază 1 și cu centrul în O(0,0) din planul complex.

Rădăcinile ecuației de mai sus sunt însă greu de calculat în practică, așa înc=t vom căuta diverse condiții suficiente de convergență.

b). Condiții de convergență pentru metoda Gauss – Seidel

Teorema 1 : Dacă are loc criteriul (9) de la metoda lui Jacobi ( , ) – criteriu de dominanță a diagonalei pe linii, atunci at=t metoda lui Jacobi, c=t și metoda lui Gauss – Seidel generează șiruri convergente pentru orice iterație inițială.

Demonstrație :

Convergența șirurilor generate de metoda Jacobi a fost demonstrată anterior.

Mai răm=ne să arătăm convergența șirurilor generate de metoda Gauss – Seidel.

Avem :

Scădem aceste relații și vom obține :

Evaluăm modulul :

Să notăm : , atunci conform acestor notații din ultima relație vom obține :

Am văzut că pentru norma duală, condiția M *1 < 1 , se obține condițiile de convergență : , (dominanța diagonalei pe linii)

atunci din ultima relație și folosind inegalitatea de indice i0 pentru care : vom obține :

Vom nota :

Atunci avem :

În aceste condiții din relația () rezultă : ceea ce ne asigură că șirul {x(k)}kN este convergent pentru orice x(0).

Observație :

Matricea M a iterației Jacobi are norma : M 1* # .

Atunci, după relația :

x(k) – x # M (x(k-1) – x) # M2 (x(k-2) – x) # … # Mk (x(0) – x) pentru un șir { (k)}kN generat de metoda Jacobi vom avea în ipoteza : , relația : .

Compar=nd cu și țin=nd seama de rezultă că pentru aceleași iterații inițiale , metoda Gauss – Seidel converge cel puțin la fel de rapid ca metoda Jacobi, eventual chiar mai rapid.

Pentru a compara riguros cele două metode din punctul de vedere al rapidității, ar fi însă necesară compararea razelor spectrale ale matricelor lor de iterații.

În anumite condiții, se poate arăta că (M) # (M)2 (unde am notat M matricea iterației Gauss – Seidel și M matricea iterației Jacobi); deci rata de convergență a metodei Gauss – Seidel este dublul ratei de convergență a celeilalte metode.

Consecința 1 : Rezultatul teoremei răm=ne valabil, dacă se înlocuiește ipoteza (9) cu criteriul care folosește normele ponderate.

Demonstrație :

Conform teoremei anterioare avem :

pi , pi > 0

()

Folosind următoarele ipoteze :

– ,

– norma ponderată :

– inegalitatea de indice io pentru care avem :

– notațiile:

atunci relația () devine :

Notăm :

Atunci avem :

În aceste condiții din relația () rezultă : , ceea ce ne asigură că șirul {x(k)}kN este convergent pentru orice x(0).

Consecința 2 : De asemenea, convergența șirurilor generate de cele două metode are loc, dacă matricea A satisface condiția de dominanță a diagonalei pe coloane (adică criteriul de dominanță a diagonalei pe coloane (16), de la metoda Jacobi).

Demonstrație :

Aplicăm teorema precedentă unei matrici H # (hij), cu hij # an^1-i, n^1-i pentru toți i și j.

Condițiile de convergență : ( , ) – criteriu de dominanță a diagonalei pe linii – pentru matricea H este echivalent cu criteriu de dominanță pe coloane : , pentru matricea A.

Conform teoremei 1 avem că : toate rădăcinile ecuațiilor :

() și () sunt subunitare în modul.

Observăm că determinanții din () și (), printr-un șir de permutări de linii și coloane, se reduc la determinanți de forma :

(’) și (’).

Deci are loc convergența celor două metode și pentru matricea A.

Luăm în condițiile de convergență : , , și vom obține : , (din pi >0 ci >0 ).

Vom înmulți inegalitatea () din teorema 1 cu și rațion=nd analog ca și la consecința 1, vom obține valabilitatea criteriului din enunț.

Definiție : Spunem că matricea A este ireductibilă, dacă există indicii i1, i2, …, im, jm^1, jm^2, …jn care formează o permutare a indicilor 1, 2, …, n așa înc=t :

Observație : 1). Această afirmație arată că există matricile de permutări P și Q așa înc=t să aibă loc :

, unde matricea A1 este de tipul (m x m), A2 este de tipul (n-m)x(n-m) și O este matricea nulă de tipul (n-m)xm.

2). Conform observației 1., un sistem de ecuații liniare cu o matrice de această formă se reduce la două sisteme de matrici A3 și A1, decide ordin mai mic.

Definiție : Spunem că matricea este ireductibilă dacă matricea nu este reductibilă.

Teorema 2 : Fie matricea A ireductibilă. Atunci, pentru ca metodele lui Jacobi și Gauss – Seidel să genereze șiruri convergente pentru orice iterație inițială, este suficient ca una din condițiile următoare să fie îndeplinite :

(5) , , () io , – criteriu de dominanță a

diagonalei pe linii

sau

(6) , , () io , – criteriu de dominanță a

diagonalei pe coloane

Demonstrație :

Să demonstrăm mai înt=i valabilitatea criteriul (5).

Convergența metodei Jacobi.

Vom arăta că rădăcinile ecuației (7) sunt subunitare.

Fie o rădăcină a acestei ecuații.

Observăm că un sistem omogen cu determinantul coeficienților dat de (7) are soluții nebale de forma : z # (z1, z2, …, zn)T 0. Deci are loc : .

Să notăm : , atunci pentru i # s, conform relației precedente, vom avea :

Să presupunem că are loc .

Atunci, avem :

Relația (8) ne arată că avem două cazuri :

Caz 1 : asj # 0 pentru toți s pentru care și pentru toți j pentru care . Dar acest lucru ne arată că matricea A este reductibilă , ceea ce este în contradicție cu ipoteza, deoarece am presupus că z este diferit de zero.

Caz 2 : toate componentele lui z au același modul egal cu . Dar acest lucru ne arată că :

ceea ce este absurd.

Conform celor două cazuri avem că : presupunerea că este falsă.

Deci , ceea ce ne asigură convergența metodei Jacobi.

Convergența metodei Gauss – Seidel.

Fie o rădăcină a ecuației : . Atunci există o soluție nebanală de forma z # (z1, z2, …, zn)T 0 așa înc=t să aibă loc :

Conform notațiilor de mai sus avem :

și deoarece suntem în ipoteza (5), obținem :

Să presupunem că . Atunci, avem : . Dar am văzut aneterior că aceasta conduce la o contradicție.

Deci , ceea ce ne asigură convergența metodei Gauss – Seidel.

Un raționament similar cu cel făcut în consecința 2 ne va asigura valabilitatea criteriului (6).

Un criteriu de convergență corespunzător normei euclidiene este dat de următoarea teoremă :

Teorema 3 : Dacă este satisfăcută inegalitatea (11) (de la metoda Jacobi), atunci șirurile generate de metoda lui Jacobi și de cea a lui Gauss – Seidel converg, pentru orice iterație inițială.

Demonstrație :

Convergența șirurilor generate de metoda Jacobi a fost dovedită anterior.

Să dovedim acum convergența șirurilor generate de metoda Gauss – Seidel.

Conform teoremei 1, avem :

(9)

Vom ridica ambii membri ai relației (9) la pătrat și obținem :

Aplicăm inegalitatea Cauchy – Schwartz, , ultimei relații și obținem :

(10)

Să notăm : și să sumăm după i de la 1 la n relația (10) :

În ultima relație vom schimba ordinea de sumare în membrul drept (i j și j i) și obținem :

(11)

Să facem notațiile :

Conform acestor notații din (11) rezultă :

În membrul st=ng facem i j și vom obține :

Să observăm că :

Dacă notăm norma ponderată cu ponderile pj # 1-j atunci din ultima relație obținem :

ceea ce asigură convergența șirului {x(k)}kN generat de metoda Gauss – Seidel.

3). Metode iterative pentru sisteme liniare algebrice bazate pe minimizarea unei forme pătratice.

În acest paragraf vom formula metode iterative bazate pe minimizarea unei forme pătratice al cărei unic minim este soluția sistemului Ax # b.

Observație : În cazul c=nd matricea A nu este simetrică și pozitiv definită, pentru construcția metodelor iterative, este indicat să se construiască forma :

(1) R # (b – Ax)T∙ (b – Ax) # rT ∙ r

Dacă A mxn(R) atunci (b – Ax)T nxm(R) și (b – Ax) mxn(R) deci produsul lor (b – Ax)T∙ (b – Ax) nxn(R) . Avem că R nxn(R), deci aceasta este o formă pătratică, precum și pozitiv definită (acest lucru se obsrvă din formularea lui R).

Am obținut că matricea R este o matrice pătratică și pozitiv definită care are un minim unic c=nd x este soluția sistemului Ax # b.

Dar metodele de minimizare ale (1) necesită mari eforturi de calcul. Datorită acestui lucru, vom considera în continuare că matricea A este simetrică și pozitiv definită.

Fie un sitem algebric de ecuații liniare de forma :

Ax # b, unde ; ; .

Din (1) se obține : Ax – b # 0, relație pe care o vom înmulți la st=nga cu xT # tx.

Vom obține relația : txAx – txb # 0.

Pentru matricea A pozitiv definită, vom folosi forma :

(2)

Se observă, ușor, că matricea Q este și ea pozitiv definită (din modul de definire a lui Q).

Vom vedea că Q are un unic minim atins în soluția sistemului (1).

Pentru aceasta vom considera doi vectori x,y cu x y.

Avem :

Cum , atunci vom avea :

(Am folosit aici că Ax – b # 0)

Fie x – soluție pentru sistemul (1). Deoarece , ceea ce înseamnă că , avem că Q are un minim în x.

Am arătat că Q admite minim în x – soluția acestui sistem. În continuare vom demonstra că acest minim este unic.

Demonstrația acestui lucru o vom face prin reducere la absurd.

Presupunem că Q admite un alt minim în x’ așa înc=t x x’; acest lucru arată că x’ este soluție a sistemului Ax # b.

Vom avea : Ax’ # b x’ # A-1 b.

Dar și x este soluție a sistemului (1). Deci vom avea : x # A-1b.

Din cele două rezultate vom obține că : x # x’, ceea ce reprezintă o contradicție (în ipoteză am presupus că x x’).

Deci Q admite un unic minim x – soluție a sistemului (1).

În continuare, vom construi un șir de vectori (xi) i 1 așa înc=t să avem Q(xi^1) < Q(xi), cu Q (xi) 0, .

Deoarece Q (xi) 0, și Q admite un unic minim, atunci un astfel de șir converge la soluția sistemului Ax # b.

Tehnica de minimizare a lui Q este următoarea :

Pasul 1 : se alege un vector de start x(1) ;

Pasul 2 : în direcția dată v1 , la distanța 1 , se consideră vectorul x(2) așa înc=t să avem :

x(2) # x(1) ^ 1v1 ;

Proced=nd inductiv și pentru următoarele pasuri vom obține pentru n :

Pasul n : (3) x(n) # x(n-1) ^ nvn .

Această tehnică de minimizare a formei Q se poate observa (reprezentată în cazul R2 ) pe figura de mai jos (s-a luat cazul i # 2) :

Observații : 1). Pentru forma Q dată de formula (2), problema de căutare a minimului este o problemă ce se rezolvă explicit.

2). Evident că sunt posibile o infinitate de proceduri de tipul (3).

În cele ce urmează vom considera două tehnici de alegere a distanței și direcției la fiecare pas:

a). metoda celei mai rapide descreșteri (metoda gradientului);

b). metoda gradientului conjugat.

a). Metoda gradientului

(metoda celei mai rapide descreșteri)

Baza acestei metode este ideea naturală de a alege vi în direcția celei mai mari rate de schimbare a lui Q, anume în direcția gradientului și, de asemenea, de a alege i așa înc=t Q să fie minimă.

Cu ajutorul relației (2) , vom calcula gradientul lui Q pentru x # x(i) :

Conform definiției gradientului : , vom avea :

(aici am sumat după i)

Am obținut următoarea relație : , pentru și .

Deci am obținut relația : ,

Această ultimă relație ne arată că : , .

Acum, pentru a găsi i , vom calcula Q ( x(i) – I ri ). În relația (2) înlocuim x cu x(i) – i ri și vom obține :

Cum Ax(i) – b # – ri , atunci vom avea că Ax(i) # b – ri .

Atunci relația de mai sus devine :

Cum b # Ax(i) ^ ri , vom avea :

Deci am obținut relația :

Derivăm parțial în raport cu i ultima relație și obținem :

Egal=nd ultima relație cu 0 (s-a ales i așa înc=t Q să fie minimă), vom obține :

Ultima relație obținută și relația (3) ne va da metoda celei mai rapide cobor=ri (metoda gradientului) :

(4)

Pentru a demonstra că această metodă este convergentă nu trebuie dec=t să arătăm că Q (x(i^1) ) < Q (x(i) ). Aceasta este echivalent cu a arăta că Q (x(i^1) ) – Q (x(i) ) < 0.

Cum x(i^1) # x(i) ^ ivi și not=nd ivi # y, vom obține că x(i^1) # x(i) ^ y.

Trebuie să arătăm că Q (x(i) ^ y) – Q (x(i) ) < 0. Dar acest lucru am arătat anterior.

Observație : Pentru sisteme rău condiționate convergența este înceată.

b). Metoda gradientului conjugat

Procedeul prezentat anterior, adică lucrul după direcția gradientului lui Q, pare o alegere logică. Totuși, dacă alegem cu mai multă grijă direcția vi , se poate garanta convergența într-un număr finit de pași.

Observație : În practică, rotunjirile pot împiedica atingerea convergenței.

Fie – o bază a spațiului Rn.

Fie x – soluția sistemului (1) Ax # b.

Fie x(1) – o aproximație inițială a sistemului (1).

Vom încerca minimizarea lui Q proced=nd în modul următor :

Pasul 1 : – alegem x(1) o aproximație inițială.

– luăm un v1 la o distanță fixată i așa înc=t Q să fie minimă.

Vom avea deci :

x # x(1) ^ iv1 .

Pasul 2 : x # x(1) ^ iv1 ^ iv2 , unde x(2) # x(1) ^ iv1

Proced=nd inductiv vom obține :

Pasul n : x # x(1) ^ iv1 ^ iv2 ^ … ^ ivn sau (5) (sumat după j)

Din (5) putem obține relația (5’)

În cazul în care vectorii vj sunt ortogonali (baza este ortonormată), atunci înmulțim relațiile (5’) la st=nga cu tvj și vom obține :

,

Să mai presupunem că tvjAvj # 0 (), caz în care vj se numesc A–ortogonali sau A–conjugați, vom obține relațiile :

(7)

Cum A (x – x(1)) # Ax – Ax(1) # r1, vom obține :

(7’) ,

Relațiile (7’) ne va da calculul lui j .

Observație : În cazul c=nd tvjAvj 0, adică A – pozitiv definită, atunci pentru calculul lui j se va folosi relațiile (6).

Observație : Metoda gradientului conjugat are ideea de bază : generarea unei secvențe de vectori vj – A ortogonali (), calculul lui j după relațiile (7’) (sau (6) (în cazul c=nd matricea A este pozitiv efinită) și aflarea lui x.

În scopul generării unei mulțimi de vectori vj – A ortogonali, să presupunem că avem o mulțime de vectori u1 , u2 , … , un , vectori care formează o bază.

Folosim procedura de ortogonalizare Gram – Schmidt pentru a calcula vectorii vn :

Notăm : , avem :

(8)

Deoarece vectorii ui sunt liniari independenți la fel sunt și vectorii vi .

Trebuie să calculăm componentele șirului .

Una dintre soluții ar fi ca vectorii să se aleagă drept vectorii bazei canonice. Dacă se procedează în acest mod vom obține o procedură echivalentă (funcțional) cu eliminarea Gauss. Totuși această alegere nu este convenabilă. Vom încerca să construim un algoritm mai convenabil.

Definim , , cu j dați de relațiile (7’).

}in=nd cont de relațiile (3), avem : (9) , .

Am văzut că reziduu este dat de ri # b – Ax(i), (relație ce rezultă din Ax(i) – b # – ri).

Baza metodei gradientului conjugat este de a alege vectorii ui sub forma : ui # ri ,

Observație : 1). {irul – formează o bază ortogonală pentru spațiul Rn.

2). Q (x(i^1)) < Q (x(i)) (demonstrația acestei relații se va face analog ca în paragraful a.)

Vom determina vectorii și printr-un proces iterativ simultan :

Pasul 1 : avem : r1 # b – Ax(1)

Din relațiile (8) avem : v1 # r1 .

Pasul 2 : din relațiile (7’) și (9) vom obține :

Cum r2 # b – Ax(2) # r1 – 1Av1 , din relațiile (8) vom avea : , cu

Am obținut la acest pas :

Proced=nd inductiv vom obține :

Pasul n : pentru acest pas vom obține ecuațiile algoritmului metodei gradientului conjugat :

(10) ,

Observații : 1). Aceste relații (10) împreună cu relațiile (5) vom obține xn^1 # x. Aceasta ne arată că metoda gradientului conjugat este o metodă iterativă finită, convergentă în n – pași.

2). Datorită acestui lucru ea este analoagă metodelor directe.

3). Această metodă poate fi descrisă și organizată mai ușor din punct de vedere al implementării pe calculator.

4). Datorită acumulărilor erorilor, nu se poate obține convergența chiar în n – pași.

O caracteristică a acestei metode este că dacă se continuă procesul iterativ dat de relațiile (10) pentru valoarea i > n, atunci :

, (cu ) p=nă c=nd erorile de rotunjire împiedică aceasta.

Capitolul al V lea . Aplicații practice.

În acest capitol vom prezenta un exemplu de aplicație la metodele prezentate – Ecuația eliptică și condiții la frontieră .

Exemplu 1.

Se consideră o placă de oțel dreptunghiulară care are una dintre laturile mici menținută la temperatura de 1000 C, iar celelalte trei laturi la temperatura de 00 C. Să se studieze distribuția temperaturii în interiorul plăcii.

Prezentăm grafic datele problemei :

Problema Dirichlet a unei ecuații eliptice (ecuația căldurii este ecuație eliptică) – în general :

Fie – mulțime conexă și deschisă. Avem aplicația u definită pe x (0,T), de clasă C2(), dată prin :

(1) – ecuația propagării căldurii

(2)

(3) – condiție inițială , unde

Problema la limită pentru ecuația căldurii :

, unde u0(x) – temperatura.

Pentru problema noastră avem :

Ecuația Laplace în două dimensiuni –

Condițiile Dirichlet :

, unde (conform condițiilor din problemă) :

Această problemă poate fi rezolvată analitic scriind u (x,y) ca o serie Fourier infinită, dar aici vom ilustra o tehnică aparțin=nd analizei numerice.

Aproximarea derivatelor (aproximarea funcției f(x) cu ajutorul formulei Taylor cu rest Lagrange de ordin 3) :

(1)

cu

(2)

cu

Adunăm cele două relații (1) și (2) :

și not=nd : vom obține :

(3)

Din relația (3) (3’)

Considerăm o rețea de forma :

Vom aplica formula (3’) pentru următoarele date : și vom obți-ne :

Vom aplica formula (3’) pentru următoarele date : și vom obține :

Înlocuind ultimele două relații în ecuația Laplace vom obține ecuația :

Considerăm x # y # h și obținem ecuațiile simplificate :

(5)

Sau pe scurt :

(5’)

Dacă considerăm placa de oțel sub formă de rețea (figură) atunci conform formulelor (5’) avem :

– din rețeaua : vom obține ecuația :

– din rețeaua : vom obține ecuația :

– din rețeaua : vom obține ecuația :

– din rețeaua : vom obține ecuația :

– din rețeaua : vom obține ecuația :

– din rețeaua : vom obține ecuația :

Deci, distribuția temperaturii în interiorul plăcii este dată de :

Exemplu 2.

Fie problema : .

Unde .

Să se studieze distribuția temperaturii în interiorul domeniului (a,0)x(0,b).

Să rezolvăm această problemă folosind aceeași tehnică ca în exemplu precedent.

Să prezentăm grafic datele problemei :

Ecuația Laplace în două dimensiuni :

Din

Pentru problema noastră avem – condițiile Dirichlet :

, unde :

Proced=nd analog ca la exemplu anterior avem :

Să considerăm o rețea de forma :

Vom aplica formula de mai sus pentru următoarele date : și vom obți-ne :

Vom aplica aceeași formulă de mai sus pentru următoarele date : și vom obține :

Înlocuind ultimele două relații în ecuația Laplace vom obține ecuația :

Considerăm x # y # h și obținem ecuațiile simplificate :

sau pe scurt :

Proced=nd analog ca în exemplu precedent vom obține ecuațiile care dau distribuția temperaturii în interiorul domeniului considerat (a,0)x(0,b) :

Capitolul al VI lea . PROGRAME TURBO PASCAL

În acest capitol vom prezenta două programe TurboPascal, programe pentru problemele dezvoltate în capitolul VI, programe referitoare la metoda Jacobi și metoda Gauss – Seidel.

Program mJacobi;

Uses crt ;

const c # 10 ;

type vector # array ș 1..c ț of real ; { iter # numărul de iterații}

matrice # array ș1..c,1..cț of real ; {itmax # numărul maxim de

var a : matrice ; iterații}

b,x,y : vector ; {eps # precizia}

n, i, j : 1 .. c ; {norma # distanța dintre

eps, temp, norma, s : real ; două iterații consecutive}

iter,itmax, flag,k : integer ;

begin

write (‘Introduceți dimensiunea sistemului, n <# 10, n # ‘ );

readln (n) ;

writeln ;

writeln (‘Introduceți elementele matricei asociate sistemului și elementele matricei coloană ’ ) ;

for i : # 1 to n do

begin

writeln (‘Ecuația ‘, i ) ;

for j : # 1 to n do

begin

write ( ‘ ‘ , j : 2, ‘ : ‘ ) ;

readln ( a ș i , j ț ) ;

end ;

write (‘ b ‘, i , ‘ # ‘ );

readln ( b ș i ț ) ;

end ;

writeln ;

write (‘Dați precizia # ‘ ) ;

readln (eps) ; writeln ;

write (‘Dați numărul de iterații # ‘) ;

readln ( itmax) ; writeln ;

{ Pas0 : Inițializări }

flag : # 0 ;

for i : # 1 to n do

x ș i ț : # 0 ;

iter : # 1 ;

repeat

{Pas1 : Se calculează noua iterată }

for i : # 1 to n do

begin

s : # 0 ;

for j : # 1 to n do

if j < > i then s : = s + a ș i , j ț x ș j ț ;

y ș i ț : = ( b ș i ț – s ) / a ș i , i ț ;

end ;

readln ;

{ Pas2 : Se calculează distanța dintre iteratele consecutive }

norma : # 0 ;

for i : # 1 to n do

begin

temp : # abs (y ș i ț – x ș i ț ) ;

if norma < temp then norma : # temp ;

end;

{ Pas3 : Testele de oprire }

if norma < eps then flag : # 1

else

begin

iter : # iter ^1 ;

if iter < # itmax then x : # y

else flag : # 99 ;

end ;

until flag < >0 ;

{ Afișăm rezultatul }

if flag # 1 then

begin

writeln (‘ Soluția sistemului este :‘ ) ;

for k : # 1 to 18 do

write ( ‘ # ‘ ) ;

writeln ;

writeln ( ‘ S-au efectuat ‘, iter : 4, ‘ iterații ’ );

writeln ;

for i : # 1 n do

writeln ( ‘ ‘, i : 2, ‘ : ‘, y ș i ț ) ;

end;

else writeln ( ‘ Depășirea numărului de iterații ‘, itmax : 4) ;

end.

Program metodaGS ;

const c # 10 ;

type vector # array ș 1..c ț of real ; { iter # numărul de iterații}

matrice # array ș1..c,1..cț of real ; {itmax # numărul maxim de

var a : matrice ; iterații}

b,x,y : vector ; {eps # precizia}

n, i, j : 1 .. c ; {norma # distanța dintre

eps, temp, norma, s : real ; două iterații consecutive}

iter,itmax, flag,k : integer ;

begin

write (‘Introduceți dimensiunea sistemului, n <# 10, n # ‘ );

readln (n) ;

writeln ;

writeln (‘Introduceți elementele matricei asociate sistemului A și elementele matricei coloană B’ ) ;

for i : # 1 to n do

begin

writeln (‘Ecuația ‘, i ) ;

for j : # 1 to n do

begin

write ( ‘ ‘ , j : 2, ‘ : ‘ ) ;

readln ( a ș i , j ț ) ;

end ;

write (‘ b ‘, i , ‘ # ‘) ;

readln ( b ș i ț );

end ;

writeln ;

write (‘Dați precizia # ‘ ) ;

readln (eps) ; writeln ;

write (‘Dați numărul de iterații # ‘) ;

readln ( itmax) ; writeln ;

{ Pas0 : Inițializări }

flag : # 0 ;

for i : # 1 to n do

x ș i ț : # 0 ;

iter : # 1 ;

repeat

{Pas1 : Se calculează noua iterată }

for i : # 1 to n do

begin

s : # 0 ;

for j : # 1 to i-1 do

s : = s + a ș i , j ț y ș j ț ;

for j : # 1^1 to n do

s : = s + a ș i , j ț x ș j ț ;

y ș i ț : = ( b ș i ț – s ) / a ș i , i ț ;

end ;

readln ;

{ Pas2 : Se calculează distanța dintre iteratele consecutive }

norma : # 0 ;

for i : # 1 to n do

begin

temp : # abs (y ș i ț – x ș i ț ) ;

if norma < temp then norma : # temp ;

end;

{ Pas3 : Testele de oprire }

if norma < eps then flag : # 1

else

begin

iter : # iter ^1 ;

if iter < # itmax then x : # y

else flag : # 99 ;

end ;

until flag < >0 ; { Afișăm rezultatul }

if flag # 1 then

begin

writeln (‘ Soluția sistemului este‘ ) ;

for k : # 1 to 18 do

write ( ‘ # ‘ ) ;

writeln ;

writeln ( ‘ S-au efectuat ‘, iter : 4, ‘ iterații ’ );

writeln ;

for i : # 1 n do

writeln ( ‘ ‘, i : 2, ‘ : ‘, y ș i ț ) ;

end;

else writeln ( ‘ Depășirea numărului de iterații ‘, itmax : 4) ;

end.

Bibliografie

1. A.G.Kurosh – “Course of Higher Algebra”, 1972 (chapter 8)

2. Bakhvalov, N.S. – “Méthodes numériques”, Edition Mir (de

Moscou), 1976.

3. Bucur, C.M. – “Metode numerice”, Editura Didactică și

Pedagogică, București, 1986.

4. Bucur G. – “Calcul numeric”, Editura Didactică și Pedagogică,

Poppeea, C. București, 1983.

Simion, G.

5. C. Ilioi – “Analiză numerică” (curs universitar), vol. I,

Editura Univ. “Al. I. Cuza”, Iași, 1990.

6. Constantin Chiruță – “Analiza numerică și programare”, Editura Univ.

“Al. I. Cuza”, Iași, 1999.

7. Dan Larionescu – “Metode numerice”, Editura Tehnică,

București, 1989.

8. Demidovich, – “Computational mathematics”, R.W.S. – Kent –

BP & Maron, J. A. Publishing Company, Boston,

1985.

9. Dodescu, Gh. – “Metode de calcul numeric”, Ed. Tehnică,

Toma, M. București, 1981.

10. I. Ichim, – “Metode de aproximare numerică”,

G. Marinescu Editura Academiei Rom=ne, București, 1986.

11. Gheorghe Coman – “Analiză numerică”, Editura Libris, Cluj, 1995.

12. Gh. Dodescu – “Metode numerice în algebră”, Editura Tehnică,

București, 1979.

13. Marciuk, Gourij Ivanovici – “Méthodes de cacul numérique”, Édition Mir,

Moscou, 1980.

14. Mihai Postolache – “Metode numerice”, Editura Sirius, București,

1994.

15. Mihu, Cerchez – “Metode numerice în algebra liniară”,

Editura Tehnică, București, 1977.

16. M. Toma, i. Odăgescu – “Metode numerice și subrutine”, Editura Tehnică,

București, 1980.

17. Vilenkin, N. Ja – “Method of succesive approximations”,

Mir Publishers, Moscow, 1979.

18. V.I. Smirnov – “Course oh Higher Mathematics”, 1933.

Cuprins

CAPITOLUL I. Introducere

CAPITOLUL II. pregătitoare

1. Vectori și matrici. Clase speciale de matricipag 4

2. Norme vectoriale și norme matriceale

3. {iruri și serii de puteri matriceale

CAPITOLUL III. Rezolvarea sistemelor algebrice de ecuații liniare

1. Noțiuni generale privind rezolvarea sistemelor liniare. Prezentare general

2. Metode directe. Prezentare general

3. Metode indirecte. Prezentare general

CAPITOLUL IV. Metode iterative de rezolvare a sistemelor algebrice liniare

Partea A.

1. Formule de iterare

2. Convergența șirurilor de iterare

3. Propagarea erorilor

4. Metode iterative staționare și nestaționare

5. Schema generală a unei metode iterative

de rezolvare a sistemelor liniare algebrice

Partea B. Exemple de metode iterative

1. Metoda Jacobi

a. Algoritmul metodei Jacobi

b. Condiții de convergență a metodei Jacobi

2. Metoda Gauss – Seidel

a. Algoritmul metodei Gauss – Seidel

b. Condiții de convergență a metodei Gauss – Seidel

3. Metode iterative pentru sisteme liniare

bazate pe minimizarea unei forme pătratice

a. Metoda gradientului (metoda celei mai rapide descreșteri

b. Metoda gradientului conjugat

CAPITOLUL V. Aplicații practice

CAPITOLUL VI. Programe TurboPascal

Bibliografie

1. A.G.Kurosh – “Course of Higher Algebra”, 1972 (chapter 8)

2. Bakhvalov, N.S. – “Méthodes numériques”, Edition Mir (de

Moscou), 1976.

3. Bucur, C.M. – “Metode numerice”, Editura Didactic\ [i

Pedagogic\, Bucure[ti, 1986.

4. Bucur G. – “Calcul numeric”, Editura Didactic\ [i Pedagogic\,

Poppeea, C. Bucure[ti, 1983.

Simion, G.

5. C. Ilioi – “Analiz\ numeric\” (curs universitar), vol. I,

Editura Univ. “Al. I. Cuza”, Ia[i, 1990.

6. Constantin Chiru]\ – “Analiza numeric\ [i programare”, Editura Univ.

“Al. I. Cuza”, Ia[i, 1999.

7. Dan Larionescu – “Metode numerice”, Editura Tehnic\,

Bucure[ti, 1989.

8. Demidovich, – “Computational mathematics”, R.W.S. – Kent –

BP & Maron, J. A. Publishing Company, Boston,

1985.

9. Dodescu, Gh. – “Metode de calcul numeric”, Ed. Tehnic\,

Toma, M. Bucure[ti, 1981.

10. I. Ichim, – “Metode de aproximare numeric\”,

G. Marinescu Editura Academiei Rom=ne, Bucure[ti, 1986.

11. Gheorghe Coman – “Analiz\ numeric\”, Editura Libris, Cluj, 1995.

12. Gh. Dodescu – “Metode numerice `n algebr\”, Editura Tehnic\,

Bucure[ti, 1979.

13. Marciuk, Gourij Ivanovici – “Méthodes de cacul numérique”, Édition Mir,

Moscou, 1980.

14. Mihai Postolache – “Metode numerice”, Editura Sirius, Bucure[ti,

1994.

15. Mihu, Cerchez – “Metode numerice `n algebra liniar\”,

Editura Tehnic\, Bucure[ti, 1977.

16. M. Toma, i. Od\gescu – “Metode numerice [i subrutine”, Editura Tehnic\,

Bucure[ti, 1980.

17. Vilenkin, N. Ja – “Method of succesive approximations”,

Mir Publishers, Moscow, 1979.

18. V.I. Smirnov – “Course oh Higher Mathematics”, 1933.