Metode Numerice de Rezolvare a Sistemelor de Ecautii Liniare

Această lucrare cuprinde câteva metode numerice de rezolvare a sistemelor de ecații liniare cum ar fi:

-Metode directe( Metoda lui Gauss, Metoda Radăcinii pătrate,

Metoda Khaletski, Metoda lui Ritz) -Metode iterative( Metoda lui Jacobi, Metoda Gauss-Seidel, Metoda Relaxării,

Metoda Relaxării successive)

-Pseudoinversa unei aplicații liniare

Rezolvarea sistemelor de ecuații algebrice liniare este una din cele mai frecvente aplicații care apar în calculul numeric.

Metodele directe permit obținerea soluției exacte a sistemului considerat, făcând abstracție de erorile de rotunjire, folosind un număr finit de operații. Eroarea care se introduce prin rotunjirea calculelor devine destul de largă pentru sisteme unde ordinul n este destul de mare.

Metodele iterative se caracterizează prin faptul că soluția sistemului considerat se obține ca o limită a unui șir de vectori ce reprezintă soluția pentru diversele iterații efectuate. În cadrul acestor metode se pune problema de alegere a acelei metode care e convenabilă din punct de vedere al vitezei de convergență. Aceste metode pornesc de la o valoare inițială pentru vectorul soluție și cu ajutorul unui algoritm de calcul iterativ se determină un șir de aproximații succesive pentru vectorul soluție, , care trebuie să conveargă către soluția exactă a sistemului. În cazul în care procedeul converge suficient de rapid procedeul de calcul a sfârșit foarte repede și se obține o aproximație foarte bună. Unul din avantajele principale ale metodelor iterative este faptul că eroarea de rotunjire poate fi eliminată.

Vom nota cu corpul , cu spațiul liniar al matricilor cu r linii, m coloane și elemente din , cu spațiul liniar al operatorilor liniari definiți pe cu valori în . Știm că aplicația

unde

(1.1.)

este un izomorfism de spații liniare între și .

Vom considera și deci relația (1.1.) devine

(1.2.)

În particular, notăm cu algebra , cu algebra iar aplicația S definită prin (1.1.) este atunci un izomorfism de algebre între și .

Dacă , cu m<n, vom scrie .

Fie și pentru notăm

,

Folosind notația

și cu notațiile corespunzătoare pentru avem:

Definiție. Matricea este subdiagonală dacă pentru , iar în cazul în care pentru vom spune că este supradiagonală. Dacă pentru , vom spune că matricea este diagonală.

Definiție. O matrice este nesingulară dacă determinantul ei este nenul și deci matricea este inversabilă în .

Vom presupune cunoscut faptul că pe orice două norme sunt echivalente. Vom nota cu urmatoarele norme uzuale pe :

Dacă p este o normă pe iar vom nota norma lui A în spațiul adică

(1.3.)

Pentru normele uzuale vom nota corespunzător

. Dacă atunci

Dacă este un spațiu Hilbert și este spațiul operatorilor lininari și continui definiți pe și cu valori în atunci pentru vom nota cu adjunctul lui A adică operatorul din pentru care . Dacă se spune că acest operator este autoadjunct.

Pe se consideră structura de spațiu Hilbert generată de produsul scalar .

Teoremă. Adjunctul unui operator are următoarele proprietăți:

dacă există, atunci și există și

Definiție. Un operator autoadjunct spunem că este pozitiv dacă oricare ar fi .

Definiție. Pentru ca un operator liniar pe spațiul Hilbert ( cu valori în același spațiu) să fie continuu, este necesar și suficient să existe un număr așa ca

oricare ar fi .

Definiție. Pentru avem , deci și vom spune atunci că este simetrică.

În cazul general, în care este un spațiu Banach, iar , spectrul său, notat ; aceasta fiind complementara mulțimii numerelor pentru care . Dacă este un spațiu Banach, spectrul lui A este o mulțime compactă, nevidă. Numărul se numește număr propriu pentru A dacă există astfel încât . În acest caz elementul se numește vector propriu corespunzător numărului . Orice număr propriu aparține spectrului. Dacă spectrul este format numai din numere proprii. Dacă este număr propriu, subspațiul se numește subspațiu propriu corespunzător lui . este un subspațiu invariant pentru și îl vom nota . Dacă este un spațiu Hilbert, iar atunci spectrul este inclus în . Dacă vom spune despre că este unitar.

Pentru , matricea diagonală

o vom nota .

Teorema 1. Matricea este simetrică dacă și numai dacă există o matrice unitară și astfel încât .

Propoziția 1. O matrice este echivalentă cu o matrice diagonală dacă și numai dacă are vectori proprii liniar independenți.

Propoziția 2. Dacă are numere proprii distincte, atunci este echivalentă cu o matrice diagonală.

Teorema 2.( forma canonică Jordan). Fie ( . Atunci este echivalentă cu o matrice de forma

unde este o matrice , sau

unde sunt numere proprii pentru .

Matricile se numesc blocuri Jordan.

Deoarece unele din afirmațiile care urmează ar putea fi inconsistente datorită faptului că spectrul unor operatori din este vid, vom avea mereu în vedere prelungirea naturală a unor asemenea operatori la și vom lucra cu spectrul acestei extensii, în definitiv cu .

Numim rază spectrală a operatorului numărul

Vom considera în continuare pe algebra o normă de algebră , deci o normă pentru care . Reamintim că dacă este o norma pe atunci

este o norma de algebră pe .

Avem

(1.4.)

Propoziția 3. Dacă atunci

Demonstrație. Matricea fiind simetrică, fie și , , . Atunci

deci

și deci

(1.5.)

Fie apoi, ca în teorema 1, matricea unitară și așa ca , unde . Avem

.

Cum constituie spectrul lui deducem

și deci

Din inegalitatea precedentă și din (1.5.) propoziția este demonstrată.

Propoziția 4. Dacă , atunci .

Demonstrație. Fără o notație specială, vom avea mereu în vedere extensia lui A la . Fie așa ca . Din rezultă și atunci , de unde rezultă . Dacă , din inegalitatea precedentă rezultă că . Reciproc, dacă , adică dacă , există N și așa ca . Atunci , deci , ceea ce încheie demonstrația.

Definiție. Se spune că vectorul este ortogonal pe vectorul și se notează dacă .

Observație. Dacă și rezultă , .

Dacă avem ceea ce este echivalent cu 0.

Dacă .

Observație. Dacă avem implică . Deci putem scrie că , unde este subspațiu vectorial( este complementul ortogonal al lui ).

Dacă .

Fie , un subspațiu vectorial cu proprietatea ( astfel încât . Descompunerea precedentă este unică).

Avem .

Observație. Fie subspatiu, rezultă .

Corespondența , este un operator liniar. Operatorul definit de egalitatea ( respectiv ) se numește operator de proiecție sau proiector.

Definiție. Fie . Proiectorul este un operator autoadjunct pozitiv, mai precis un operator autoadjunct care se bucură de următoarele proprietăți:

0

Definiție. Elementele ale spațiului liniar se zic că sunt liniar independente dacă există n numere dintre care cel puțin unul nenul , așa ca

0 (1.6.)

Dacă elementele nu sunt liniar dependente, atunci ele se zic liniar independente.

Cu alte cuvinte, elementele sunt liniar independente dacă din (1.6.) rezultă 0.

Se verifică imediat că:

Pentru ca n elemente să fie liniar dependente, este necesar și suficient ca cel puțin unul dintre ele să fie o combinație liniară de celelalte.

Este evident că dacă printre elementele se află elementul nul, atunci elementele sunt liniar dependente.

Definiție. O mulțime (nevidă) spunem că este liniar independentă dacă orice n elemente (distincte) ale lui sunt liniar independente.

Din această definiție rezultă că mulțimea nevidă este liniar independentă dacă și numai dacă nici un element al său nu este o combinație liniară de alte elemente ale sale.

O mulțime liniar independentă nu conține elementul nul.

Definiție. Se numește bază algebrică o mulțime liniar independentă maximală.

Cu alte cuvinte, mulțimea liniar independentă este o bază dacă pentru orice mulțime liniar independentă , din , rezultă .

Definiție. O bază se numește ortonormată dacă vectorii bazei sunt unitari și ortogonali doi câte doi.

Teorema( Sylvester). Matricea este pozitiv definită dacă toți determinanții principali ai matricei sunt strict pozitivi.

Aceste sisteme au forma

(2.1.)

unde .

Ne vom situa în ipoteza că sistemul precedent are soluție unică, deci că det A≠0.

Metoda lui Gauss

Diversele variante ale acestei metode realizează transformarea sistemului (2.1.) într-un sistem echivalent, având matrice triunghiulară. Transformarea se face prin operații liniare, prin care se obține eliminarea succesivă a necunoscutelor. În procedeul pe care îl vom prezenta sunt necesare anumite permutări de linii și coloane, operații numite de pivotare.

Dacă fixăm prima coloană în (2.1.), putem face o permutare a două linii, astfel încât, elementul de modul maxim din prima coloană să fie în prima linie. În acest caz, elementul considerat este numit pivot, iar prima linie este numită linie pivot. Sistemul devine

(2.2.)

unde . Elementul este diferit de zero căci det A0. Împărțind prima linie a sistemului (2.2.) cu ea se transformă în:

(2.3.)

unde

Pentru , înmulțim ecuația (2.3.) cu și o scădem din linia i a sistemului (2.2.). Obținem astfel sistemul:

(2.4.)

unde

Determinantul matricii sistemului format cu ultimele n-1 linii ale sistemului (2.4.) este nenul, căci în caz contrar det A=0.

Procedeul aplicat sistemului (2.1.) se poate aplica sistemului redus:

sistem care conține numai necunoscutele și care, cu notațiile,

se scrie

După p-1 pași ajungem la sistemul

nali doi câte doi.

Teorema( Sylvester). Matricea este pozitiv definită dacă toți determinanții principali ai matricei sunt strict pozitivi.

Aceste sisteme au forma

(2.1.)

unde .

Ne vom situa în ipoteza că sistemul precedent are soluție unică, deci că det A≠0.

Metoda lui Gauss

Diversele variante ale acestei metode realizează transformarea sistemului (2.1.) într-un sistem echivalent, având matrice triunghiulară. Transformarea se face prin operații liniare, prin care se obține eliminarea succesivă a necunoscutelor. În procedeul pe care îl vom prezenta sunt necesare anumite permutări de linii și coloane, operații numite de pivotare.

Dacă fixăm prima coloană în (2.1.), putem face o permutare a două linii, astfel încât, elementul de modul maxim din prima coloană să fie în prima linie. În acest caz, elementul considerat este numit pivot, iar prima linie este numită linie pivot. Sistemul devine

(2.2.)

unde . Elementul este diferit de zero căci det A0. Împărțind prima linie a sistemului (2.2.) cu ea se transformă în:

(2.3.)

unde

Pentru , înmulțim ecuația (2.3.) cu și o scădem din linia i a sistemului (2.2.). Obținem astfel sistemul:

(2.4.)

unde

Determinantul matricii sistemului format cu ultimele n-1 linii ale sistemului (2.4.) este nenul, căci în caz contrar det A=0.

Procedeul aplicat sistemului (2.1.) se poate aplica sistemului redus:

sistem care conține numai necunoscutele și care, cu notațiile,

se scrie

După p-1 pași ajungem la sistemul

Sistemului

scris prescurtat îi aplicăm apoi tratamentul din prima etapă.

După n pași se ajunge la sistemul

(2.5.)

Sistemul (2.5.) este atunci echivalent cu (2.1.) și din , prin înlocuire, obținem succesiv valorile necunoscutelor :

Observația 1. Deoarece în procedeul prezentat am permutat numai linii vom spune că rezolvarea s-a facut prin pivotare parțială. Pentru eliminarea la fiecare pas a necunoscutei cu coeficientul de modul maxim în matricea sunt necesare de obicei permutări de linii si coloane și spunem în acest caz că rezolvarea s-a facut prin pivotare completă.

Obsrevația 2. Dacă în procedeul prezentat alegerea lui (respectiv a lui se face impunând doar condiția , ajungem de asemenea la un sistem de forma (2.5.). Se poate însă întâmpla ca datorită unor asemenea coeficienți nenuli, dar mici în modul, să ajungem, prin contribuția erorilor, la o soluție cu eroare mare sau chiar la un sistem neechivalent cu sistemul inițial, spre exemplu la un sistem nedeterminat sau incompatibil.

Aplicație. Să se rezolve cu ajutorul metodei lui Gauss sistemul:

Rezolvare:

Matricea A și vectorul termenilor liberi b este:

,

Calculând, constatăm că det A0.

Deoarece este elementul în modul maxim din prima coloana a matricei A, nu este necesară o permutare a ecațiilor sistemului.

Regulile algoritmului indicat de metoda lui Gauss sunt:

-împărțim prima ecuație din sistem cu :

(*)

-ecuația (*) înmulțită cu o scădem din ecuația a doua a sistemului:

-ecuația (*) înmulțită cu o scădem din ecuația a treia a sistemului:

-ecuația (*) înmulțită cu o scădem din ecuația a treia a sistemului:

Se obține sistemul echivalent cu matricea coeficienților și vectorul termenilor liberi :

și

Pentru sistemul cu matricea coeficienților

și vectorul termenilor liberi .

observăm că elementul este în modul maxim pe coloana lui și deci nu este necesară o permutare a ecuațiilor din acest sistem.

Procedând ca mai sus,obținem sistemul cu matricea coeficienților și vectorul termenilor liberi .

, ,

adăugând prima ecuație de la pasul anterior care nu se mai modifică.

În continuare, observăm că nu este elementul în modul maxim pe coloana lui și de aceea vom schimba ultimele două ecuații între ele

,

iar apoi aplicăm rețeta cunoscută și obținem:

, .

Analog se obține

, .

Rezolvând sistemul , adică

începând cu ultima ecuație găsim

Se obține soluția:

.

Metoda rădăcinii pătrate

Cu notațiile și definițiile din capitolul precedent matricea simetrică este strict pozitivă( sau pozitiv definită) dacă deci dacă

Teoremă. Matricea A este simetrică și strict pozitivă dacă și numai dacă există matricea nesingulară subdiagonală astfel încât .

Demonstrație. Dacă , unde B este nesingulară, atunci este evident simetrică. Avem pentru, căci este de asemenea nesingulară.

Afirmația reciprocă o vom demonstra prin inducție după dimensiunea spațiului.

Pentru n=1, , , din ipoteză, , adică , de unde . Vom lua atunci și deci .

Să presupunem afirmația adevărată pentru n-1 și fie o matrice simetrică și strict pozitivă. Scriem

unde

.

Deoarece este simetrică, rezultă că este simetrică. Tot de aici rezultă că . Vom folosi în continuare notațiile:

, .

Fie și . Avem

și deci este strict pozitivă. Conform ipotezei, o matrice nesingulară, subdiagonală, astfel încât . Vom arăta că descompunerea căutată se poate realiza cu o matrice B de forma

unde

0= și .

Trebuie deci să avem

și deci

Este deci suficient să luăm

de unde

Rămâne să arătăm că

și atunci

Fie pentru aceasta

unde . Avem atunci

ceea ce încheie demonstrația.

Din demonstrație reținem că putem alege .

Teorema precedentă constituie fundamentul teoretic pentru metoda rădăcinii pătrate, folosită pentru rezolvarea unor sisteme de forma , unde este o matrice simetrică și strict pozitivă. În această metodă, după descompunerea cu nesingulară, subdiagonală, sistemul este echivalent cu

(2.6.)

sistem a cărui rezolvare este simplă deoarece este o matrice nesingulară supradiagonală.

Vom prezenta în continuare formulele prin care se ajunge la rezolvarea sistemului (2.6.) precum și cele ce realizează rezolvarea acestui ultim sistem.

Dacă , deoarece pentru avem rezultă că

(2.7.)

de unde obținem succesiv

Pentru trecerea de la coloana la coloana scriem pentru început

Rezultă de aici

Continuăm cu determinarea celorlalte elemente de pe coloana .

Tot din (2.6.), pentru :

și deci

Pentru rezolvarea sistemului (2.6.), notând avem de unde deducem:

Analog se rezolvă apoi sistemul :

Aplicație. Să se rezolve cu ajutorul metodei radăcinii pătrate sistemul

Rezolvare. Matricea A și vectorul termenilor liberi sunt

și

Matricea A este simetrică. Pentru a verifica dacă A este pozitiv definită va trebui să găsim toți determinanții principali ai matricei A strict pozitivi:

.

Deci A este pozitiv definită.

Există deci o descompunere a matricei A în produsul , unde este o matrice subdiagonală.

Conform formulelor de mai sus avem

elemente ce reprezintă prima coloană din matricea .

În continuare obținem elementele coloanei a doua din matricea :

Analog obținem

Rezolvăm sistemul , unde

,

și obținem

Cu vectorul

se va rezolva sistemul , unde

astfel:

Vectorul soluție este:

.

Metoda lui Khaletski

Pentru comoditatea raționamentului scrierea sistemului de ecuații are forma matriceală

(2.1.)

unde o matrice pătratică și , sunt vectorii coloană. Matricea are forma unui produs dintre o matrice triunghiulară subdiagonală și o matrice triunghiulară supradiagonală ale cărei elemente de pe diagonală sunt 1. Cu alte cuvinte

(2.8.)

unde

și .

Elementele și se obțin astfel

(2.9.)

și

(2.10.)

Vectorul căutat este calculat din sistemul de ecuații

(2.11.)

Matricile și fiind triunghiulare, sistemele (2.11.) se rezolvă ușor

(2.12.)

și

(2.13.)

Formulele (2.12.) au avantajul de a calcula numerele în același timp cu coeficienții . Această metodă este cunoscută sub numele de metoda lui Khaletski. Metoda face apel la calculul unor sume.

Remarcăm că dacă este simetrică, rezultă

.

Metoda Khaletski este cunoscută pentru ușurința aplicării calculului electronioc, în cazul folosirii computerelor deoarece elimină înregistrarea unor rezultate intermediare.

Factorizarea din metodă este posibilă dacă toți minorii principali sunt nenuli ( [5] ).

Aplicație. Folosind metoda Khaletski să se rezolve sistemul

Rezolvare.

Matricea sistemului și vectorul termenilor liberi sunt

și .

Calculăm elementele și astfel

Aplicăm aceeași rețetă și găsim

Obținem, deci

și

Rezolvând sistemul se obține

Din sistemul obținem succesiv valorile necunoscutelor începând cu ultima

Vectorul soluție este

.

Metoda lui Ritz de inversare a unei matrici

simetrice și pozitiv definite

Propoziția 1. Fie simetrică, pozitiv definită și

-ortogonali și nenuli. Fie ,

Amintim că și că sunt ortogonali dacă .

Atunci și .

Demonstrație. Fie ,

. Pentru avem , de unde rezultă că 0, . Deoarece sunt ortogonali și nuli formează bază. Din ultimele două relații rezultă 0 ceea ce este echivalent cu .

Propoziția 2( Algoritm pentru ). Fie

și astfel încâ 0. Pentru se alege astfel îcât 0 și . Fie . Atunci .

Demonstrație. Pentru vom alege astfel încât 0 și . Să presupunem că am construit

astfel încâtsă fie ortogonali și nenuli. Se poate construi și fie ca în enunț astfel încât 0 și , . Atunci .

Vom arăta că sunt ortogonali. Arătăm doar că pentru . Deci

. Rezultă că sunt ortogonali. Atunci sunt ortogonali, nenuli și deci .

Aplicație. Folosind metoda lui Ritz să se inverseze următoarea matrice

.

Rezolvare. Observăm că matricea este simetrică și ne rămâne de arătat că este pozitiv definită, adică toți determinanții principali ai matricei trebuie să fie strict pozitivi:

Deci pozitiv definită.

Se aleg vectorii coloană formați din coloanele lui :

, , , .

Fie . Fie astfel încât

0. Observăm că 0 deci algoritmul poate continua și vom calcula . Atunci , unde ,

=, ,

. Deci

.

Fie . Verificăm dacă algoritmul poate continua, adică dacă 0. Făcând calculele obținem

, și

deci 0. Luăm

și , unde

,

, ,

. Obținem

.

Fie , verificăm 0 și găsim

, ,

deci 0 rezultă, algoritmul poate continua și luăm

și , unde

,

, ,

. Deci

.

Fie . Verificăm condiția 0. Calculăm

și obținem ,

. Deci condiția este îndeplinită și luăm

iar , unde

,

, ,

. Obținem

și

.

O metodă iterativă convergentă de rezolvare a unui sistem de ecuații liniare se caracterizează prin construcția unui șir de vectori din , convergent către soluția a sistemului. Vom calcula și diferența numită eroare la pasul n și vom evalua această eroare cu una din normele uzua-le pe .

Metoda lui Jacobi

Fie sistemul

(3.1.)

unde , I fiind operatorul identitate, iar .

Fie și

(3.2.)

Referindu-ne la acest șir vom spune că este șirul Jacobi asociat sistemului (3.2.).

Teorema 1. Sistemul (3.1.) are pentru orice soluția unică z și pentru orice șirul (3.2.) converge către z dacă și numai dacă

Demonstrație. Dacă sistemul are pentru orice b soluție unică, atunci dacă și numai dacă x=0 și tot din ipoteză, pentru orice , șirul (3.2.) converge la zero.Pe de altă parte, pentru b=0, ecuația (3.1.) devine

Avem și deci adică (se folosește faptul că convergența punctuală este echivalentă cu convergența în orice normă pe acest spațiu).

Reciproc, să presupunem că Dacă atunci de unde rezultă . Conform ipotezei avem atunci y=0, adică este injectiv de unde rezultă că este și surjectiv, deci bijectiv, ceea ce arată că sistemul (3.1.) are pentru orice b soluția unică z adică . Observând că , rezultă și din ipoteză deducem atunci

Evaluarea erorii la metoda lui Jacobi

În evaluările care urmează vom considera o normă pe , iar

pe norma indusă, adică dacă p este o normă pe , iar , vom nota norma lui A în spațiul adică.Vom folosi notațiile din teorema precedentă și ne vom situa în ipoteza

Avem

și deci

, (3.3.)

Rezultă evaluarea

, (3.4.)

Inegalitatea precedentă constituie o formulă de evaluare aposteriori a erorii la metoda Jacobi.

Din (3.3.) deducem apoi

de unde rezultă

, (3.5.)

ceea ce constitue o formulă de evaluare apriori a erorii la metoda lui Jacobi.

În practică metoda se aplică în ipoteza , de unde rezultă , . Din (3.3.) deducem atunci

(3.6.)

Prin diverse metode, dintr-un sistem de forma , se poate obține un sistem echivalent de forma (3.1.), transformarea fiind adecvată dacă se ajunge la un sistem pentru care metoda Jacobi este convergentă. Fie deci un sistem de forma unde .

Observația 1. Să presupunem că

(3.7.)

Fie matricea diagonală . Rezultă că există astfel încât sistemul este echivalent cu

sistemul , ceea ce se mai poate scrie sau unde și sau, scris pe componente, pentru si iar .

Avem

de unde rezultă că se poate aplica teorema 1 și formula (3.6.) (în raport cu ).

Observația 1.1. Soluția sistemului este aproximată cu . Șirul iteratelor, scris pe componente, este

.

Observația 2. Să presupunem că

, (3.8.)

atunci și deci există . Fie Sistemul este echivalent cu , unde = sau, scris pe componente pentru si . Deoarece

sistemului i se poate aplica metoda Jacobi și formula (3.6.), în raport cu . Fie z soluția sistemului . Dacă w este soluția unică a sistemului rezultă . Sistemul are soluția . Notând cu șirul Jacobi asociat sistemului , atunci șirul , converge către z. Avem

În situațiile prezentate în observațiile precedente se spune că matricea A este diagonal dominantă.

Metoda Jacobi conduce deci la un proces iterativ convergent în cazul sistemelor de ecuații liniare cu matrice diagonal dominantă.

Observatia 3. Fiind dat sistemul , cu det A0, dacă elementele de pe diagonala matricei A nu sunt diferite de zero, atunci se pot face combinații liniare de ecuațiile sistemului, în așa fel încât să obținem un sistem echivalent cu cel dat, care să aibă matricea cu elementele de pe diagonală diferite de zero.

Aplicație. Folosind metoda Jacobi, să se aproximeze soluția pentru următoarele sisteme, de forma , astfel încât eroarea să fie mai mică decât :

a) b)

Rezolvare:

a) Matricea A și vectorul termenilor liberi a este:

Elementele matricei A satisfac condițiile (3.7) și anume

deci, matricea A este dominant diagonală pe linii.

Transformăm sistemul inițial într-un sistem echivalent de forma cu și unde

, ,

și .

Obținem

, .

Pentru acest sistem este îndeplinită condiția

. Avem

Deci metoda converge iar q=0,806452.

Șirul iteratelor este convergent la soluția sistemului , oricare ar fi .

Scris pe componente acest șir arată astfel

Pentru x(0) ={ 0, 0, 0 } rezultatele iteratelor sunt:

x(0) = { 0 0 0 }

x(1) = { 3.493548 3.68 4.071429 }

x(2) = { 0.399539 0.769585 2.801536 }

x(3) = { 2.217447 2.879969 3.884683 }

x(4) = { 0.846891 1.572595 3.200611 }

x(5) = { 1.70016 2.531743 3.682574 }

x(6) = { 1.080584 1.923389 3.365231 }

x(7) = { 1.477318 2.358603 3.585172 }

x(8) = { 1.195782 2.076575 3.4389 }

x(9) = { 1.379431 2.274751 3.539507 }

x(10) = { 1.251086 2.14444 3.472189 }

x(11) = { 1.335855 2.234911 3.518261 }

x(12) = { 1.277217 2.174835 3.487307 }

x(13) = { 1.316271 2.216208 3.50842 }

x(14) = { 1.289441 2.188553 3.494196 }

x(15) = { 1.307411 2.207496 3.503877 }

x(16) = { 1.295122 2.194778 3.497343 }

x(17) = { 1.303384 2.203458 3.501783 }

x(18) = { 1.297752 2.197613 3.498783 }

x(19) = { 1.301548 2.201593 3.500819 }

x(20) = { 1.298965 2.198908 3.499442 }

x(21) = { 1.300709 2.200733 3.500376 }

x(22) = { 1.299524 2.1995 3.499744 }

x(23) = { 1.300325 2.200337 3.500173 }

x(24) = { 1.299781 2.199771 3.499883 }

x(25) = { 1.300149 2.200155 3.500079 }

x(26) = { 1.2999 2.199895 3.499946 }

x(27) = { 1.300068 2.200071 3.500036 }

x(28) = { 1.299954 2.199952 3.499975 }

x(29) = { 1.300031 2.200033 3.500017 }

x(30) = { 1.299979 2.199978 3.499989 }

x(31) = { 1.300014 2.200015 3.500008 }

x(32) = { 1.29999 2.19999 3.499995 }

x(33) = { 1.300007 2.200007 3.500004 }

Pentru evaluarea erorii vom folosi formula (3.6) ( în raport cu ) și obținem:

b) Matricea A și vectorul termenilor liberi a este:

Elementele matricei A satisfac condițiile (3.8.) și anume

deci, matricea A este dominant diagonală pe coloane.

Fie

.

Făcând schimbarea sistemul dat devine , cu

. Luând se obține sistemul echivalent

cu

Pentru acest sistem este îndeplinită condiția . Avem Deci metoda converge iar .

Șirul iteratelor este convergent la soluția sistemului , oricare ar fi .

Scris pe componente, acest șir arată astfel:

Pentru y ( 0 ) = { 0, 0, 0, 0 } rezultatele iteratelor sunt:

y(0) = { 0 0 0 0 }

y(1 ) = { 1 1 1 1 }

y(2 ) = { 0.5 0.6 0.6 0.7 }

y(3) = { 0.69 0.77 0.77 0.83 }

y(4) = { 0.609 0.702 0.702 0.777 }

y(5) = { 0.6415 0.7303 0.7303 0.7987 }

y(6) = { 0.62801 0.7188 0.7188 0.78979 }

y(7) = { 0.633501 0.723539 0.723539 0.793439 }

y(8) = { 0.63124 0.721602 0.721602 0.791942 }

y(9) = { 0.632165 0.722397 0.722397 0.792556 }

y(10) = { 0.631785 0.722072 0.722072 0.792304 }

y(11) = { 0.631941 0.722205 0.722205 0.792407 }

y(12) = { 0.631877 0.722151 0.722151 0.792365 }

y(13) = { 0.631903 0.722173 0.722173 0.792382 }

Pentru evaluarea erorii vom folosi formula (3.6.) (în raport cu ) și obtinem:

Soluția sistemului dat se obține ținând seama de substituția făcută. Deci soluția finală este:

x={ 0.0631903 0.0722173 0.0722173 0.0792382 }

Metoda Gauss-Seidel

În această metodă, care este o rafinare a metodei lui Jacobi, pentru rezolvarea unui sistem de forma (3.1.), în construcția iteratei se folosesc componentele lui anterior calculate.

Pentru se notează

și

unde R=B-L. Sistemul (3.1.), adică, este echivalent cu sistemul

. Calculând det(I-L) observăm că este diferit de zero, de unde rezultă că există astfel încât este echivalent cu

. Aplicând obținem . Deci sistemul (3.1.) este identic cu sistemul

(3.9.)

unde

C=, c=.

Metoda Gauss-Seidel constă în aproximarea soluției sistemului (3.1.) cu ajutorul șirului Jacobi asociat sistemului (3.9.). Pentru considerăm deci adică .Aplicăm (I-L) și obținem

, ceea ce se mai poate scrie:

(3.10.)

Componentele șirului precedent sunt date prin:

(3.11.)

Vom nota

, și .

În cazul în care șirul iteratelor este:

.

Teorema 2. Dacă q<1 atunci sistemul (3.1.) are pentru orice soluție unică z, pentru orice șirul construit prin (3.10.) converge către această soluție și au loc evaluările:

(3.12.)

Demonstrație. Vom arăta că și vom aplica apoi teorema 1. Pentru , avem:

, (3.13.)

Într-adevar, care este identic cu sau , ceea ce pe componente se scrie:

,

Pentru avem , de unde deducem

Să presupunem inegalitațile (3.13.) adevărate pentru unde . Atunci

Din (3.13.) obținem sau și deci . Deoarece se poate aplica teorema 1.Rezultă că sistmul (3.1.) are pentru orice soluție unică z, șirul (3.10) converge la această soluție și au loc evaluările (3.12.).

Vom nota cu șirul Jacobi asociat sistemului (3.1.) și cu șirul (3.10.), pe care îl vom numi șirul Gauss-Seidel corespunzător sistemului (3.1.).

Numărul q este cel din teorema 2.

Teorema 3. Dacă

și (3.14.)

atunci (3.1.) are pentru orice soluție unică z, șirurile și converg către această soluție și au loc:

Demonstrație. Vom arăta pentru început că . Prin ipoteză avem ceea ce este echivalent cu și vom demonstra prin inducție că Presupunând afirmația adevărată până la k evaluăm . Fie . Dacă atunci (conform cu (3.14.)). În caz contrar, există astfel încât . Conform ipotezei de inducție avem și deci din prima dintre formulele (3.14.) rezultă

Avem deci . Conform teoremei 2 sistemul (3.1.) are pentru orice soluție unică z și avem

Vom arăta în continuare că în ipotezele considerate converge și șirul Jacobi asociat sistemului (3.1.). Din prima dintre relațiile (3.14.) deducem că și de aici

N (3.15.)

Fie si , .

Avem

, , (3.16.)

Într-adevăr, dacă în (3.16.) considerăm n=1 va rezulta k=1. Conform ipotezei avem ceea ce înseamnă că relațiile (3.16.) sunt adevărate pentru n=1. Să presupunem (3.16.) adevărate pentru n<m. Pentru avem

Folosind ipoteza de inducție și (3.15.) obținem Pentru n=m, din (3.16.) deducem atunci

sau

și deci

(3.17.)

Pentru N, , notăm cu partea întreagă a lui și fie N, astfel încât . Folosind (3.15.) și (3.17.) avem

Reținem deci că

(3.18.)

Reamintim că de unde obținem

Folosind și (3.18.) obținem

(3.19.)

ceea ce încheie demonstrația.

Din (3.15.) rezultă că evaluarea aposteriori din (3.19.) are loc sub forma

Observație. Sistemul Ax=a este echivalent cu sistemul (3.1.), unde si b=. Atunci relațiile (3.14.) pentru sistemul (3.1.) devin:

și

sau

și

Aplicație. Folosind metoda Gauss-Seidel, să se rezolve cu o eroare mai mică decat sistemul

Rezolvare. Matricea A și vectorul termenilor liberi asociați sistemului este:

și .

Elementele matricei A satisfac condițiile și anume

Transformăm sistemul dat într-un sistem echivalent de forma , cu , unde

Obținem

și , pentru care sunt verificate condițiile ( 3.14. ), adică

și deci șirul Gauss-Seidel converge la soluția sistemului dat. Calculăm :

și obținem q=0,539216.

Deoarece q=0,539216<1, șirul iteratelor Gauss-Seidel,

,

converge la soluția sistemului dat, oricare ar fi .

Atunci avem

,

și șirul iteratelor Gauss-Seidel, scris pe componente este

Luând x(0) ={ 0, 0, 0 }, primele iterate calculate cu 3 zecimale exacte sunt:

x(1) = { 0.504902 1.559301 2.582254 }

x(2) = { 1.64657 2.316311 2.905724 }

x(3) = { 1.92725 2.461089 2.980467 }

x(4) = { 1.984718 2.491862 2.995901 }

x(5) = { 1.9968 2.498295 2.999142 }

x(6) = { 1.99933 2.499643 2.99982 }

x(7) = { 1.99986 2.499925 2.999962 }

Pentru evaluarea erorii folosim formula

Metode de relaxare

În metodele iterative studiate până acum în acest capitol am aproximat soluția unui sistem de ecuații liniare prin intermediul unui șir , convergent la soluție, șir care se construia formal prin unde la metoda Jacobi și la metoda Gauss-Seidel. Privind pe ca o “corecție” efectuată asupra lui pentru a obține pe , metodele de relaxare constau într-o “supracorectare” în care trecerea de la la se face prin unde iar are semnificația de mai sus. Cele precedente constituie o încercare de acreditare a ideii că procedeele de relaxare sunt naturale, existând însă și alte modalități de prezentare a apariției lor, spre exemplu ca rezultat al unor fenomene de optimizare.

Metoda relaxării simultane

Vom considera sistemul

(3.20.)

Presupunând că elementele de pe diagonala matricii sunt nenule, sistemul (3.20.) se aduce la forma echivalentă

(3.21.)

unde și , , . Pentru notăm , și sistemul (3.21.) este echivalent cu sistemul

(3.22.)

Șirul Jacobi asociat acestui ultim sistem, se scrie sub forma echivalentă

(3.23.)

Deoarece șirul Jacobi asociat sistemului (3.21.) este

rezultă că șirul (3.23.) realizează supracorectarea despre care spunem că este caracteristică unei metode de relaxare.

Observație. Șirul iteratelor, scris pe componente, pentru metoda relaxării simultane este:

.

Fie simetrică și pozitiv definită. Atunci de unde rezultă că este un produs scalar pe . Vom nota cu norma generată de acest produs scalar pe ( respectiv norma operatorială indusa pe ).

Observație. Observând că o matrice -hermitică considerăm . Atunci există 0, astfel încât , de unde rezultă . Deci . Știm că

. Rezultă deci că.

Rezultă deci că dacă este spectrul matricii atunci iar este spectrul matricii . Vom presu-pune că .

Teorema 4. Dacă este simetrică și pozitiv definită atunci, pentru orice șirul (3.23.) converge la soluția z a sistemului (3.20.) dacă și numai dacă . Avem

unde .

Demonstrație. Dacă , spre exemplu dacă avem . Atunci șirul nu converge la zero și din teorema 1 rezultă că șirul (3.23.) nu converge întotdeauna către soluția sistemului (3.20.).

Dacă apoi , atunci , deci

ceea ce arată că . Conform teoremei 1 șirul (3.23.) converge la soluția z a sistemului considerat.

Evaluarea erorii va rezulta din (3.6.) observând că .Într-adevar, deoarece este -simetrică, din propoziția 3, scrisă pentru avem

.

ceea ce încheie demonstrația.

Deoarece pentru orice norma operatorială avem

, rezultă că pentru a avea este suficient să luăm .

Teoremă( Optimizarea parametrului de relaxare). Dacă , atunci . Deci este valoarea minimă pe care o poate lua parametrul .

Aplicație. Folosind metoda Relaxării simultane să se rezolve cu o eroare de sistemul

Rezolvare.

Matricea A și vectorul termenilor liberi a este:

, .

A este simetrică și verificăm dacă este pozitiv definită, adică dacă determinanții principali ai matricei A sunt strict pozitivi:

.

Deci matricea A este pozitiv definită.

Alegem astfel încât , unde iar . Avem deci de unde rezultă că

Deoarece , adica vom considera =0,9( pentru acest sigma metoda converge mai rapid ).

Fie , , unde

și .

Obținem și .

Pentru acest sistem avem , deci

Șirul iteratelor este

Primele iterate calculate pentru x(0) ={ 0, 0, 0 } sunt:

x(1) = { 9.72 -0.54 -1.44 }

x(2) = { 10.7244 1.2852 -2.4102 }

x(3) = { 11.2407 1.73583 -2.76188 }

x(4) = { 11.4051 1.90548 -2.88408 }

x(5) = { 11.4631 1.96303 -2.92636 }

x(6) = { 11.483 1.98303 -2.94098 }

x(7) = { 11.4899 1.98994 -2.94604 }

x(8) = { 11.4923 1.99233 -2.94779 }

x(9) = { 11.4932 1.99315 -2.9484 }

Pentru evaluarea erorii vom folosi formula (3.6.)(în raport cu ) și obținem:

Metoda relaxării successive

Vom considera în continuare sistemul care, în ipoteza că elementele de pe diagonala matricii A sunt nenule, se scrie sub forma echivalentă . Notând

și șirul Gauss-Seidel asociat sistemului este

Supracorectarea specifică unei metode de relaxare se realizează considerând șirul

care se poate scrie

sau

(3.24.)

Notând

și

șirul precedent se construiește prin

(3.25.)

Șirul (3.25.) este șirul Jacobi asociat sistemului, sistem echivalent cu (3.20.). Metoda relaxării successive constă în aproximarea soluției sistemului (3.20.) cu ajutorul șirului (3.25.).

Lemă. Dacă matricea A are numai elemente nenule pe diagonală atunci .

Demonstrație. Fie sistemul rădăcinilor ecuației . Se știe că avem

(3.26.)

Pe de altă parte

Din (3.26.) rezultă că și deci

La fel ca în teorema precedentă, dacă matricea A este simetrică și pozitiv definită atunci aplicația este un produs scalar pe și notăm cu norma corespunzatoare pe ( respectiv norma operatorială indusă pe ).

Teorema 5. Fie A o matrice simetrică și pozitiv definită. Pentru orice iterație inițială , șirul (3.25.) converge la soluția z a sistemului (3.20) dacă și numai dacă . Avem atunci

unde .

Demonstrație. Dacă atunci din lema precedentă rezultă că . Conform propoziției 4 șirul nu converge la zero și atunci conform teoremei 1 șirul (3.25.) nu converge întotdeauna la soluția sistemului considerat.

Vom arăta că pentru matricea este pozitiv definită de unde va rezulta că . Fie

și

Atunci

Din rezultă deci că

și prin urmare

Deoarece , iar , matricea este pozitiv definită. Matricea fiind nesingulară rezultă atunci că , și deci , este pozitiv definită.Aceasta înseamnă că

ceea ce se poate scrie:

(3.27.)

Deoarece este compactă, iar este continuă, din (3.27.) rezultă că . Din teorema 1 rezultă atunci că șirul (3.25.) converge către soluția z a sistemului (3.20.).

Observația 1. Pentru , metoda se numește subrelaxare iar pentru suprarelaxare. Pentru se obține metoda Gauss-Seidel.

Observația 2. Pe componente, șirul (3.25.) este

Aplicație. Să se rezolve cu ajutorul metodei Relaxării succesive sistemul

Rezolvare.

Matricea și vectorul termenilor liberi este

și .

Observăm că matricea A este simetrică. Rămâne să verificăm dacă A este pozitiv definită, adică toți determinanții principali ai matricei să fie strict pozitivi:

Deci A este pozitiv definită.

Așa cum am văzut mai sus avem

, ,

,

Pentru șirul iteratelor scris pe componente este:

Pentru x ( 0 ) = { 0, 0, 0, 0 } primele iterate sunt

x ( 1 ) = { 4.74283 -0.427668 1.04296 0.989077 }

x ( 2 ) = { 4.33056 -0.972266 1.48426 1.16784 }

x ( 3 ) = { 4.06861 -1.2072 1.74225 1.22887 }

x ( 4 ) = { 3.91474 -1.33187 1.89225 1.2593 }

x ( 5 ) = { 3.82488 -1.40263 1.97935 1.27642 }

x ( 6 ) = { 3.77257 -1.44349 2.02993 1.2863 }

x ( 7 ) = { 3.74216 -1.46719 2.0593 1.29204 }

x ( 8 ) = { 3.72449 -1.48094 2.07636 1.29538 }

x ( 9 ) = { 3.71422 -1.48893 2.08627 1.29732 }

x ( 10 ) = { 3.70826 -1.49357 2.09203 1.29844 }

x ( 11 ) = { 3.7048 -1.49627 2.09537 1.29909 }

x ( 12 ) = { 3.70279 -1.49783 2.09731 1.29947 }

x ( 13 ) = { 3.70162 -1.49874 2.09844 1.29969 }

x ( 14 ) = { 3.70094 -1.49927 2.09909 1.29982 }

x ( 15 ) = { 3.70055 -1.49957 2.09947 1.2999 }

x ( 16 ) = { 3.70032 -1.49975 2.09969 1.29994 }

x ( 17 ) = { 3.70018 -1.49986 2.09982 1.29997 }

x ( 19 ) = { 3.70011 -1.49992 2.0999 1.29998 }

x (19 ) = { 3.70006 -1.49995 2.09994 1.29999 }

x ( 20 ) = { 3.70004 -1.49997 2.09997 1.29999 }

x ( 21 ) = { 3.70002 -1.49998 2.09998 1.3 }

x ( 22 ) = { 3.70001 -1.49999 2.09999 1.3 }

Pentru șirul iteratelor scris pe componente este

Pentru x ( 0 ) = { 0, 0, 0, 0 } primele iterate sunt

x ( 1 ) = { 6.32377 -0.771945 0.796784 1.01338 }

x ( 2 ) = { 4.13341 -0.923111 1.74546 1.23128 }

x ( 3 ) = { 3.77932 -1.41511 2.07691 1.28788 }

x ( 4 ) = { 3.6875 -1.50068 2.11412 1.30372 }

x ( 5 ) = { 3.6901 -1.50711 2.10718 1.30154 }

x ( 6 ) = { 3.69745 -1.50254 2.10156 1.3004 }

x ( 7 ) = { 3.69973 -1.5004 2.10006 1.30002 }

x ( 8 ) = { 3.70009 -1.49997 2.09991 1.29998 }

x ( 9 ) = { 3.70005 -1.49996 2.09996 1.29999 }

x ( 10 ) = { 3.70001 -1.49999 2.09999 1.3 }

x ( 11 ) = { 3.7 -1.5 2.1 1.3 }

x ( 12 ) = { 3.7 -1.5 2.1 1.3 }

Fie operatorul liniar generat de matricea , , unde dacă și atunci .

Se notează operatorul generat de matricea . se numește adjunctul lui și are proprietatea .

Observație. Fie o aplicație liniară. Se notează

0,

Propoziția 1. Fie operator liniar. Atunci

. (4.1.)

Demonstrație. Fie , . Fie astfel încât 0 Deci rezultă . Am demonstrat

Incluziunea este echivalentă cu ceea ce este identic cu Fie Atunci rezultă că Deci are loc (4.1.).

Corolar. Cu notațiile precedente avem:

Demonstrație. Formulele rezultă din (4.1.) deoarece .

Propoziția 2. Fie operator liniar. Notăm Atunci este o bijecție liniară.

Demonstrație. Pentru a arata că este o bijecție trebuie să arătăm că este injectiv și surjectiv.

Din definiția injectivității rezultă că trebuie să arătăm că 0. Fie 0 și 0. Deci este injectiv.

Trebuie să demonstrăm că surjectiv . Fie

(Am folosit faptul că

Deci ). Deci este

surjectiv.

Definiție. În condițiile prezentate mai sus există atunci , . Notăm și îl numim pseudoinversul lui T( sau pseudoinversa matricii ,unde este proiectorul generat de ).

Observatie. Fie proiectorul generat de . Atunci avem:

0

Observatie. Dacă este proiectorul generat de avem:

0

Observație. Dacă și 0 rezultă că există și

Teorema( Penrose). Operatorul este soluția unică a următorului sistem de ecuații operatoriale:

(4.2.)

(4.3.)

(4.4.)

unde este operator liniar.

Demonstrație. Fie Vom arăta că sunt verificate ecuațiile (4.2.), (4.3.), (4.4.).

Vrem să arătăm că care este echivalent cu Vom distinge două cazuri:

Cazul I. Luăm 00. Deci

Cazul II. Fie .

Deci . Din ,folosind liniaritatea operatorilor care intervin rezulta că

Să demonstrăm (4.3.) pentru . Avem într-adevăr: .

Apoi

Reciproc, fie operator liniar care verifică (4.2.), (4.3.), (4.4.). Din rezultă

Din (4.3.) rezultă

Teorema( Moore). Fie și Sunt echivalente:

0

Dacă atunci verifică și dacă verifică atunci .

Demonstrație. Presupunem că verifică : . Dar este unicul element din astfel încât . Deci ( și reciproc). Prin urmare relația este echivalentă cu .

Relația este echivalentă cu și atunci avem

, deci relația este echivalentă cu .

Fie apoi . Aplicăm și obținem . În continuare vom considera un element care verifică Știm că și ( am folosit , , ). Avem deci , , si deci .

Teorema 1. Fie bază algebrică în și liniar:

și

0.

Atunci .

Demonstrație. Deoarece este o bijecție liniară iar este bază algebrică în , rezultă că este bază în . Vom arăta că .

Dacă rezultă , , deci . Dacă adică 0, deci 0 și deci și din rezultă .

Arătăm că . Fie deoarece este bază în , rezultă . Deci .

În cazul în care avem 0 și rezultă 0=. Folosind și rezultă că .

Analog se arată că

Atunci .

Fie . Notăm aplicația liniară, .

Teorema 2. Fie liniară. Presupunem că dim . Fie și astfel încât dacă și . Atunci

.

Demonstrație. Din ipoteză va rezulta că este bază în . Trebuie să arătăm că este mulțime liniar independentă. Presupunem că 0. Înmulțim scalar cu și obținem 0 =

. Deci este bază în .

Fie atunci,

.

Dacă , deoarece rezultă atunci că ceea ce este echivalent cu . Observăm că 0 deci

{0}. Deci .

Deoarece , atunci și fie .. Avem deci , adica este simetrică. Apoi . Dacă implică 0 rezultă , dar cum , rezultă deci că ceea ce înseamnă că . Atunci este pozitiv definită. Avem simetrică și pozitiv definită, rezultă că există o bază -ortonormată în , pentru ( cu alte cuvinte există o bază în spațiu în care matricea este diagonală) și . Din rezultă

pentru

sau

pentru (4.5.) Din rezultă și deci . Deci

0 (4.6.)

Fie , , bază în , , atunci (4.5.) și (4.6.) devin

respectiv

Conform cu teorema 2 avem .

Are loc deci

Teorema 3. Dacă este o bază -ortogonală în atunci

.

Teorema 4( Algoritm de calcul pentru ). Fie liniar. Fie sistemele , și , . Fie operatori liniari construiți astfel: 0. Presupunem că am construit elemente din care sunt ortogonali și nenuli. Fie .

Fie astfel încât 0 și . Fie . Atunci .

Demonstrație. Fie ca în ipoteză 0. Se alege astfel încât 0. Se alege deci . Apoi . Presupunem apoi că am construit care sunt -ortogonali nenuli. Fie ca în enunț , 0 și . Vom arăta că sunt -ortogonali. Ca să fie tot sistemul -ortogonal avem de arătat că care este echivalent cu

Prin inducție deci sunt -ortogonali, rezultă conform teoremei 3 că .

Observație. Algoritmul poate continua atâta timp cât 0. Deci algoritmul se oprește atunci când 0 .

Observație. O bază în este . Fie , . Aplicăm și obținem . Deci formează sistem de generatori în . Algoritmul se oprește dacă . Fie , atunci algoritmul se oprește dacă 0 . Dar sunt coloanele lui .

Observație. În aplicații se consideră vectorii coloană

și .

Aplicație. Fie . Să se calculeze cu algoritmul de mai înainte .

Rezolvare.

Fie , rezultă . Se aleg vectorii coloană formați din coloanele lui :

, , .

Fie matricea . Fie astfel încât 0. Observăm că 0 deci algoritmul poate continua și vom calcula .Atunci , unde și

iar . Deci

.

Fie . Verificăm dacă algoritmul poate continua, adică dacă 0. Făcând calculele obținem

și

iar 0. Luăm

și , unde

și iar

. Obținem .

Fie astfel încât 0. Făcând calculul observăm că 0, deci algoritmul se oprește și

.

Metode directe

Metoda lui Gauss

#include<iomanip.h>

#include<iostream.h>

#include<math.h>

#include<conio.h>

#define N 50

//funcție pentru afișarea matricii și a vectorului

void afis_mat(int k,long double A[N][N],int n)

{

cout<<"\n Matricea A("<<k<<") este";

for(int j,i=0;i<n;i++)

{

cout<<"\n";

for(j=0;j<n;j++)

{

cout<<setw(15)<<A[i][j];

}

}

cout<<"\n Vectorul termenilor liberi b("<<k<<") este";

for(i=0;i<n;i++)

{

cout<<"\n "<<setw(15)<<A[i][n];

}

}

//funcție pentru aducerea maximului de pe coloana k pe linia k

void max_kk(long double A[N][N],int k,int n)

{

//det max pe coloana k

long double max=fabs(A[k][k]);

int imax=k;//pp că max este pe linia k

for(int i=k+1;i<n;i++)

{

if(fabs(A[i][k])>max)

{

imax=i;

max=fabs(A[i][k]);

}

}

cout<<"\n imax="<<imax<<" max="<<max;

//schimbarea liniilor k cu imax

long double t;

if(imax>k)

{

for(int j=k;j<=n;j++)

{

t=A[k][j];

A[k][j]=A[imax][j];

A[imax][j]=t;

}

}

}

//funcție care împarte linia k la a[k][k]

void impart(long double A[N][N],int n,int k)

{

long double m=A[k][k];

for(int j=k;j<=n;j++)

{

A[k][j]=A[k][j]/m;

}

}

//funcție care face zero sub a[k][k]

void zero(long double A[N][N],int n,int k)

{

for(int j,i=k+1;i<n;i++)

{

for(j=k+1;j<=n;j++)

{

A[i][j]=A[i][j]-A[k][j]*A[i][k];

}

}

for(i=k+1;i<n;i++)

{

A[i][k]=0;

}

}

//funcția Gauss

void gauss(long double A[N][N],long double x[N],int n)

{

for(int k=0;k<n;k++)

{

cout<<"\n Pasul "<<k+1;

max_kk(A,k,n);

if(A[k][k]==0)

{

cout<<"\n Determinantul sistemului este nul ";

return;

}

împart(A,n,k);

zero(A,n,k);

afis_mat(k+1,A,n);

}

x[n-1]=A[n-1][n];

for(int j,i=n-2;i>=0;i–)

{

long double s=0;

for(j=i+1;j<n;j++)

{

s=s+A[i][j]*x[j];

}

x[i]=A[i][n]-s;

}

cout<<"\n Vectorul soluție este ";

for(i=0;i<n;i++)

{

cout<<"\n"<<setw(10)<<x[i];

}

}

void main(void)

{

int n=4;

long double A[N][N]={{8.3,2.62,4.1,1.9,-10,65},

{3.92,8.45,7.78,2.46,12,21},

{3.77,7.21,8.04,2.28,15.45},

{2.21,3.65,1.69,6.99,-8.35}

};

long double x[N];

afis_mat(0,A,n);

gauss(A,x,n);

cout<<”\n”;

}

Metoda rădăcinii pătrate

#include<iomanip.h>

#include<iostream.h>

#include<math.h>

#include<conio.h>

#define N 30

//funcție pentru afișarea determinantului

void afis_det(long double A[N][N],int n)

{

cout<<"\n Determinantul este";

for(int j,i=0;i<n;i++)

{

cout<<"\n";

for(j=0;j<n;j++)

{

cout<<setw(15)<<A[i][j];

}

}

}

//funcție pentru afișarea vectorului

void afis_v(char *nume,long double a[N],int m)

{

cout<<"\n Vectorul "<<nume<<" este";

for(int i=0;i<m;i++)

{

cout<<"\n "<<setw(15)<<a[i];

}

}

//funcție pentru afișarea matricii

void afis_mat(char *nume,long double A[N][N],int n)

{

cout<<"\n Matricea "<<nume<<" este ";

for(int j,i=0;i<n;i++)

{

cout<<"\n";

for(j=0;j<n;j++)

{

cout<<setw(15)<<A[i][j];

}

}

}

//functie care verifica daca o matrice este simetrică

int simetrica(long double A[N][N],int n)

{

for(int j,i=0;i<n;i++)

{

for(j=0;j<n;j++)

{

if(A[i][j]!=A[j][i])

{

return 0;

}

}

}

return 1;

}

//funcție pentru aducerea maximului de pe coloana k pe linia k

int max_kk(long double A[N][N],int k,int n)

{

int ns=0;

//determinăm max pe coloana k

long double max=fabs(A[k][k]);

int imax=k;//presupunem că max este pe linia k

for(int i=k+1;i<n;i++)

{

if(fabs(A[i][k])>max)

{

imax=i;

max=fabs(A[i][k]);

}

}

//schimbarea liniilor k cu imax

long double t;

if(imax>k)

{

ns++;

for(int j=k;j<=n;j++)

{

t=A[k][j];

A[k][j]=A[imax][j];

A[imax][j]=t;

}

}

return ns;

}

//funcție care împarte linia k la akk

long double impart(long double A[N][N],int n,int k)

{

long double m=A[k][k];

for(int j=k;j<=n;j++)

{

A[k][j]=A[k][j]/m;

}

return m;

}

//funcție care face zero sub a[k][k]

void zero(long double A[N][N],int n,int k)

{

for(int j,i=k+1;i<n;i++)

{

for(j=k+1;j<=n;j++)

{

A[i][j]=A[i][j]-A[k][j]*A[i][k];

}

}

for(i=k+1;i<n;i++)

{

A[i][k]=0;

}

}

//funcție care calculează determinantul

long double det_gauss(long double A[N][N],int n)

{

long double d=1;

for(int k=0;k<n;k++)

{

int ns=max_kk(A,k,n);

if(A[k][k]==0)

{

cout<<"\n Determinantul sistemului este nul ";

return 0;

}

d=d*impart(A,n,k)*pow(-1,ns);

zero(A,n,k);

}

return d;

}

//funcție care calculează matricea B

void calc_B(long double A[N][N],long double B[N][N],int m)

{

//i<j rezulta B[i][j]=0

for(int j,i=0;i<m;i++)

{

for(j=0;j<m;j++)

{

if(i<j)

{

B[i][j]=0;

}

}

}

//calculăm celelalte elememte din B

//calculăm B[i][i]

long double s;

for(i=0;i<m;i++)

{

s=0;

for(int k=0;k<=i-1;k++)

{

s=s+fabs(pow(B[i][k],2));

}

B[i][i]=sqrt(A[i][i]-s);

//calculăm B[j][i]

for(j=i+1;j<m;j++)

{

s=0;

for(k=0;k<=i-1;k++)

{

s=s+B[j][k]*B[i][k];

}

B[j][i]=1/B[i][i]*(A[j][i]-s);

}

}

}

//funcție care calculează vectorul c

void calc_c(long double B[N][N],long double b[N],long double c[N],

int m)

{

long double s;

for(int i=0;i<m;i++)

{

s=0;

for(int k=0;k<=i-1;k++)

{

s=s+B[i][k]*c[k];

}

c[i]=1/B[i][i]*(b[i]-s);

}

}

//funcție care calculează vectorul soluție x

void calc_x(long double B[N][N],long double c[N],long double x[N],

int m)

{

long double s;

for(int k,i=m-1;i>=0;i–)

{

s=0;

for(k=i+1;k<m;k++)

{

s=s+B[k][i]*x[k];

}

x[i]=1/B[i][i]*(c[i]-s);

}

}

void main(void)

{

int m=4;

long double A[N][N]={{1,0.42,0.54,0.66},

{0.42,1,0.32,0.44},

{0.54,0.32,1,0.22},

{0.66,0.44,0.22,1}

};

long double b[N]={0.3,0.5,0.7,0.9};

long double v_det,Det[N][N],B[N][N],c[N],x[N];

afis_mat("A",A,m);

int poz_def=1;

for(int k=1;k<=m;k++)

{

for(int i=0;i<k;i++)

{

for(int j=0;j<k;j++)

{

Det[i][j]=A[i][j];

}

}

afis_det(Det,k);

v_det=det_gauss(Det,k);

cout<<"\n Valoarea determinantului este:"<<v_det;

if(v_det<=0)

{

poz_def=0;

}

}

if(simetrica(A,m)==0)

{

cout<<"\n A nu este simetrica";

}

else

{

cout<<"\n A este simetrica";

if(poz_def==0)

{

cout<<"\n A nu este pozitiv definita";

}

else

{

cout<<"\n A este pozitiv definita";

}

calc_B(A,B,m);

afis_mat("B",B,m);

calc_c(B,b,c,m);

afis_v("c",c,m);

calc_x(B,c,x,m);

afis_v("x",x,m);

}

cout<<"\n";

}

Metoda lui Khaletski

#include<iomanip.h>

#include<iostream.h>

#include<math.h>

#include<conio.h>

#define N 30

//funcție pentru afișarea vectorului

void afis_v(char *nume,long double a[N],int m)

{

cout<<"\n Vectorul "<<nume<<" este";

for(int i=0;i<m;i++)

{

cout<<"\n "<<setw(15)<<a[i];

}

}

//funcție pentru afișarea matricii

void afis_mat(char *nume,long double A[N][N],int n)

{

cout<<"\n Matricea "<<nume<<" este ";

for(int j,i=0;i<n;i++)

{

cout<<"\n";

for(j=0;j<n;j++)

{

cout<<setw(15)<<A[i][j];

}

}

}

//funcție care verifică dacă o matrice este simetrică

int simetrica(long double A[N][N],int n)

{

for(int j,i=0;i<n;i++)

{

for(j=0;j<n;j++)

{

if(A[i][j]!=A[j][i])

{

return 0;

}

}

}

return 1;

}

//funcție care calculează matricea B și matricea C

void calc_BC(long double A[N][N],long double B[N][N],

long double C[N][N],int m)

{

for(int j,i=0;i<m;i++)

{

for(j=0;j<m;j++)

{

if(i<j)

{

B[i][j]=0;

}

if(i==j)

{

C[i][j]=1;

}

if(i>j)

{

C[i][j]=0;

}

}

}

long double s;

int k,p;

for(p=0;p<m;p++)//parcurgem pe diagonală

{

//calculam coloana p din B

for(i=p;i<m;i++)

{

s=0;

for(k=0;k<=p-1;k++)

{

s=s+B[i][k]*C[k][p];

}

B[i][p]=A[i][p]-s;

}

for(j=p+1;j<m;j++)

{

s=0;

for(k=0;k<=p-1;k++)

{

s=s+B[p][k]*C[k][j];

}

C[p][j]=(1/B[p][p])*(A[p][j]-s);

}

}

}

//funcție care calculează vectorul y

void calc_y(long double B[N][N],long double y[N],long double b[N],

int m)

{

long double s;

for(int k,i=0;i<m;i++)

{

s=0;

for(k=0;k<=i-1;k++)

{

s=s+B[i][k]*y[k];

}

y[i]=(1/B[i][i])*(b[i]-s);

}

}

//funcție care calculează vectorul soluție

void calc_x(long double C[N][N],long double y[N],long double x[N],

int m)

{

long double s;

for(int k,i=m-1;i>=0;i–)

{

s=0;

for(k=i+1;k<m;k++)

{

s=s+C[i][k]*x[k];

}

x[i]=y[i]-s;

}

}

void main(void)

{

int m=3;

long double B[N][N],C[N][N],y[N],x[N];

long double A[N][N]={{1,4,5},

{4,20,32},

{5,32,70}

};

long double b[N]={10,56,107};

afis_mat("A",A,m);

afis_v("b",b,m);

calc_BC(A,B,C,m);

afis_mat("B",B,m);

afis_mat("C",C,m);

calc_y(B,y,b,m);

afis_v("y",y,m);

calc_x(C,y,x,m);

afis_v("x",x,m);

cout<<"\n";

}

Metoda lui Ritz de inversare a unei matrici simetrice și pozitiv definite

#include<iomanip.h>

#include<iostream.h>

#include<math.h>

#include<stdio.h>

#include<string.h>

#define N 20

//funcție pentru afișarea unei matrici

void afis_m(char *nume,long double A[N][N],int m)

{

cout<<"\n Matricea "<<nume<<" este";

for(int j,i=0;i<m;i++)

{

cout<<"\n";

for(j=0;j<m;j++)

{

cout<<setw(15)<<setiosflags(ios::fixed)<<A[i][j];

}

}

}

//funcție pentru afișarea unui determinant

void afis_det(long double A[N][N],int n)

{

cout<<"\n Determinantul este";

for(int j,i=0;i<n;i++)

{

cout<<"\n";

for(j=0;j<n;j++)

{

cout<<setw(15)<<A[i][j];

}

}

}

//funcție pentru afișarea matricei inițiale

void H_0(long double H0[N][N],int m)

{

for(int j,i=0;i<m;i++)

{

for(j=0;j<m;j++)

{

H0[i][j]=0;

}

}

}

//funcție de afișare a vectorului y

void afis_y(long double A[N][N],int m)

{

for(int i,k=1;k<=m;k++)

{

cout<<"\n y"<<k;

for(i=0;i<m;i++)

{

cout<<"\n"<<A[i][k-1];

}

}

}

//funcție pentru calculul lui

void calc_xk(long double yk[N],long double Hk_1[N][N],

long double A[N][N],long double xk[N],int m)

{

long double P[N][N];

for(int k,j,i=0;i<m;i++)

{

for(k=0;k<m;k++)

{

P[i][k]=0;

for(j=0;j<m;j++)

{

P[i][k]=P[i][k]+Hk_1[i][j]*A[j][k];

}

}

}

afis_m("P",P,m);

long double V[N];

for(i=0;i<m;i++)

{

V[i]=0;

for(j=0;j<m;j++)

{

V[i]=V[i]+P[i][j]*yk[j];

}

}

cout<<"\n V=";

for(i=0;i<m;i++)

{

cout<<setw(14)<<V[i];

}

for(i=0;i<m;i++)

{

xk[i]=yk[i]-V[i];

}

cout<<"\n xk=";

for(i=0;i<m;i++)

{

cout<<setw(14)<<xk[i];

}

}

//funcție pentru calculul lui

void calc_Hk(long double Hk[N][N],long double Hk_1[N][N],

long double xk[N],long double A[N][N],int m)

{

long double Axk[N];

for(int j,i=0;i<m;i++)

{

Axk[i]=0;

for(j=0;j<m;j++)

{

Axk[i]=Axk[i]+A[i][j]*xk[j];

}

}

cout<<"\n Axk=";

for(i=0;i<m;i++)

{

cout<<setw(14)<<Axk[i];

}

long double Axk_xk=0;

for(i=0;i<m;i++)

{

Axk_xk=Axk_xk+Axk[i]*xk[i];

}

cout<<"\n Axk_xk="<<Axk_xk;

//wk*traspTwk

long double xk_tranxk[N][N];

int k;

for(k,i=0;i<m;i++)

{

for(k=0;k<m;k++)

{

xk_tranxk[i][k]=xk[i]*xk[k];

}

}

afis_m("xk_tranxk",xk_tranxk,m);

long double impart[N][N];

for(i=0;i<m;i++)

{

for(j=0;j<m;j++)

{

impart[i][j]=xk_tranxk[i][j]/Axk_xk;

}

}

for(i=0;i<m;i++)

{

for(j=0;j<m;j++)

{

Hk[i][j]=Hk_1[i][j]+impart[i][j];

}

}

}

//funcție pentru verificarea condiției 0

int alg_cont(long double yk[N],long double Hk_1[N][N],

long double A[N][N],int m)

{

long double Hk_1A[N][N];

for(int k,j,i=0;i<m;i++)

{

for(k=0;k<m;k++)

{

Hk_1A[i][k]=0;

for(j=0;j<m;j++)

{

Hk_1A[i][k]=Hk_1A[i][k]+Hk_1[i][j]*A[j][k];

}

}

}

long double Hk_1Ay[N];

for(i=0;i<m;i++)

{

Hk_1Ay[i]=0;

for(j=0;j<m;j++)

{

Hk_1Ay[i]=Hk_1Ay[i]+Hk_1A[i][j]*yk[j];

}

}

for(i=0;i<m;i++)

{

if(fabs(yk[i]-Hk_1Ay[i])>0.0000000000001)

{

return 1;

}

}

return 0;

}

//funcție pentru aducerea maximului de pe coloana k pe linia k

int max_kk(long double A[N][N],int k,int n)

{

int ns=0;

//determinăm max pe coloana k

long double max=fabs(A[k][k]);

int imax=k;//presupunem că max este pe linia k

for(int i=k+1;i<n;i++)

{

if(fabs(A[i][k])>max)

{

imax=i;

max=fabs(A[i][k]);

}

}

long double t;

if(imax>k)

{

ns++;

for(int j=k;j<=n;j++)

{

t=A[k][j];

A[k][j]=A[imax][j];

A[imax][j]=t;

}

}

return ns;

}

//funcție care împarte linia k la a[k][k]

long double impart(long double A[N][N],int n,int k)

{

long double m=A[k][k];

for(int j=k;j<=n;j++)

{

A[k][j]=A[k][j]/m;

}

return m;

}

//funcție care face zero sub a[k][k]

void zero(long double A[N][N],int n,int k)

{

for(int j,i=k+1;i<n;i++)

{

for(j=k+1;j<=n;j++)

{

A[i][j]=A[i][j]-A[k][j]*A[i][k];

}

}

for(i=k+1;i<n;i++)

{

A[i][k]=0;

}

}

//funcție care calculează determinantul cu ajutorul metodei lui Gauss

long double det_gauss(long double A[N][N],int n)

{

long double d=1;

for(int k=0;k<n;k++)

{

int ns=max_kk(A,k,n);

if(A[k][k]==0)

{

cout<<"\n Determinantul sistemului este nul ";

return 0;

}

d=d*impart(A,n,k)*pow(-1,ns);

zero(A,n,k);

}

return d;

}

//funcție care verifică dacă o matrice este simetrică

int simetrica(long double A[N][N],int n)

{

for(int j,i=0;i<n;i++)

{

for(j=0;j<n;j++)

{

if(A[i][j]!=A[j][i])

{

return 0;

}

}

}

return 1;

}

void main(void)

{

int m=4;

long double H0[N][N],H1[N][N],xk[N],yk[N],v_det,Det[N][N];

long double A[N][N]={{50,2,1,1},

{2,50,1,1},

{1,1,50,-1},

{1,1,-1,50}

};

afis_m("A",A,m);

int poz_def=1;

int i,j,k;

for( k=1;k<=m;k++)

{

for(i=0;i<k;i++)

{

for(j=0;j<k;j++)

{

Det[i][j]=A[i][j];

}

}

afis_det(Det,k);

v_det=det_gauss(Det,k);

cout<<"\n Valoarea determinantului este:"<<v_det;

if(v_det<=0)

{

poz_def=0;

}

}

if(poz_def==0)

{

cout<<"\n A nu este pozitiv definită";

}

else

{

cout<<"\n A este pozitiv definită";

if(simetrica(A,m)==0)

{

cout<<"\n A nu este simetrică";

}

else

{

cout<<"\n A este simetrică";

H_0(H0,m);

afis_m("H0",H0,m);

afis_y(A,m);

for(k=1; ;k++)

{

for(i=0;i<m;i++)

{

yk[i]=A[i][k-1];

}

if(!alg_cont(yk,H0,A,m))

{

cout<<"\n Algoritmul se oprește ";

break;

}

calc_xk(yk,H0,A,xk,m);

calc_Hk(H1,H0,xk,A,m);

char nume[5];

sprintf(nume,"H%d",k);

afis_m(nume,H1,m);

for(i=0;i<m;i++)

{

for(j=0;j<m;j++)

{

H0[i][j]=H1[i][j];

}

}

}

}

}

cout<<"\n";

}

Metode iterative

Metoda lui Jacobi pentru matrici dominant diagonale pe linie

#include<stdio.h>

#include<math.h>

#include<conio.h>

#include<iostream.h>

#include<iomanip.h>

#define N 20

//funcție pentru afișarea vectorului

void afis_v(char *nume,long double a[N],int m)

{

cout<<"\n Vectorul "<<nume<<" este";

for(int i=0;i<m;i++)

{

cout<<"\n "<<setw(15)<<a[i];

}

}

//funcție pentru afișarea matricii

void afis_m(char *nume,long double A[N][N],int m)

{

printf("\n Matricea %s este :",nume);

for(int j,i=0;i<m;i++)

{

printf("\n ");

for(j=0;j<m;j++)

{

printf("%11.6Lf",A[i][j]);

}

}

}

//funcție care calculează produsul dintre două matrici

void produs_m(long double D[N][N],long double A[N][N],

long double D_1A[N][N],int ld,int cd,int ca)

{

for(int j,k,i=0;i<ld;i++)

{

for(k=0;k<ca;k++)

{

D_1A[i][k]=0;

for(j=0;j<cd;j++)

{

D_1A[i][k]=D_1A[i][k]+D[i][j]*A[j][k];

}

}

}

}

//funcție care afișează matricea D

void matr_D(long double D[N][N],long double A[N][N],int n)

{

for(int j,i=0;i<n;i++)

{

for(j=0;j<n;j++)

{

if(i==j)

{

D[i][j]=A[i][j];

}

else

{

D[i][j]=0;

}

}

}

}

//funcție care calculează matricea

void matr_D_1(long double D[N][N],int n)

{

for(int i=0;i<n;i++)

{

D[i][i]=1/D[i][i];

}

}

//funcție care afișează matricea identitate I

void matr_I(long double I[N][N],int n)

{

for(int j,i=0;i<n;i++)

{

for(int j=0;j<n;j++)

{

if(i==j)

{

I[i][j]=1;

}

else

{

I[i][j]=0;

}

}

}

}

//funcție care verifică dacă o matrice este dominant diagonală pe linie

int dominant_diag_linie(long double A[N][N],int m)

{

long double s;

for(int j,i=0;i<m;i++)

{

for(s=0,j=0;j<m;j++)

{

if(i==j)

{

continue;

}

s=s+fabs(A[i][j]);

}

if(fabs(A[i][i])<=s)

{

return 0;

}

}

return 1;

}

//funcție care returnează valoarea lui q

long double ret_q(long double B[N][N],int n)

{

long double suma_linie[N];

for(int j,i=0;i<n;i++)

{

suma_linie[i]=0;

for(j=0;j<n;j++)

{

suma_linie[i]=suma_linie[i]+fabs(B[i][j]);

}

}

long double max=suma_linie[0];

for(i=1;i<n;i++)

{

if(suma_linie[i]>max)

{

max=suma_linie[i];

}

}

return max;

}

//funcție pentru diferența a două matrici

void diferenta(long double I[N][N],long double D_1A[N][N],

long double B[N][N],int n)

{

for(int j,i=0;i<n;i++)

{

for( j=0;j<n;j++)

{

B[i][j]=I[i][j]-D_1A[i][j];

}

}

}

//funcție care împarte termenul liber la a[i][i]

void term_liber(long double A[N][N],long double b[N],

long double a[N],int n)

{

for(int i=0;i<n;i++)

{

b[i]=a[i]/A[i][i];

}

}

//funcție care calculează iterațiile

void iteratii(long double B[N][N],long double xn_1[N],int n,

long double b[N],long double xn[N])

{

long double s;

for(int j,i=0;i<n;i++)

{

s=0;

for(j=0;j<n;j++)

{

s=s+B[i][j]*xn_1[j];

}

xn[i]=s+b[i];

}

}

//funcție care calculează

long double norma1(long double xn_1[N],long double xn[N],int n)

{

long double norm[N];

for(int i=0;i<n;i++)

{

norm[i]=fabs(xn_1[i]-xn[i]);

}

long double max=norm[0];

for(i=1;i<n;i++)

{

if(norm[i]>max)

{

max=norm[i];

}

}

return max;

}

void main(void)

{

long double b[N],x[N],epsilon=0.0001;

int n=3;

long double A[N][N]={{3.1,1.5,1},

{1.5,2.5,0.5},

{1,0.5,4.2}

};

long double a[N]={10.83,9.2,17.1};

long double x0[N]={0,0,0,0};

afis_m("A",A,n);

afis_v("a",a,n);

afis_v("x0",x0,n);

if(dominant_diag_linie(A,n))

{

cout<<"\n Matricea A este dominant diagonală pe linie ";

}

else

{

cout<<"\n Matricea nu este dominant diagonală pe linie";

}

long double B[N][N],I[N][N],D[N][N],D_1A[N][N];

matr_D(D,A,n);

afis_m("D",D,n);

matr_I(I,n);

afis_m("I",I,n);

matr_D_1(D,n);

afis_m("D-1",D,n);

produs_m(D,A,D_1A,n,n,n);

afis_m("D-1*A",D_1A,n);

diferenta(I,D_1A,B,n);

afis_m("B=I-D-1*A",B,n);

term_liber(A,b,a,n);

afis_v("b",b,n);

long double q=ret_q(B,n);

cout<<"\n q="<<q;

if(q<1)

{

cout<<"\n Metoda converge";

}

else

{

cout<<"\n Metoda nu converge";

}

if(dominant_diag_linie(A,n)&&q<1)

{

cout<<"\n Epsilon este "<<epsilon;

for(int k=0; ;k++)

{

iteratii(B,x0,n,b,x);

cout<<"\n Iteratia "<<k+1;

cout<<"\n x={" ;

for(int i=0;i<n;i++)

{

cout<<setw(15)<<x[i];

}

cout<<"}";

if(q/(1-q)*norma1(x0,x,n)<epsilon)

{

break;

}

for( i=0;i<n;i++)

{

x0[i]=x[i];

}

}

cout<<"\n Numarul de iteratii "<<k+1;

}

cout<<”\n”;

}

Metoda lui Jacobi pentru matrici dominant diagonale pe coloană

#include<stdio.h>

#include<math.h>

#include<conio.h>

#include<iostream.h>

#include<iomanip.h>

#define N 30

//funcție pentru afișarea matricii

void afis_m(char *nume,long double A[N][N],int m)

{

printf("\n Matricea %s este :",nume);

for(int j,i=0;i<m;i++)

{

printf("\n ");

for(j=0;j<m;j++)

{

printf("%11.6Lf",A[i][j]);

}

}

}

//funcție pentru afișarea unui vector

void afis_v(char *nume,long double v[N],int m)

{

printf("\n Vectorul %s este : ",nume);

for(int i=0;i<m;i++)

{

printf("\n %11.6Lf",v[i]);

}

}

//funcție care calculează produsul dintre două matrici

void produs_m(long double A[N][N],long double B[N][N],

long double C[N][N],int la,int ca,int cb)

{

for(int j,k,i=0;i<la;i++)

{

for(k=0;k<cb;k++)

{

C[i][k]=0;

for(j=0;j<ca;j++)

{

C[i][k]=C[i][k]+A[i][j]*B[j][k];

}

}

}

}

//funcție care afișează matricea D

void matr_D(long double D[N][N],long double A[N][N],int n)

{

for(int j,i=0;i<n;i++)

{

for(j=0;j<n;j++)

{

if(i==j)

{

D[i][j]=A[i][j];

}

else

{

D[i][j]=0;

}

}

}

}

//funcție care calculează matricea

void matr_D_1(long double D[N][N],int n)

{

for(int i=0;i<n;i++)

{

D[i][i]=1/D[i][i];

}

}

//funcție care afișează matricea unitate I

void matr_I(long double I[N][N],int n)

{

for(int j,i=0;i<n;i++)

{

for(int j=0;j<n;j++)

{

if(i==j)

{

I[i][j]=1;

}

else

{

I[i][j]=0;

}

}

}

}

//funcție care verifică dacă o funcție este dominant diagonală pe coloană

int dominant_diag_col(long double A[N][N],int m)

{

long double s;

for(int i,j=0;j<m;j++)

{

for(s=0,i=0;i<m;i++)

{

if(i==j)

{

continue;

}

s=s+fabs(A[i][j]);

}

if(fabs(A[j][j])<=s)

{

return 0;

}

}

return 1;

}

//funcție care returneză valoarea lui q

long double ret_q(long double C[N][N],int n)

{

long double suma_col[N];

for(int i,j=0;j<n;j++)

{

suma_col[j]=0;

for(i=0;i<n;i++)

{

suma_col[j]=suma_col[j]+fabs(C[i][j]);

}

}

long double max=suma_col[0];

for(i=1;i<n;i++)

{

if(suma_col[i]>max)

{

max=suma_col[i];

}

}

return max;

}

//funcție pentru diferența a două matrici

void diferenta(long double A[N][N],long double B[N][N],

long double C[N][N],int n)

{

for(int j,i=0;i<n;i++)

{

for( j=0;j<n;j++)

{

C[i][j]=A[i][j]-B[i][j];

}

}

}

//funcție care calculează iterațiile

void iteratii(long double C[N][N],long double yn_1[N],int n,

long double a[N],long double yn[N])

{

long double s;

for(int j,i=0;i<n;i++)

{

s=0;

for(j=0;j<n;j++)

{

s=s+C[i][j]*yn_1[j];

}

yn[i]=s+a[i];

}

}

// funcție care calculează

long double norma1(long double yn_1[N],long double yn[N],int n)

{

long double s=fabs(yn_1[0]-yn[0]);

for(int i=1;i<n;i++)

{

s=s+fabs(yn_1[i]-yn[i]);

}

return s;

}

void main(void)

{

long double y[N],epsilon=0.0001;

int n=4;

long double A[N][N]={{10,2,2,1},

{2,10,1,1},

{2,1,10,1},

{1,1,1,10}

};

long double a[N]={1,1,1,1};

long double y0[N]={0,0,0,0};

afis_m("A",A,n);

afis_v("a",a,n);

afis_v("y0",y0,n);

if(dominant_diag_col(A,n))

{

cout<<"\n Matricea A este dominant diagonală pe coloană ";

}

else

{

cout<<"\n Matricea nu este dominant diagonală pe coloană";

}

long double C[N][N],I[N][N],D[N][N],P[N][N];

matr_D(D,A,n);

afis_m("D",D,n);

matr_I(I,n);

afis_m("I",I,n);

matr_D_1(D,n);

afis_m("D-1",D,n);

produs_m(A,D,P,n,n,n);

afis_m("A*D-1",P,n);

diferenta(I,P,C,n);

afis_m("C=I-A*D-1",C,n);

long double q=ret_q(C,n);

cout<<"\n q="<<q;

if(q<1)

{

cout<<"\n Metoda converge";

}

else

{

cout<<"\n Metoda nu converge";

}

if(dominant_diag_col(A,n)&&q<1)

{

cout<<"\n Epsilon este "<<epsilon;

for(int k=0; ;k++)

{

iteratii(C,y0,n,a,y);

cout<<"\n Iteratia "<<k+1;

cout<<"\n y={" ;

for(int i=0;i<n;i++)

{

cout<<setw(15)<<y[i];

}

cout<<"}";

if(q/(1-q)*norma1(y0,y,n)<epsilon)

{

break;

}

for( i=0;i<n;i++)

{

y0[i]=y[i];

}

}

cout<<"\n Numarul de iterații "<<k+1;

long double x_final[N];

for(int i=0;i<n;i++)

{

x_final[i]=y[i]/A[i][i];

}

cout<<"\n X_final={";

for(i=0;i<n;i++)

{

cout<<setw(10)<<x_final[i];

}

cout<<"}";

cout<<”\n”;

}

Metoda Gauss-Seidel

#include<stdio.h>

#include<math.h>

#include<conio.h>

#include<iostream.h>

#include<iomanip.h>

#define N 20

//funcție pentru afișarea unui vector

void afis_v(char *nume,long double a[N],int m)

{

cout<<"\n Vectorul "<<nume<<" este";

for(int i=0;i<m;i++)

{

cout<<"\n "<<setw(15)<<a[i];

}

}

//funcție pentru afișarea unei matrici

void afis_m(char *nume,long double A[N][N],int m)

{

printf("\n Matricea %s este :",nume);

for(int j,i=0;i<m;i++)

{

printf("\n ");

for(j=0;j<m;j++)

{

printf("%11.6Lf",A[i][j]);

}

}

}

//funcție care calculează produsul dintre două matrici

void produs_m(long double D[N][N],long double A[N][N],

long double D_1A[N][N],int ld,int cd,int ca)

{

for(int j,k,i=0;i<ld;i++)

{

for(k=0;k<ca;k++)

{

D_1A[i][k]=0;

for(j=0;j<cd;j++)

{

D_1A[i][k]=D_1A[i][k]+D[i][j]*A[j][k];

}

}

}

}

//funcție care afișează matricea D

void matr_D(long double D[N][N],long double A[N][N],int n)

{

for(int j,i=0;i<n;i++)

{

for(j=0;j<n;j++)

{

if(i==j)

{

D[i][j]=A[i][j];

}

else

{

D[i][j]=0;

}

}

}

}

//funcție care calculează matricea

void matr_D_1(long double D[N][N],int n)

{

for(int i=0;i<n;i++)

{

D[i][i]=1/D[i][i];

}

}

//funcție care afișează matricea unitate I

void matr_I(long double I[N][N],int n)

{

for(int j,i=0;i<n;i++)

{

for(int j=0;j<n;j++)

{

if(i==j)

{

I[i][j]=1;

}

else

{

I[i][j]=0;

}

}

}

}

//funcție pentru diferența a două matrici

void diferenta(long double I[N][N],long double D_1A[N][N],

long double B[N][N],int n)

{

for(int j,i=0;i<n;i++)

{

for( j=0;j<n;j++)

{

B[i][j]=I[i][j]-D_1A[i][j];

}

}

}

//funcție care verifică dacă sunt îndeplinite condițiile pentru A

int conditii_A(long double A[N][N],int m)

{

long double s1,s2;

for(int j,k=0;k<m;k++)

{

s1=0;

for(j=0;j<m;j++)

{

if(k==j)

{

continue;

}

s1=s1+fabs(A[k][j]);

}

if(fabs(A[k][k])<s1)

{

return 0;

}

s2=0;

for(j=k+1;j<m;j++)

{

s2=s2+fabs(A[k][j]);

}

if(s2>=fabs(A[k][k]))

{

return 0;

}

}

return 1;

}

//funcție care verifică dacă sunt îndeplinite condițiile pentru B

int conditii_B(long double B[N][N],int m)

{

long double s1,s2;

for(int j,k=0;k<m;k++)

{

s1=0;

for(j=0;j<m;j++)

{

s1=s1+fabs(B[k][j]);

}

s1=int(s1*pow(10,6))/pow(10,6);

if(s1>1)

{

return 0;

}

s2=0;

for(j=k;j<m;j++)

{

s2=s2+fabs(B[k][j]);

}

s2=int(s2*pow(10,6))/pow(10,6);

if(s2>=1)

{

return 0;

}

}

return 1;

}

//funcție care afișează matricea L

void mat_L(long double B[N][N],long double L[N][N],int m)

{

for(int j,k=0;k<m;k++)

{

for(j=0;j<m;j++)

{

if(k>j)

{

L[k][j]=B[k][j];

}

else

{

L[k][j]=0;

}

}

}

}

//funcție care afișează matricea R

void mat_R(long double B[N][N],long double R[N][N],int m)

{

for(int j,k=0;k<m;k++)

{

for(j=0;j<m;j++)

{

if(k<=j)

{

R[k][j]=B[k][j];

}

else

{

R[k][j]=0;

}

}

}

}

//funcție care calculează iterațiile

void iteratii(long double B[N][N],long double xn[N],

long double xn_1[N], long double b[N],int m)

{

long double s1,s2;

for(int k=0;k<m;k++)

{

s1=0,s2=0;

for(int j=0;j<=k-1;j++)

{

s1=s1+B[k][j]*xn[j];

}

for(j=k;j<m;j++)

{

s2=s2+B[k][j]*xn_1[j];

}

xn[k]=s1+s2+b[k];

}

}

//funcție care împarte termenul liber la a[k][k]

void termen_liber(long double A[N][N],long double a[N],

long double b[N],int m)

{

for(int k=0;k<m;k++)

{

b[k]=a[k]/A[k][k];

}

}

//funcție care returnează valoarea lui q

long double ret_q(long double B[N][N],int m,long double q_v[N])

{

long double s1,s2;

for(int j,k=0;k<m;k++)

{

s1=0,s2=0;

for(j=0;j<=k-1;j++)

{

s1=s1+fabs(B[k][j])*q_v[j];

}

for(j=k;j<m;j++)

{

s2=s2+fabs(B[k][j]);

}

q_v[k]=s1+s2;

}

long double max=q_v[0];

for(k=1;k<m;k++)

{

if(q_v[k]>max)

{

max=q_v[k];

}

}

return max;

}

//funcție care calculează

long double norma1(long double xn_1[N],long double xn[N],int n)

{

long double norm[N];

for(int i=0;i<n;i++)

{

norm[i]=fabs(xn_1[i]-xn[i]);

}

long double max=norm[0];

for(i=1;i<n;i++)

{

if(norm[i]>max)

{

max=norm[i];

}

}

return max;

}

void main(void)

{

int m=3;

long double D[N][N],I[N][N],D_1A[N][N],B[N][N],L[N][N],R[N][N];

long double A[N][N]={{1.02,-0.25,-0.3},

{-0.41,1.13,-0.15},

{-0.25,-0.14,1.21}

};

long double a[N]={0.515,1.555,2.78};

long double x[N],b[N];

cout<<"\n Consideram sistemul Ax=a ";

afis_m(" A",A,m);

afis_v(" a",a,m);

if(conditii_A(A,m))

{

cout<<"\n Conditii indeplinite pentru A";

}

else

{

cout<<"\n Nu sunt indeplinite conditiile pentru A";

}

if(conditii_A(A,m))

{

matr_D(D,A,m);

afis_m("D",D,m);

matr_D_1(D,m);

afis_m("D_1",D,m);

matr_I(I,m);

afis_m("I",I,m);

produs_m(D,A,D_1A,m,m,m);

afis_m("D_1*A",D_1A,m);

diferenta(I,D_1A,B,m);

afis_m("B=I-D_1*A",B,m);

termen_liber(A,a,b,m);

afis_v("b",b,m);

if(conditii_B(B,m))

{

cout<<"\n Conditii indeplinite pentru B";

}

else

{

cout<<"\n Conditii neeindeplinite pentru B";

}

long double q,q_v[N];

q=ret_q(B,m,q_v);

afis_v("q_v",q_v,m);

cout<<"\n q="<<q;

if(q<1)

{

cout<<"\n Metoda converge deoarece q="<<q<<" < 1";

}

else

{

cout<<"\n Metoda nu converge";

}

mat_L(B,L,m);

afis_m("L",L,m);

mat_R(B,R,m);

afis_m("R",R,m);

long double x0[N]={0,0,0,0};

long double epsilon=0.001;

cout<<"\n Epsilon este "<<epsilon;

afis_v("x0",x0,m);

for(int k=0; ;k++)

{

iteratii(B,x,x0,b,m);

cout<<"\n Iteratia "<<k+1;

cout<<"\n x={" ;

for(int i=0;i<m;i++)

{

cout<<setw(15)<<x[i];

}

cout<<"}";

if(q/(1-q)*norma1(x0,x,m)<epsilon)

{

break;

}

for( i=0;i<m;i++)

{

x0[i]=x[i];

}

}

cout<<"\n Nr iteratii "<<k+1;

}

cout<<”\n”;

}

Metoda relaxării simultane

#include<iomanip.h>

#include<iostream.h>

#include<math.h>

#include<conio.h>

#define N 30

//funcție pentru afișarea determinantului

void afis_det(long double A[N][N],int n)

{

cout<<"\n Determinantul este";

for(int j,i=0;i<n;i++)

{

cout<<"\n";

for(j=0;j<n;j++)

{

cout<<setw(15)<<A[i][j];

}

}

}

//funcție pentru afișarea unui vector

void afis_v(char *nume,long double a[N],int m)

{

cout<<"\n Vectorul "<<nume<<" este";

for(int i=0;i<m;i++)

{

cout<<"\n "<<setw(15)<<a[i];

}

}

//funcție care calculează produsul dintre două matrici

void produs_m(long double D[N][N],long double A[N][N],

long double D_1A[N][N],int ld,int cd,int ca)

{

for(int j,k,i=0;i<ld;i++)

{

for(k=0;k<ca;k++)

{

D_1A[i][k]=0;

for(j=0;j<cd;j++)

{

D_1A[i][k]=D_1A[i][k]+D[i][j]*A[j][k];

}

}

}

}

//funcție care afișează matricea D

void matr_D(long double D[N][N],long double A[N][N],int n)

{

for(int j,i=0;i<n;i++)

{

for(j=0;j<n;j++)

{

if(i==j)

{

D[i][j]=A[i][j];

}

else

{

D[i][j]=0;

}

}

}

}

//funcție care calculează matricea

void matr_D_1(long double D[N][N],int n)

{

for(int i=0;i<n;i++)

{

D[i][i]=1/D[i][i];

}

}

//funcție care afișează matricea unitate I

void matr_I(long double I[N][N],int n)

{

for(int j,i=0;i<n;i++)

{

for(int j=0;j<n;j++)

{

if(i==j)

{

I[i][j]=1;

}

else

{

I[i][j]=0;

}

}

}

}

//funcție care returnează valoarea lui lambda

long double lambda1(long double suma_linie[N],

long double D_1A[N][N],int n)

{

for(int j,i=0;i<n;i++)

{

suma_linie[i]=0;

for(j=0;j<n;j++)

{

suma_linie[i]=suma_linie[i]+fabs(D_1A[i][j]);

}

}

long double max=suma_linie[0];

for(i=1;i<n;i++)

{

if(suma_linie[i]>max)

{

max=suma_linie[i];

}

}

return max;

}

//funcție pentru diferența a doua matrici

void diferenta(long double I[N][N],long double D_1A[N][N],

long double B[N][N],int n)

{

for(int j,i=0;i<n;i++)

{

for( j=0;j<n;j++)

{

B[i][j]=I[i][j]-D_1A[i][j];

}

}

}

//funcție care împarte termenul liber la a[i][i]

void term_liber(long double A[N][N],long double b_sigma[N],

long double a[N],long double sigma,int n)

{

for(int i=0;i<n;i++)

{

b_sigma[i]=sigma*(a[i]/A[i][i]);

}

}

//funcție care verifică dacă o matrice este simetrică

int simetrica(long double A[N][N],int n)

{

for(int j,i=0;i<n;i++)

{

for(j=0;j<n;j++)

{

if(A[i][j]!=A[j][i])

{

return 0;

}

}

}

return 1;

}

//funcție care returnează valoarea lui q

long double ret_q(long double s_l[N],long double B_sigma[N][N],int m)

{

for(int j,i=0;i<m;i++)

{

s_l[i]=0;

for(j=0;j<m;j++)

{

s_l[i]=s_l[i]+fabs(B_sigma[i][j]);

}

}

long double max=s_l[0];

for(i=0;i<m;i++)

{

if(s_l[i]>max)

{

max=s_l[i];

}

}

return max;

}

//funcție pentru aducerea maximului de pe coloana k pe linia k

int max_kk(long double A[N][N],int k,int n)

{

int ns=0;

//deterinam maximul pe coloana k

long double max=fabs(A[k][k]);

int imax=k;//presupunem ca max este pe linia k

for(int i=k+1;i<n;i++)

{

if(fabs(A[i][k])>max)

{

imax=i;

max=fabs(A[i][k]);

}

}

long double t;

if(imax>k)

{

ns++;

for(int j=k;j<=n;j++)

{

t=A[k][j];

A[k][j]=A[imax][j];

A[imax][j]=t;

}

}

return ns;

}

//funcție care împarte linia k la a[k][k]

long double impart(long double A[N][N],int n,int k)

{

long double m=A[k][k];

for(int j=k;j<=n;j++)

{

A[k][j]=A[k][j]/m;

}

return m;

}

//funcție care face zero sub a[k][k]

void zero(long double A[N][N],int n,int k)

{

for(int j,i=k+1;i<n;i++)

{

for(j=k+1;j<=n;j++)

{

A[i][j]=A[i][j]-A[k][j]*A[i][k];

}

}

for(i=k+1;i<n;i++)

{

A[i][k]=0;

}

}

//funcție care calculează determinantul cu ajutorul metodei lui Gauss

long double det_gauss(long double A[N][N],int n)

{

long double d=1;

for(int k=0;k<n;k++)

{

int ns=max_kk(A,k,n);

if(A[k][k]==0)

{

cout<<"\n Determinantul sistemului este nul ";

return 0;

}

d=d*impart(A,n,k)*pow(-1,ns);

zero(A,n,k);

}

return d;

}

//funcție pentru afișarea unei matrici

void afis_mat(char *nume,long double A[N][N],int n)

{

cout<<"\n Matricea "<<nume<<" este ";

for(int j,i=0;i<n;i++)

{

cout<<"\n";

for(j=0;j<n;j++)

{

cout<<setw(15)<<A[i][j];

}

}

}

//funcție care calculează matricea

void matr_sigma_D_1A(long double sigma,long double D_1A[N][N],

int m)

{

for(int j,i=0;i<m;i++)

{

for(j=0;j<m;j++)

{

D_1A[i][j]=sigma*D_1A[i][j];

}

}

}

//funcție care calculează iterațiile

void iteratie(long double B_sigma[N][N],long double b_sigma[N],

long double xn[N],long double xn_1[N],int m)

{

long double s;

for(int j,i=0;i<m;i++)

{

s=0;

for(j=0;j<m;j++)

{

s=s+B_sigma[i][j]*xn_1[j];

}

xn[i]=s+b_sigma[i];

}

}

//funcție care calculează

long double norma(long double xn_1[N],long double xn[N],int m)

{

long double norm[N];

for(int i=0;i<m;i++)

{

norm[i]=fabs(xn_1[i]-xn[i]);

}

long double max=norm[0];

for(i=1;i<m;i++)

{

if(norm[i]>max)

{

max=norm[i];

}

}

return max;

}

void main(void)

{

long double lambda,suma_linie[N],x[N];

int m=3;

long double A[N][N]={{10,-2,1},

{-2,10,1},

{1,1,10}

};

long double a[N]={108,-6,-16,5};

long double x0[N]={0,0,0,0};

long double v_det,D[N][N],Det[N][N],I[N][N],D_1A[N][N],

B_sigma[N][N],b_sigma[N];

long double sigma=0.9;

afis_mat("A",A,m);

afis_v("a",a,m);

int i,j;

int poz_def=1;

for(int k=1;k<=m;k++)

{

for(i=0;i<k;i++)

{

for(j=0;j<k;j++)

{

Det[i][j]=A[i][j];

}

}

afis_det(Det,k);

v_det=det_gauss(Det,k);

cout<<"\n Valoarea determinantului este:"<<v_det;

if(v_det<=0)

{

poz_def=0;

}

}

if(simetrica(A,m)==0)

{

cout<<"\n A nu este simetrică";

}

else

{

cout<<"\n A este simetrică";

if(poz_def==0)

{

cout<<"\n A nu este pozitiv definită";

}

else

{

cout<<"\n A este pozitiv definită";

cout<<"\n m="<<m;

matr_D(D,A,m);

afis_mat("D",D,m);

matr_D_1(D,m);

afis_mat("D_1",D,m);

produs_m(D,A,D_1A,m,m,m);

afis_mat("D_1A",D_1A,m);

lambda=lambda1(suma_linie,D_1A,m);

afis_v("suma_linie",suma_linie,m);

cout<<"\n Lambda="<<lambda;

cout<<"\n sigma="<<sigma;

matr_sigma_D_1A(sigma,D_1A,m);

afis_mat("sigma*D_1A",D_1A,m);

matr_I(I,m);

afis_mat("I",I,m);

diferenta(I,D_1A,B_sigma,m);

afis_mat("B_sigma",B_sigma,m);

term_liber(A,b_sigma,a,sigma,m);

afis_v("b_sigma",b_sigma,m);

long double s_l[N],q=ret_q(s_l,B_sigma,m);

afis_v("s_l",s_l,m);

cout<<"\n q="<<q;

long double epsilon=0.001;

cout<<"\n epsilon="<<epsilon;

if(q<1)

{

for(int k=0; ;k++)

{

iteratie(B_sigma,b_sigma,x,x0,m);

cout<<"\n Iteratia "<<k+1<<"\n x={";

for(int i=0;i<m;i++)

{

cout<< setw(10)<<x[i];

}

cout<<"}";

if(q/(1-q)*norma(x0,x,m)<epsilon)

{

break;

}

for(i=0;i<m;i++)

{

x0[i]=x[i];

}

}

}

else

{

cout<<"\n Metoda nu converge deoarece "<<q<<">=1" ;

}

}

}

cout<<"\n";

}

Metoda relaxării succesive

#include<iostream.h>

#include<iomanip.h>

#include<math.h>

#define N 10

//funcție pentru afișarea unei matrici

void afis_m(char *nume,long double A[N][N],int m)

{

cout<<"\n Matricea "<<nume<<" este";

for(int j,i=0;i<m;i++)

{

cout<<"\n";

for(j=0;j<m;j++)

{

cout<<setw(15)<<A[i][j];

}

}

}

//funcție pentru afișarea unui vector

void afis_v(char *nume,long double a[N],int m)

{

cout<<"\n Vectorul "<<nume<<" este";

for(int i=0;i<m;i++)

{

cout<<"\n"<<setw(15)<<a[i];

}

}

//funcție pentru afișarea matricei D

void m_D(long double A[N][N],long double D[N][N],int m)

{

for(int j,i=0;i<m;i++)

{

for(j=0;j<m;j++)

{

if(i==j)

{

D[i][j]=A[i][j];

}

else

{

D[i][j]=0;

}

}

}

}

//funcție care calculează matricea

void m_D_1(long double D[N][N],int m)

{

for(int j,i=0;i<m;i++)

{

for(j=0;j<m;j++)

{

if(i==j)

{

D[i][j]=1/D[i][j];

}

else

{

D[i][j]=0;

}

}

}

}

//funcție care calculează matricea L

void m_L(long double A[N][N],long double L[N][N],int m)

{

for(int j,i=0;i<m;i++)

{

for(j=0;j<m;j++)

{

if(i>j)

{

L[i][j]=A[i][j];

}

else

{

L[i][j]=0;

}

}

}

}

//funcție care calculează diferența dintre două matrici

void diferenta(long double A[N][N],long double B[N][N],

long double C[N][N],int m)

{

for(int j,i=0;i<m;i++)

{

for(j=0;j<m;j++)

{

C[i][j]=A[i][j]-B[i][j];

}

}

}

//funcție care calculează matricea R

void m_R(long double A[N][N],long double D[N][N],

long double L[N][N],long double R[N][N],int m)

{

long double R1[N][N];

diferenta(A,D,R1,m);

diferenta(R1,L,R,m);

}

//funcție care calculează iterațiile

void iteratii(long double xn[N],long double xn_1[N],

long double A[N][N],long double sigma,long double a[N],int m)

{

long double s1,s2;

for(int j,i=0;i<m;i++)

{

s1=0;

for(j=0;j<=i-1;j++)

{

s1=s1+A[i][j]*xn[j];

}

s2=0;

for(j=i+1;j<m;j++)

{

s2=s2+A[i][j]*xn_1[j];

}

xn[i]=-(sigma/A[i][i])*s1+(1-sigma)*xn_1[i]-(sigma/A[i][i])*s2+

+(sigma/A[i][i])*a[i];

}

}

//funcție care verifică dacă o funcție este simetrică

int simetrica(long double A[N][N],int n)

{

for(int j,i=0;i<n;i++)

{

for(j=0;j<n;j++)

{

if(A[i][j]!=A[j][i])

{

return 0;

}

}

}

return 1;

}

//funcție pentru aducerea maximului de pe coloana k pe linia k

int max_kk(long double A[N][N],int k,int n)

{

int ns=0;

//determinăm maximul pe coloana k

long double max=fabs(A[k][k]);

int imax=k;//pp ca max este pe linia k

for(int i=k+1;i<n;i++)

{

if(fabs(A[i][k])>max)

{

imax=i;

max=fabs(A[i][k]);

}

}

//schimbarea liniilor k cu imax

long double t;

if(imax>k)

{

ns++;

for(int j=k;j<=n;j++)

{

t=A[k][j];

A[k][j]=A[imax][j];

A[imax][j]=t;

}

}

return ns;

}

//funcție care împarte linia k la a[k][k]

long double impart(long double A[N][N],int n,int k)

{

long double m=A[k][k];

for(int j=k;j<=n;j++)

{

A[k][j]=A[k][j]/m;

}

return m;

}

//funcție care face zero sub a[k][k]

void zero(long double A[N][N],int n,int k)

{

for(int j,i=k+1;i<n;i++)

{

for(j=k+1;j<=n;j++)

{

A[i][j]=A[i][j]-A[k][j]*A[i][k];

}

}

for(i=k+1;i<n;i++)

{

A[i][k]=0;

}

}

//funcție care calculează determinantul cu ajutorul metodei lui Gauss

long double det_gauss(long double A[N][N],int n)

{

long double d=1;

for(int k=0;k<n;k++)

{

int ns=max_kk(A,k,n);

if(A[k][k]==0)

{

cout<<"\n Determinantul sistemului este nul ";

return 0;

}

d=d*impart(A,n,k)*pow(-1,ns);

zero(A,n,k);

}

return d;

}

//funcție care calculează

long double norma(long double xn_1[N],long double xn[N],int m)

{

long double norm[N];

for(int i=0;i<m;i++)

{

norm[i]=fabs(xn_1[i]-xn[i]);

}

long double max=norm[0];

for(i=1;i<m;i++)

{

if(norm[i]>max)

{

max=norm[i];

}

}

return max;

}

void main(void)

{

int m=4;

long double D[N][N],L[N][N],R[N][N],Det[N][N];

long double A[N][N]={{2.12,0.42,1.34,0.88},

{0.42,3.95,1.87,0.43},

{1.34,1.87,2.98,0.46},

{0.88,0.43,0.46,4.44}};

long double a[N]={11.172,0.115,9.009,9.349};

long double x0[N]={0,0,0,0},x[N];

long double v_det;

afis_m("A",A,m);

afis_v("a",a,m);

int poz_def=1;

int i,j;

for(int k=1;k<=m;k++)

{

for(i=0;i<k;i++)

{

for(j=0;j<k;j++)

{

Det[i][j]=A[i][j];

}

}

afis_m("Det",Det,k);

v_det=det_gauss(Det,k);

cout<<"\n Valoarea determinantului este:"<<v_det;

if(v_det<=0)

{

poz_def=0;

}

}

if(simetrica(A,m)==0)

{

cout<<"\n A nu este simetrică";

}

else

{

cout<<"\n A este simetrică";

if(poz_def==0)

{

cout<<"\n A nu este pozitiv definită";

}

else

{

cout<<"\n A este pozitiv definită";

m_D(A,D,m);

afis_m("D",D,m);

m_L(A,L,m);

afis_m("L",L,m);

m_R(A,D,L,R,m);

afis_m("R",R,m);

long double epsilon=0.00001;

cout<<"\n epsilon="<<epsilon;

long double sigma=0.9;

cout<<"\n sigma="<<sigma;

for(k=1; ;k++)

{

iteratii(x,x0,A,sigma,a,m);

cout<<"\n Iteratia "<<k<<"\n x={";

for(i=0;i<m;i++)

{

cout<< setw(15)<<x[i];

}

cout<<"}";

if(norma(x0,x,m)<epsilon)

{

break;

}

for(i=0;i<m;i++)

{

x0[i]=x[i];

}

}

}

}

cout<<"\n";

}

Pseudoinversa unei aplicații liniare

#include<iomanip.h>

#include<iostream.h>

#include<math.h>

#include<stdio.h>

#include<string.h>

#define N 20

//funcție pentru afișarea unei matrici

void afis_m(char *nume,long double T[N][N],int l,int c)

{

cout<<"\n Matricea "<<nume<<" este";

for(int j,i=0;i<l;i++)

{

cout<<"\n";

for(j=0;j<c;j++)

{

cout<<setw(15)<<setiosflags(ios::fixed)<<T[i][j];

}

}

}

//funcție care calculează transpuusa unei matrici

void transpus(long double T[N][N],long double Tran[N][N],int l,int c)

{

for(int j,i=0;i<c;i++)

{

for(j=0;j<l;j++)

{

Tran[i][j]=T[j][i];

}

}

}

//funcție care afișează matricea inițială

void H_0(long double H0[N][N],int l,int c)

{

for(int j,i=0;i<c;i++)

{

for(j=0;j<l;j++)

{

H0[i][j]=0;

}

}

}

//funcție care afișează vectorii y

void afis_y(long double Tran[N][N],int l,int c)

{

for(int i,k=1;k<=l;k++)

{

cout<<"\n y"<<k;

for(i=0;i<c;i++)

{

cout<<"\n"<<Tran[i][k-1];

}

}

}

//funcție care calculează

void calc_Wk(long double yk[N],long double Hk_1[N][N],

long double T[N][N],long double wk[N],int l,int c)

{

long double P[N][N];

for(int k,j,i=0;i<c;i++)

{

for(k=0;k<c;k++)

{

P[i][k]=0;

for(j=0;j<l;j++)

{

P[i][k]=P[i][k]+Hk_1[i][j]*T[j][k];

}

}

}

afis_m("P",P,c,c);

long double V[N];

for(i=0;i<c;i++)

{

V[i]=0;

for(j=0;j<c;j++)

{

V[i]=V[i]+P[i][j]*yk[j];

}

}

cout<<"\n V=";

for(i=0;i<c;i++)

{

cout<<setw(14)<<V[i];

}

for(i=0;i<c;i++)

{

wk[i]=yk[i]-V[i];

}

cout<<"\n w=";

for(i=0;i<c;i++)

{

cout<<setw(14)<<wk[i];

}

}

//funcție care calculează

void calc_Hk(long double Hk[N][N],long double Hk_1[N][N],

long double wk[N],long double T[N][N],int l,int c)

{

long double Twk[N];

for(int j,i=0;i<l;i++)

{

Twk[i]=0;

for(j=0;j<c;j++)

{

Twk[i]=Twk[i]+T[i][j]*wk[j];

}

}

cout<<"\n Twk=";

for(i=0;i<l;i++)

{

cout<<setw(14)<<Twk[i];

}

long double normaTwk=0;

for(i=0;i<l;i++)

{

normaTwk=normaTwk+pow(Twk[i],2);

}

cout<<"\n norma="<<normaTwk;

//wk*traspTwk

long double wktranTwk[N][N];

int k;

for(k,i=0;i<c;i++)

{

for(k=0;k<l;k++)

{

wktranTwk[i][k]=wk[i]*Twk[k];

}

}

afis_m("wktranTwk",wktranTwk,c,l);

long double wktranTwkimpart[N][N];

for(i=0;i<c;i++)

{

for(j=0;j<l;j++)

{

wktranTwkimpart[i][j]=wktranTwk[i][j]/normaTwk;

}

}

for(i=0;i<c;i++)

{

for(j=0;j<l;j++)

{

Hk[i][j]=Hk_1[i][j]+wktranTwkimpart[i][j];

}

}

}

//funcție pentru verificarea condiției 0

int alg_cont(long double yk[N],long double Hk_1[N][N],

long double T[N][N],int l,int c)

{

long double Hk_1T[N][N];

for(int k,j,i=0;i<c;i++)

{

for(k=0;k<c;k++)

{

Hk_1T[i][k]=0;

for(j=0;j<l;j++)

{

Hk_1T[i][k]=Hk_1T[i][k]+Hk_1[i][j]*T[j][k];

}

}

}

long double Hk_1Ty[N];

for(i=0;i<c;i++)

{

Hk_1Ty[i]=0;

for(j=0;j<c;j++)

{

Hk_1Ty[i]=Hk_1Ty[i]+Hk_1T[i][j]*yk[j];

}

}

for(i=0;i<c;i++)

{

if(fabs(yk[i]-Hk_1Ty[i])>0.000000001)

{

return 1;

}

}

return 0;

}

void main(void)

{

long double yk[N],wk[N];

int l=3,c=4;

long double Tran[N][N];// c linii si l coloane

long double H1[N][N];

long double H0[N][N];

long double T[N][N]={{1,-1,2,0},{-1,2,-3,1},{0,1,-1,1}};

afis_m("T",T,l,c);

transpus(T,Tran,l,c);

afis_m("Tran",Tran,c,l);

H_0(H0,l,c);

afis_m("H0",H0,c,l);

afis_y(Tran,l,c);

for(int j,i,k=1;;k++)

{

//construim yk

for( i=0;i<c;i++)

{

yk[i]=Tran[i][k-1];

}

if(!alg_cont(yk,H0,T,l,c))

{

cout<<"\n Algoritmul se oprește ";

break;

}

calc_Wk(yk,H0,T,wk,l,c);

calc_Hk(H1,H0,wk,T,l,c);

char nume[5];

sprintf(nume,"H%d",k);

afis_m(nume,H1,c,l);

for(i=0;i<c;i++)

{

for(j=0;j<l;j++)

{

H0[i][j]=H1[i][j];

}

}

}

cout<<"\n";

}

1. Cristescu R., Analiză funcțională, Editura Didactică și Pedagogică București, 1965.

2. Demidovici B., Maron I., Elémente de calcul numérique, Éditions Mir Moscou, 1973.

3. Grigore Gh., Lecții de analiză numerică, Editura Universității din București, 1990.

4. Marinescu Gh., Grigore Gh., Jambor C., Mazilu P., Rizzoli I., Ștefan C., Probleme de analiză numerică, Editura Didactică și Pedagogică București, 1978.

5. V. Voievodine, Principes numériques d ΄algèbre linéaire, Edition Mir Moscv 1980.