Algebra Matriceala

CUPRINS

Capitolul 1. – Elemente de algebră matriceală

1.1. Prezentarea problemelor care utilizează algebra matriceală

1.2.1. Definirea inversei

1.2. Inversarea matricelor

1.2.2. Metode de calcul a inversei

1.3. Factorizarea matricelor

1.3.1. Considerații preliminare

1.3.2. Factorizarea LR

Capitolul 2. – Analiza circuitelor rezistive liniare

2.1. Metode și algoritmi de calcul

2.1.1.Metoda teoremelor lui Kirchhoff

2.1.2.Metoda curenților de buclă

2.1.3.Metoda potențialelor nodurilor

Capitolul 3. – Soluționarea numerică a sistemelor de ecuații liniare

3.1. Prezentarea problemei

3.2. Metode directe de soluționare

3.2.1. Rezolvarea unui sistem n ecuații liniare cu n necunoscute, folosind regula lui Cramer

3.2.2. Metoda inversării matriceale

3.2.3. Metoda de eliminare a lui Gauss

3.2.4. Metode bazate pe factorizarea LR a matricei coeficienților

3.3. Metode iterative de soluționare

3.3.1. Metoda aproximațiilor succesive, versiunea Jacobi

3.3.2. Metoda aproximațiilor succesive, versiunea Gauss – Seidel

3.3.3. Metoda gradientului (celei mai abrupte coborâri)

Bibliografie

Capitolul 1. – Elemente de algebră matriceală

1.1. Prezentarea problemelor care utilizează algebra matriceală

Utilizarea algebrei matriceale s-a extins considerabil în ultimii 25 – 30 de ani, paralel cu evoluția tehnicii de calcul și cu folosirea ei din ce în ce mai accentuată la soluționarea problemelor științifice și tehnice. Introducerea notației matriceale conduce la o formulare simplă și concisă a unor aplicații deosebit de complexe, mai ales în situațiile în care modelul matematic conține sisteme de ecuații liniare de dimensiuni mari. Exemplele sunt numeroase, atât din domeniul tehnicii în general, cât și din cel al aplicațiilor electrotehnice în particular: modelarea numerică a funcționării mașinilor electrice, hidraulice sau termice; analiza, optimizarea și estimarea regimurilor de funcționare a sistemelor electroenergetice; modelarea numerică a funcționării aparatelor și echipamentelor electrice. De altfel, pentru orice problemă se poate elabora în primă instanță un model liniar sau, dacă modelul este neliniar, el poate fi liniarizat în primă aproximație, o singură dată sau la fiecare pas al unui proces iterativ de soluționare.

În acest context, prezentul capitol se referă la inversarea matricelor și la factorizarea matricelor. Ele vor fi privite atât ca operații de sine stătătoare, cât și ca elemente pregătitoare pentru capitolele următoare, referitoare la soluționarea numerică a sistemelor de ecuații liniare și la determinarea valorilor proprii și a vectorilor proprii [7]. Se vor prezenta atât metodele generale de inversare și de factorizare, cât și elementele specifice de aplicare pentru cazul matricelor simetrice, al celor triunghiulare sau al celor de mari dimensiuni.

Înainte de a trece la metodele de inversare și la cele de factorizare, se reamintește terminologia legată de algebra matriceală, utilizată în acest capitol și în cele următoare.

În spațiul vectorial n – dimensional Cn peste corpul numerelor complexe C, se definește vectorul coloană x ca fiind un ansamblu ordonat de n numere complexe, x1, x2, x3, …, xn , numite componentele sau coordonatele vectorului. Se definește de asemenea transpusa vectorului coloană x ca fiind vectorul linie xt, respectiv conjugata sa, ca fiind vectorul coloană . Definirea transpusei permite scrierea vectorului x și sub forma (1.4). Vectorul coloană care are toate componentele nule se notează cu 0.

(1.1)

(1.2)

Xt = [ x1 x2……xn ] (1.3)

(1.4)

Vectorul x care are modulul egal cu unitatea se numește vector normalizat sau vector unitate. Pentru doi vectori x și y se definește produsul scalar:

(x, y)=xy (1.5)

ca fiind suma produselor componentelor celor doi vectori. Dacă (x, y) = 0, atunci cei doi vectori se numesc ortogonali. Mulțimea finită de n vectori x1, x2, x3, …, xn este o mulțime ortogonală, dacă, pentru orice i ≠ j, (xi.xj) = 0 , respectiv ortonormală dacă, în plus, fiecare vector are modulul egal cu unitate.

Un spațiu vectorial este n – dimensional dacă și numai dacă spațiul respectiv conține n vectori liniar independenți și orice mulțime de n + 1 vectori din acel spațiu este liniar dependentă. Aceasta înseamnă ca dacă există o asemenea mulțime de vectori independenți, fiecare vector al spațiului poate fi exprimat printr-o combinație liniară a lor [11]. Fiecare mulțime de vectori independenți din acel spațiu alcătuiește o bază a spațiului respectiv.

O matrice A de dimensiune m x n este definită ca fiind un tablou dreptunghiular, având elemente reale sau complexe, dispuse în mod sistematic pe m linii și n coloane.

A = (1.6)

Transpusa unei matrice A, de dimensiune m x n, este matricea At, de dimensiune n x m, obținută prin schimbarea liniilor cu coloanele. Transpusa sumei și transpusa produsului a două matrice, A și B, satisfac relațiile (1.7) și (1.8).

(A + B)t = At + Bt (1.7)

(A B)t = Bt At (1.8)

Conjugata unei matrice A este matricea , obținută prin înlocuirea elementelor matricei originale cu conjugatele lor.

Dacă n = 1, atunci matricea este de tip coloană (dimensiunea m x 1), iar dacă m = 1, atunci este de tip linie (1xn). Dacă m = n, atunci matricea A este pătrată (de ordinul n). O matrice pătrată A este simetrică dacă At = A, adică aij = aji, i, j = 1, 2, … , n.

O matrice A pătrată este superior triunghiulară dacă pentru j i, aij = 0. Ea este inferior triunghiulară dacă pentru i j, aij = 0. Dacă elementele unei matrice reale pătrate A sunt nule, mai puțin cele de pe diagonala principală (aij =1, i = j), atunci ea este o matrice diagonală [3]. Dacă la o matrice diagonală elementele nenule au valoarea unitară (aij = 1, i = j), atunci ea este matrice unitate, fiind notată cu I.

O matrice A de forma (1.9) se numește matrice tridiagonală, iar cea de forma (1.10) este cunoscută sub numele de matrice Hessenberg.

A = (1.9)

A = (1.10)

O matrice A este bloc diagonală dacă are forma din relația (1.11), unde "blocurile" A1, A2, A3 , …, Ap sunt matrice pătrate, nu neapărat de același ordin.

A = (1.11)

O matrice A este ortogonală dacă A-1 = At , respectiv A At = I. Dacă a1 , a2 , a3 , …, an sunt vectorii coloană corespunzători coloanelor lui A , atunci ei sunt ortonormali.

O matrice A se spune că este hermitiană dacă este satisfăcută condiția A = ()t. O matrice reală simetrică îndeplinește condiția de a fi hermitiană [8].

O matrice reală pătrată A este singulară dacă determinantul său det A= 0, în caz contrar ea este nesingulară. Ordinul lui A fiind n, se spune că ea este pozitiv definită dacă pentru x n este satisfăcută condiția (1.12):

(A x, x)t = xt A x > 0

1.2. Inversarea matricelor

1.2.1. Definirea inversei

Se consideră matricea pătrată reală nesingulară A de ordinul n din relația:

A = (1.13)

Matricea inversă A-1 a matricei A se definește ca fiind acea matrice care satisface relația:

A A-1 = A-1 A = I (1.14)

unde: I este matricea unitate de ordinul n.

Matricea inversă A-1 se poate exprima sub forma (1.15), unde A+ este matricea adjunctă, iar detA determinantul lui A. Matricea adjunctă este transpusa matricei obținute prin înlocuirea elementelor lui A cu cofactorii lor. Cofactorul lui aij este minorul acelui element, înmulțit cu (-1)i+j. Minorul lui aij este determinantul de ordinul n – 1 obținut prin înlăturarea liniei i și a coloanei j a matricei A.

A-1= (1.15)

Relația (1.15) sugerează și o modalitate de calcul a inversei, apreciază ca nepractică, din cauza volumului mare de calcule și a dificultăților de transpunere pe calculator [14]. Se vor prezenta metode mai eficiente de obținere a inversei (și a valorii determinantului).

În final se remarcă unele proprietăți ale inversei și ale operației de inversare, ca fiind de interes în cele ce urmează:

între determinantul matricei A și cel al matricei A-1 există relația (1.16);

det(A-1) = (detA)-1 (1.16)

inversa produsului a două matrice, A și B , se poate determina cu relația (1.17);

(A B)-1 = B-1 A-1 (1.17)

transpusa matricei inverse A-1 este egală cu inversa transpusei matricei A.

(A-1)t = (At)-1 (1.18)

În continuare se prezintă metodele practice de determinare a inversei și unele particularizări de interes în domeniul aplicațiilor electrotehnice, referitoare la matricole simetrice, la cele triunghiulare și la cele de mari dimensiuni, partiționate. Se menționează că metodele abordate se referă la matricele reale, dar ele pot fi aplicate și pentru matricele complexe, modificările și adaptările necesare fiind simple.

1.2.2. Metode de calcul a inversei

Între metodele de calcul a inversei, descrise în literatura de specialitate, se prezintă cea a reducerii la matricea unitate, cunoscută sub numele de metoda eliminării, versiunea Gauss – Jordan (diagonalizării) [5]. Alegerea este justificată atât de calitățile algoritmului, cât și de utilizarea ei pe scară largă în aplicațiile din domeniul aplicațiilor electrotehnice. În plus, ea permite și obținerea valorii determinantului, fără a efectua calcule suplimentare.

Metoda eliminării are la bază următoarea teoremă din algebra matriceală: dacă o matrice nesingulară A poate fi redusă la matricea unitate l prin înmulțire la stânga cu un șir de matrice, atunci matricea inversă A-1 se poale calcula prin premultiplicarea lui I cu același șir de matrice.

Calculul pornește de la matricea A0 = A , dată de relația (1.13), având elementele ai0j = aij, i = 1, 2, …, n , j = 1, 2, …., n (indicele superior se referă la numărul pasului). El se desfășoară în n etape, fiecare dintre ele implicând efectuarea unor operații pe toate liniile matricei A, prezentate intuitiv în cele ce urmează [10].

După etapa 1 se urmărește obținerea lui A-1, la care elementul diagonal a111 = 1, iar restul elementelor de pe coloana 1 sunt nule (ai11 = 0, i = 1, 2, …, n, i ≠ 1).

A1= (1.19)

Elementele lui A-1 , în ipoteza a110 0, rezultă cu relatiile (1.20) – (1.21):

a1j1 = , j=1, 2, …, n (1.20)

aij1 = aij0 – ai10 a1j1, i = 2 , 3, …, n (1.21)

După pasul 2 se urmărește obținerea matricei A2, caracterizată prin faptul că elementul diagonal a222 = 1, iar restul elementelor de pe coloana 2 sunt nule (ai22 = 0 , i = 1, 2, …, n, i ≠ 2). Evident, elementele coloanei 1 trebuie să rămână nemodificate.

A2= (1.22)

Elemcntele lui A2, în ipoteza a221 0, rezultă cu relațiile (1.23) – (1.24).

a2j2=, j=1, 2, …, n (1.23)

aij2 = aij1– ai21 a2j2, i = 1, 3, 4, …, n (1.24)

Generalizand procedura de la pașii anteriori, rezultă că după pasul k de eliminare se urmărește obținerea matricei Ak, caracterizată prin faptul că elementul diagonal akkk = 1, iar restul elementelor de pe coloana k sunt nule (aikk = 0, i = 1, 2, …, n, i ≠ k). Evident, elementele coloanelor 1, 2, …, k – 1 trebuie să rămână nemodificate.

Atiile (1.20) – (1.21):

a1j1 = , j=1, 2, …, n (1.20)

aij1 = aij0 – ai10 a1j1, i = 2 , 3, …, n (1.21)

După pasul 2 se urmărește obținerea matricei A2, caracterizată prin faptul că elementul diagonal a222 = 1, iar restul elementelor de pe coloana 2 sunt nule (ai22 = 0 , i = 1, 2, …, n, i ≠ 2). Evident, elementele coloanei 1 trebuie să rămână nemodificate.

A2= (1.22)

Elemcntele lui A2, în ipoteza a221 0, rezultă cu relațiile (1.23) – (1.24).

a2j2=, j=1, 2, …, n (1.23)

aij2 = aij1– ai21 a2j2, i = 1, 3, 4, …, n (1.24)

Generalizand procedura de la pașii anteriori, rezultă că după pasul k de eliminare se urmărește obținerea matricei Ak, caracterizată prin faptul că elementul diagonal akkk = 1, iar restul elementelor de pe coloana k sunt nule (aikk = 0, i = 1, 2, …, n, i ≠ k). Evident, elementele coloanelor 1, 2, …, k – 1 trebuie să rămână nemodificate.

A= (1.25)

Elementele lui Ak, în ipoteza akkk-1 ≠ 0, rezultă cu relațiile (1.26) – (1.27):

akjk = , j = 1, 2, …, n (1.26)

aijk = aijk-1 – aikk-1 akjk, i = 1, 2, …, n (1.27)

După ultimul pas de eliminare, k = n, se ajunge la matricea unitate:

An = I (1.28)

Pornind de la matricea D0 = 1 și efectuând același șir de operații, rezultă matricele D1, D2, …, Dn, ultima dintre ele fiind chiar A-1 (indicele superior de la matricele D se referă la etapa de calcul) [1].

După etapa 1 se obține forma intermediară D1:

D = (1.29)

ale cărei elemente, în ipoteza a110 ≠ 0, se determină cu relațiile (1.30) – (1.31).

d1j1 = , j = 1, 2, …, n (1.30)

dij1 = dij0 – ai10 d1j1, i = 2, 3, …, n; j = 1, 2, …, n (1.31)

După pasul 2 rezultă matricea D2, ale cărei elemente, în ipoteza a212 ≠ 0, se calculează cu relațiile (1.33) – (1.34) [10].

D = (1.32)

d2j2 = , j = 1, 2, …, n (1.33)

dij2 = dij1 – ai21 d2j2, i = 1, 3, …, n; j = 1, 2, …, n (1.34)

Generalizând procedura de la pașii anteriori, după pasul k se obține Dk , ale cărei elemente, în ipoteza akkk-1 ≠ 0, rezultă cu relațiile (1.36) – (1.37).

Dk = (1.35)

dkjk = , j = 1, 2, …, n (1.36)

dijk = dijk-1 – aikk-1 dkjk, i = 1, 2, …, n; i ≠ k; j = 1, 2 (1.37)

După ultimul pas, k = n, se ajunge la rezultatul urmărit:

A-1 = Dn (1.38)

Analiza celor două procese de calcul cu câte n pași (reducerea matricei A la matricea unitate I, respectiv obținerea matricei inverse A-1 pornind de la I) evidențiază următoarea concluzie: la un pas oarecare k se utilizează în exclusivitate informațiile de la pasul k – 1 [4]. În consecință, cele două procese de calcul se pot derula concomitent (în paralel), rezultând următorul algoritm practic:

se efectuează inițializările necesare

A0 = A (1.39)

D0 = I (1.40)

la un pas oarecare k, k = 1, 2, …, n, se calculează elementele matricelor Ak și Dk

akjk = , j = k + 1, k + 2, …, n (1.41)

dkjk = , j = 1, 2, …, k (1.42)

aiik = 1, i = 1, 2, …, k (1.43)

diik = 1, i = k + 1, k + 2, …, n (1.44)

aijk = aijk-1 – aikk-1 akjk; i = 1, 2, …, n, i ≠ k, j = k + 1, k + 2, …, n (1.45)

dijk = dijk-1 – aikk-1 dkjk; i = 1, 2, …, n, i ≠ k, j = 1, 2, …, k (1.46)

aijk = 0; i = 1, 2, …, n, i ≠ k, j = 1, 2, …, k (1.47)

dijk = 0; i = 1, 2, …, n, i ≠ j, j = k + 1, k + 2, …, n (1.48)

în final se obține matricea unitate, respectiv inversă, cu relațiile (1.28) – (1.38).

Concomitent cu reducerea matricei A la I, se poate calcula simplu și valoarea determinantului detA (știind că împărțirea unei limi cu un scalar determină reducerea în aceeași măsură a valorii determinantului, iar restul de combinații liniare nu o afectează).

(1.49)

În final rezultă valoarea determinantului detA al matricei A:

(1.50)

Expresia din relația (1.50) oferind și interpretarea corectă a situațiilor când un element diagonal curent are valoare nulă.

Pe baza celor prezentate se realizează codul sursă pentru această metodă cu observația că s-a utilizat câte un singur tablou pentru matricele A1, A2, …, An, respectiv pentru matricele D1, D2, …, Dn (tablourile A și D, cu pierderea informațiilor originale despre matricea A) [9]. Se prezintă și o fereastra din programul Delphi cu un exemplu numeric rezolvat.

Program Inv_Mat;

var i,j:integer; {variabile ajutatoare}

k:integer; {contorul pasilor procesului de inversare}

n:integer; {ordinul matricei}

detA:real; {valoarea determinantului matricei A}

eps:real; {pragul sub care un pivot se considera nul}

t:real; {variabila ajutatoare}

A:array[1..8,1..8]of real; {matricea A}

D:array[1..8,1..8]of real; {matricea A^-1}

{procedura pentru afisarea matricei A}

procedure VizMatA;

begin

writeln('Matricea A');

for i:=1 to n do begin

for j:=1 to n do write(A[i,j]:10:5);

if n<8 then writeln;

end;

writeln;

end;

{procedura pentru afisarea matricei A^-1}

procedure VizMatD;

begin

writeln('Matricea A^-1');

for i:=1 to n do begin

for j:=1 to n do write(D[i,j]:10:5);

if n<8 then writeln;

end;

writeln;

end;

begin

writeln('Inversarea matriceala prin reducerea la matricea unitate');

writeln('********************************************************');

{citire A}

writeln;

write('Ordinul matricei A (maxim 8): n=');

readln(n);

writeln;

for i:=1 to n do begin

writeln('Linia ',i);

for j:=1 to n do begin

write('A[',i,',',j,']=');

readln(a[i,j]);

end;

end;

writeln;

{Citire prag 0}

write('Pragul sub care un pivot se considera nul: eps=');

readln(eps);

writeln;

{Afisare A}

VizMatA;

{Init D cu I}

for i:=1 to n do begin

for j:=1 to n do d[i,j]:=0;

d[i,i]:=1;

end;

{Inversarea in n pasi, k=1..n}

detA:=1;

for k:=1 to n do begin

{test pivot}

if (abs(a[k,k])<eps then begin

write('Pivot nul la pasul ',k);

readln;

exit;

end;

detA:=detA*a[k,k];

{elementele liniei k a matricei A}

t:=1/A[k,k];

a[k,k]:=1;

if k<>n then for j:=k+1 to n do a[k,j]:=a[k,j]*t;

{elementele liniei k a matricei D}

for j:=1 to k do d[k,j]:=d[k,j]*t;

{elementele liniilor diferite de k ale matricelor A si D}

for i:=1 to n do begin

if i<>k then begin

t:=a[i,k];

a[i,k]:=0;

if k<>n then

for j:=k+1 to n do

a[i,j]:=a[i,j]-t*a[k,j];

for j:=1 to k do

d[i,j]:=d[i,j]-t*d[k,j];

end;

end;

{Vizualizarea matricelor A si D la pasul k}

writeln('Pasul=',k);

writeln;

VizMatA;

VizMatD;

writeln;

writeln('Apasa Enter');

readln;

end;

{Rezultate:}

writeln;

writeln('Rezultate finale:');

writeln('––––––');

writeln(' |A|=',detA);

writeln;

VizMatA;

VizMatD;

writeln('Apasa Enter');

readln;

end.

Figura 1.1. – Calculul inversei unei matrici. Fereastră program Delphi

Așa cum se observă din relațiile (1.24) și (1.36), dacă elementul akkk-1, numit pivot, are valoare nulă (sau valoarea sa absolută este sub un "prag de zero" prestabilit, atunci apar probleme la efectuarea operațiilor de împărțire. Existența unui pivot nul nu înseamnă neapărat că matricea este singulară. Trebuie încercate toate posibilitățile (poate fi adus în poziția pivotului oricare element aijk-1, i = k, k + 1, ..:, n, j = k, k + 1, …, n, cu schimbarea corespunzătoare a liniei k cu linia i, respectiv a coloanei k cu coloana j) și numai dacă toate duc la eșec, atunci matricea A este singulară [6]. De altfel, în scopul reducerii la minim a erorii de rotunjire, se recomandă ca la fiecare pas de inversare să se aducă pe poziția pivotului elementul de valoare absolută maximă, selectat în modul indicat anterior (procedeu numit "pivotare"). Se menționează că în aplicațiile din domeniul aplicațiilor electrotehnice, de regulă, elementele diagonale sunt "dominante" (mai mari în valoare absolută decât cele nediagonale), ceea ce face ca procesul de pivotare să nu fie necesar. Analiza evoluției matricelor de tip A, reflectată de relațiile (1.19), (1.22), (1.25), respectiv a matricelor de tip D , reflectată de relațiile (1.29), (1.32), (1.35), sugerează clar posibilitatea de a folosi un singur tablou de lucru pentru toată operația de inversare: la demararea inversării tabloul conține matricea originală A, iar în final matricea inversă A-1 (modificările în schema logică sunt minore: toate elementele de tip dij se înlocuiesc cu elemente de tip aij). O asemenea manieră de lucru conduce la importante economii de memorie în cazul matricelor de mari dimensiuni.

În final, se reamintește că în teorema care stă la baza inversării prin reducere la I era vorba despre "premultiplicare cu un șir de matrice". Analiza algoritmului de inversare, prezentat intuitiv, permite stabilirea formei acestor matrice, notate cu Mk, k = 1, 2, …. n, evidențiată în relațiile (1.50) – (1.51). Se apreciază însă că este nepractică și ineficientă utilizarea algoritmului cu înmulțiri matriceale, mai ales pentru valori mari ale lui n.

Mk = (1.51)

, i = 1, 2, …, n (1.52)

1.3. Factorizarea matricelor

1.3.1. Considerații preliminare

Prin factorizarea unei matrice pătrate A de ordinul n se înțelege exprimarea sa sub forma unui produs de alte două matrice de același ordin. De regulă, factorii sunt de diferite forme speciale, particulare, astfel alese încât să rezulte avantaje la aplicarea unor metode de inversare matriceală, de soluționare a sistemelor liniare de ecuații, de determinare a valorilor proprii și a vectorilor proprii [10].

În cadrul acestui subcapitol se prezintă următoarele metode de factorizare:

factorizarea LR (paragraful 1.3.2.), care presupune exprimarea lui A sub forma (1.31)

A = L R (1.53)

unde: L este o matrice inferior triunghiulară, iar R este o matrice superior triunghiulară;

factorizarea QR (paragraful ?); care presupune exprimarea lui A sub forma (1.54)

A = Q R (1.54)

unde: Q este o matrice ortogonală, iar R este o matrice superior triunghiulară.

Algoritmele concrete de factorizare LR, respectiv QR, diferă între ele prin unele elemente concrete, legate de forma particulară a matricelor factori [12].

1.3.2. Factorizarea LR

1.3.2.1. Forma generală

Factorizarea LR are la bază următoarea teoremă din algebra matriceală: o matrice pătrată A, de ordinul n, poate fi exprimată sub forma unui produs de doua matrice, L, inferior triunghiulară, și R , superior triunghiulară, dacă sunt satisfăcute condițiile (1.55).

(1.55)

Factorizarea LR a matricei A presupune exprimarea ei sub forma:

A = = (LR) =(1.56)

Efectuând produsul matriceal, rezultă relațiile care exprimă legătura dintre elementele matricei A și cele ale factorilor L și R:

aij = , i = 1, 2, …, n; j = 1, 2, …, n (1.57)

Ținând cont de forma particulară a factorilor L și R , relațiile anterioare devin:

aij = , i = 1, 2, …, n; j = 1, 2, …, n (1.58)

Privite prin prisma determinării elementelor factorilor L și R , elementele matricei A fiind cunoscute, relația (1.58) definește un sistem de n ecuații cu n2 + n necunoscute. Rezultă că, la modul general, factorii L și R nu sunt univoc definiți, existând o infinitate de soluții: valorile pentru n elemente se aleg arbitrar, celelalte rezultând din soluționarea sistemului definit de relația (1.58) [13]. Modul cum se aleg valorile pentru cele n elemente diferențiază între ele algoritmele particulare de factorizare LR:

algoritmul LR Doolittle, la care lii = 1, i = 1, 2, …, n, adică elementele diagonale ale matricei L au valoarea 1;

algoritmul LR Crout, la care r11 = 1, i = 1, 2, …, n, adică elementele diagonale ale matricei R au valoarea 1.

Dacă matricea A este simetrică, atunci R = Ll în acest caz numărul de ecuații al sistemului definit de relația (1.58) este (n2 + n) / 2 , similar cu numărul necunoscutelor, deci factorii sunt univoc definiți (nu este necesară fixarea valorilor unor elemente ale factorilor). Algoritmul de factorizare corespunzător poartă numele de Cholesky.

Pe baza algoritmilor de mai sus se scrie codul sursă pentru această metodă și se prezintă si o fereastra din programul Delphi cu un exemplu numeric rezolvat.

Program Fact_Mat;

type matrice=array[1..8,1..8]of real;

var i,j,k:integer; {variabile ajutatoare}

n:integer; {ordinul matricei}

eps:real; {pragul sub care un pivot se considera nul}

t:real; {variabila ajutatoare}

A:matrice; {matricea A}

L:matrice; {matricea L}

R:matrice; {matricea R}

{procedura pentru afsiarea matricelor}

procedure VizMat(mat:matrice);

begin

for i:=1 to n do begin

for j:=1 to n do write(mat[i,j]:10:5);

if n<8 then writeln;

end;

writeln;

end;

{procedura pentru afisarea liniei curente a matricei L si a matricei R}

procedure VizLin;

begin

writeln;

writeln('Linia ',i,' a matricei L');

writeln('––––––––');

writeln;

for i:=1 to n do write(r[i,j]:10:5);

if n<8 then writeln;

writeln;

writeln('Linia ',i,' a matricei R');

writeln('––––––––');

writeln;

for i:=1 to n do write(l[i,j]:10:5);

if n<8 then writeln;

writeln;

writeln('Apasa Enter');

readln;

end;

begin

writeln('Factorizarea matricelor cu algoritmul LR Doolittle');

writeln('**************************************************');

{citire A}

writeln;

write('Ordinul matricei A (maxim 8): n=');

readln(n);

writeln;

for i:=1 to n do begin

writeln('Linia ',i,);

for j:=1 to n do begin

write('A[',i,',',j,']=');

readln(a[i,j]);

end;

end;

writeln;

write('Pragul sub care un pivot se considera nul: eps=');

readln(eps);

writeln;

{afisare A}

writeln('Matricea A:');

VizMat(A);

writeln;

{init factorul R cu mat O si L cu I}

for i:=1 to n do begin

for j:=1 to n do begin

L[i,j]:=0;

R[i,j]:=0;

end;

L[i,i]:=1;

end;

{inceputul procesului de factorizare, efectuat in n pasi}

{pasul 1:}

for j:=1 to n do r[1,j]:=a[1,j];

i:=1;

vizlin;

{pasii 2..n:}

for i:=2 to n do begin

{linia i a factorului L}

for j:=1 to i-1 do begin

{verificare pivot}

if abs(r[j,j])<eps then begin

writeln('Pivot nul la pasul ',i);

writeln;

writeln('Apasa Enter');

readln;

exit;

end;

t:=a[i,j];

if j<>1 then for k:=1 to j-1 do t:=t-L[i,k]*r[k,j];

l[i,j]:=t/r[j,j];

end;

{linia i a factorului R}

for j:=i to n do begin

t:=a[i,j];

for k:=1 to i-1 do t:=t-l[i,k]*r[k,j];

r[i,j]:=t;

end;

vizlin;

end;

writeln;

writeln('Rezultate Finale:');

writeln('––––––');

writeln;

writeln('Factorul L');

VizMat(L);

writeln;

writeln('Factorul R');

VizMat(R);

writeln;

writeln('Apasa Enter');

readln;

end.

Figura 1.2. – Factoriazarea LR Doolittle. Fereastră program Delphi

1.3.2.2. Factorizarea LR Cholesky (metoda rădăcinii pătrate)

Factorizarea LR Cholesky se aplică pentru matricele simetrice, relația generală de factorizare devenind:

A = L Lt (1.59)

Analiza relațiilor generale de echivalare (1.57) conduce la concluzia că acum numărul de ecuații este (n2 + n) / 2, similar cu numărul necunoscutelor, deci factorii sunt univoc definiți.

Algoritmul metodei se poate obține și prin particularizarea corespunzătoare a relațiilor generale de factorizare (1.57) [5].

Se prezintă în continuare algoritmul metodei:

se calculează elemetele primei coloane a factorului L:

elementul diagonal cu relația (1.59);

l11 = (1.60)

elementele nediagonale cu relația (1.61);

li1 = , i = 2, 3, …, n (1.61)

se calculează elementele coloanei j, j = 2, 3, …, n , a factorului L:

elementul diagonal cu relația (1.62)

(1.62)

elementele nediagonale cu relația (1.63).

, i = j + 1, j + 2, …, n (1.63)

Rămân valabile toate comentariile finale de la algoritmul LR Doolittle. Suplimentar, apare condiția ca expresiile de sub radical din relațiile (1.60) și (1.62) să fie pozitive.

Capitolul 2. – Analiza circuitelor rezistive liniare

2.1. Metode și algoritmi de calcul

2.1.1.Metoda teoremelor lui Kirchhoff

Pentru un circuit rezistiv liniar conex invariabil în timp, cu l laturi, conținând rezistoare și surse independente de tensiune, și n noduri, aplicarea celor două teoreme ale lui Kirchhoff conduce la obținerea unui sistem complet de l ecuații, din care n- 1 ecuații de tipul primei teoreme a lui Kirchhoff și l – n +1 ecuații de tipul celei de-a doua teoreme a lui Kirchhoff, in l necunoscute curenți de laturi.

Dacă circuitul conține și surse de tensiune comandate, teorema a doua a lui Kirchhoff devine:

(2.1)

iar t.e.m. ale surselor comandate (Eck) se exprimă prin ecuațiile de comandă prelucrate în funcție de necunoscutele curenți de laturi, în cazul circuitelor care conțin și surse de curent independente și / sau comandate, numărul de necunoscute curenți e laturi este l – (lJ+ lJc). Acestea se obțin prin rezolvarea unui sistem de ecuații obținut prin aplicarea primei teoreme a lui Kirchhoff în n – 1 noduri independente și a celei de-a doua ecuații într-un număr de bucle independente redus la br= l – n + 1 – (lJ+lJc), unde lJ reprezintă numărul de laturi cu surse de curent independente, iar lJc reprezintă numărul de laturi cu surse de curent comandate. Acestui sistem i se adaugă ecuațiile de comandă ale surselor comandate prelucrate în funcție de necunoscutele curenți de laturi [8].

Observatii:

pentru a se obține numărul de bucle br, deci pentru a se obține un număr redus de ecuații ale sistemului, este necesară o alegere corespunzătoare a buclelor independente, astfel încât nici una din ele să nu treacă prin laturi cu surse independente și / sau comandate de curent, în caz contrar, numărul de necunoscute ale sistemului va fi l + lJc, din care l – lJ, vor fi necunoscute curenți de laturi, iar restul de lJ + lJc vor fi necunoscutele tensiuni la bornele surselor independente și / sau comandate de curent, ecuația generală, corespunzătoare celei de-a doua teoreme a lui Kirchhoff fiind exprimată cu relația (2.2):

(2.2)

este evident că alegerea unui număr redus de bucle br prezintă avantajul obținerii unui sistem redus de ecuații, deci a reducerii efortului de calcul; pe de altă parte relația permite scrierea sistematică a sistemului de ecuații;

odată calculați curenții din laturi, tensiunile la bornele laturilor se pot determina în modul următor:

pentru laturile fără surse de curent se aplică ecuația caracteristică sau teorema a doua a lui Kirchhoff;

pentru laturile formate din surse independente sau comandate de curent, numai cu ajutorul teoremei a doua a lui Kirchhoff;

dacă circuitul conține și surse de curent, atunci ecuația devine:

(2.3)

unde

(2.4)

este curentul de scurtcircuit injectat în nodul n j.

Algoritmul de aplicare al metodei teoremelor lui Kirchhoff

pasul 1: Se determină numărul nodurilor și al laturilor circuitului;

pasul 2: Se aleg sensuri de referință și se atașează simboluri pentru intensitățile curenților din laturi;

pasul 3: Se calculează numărul redus de bucle ale circuitului și se aleg aceste bucle stabilindu-se un sens de parcurgere pentru fiecare;

pasul 4: Se scriu ecuațiile corespunzătoare primei teoreme a lui Kirchhoff în (n – 1) noduri independente și ecuațiile corespunzătoare celei de-a doua teoreme pe cele b, = l – n + 1 – (/J, + lJc) bucle independente;

pasul 5: Se rezolvă sistemul de ecuații obținut prin completarea celui de la pasul 4 cu ecuațiile de comandă ale surselor de curent și de tensiune comandate, prelucrate în funcție de curenții laturilor, determinându-se intensitățile curenților din laturi;

pasul 6: Se validează rezultatul cu ajutorul bilanțului puterilor.

Analiza asistată de calculator a circuitelor electrice necesită formularea matriceală a ecuațiilor circuitului [5]. Pentru circuite reciproce, luând în considerare structura laturii standard prezentată în figura 2.1, formularea matriceală a teoremelor lui Kirchhoff și a ecuațiilor caracteristice (constitutive) ale laturilor conduce la ecuațiile:

Figura 2.1 – Schemă latură

(2.5)

(2.6)

(2.7)

Cele 2l ecuații (2.5), (2.6) și (2.7) determină în mod univoc curenții și tensiunile laturilor circuitului, dacă se dau valorile rezistențelor laturilor, ale t.e.m. ale surselor de tensiune și ale intensităților curenților surselor de curent [9].

Înlocuind relația (2.7) în (2.6) și cuplând apoi cu (2.5) se obține forma matriceală a ecuațiilor circuitului în curenții laturilor:

(2.8)

unde A (B) este matricea (n – l) x l (bxl) de incidență redusă laturi – noduri (laturi – bucle), Il (Ul) este vectorul (lx1) al curenților (tensiunilor) laturilor circuitului, Rl este matricea diagonală (lxl) a rezistențelor laturilor circuitului, iar El (Jl) este vectorul (lx1) al t.e.m. (curenților) surselor independente de tensiune (curent).

Sistemul de ecuații (2.8) se rezolvă în raport cu vectorul curenților laturilor Il, apoi cu ecuațiile (2.7) se determină tensiunile la bornele laturilor.

O formulare matriceală alternativă se obține dacă se face următoarea partiționare a laturilor circuitului:

lR- laturi conținând rezistoare și eventual surse de tensiune independente;

lE – laturi conținând surse ideale independente de tensiune;

lJ – laturi conținând surse independente de curent și, eventual, rezistoare.

Schema standard a laturii care va fi luată în considerare în continuare este reprezentată în figura 2.2, unde UbE(J) este tensiunea la bornele sursei de tensiune, ui,e = ee, sau tensiunea la bornele sursei de curent, echivalentă cu o t.e.m. UbJ = EJ.

Figura 2.2 – Schemă latură

În această situație ecuația (2.7) se scrie sub forma:

(2.9)

Dacă circuitul are numai laturi de tipul lR și lE, toți curenții laturilor și l-lE tensiuni la borne (pentru laturile tip lE, UE= – EE) sunt necunoscute. Înlocuind relația

(2.9) în (2.6) și cuplând apoi cu (2.5) se obține:

(2.10)

care se rezolvă în raport cu Il, apoi se determină tensiunile la bornele laturilor din sistemul (2.9).

În cazul în care circuitul conține și surse de curent, numărul necunoscutelor curenți de laturi este l-lJ. În locul curenților de laturi apar ca necunoscute tensiunile la bornele surselor de current [11]. Pe baza partiționării propuse pentru laturile circuitului și notând ll =lR + lE, vectorii Il, Ul, El capătă forma:

(2.11)

ce pune în evidență necunoscutele circuitului Il, UR, UJ, EJ.

Cu partiționarea propusă pentru laturi, matricea Rl a rezistențelor laturilor, de dimensiune (lxl) are forma:

(2.12)

Analiza circuitelor care conțin toate tipurile de laturi definite mai sus se poate face în două variante:

Varianta 1

(2.13)

Cu partiționarea laturilor în ll și lJ, ecuațiile (2.5), (2.6) și (2.9) devin

(2.14)

respectiv

(2.15)

Prelucrând aceste ecuații matriceale se obține:

(2.16)

(2.17)

(2.18)

(2.19)

Înlocuind relația (2.18) în (2.17), separând mărimile necunoscute și cuplând ecuațiile, se obține forma compactă

(2.20)

Ecuația (2.20) se numește ecuația matriceală hibridă a circuitului și are forma:

(2.21)

Sistemul ecuațiilor circuitului fiind complet și liniar independent, matricea C de dimensiune (lxl) este nesingulară, deci cele l necunoscute ale circuitului din care ll = lR + lE = l – lJ, curenți de laturi, și lJ tensiuni la bornele laturilor cu surse independente de curent, se determină din ecuația matriceală (2.20) [4].

Din ecuația (2.18) se determină apoi vectorul Ul, iar din ecuația (2.19) se calculează tensiunile la bornele surselor de curent:

(2.22)

Varianta 2

Această variantă permite obținerea numărului minim de ecuații ale circuitului în necunoscute curenți de laturi. Pentru aceasta se observă că se disting două tipuri de bucle în construcția matricei B: bucle fară surse de curent (f) și bucle cu surse de curent (c) [7]. Alegând sistemul buclelor fundamentale astfel încât fiecare buclă să conțină cel mult o sursă de curent (alegere asigurată plecând de la un arbore construit cu laturi din mulțimea ll și formând bucle prin laturile arborelui cu câte o coardă a coarborelui care conține toate lJ, cu partiționarea laturilor în ll și lJ, matricea B se scrie sub forma:

(2.23)

Cu această precizare, ecuațiile matriceale (2.5), (2.6) și (2.9) devin:

(2.24)

(2.25)

(2.26)

(2.27)

(2.28)

Înlocuind ecuația (2.27) în (2.25) și separând mărimile cunoscute, sistemul redus al ecuațiilor circuitului poate fi pus sub forma:

(2.29)

Ecuația (3.146) reprezintă ecuația matriceală redusă în variabile curenți de laturi a circuitului, de forma:

(2.30)

unde matricea C, de dimensiunea (l-lJ)x(l-lJ) este nesingulară.

Vectorul necunoscutelor curenți de laturi se obține deci cu ecuația (2.29) apoi din ecuația (2.27) se obține vectorul tensiunilor Ul și din ecuația (2.26) se calculează tensiunile la bornele laturilor cu surse de curent

(2.31)

Pentru o ordonare corespunzătoare a buclelor cu sursă de curent și a laturilor lJ matricea BcJ este chiar matricea unitate. Tensiunile la bornele surselor de curent se determină din ecuația (2.28).

2.1.2.Metoda curenților de buclă

Pentru circuitele de mari dimensiuni, sistemul de ecuații obținut prin aplicarea teoremelor lui Kirchhoff, fie în varianta hibridă cu l necunoscute, fie în varianta cu număr redus la l – lJ necunoscute, poate fi de dimensiuni prea mari. Apare deci necesitatea utilizării unor metode alternative de analiză a circuitelor electrice, care să reducă numărul ecuațiilor ce descriu funcționarea circuitului, respectiv numărul variabilelor independente [2].

Una din aceste metode este metoda curenților de buclă, care asociază circuitului un nou set de necunoscute – curenții de bucle Ib, în număr de l – n + l, introduse astfel încât să verifice prima teoremă a lui Kirchhoff. Prin urmare, curenții laturilor se exprimă ca sumă algebrică a curenților de buclă ce trec prin latura respectivă (fig. 2.3):

Figura 2.3 – Latura unui circuit

(2.32)

Noile necunoscute se determină cu ajutorul teoremei a doua a lui Kirchhoff, în funcție de (2.32). Se obține forma compactă a ecuațiilor circuitului în necunoscute curenți de bucle:

(2.33)

Dacă în circuit există surse de tensiune și / sau de curent comandate, sistemul de ecuații (2.33) se completează cu ecuațiile de comandă exprimate in funcție de necunoscutele metodei, Ibg.

Numărul de variabile independente introdus de metoda curenților de buclă este b = l – n + 1. Pentru circuitele fără surse de curent aceste necunoscute se determină prin rezolvarea sistemului (2.33).

În cazul când circuitul conține surse de curent independente și / sau comandate, unor variabile li se pot atribui valorile curenților acestor surse [3]. Pentru aceasta, asocierea variabilelor Ib, cu cele l – n + 1 bucle independente ale circuitului se face astfel încât prin fiecare latură cu sursă de curent J sau Jc să treacă un singur curent de buclă și numai unul. Conform relației (2.32) acest curent de buclă va avea valoarea curentului sursei:

și / sau

(2.34)

Pentru restul variabilelor, în număr de l – n + 1 – (lJ + lJc), se aplică relația (2.33) într-un număr redus de bucle, deci h = 1, …, br. Sistemului obținut din ecuațiile (2.33) și (2.34) i se adaugă ecuațiile de comandă ale surselor comandate, exprimate în funcție de necunoscutele metodei, Ibg.

Observație: pentru circuitele care conțin surse de curent independente J și / sau comandate Jc, metoda curenților de buclă permite o reducere a numărului de ecuații ale sistemului (2.33) cu numărul total al acestor surse de curent, în condițiile alegerii corespunzătoare a buclelor.

Algoritmul de aplicare al metodei curenților de buclă

pasul 1: se determină numărul de noduri, numărul de laturi și numărul surselor de curent;

pasul 2: se determină numărul variabilelor independente introduse de metodă: l- n+1;

pasul 3: se aleg (lJ+lJc) bucle care să conțină câte o singură sursă de curent independentă sau comandată și li se atașează câte un curent de buclă al cărui sens va fi același cu al sursei de curent; curenții acestor bucle vor fi exprimați cu relațiile (2.34);

pasul 4: pentru restul de bucle independente, în număr redus la br = l – n +1 – (lJ + lJc ), se atribuie tot atâtea variabile curenți de buclă cu sensuri oarecare, fiecare reprezentând și sensul de parcurgere al buclei respective;

pasul 5: se aplică teorema a doua a lui Kirchhoff în aceste bucle, ținând seama că în membrul stâng al ecuației (2.33) pot apare căderi de tensiune determinate de variabile exprimate cu relațiile (2.34);

pasul 6: se atașează sistemului obținut cu ecuațiile.(2.33) și (2.34) ecuațiile de comandă ale surselor comandate, exprimate în funcție de variabilele metodei;

pasul 7: se rezolvă sistemul astfel obținut în variabile curenți de buclă;

pasul 8: se determină curenții laturilor cu ecuația (2.32);

pasul 9: tensiunile la bornele laturilor se determină în modul descris în §2.1.1

pasul 10: se verifică bilanțul puterilor.

Formularea matriceală a metodei curenților de buclă în cazul circuitelor reciproce depinde de tipul adoptat pentru latura standard [9].

Notând cu h vectorul (bx1) al curenților de buclă asociați celor b bucle independente, b = l – n + 1, curenții laturilor circuitului se pot exprima în funcție de curenții de bucle cu relația:

(2.35)

Pentru latura reprezentată în figura 2.1, înlocuind relația (2.35) în (2.7) și relația astfel obținută în (2.6), se obține ecuația matriceală a curenților de buclă sub forma

(2.36)

sau

(2.37)

cu

(2.38)

si

(2.39)

unde matricea Rb de dimensiune (bxb) este matricea rezistențelor proprii ale buclelor, iar Eb de dimensiune (bx1), este vectorul t.e.m. ale buclelor.

Rezolvând ecuația (2.37) se obține vectorul necunoscutelor curenți de bucle, cu ajutorul cărora, aplicând relația (2.35), se determină curenții din laturile circuitului.

Dacă se adoptă latura standard din figura 2.2, se folosește ecuația caracteristică a laturii în varianta (2.9) și dacă circuitul conține numai laturi tip lR și lE, în locul relației (2.36) se obține ecuația

(2.40)

unde

(2.41)

În cazul când circuitul conține și laturi tip lJ, se poate reduce sistemul (2.37) la un număr de ecuații corespunzător unui număr redus de bucle br= l – n + 1 – lJ.

Pentru aceasta, toate sursele de curent din graful circuitului vor fi introduse în coarbore. Cum sistemul buclelor fundamentale (în număr de l – n + l) se formează cu laturi ale arborelui și cu câte o coardă, fiecare curent de buclă va fi asociat unei coarde având același sens cu sensul ei de orientare. O parte din aceste coarde sunt laturile tip lJ, iar restul, laturi tip lR. Se asigură astfel condiția ca fiecare buclă să conțină cel mult o sursă de curent. Prin urmare, matricea B se poate pune sub forma (2.23) [6].

Partiționând curenții laturilor și buclelor după l1 și lJ, ecuația (2.35) devine

(2.42)

din care se obtin ecuatiile

(2.43)

(2.44)

Din ecuația (2.44) rezultă:

(2.45)

cu semnificația: curenții buclelor care conțin o sursă de curent (sunt formate cu o latură coardă de tip lJ) au valoarea egală cu a curenților surselor de curent corespunzătoare.

Exprimând matricea Rl cu relația (2.12) și matricea B cu relația (2.23), ecuația (2.40) se scrie ca

(2.46)

din care rezultă subsistemele matriceale de ecuații:

(2.47)

(2.48)

Ținând seama de relația (2.45), din subsistemul (2.47), cuprinzând l – n + l – lJ ecuații, se determină vectorul necunoscutelor Ib1. Utilizând apoi relația (2.43) se obține vectorul necunoscutelor curenți de laturi tip l1, I1. Din subsistemul (2.48) se obține vectorul tensiunilor la bornele surselor de curent UbJ = EJ, cu relația:

(2.49)

Utilizând apoi relația (2.9) se poate calcula vectorul necunoscutelor tensiuni la bornele laturilor circuitului.

2.1.3.Metoda potențialelor nodurilor

Această metodă asociază circuitului setul de necunoscute potențiale ale nodurilor, Vn-1, în număr de n- l, introduse astfel încât să satisfacă a doua teoremă a lui Kirchhoff. Ca urmare, tensiunile laturilor [9] se exprimă ca sumă algebrică a potențialelor adiacente laturii respective (fig. 2.4):

Figura 2.4

(2.50)

Unul din cele n noduri ale circuitului este ales ca nod de referință, având potențialul nul.

Noile necunoscute se determină cu ajutorul primei teoreme a lui Kirchhoff sub forma obținută prin substituirea relației (2.50) în , reprezentând forma compactă a ecuațiilor circuitului în variabile potențiale la noduri:

(2.51)

unde:

Dacă circuitul conține surse de tensiune și / sau de curent comandate, sistemul de ecuații (2.5) se completează cu ecuațiile de comandă exprimate în funcție de necunoscutele metodei, Vj.

Numărul de variabile independente introdus de această metodă este n- 1

Pentm circuitele fără surse ideale de tensiune independente sau comandate E', respectiv E'c, necunoscutele potențiale la noduri se determină rezolvând sistemul (2.5) format din n – 1 ecuații independente [6].

În cazul circuitelor care conțin surse ideale de tensiune, potențialele nodurilor i și j la care este conectată o astfel de latură (fig. 2.5a si 2.5b) se exprima cu relatia (2.52), respectiv (2.53).

Figura 2.5a

Figura 2.5b

(2.52) (2.53)

de unde:

(2.54) (2.55)

Rezultă deci că pentru (lEi + lEic ) necunoscute se pot formula ecuații de tipul

(2.56)

și / sau

Pentru restul necunoscutelor ar trebui să se aplice ecuațiile (2.51) într-un număr redus de noduri, nr = n-1-(lEi+ lEic), deci i = 1,nr.

Observații:

ecuația (2.5) nu se poate aplica într-un nod în care este incidență o latură cu sursă ideală de tensiune, deoarece curentul de scurtcircuit al acestei surse este infinit (rezistența ei internă este zero), în acest caz se poate recurge la următoarea tehnică: se alege o suprafață închisă ∑, care să cuprindă în interior latura ij ce conține sursa ideală de tensiune, sau, dacă este cazul, toate laturile conectate în paralel între nodurile i și j, pe care se scrie apoi prima teoremă a lui Kirchhoff. Se aplică apoi sistemul (2.51) în nr = n – l – (lEi + lEic), noduri și suprafețe ∑, adică pentru i = 1, nr;

rezultă că pentru circuitele care conțin surse ideale de tensiune independente și / sau comandate, (Ei și/sau Eic), metoda potențialelor nodurilor permite o reducere a numărului de ecuații de forma (2.51) cu numărul total al acestor surse de tensiune.

Algoritmul de aplicare al metodei potențialelor nodurilor

pasul 1: se determină numărul de noduri ale circuitului;

pasul 2: se alege un nod j de referință al cărui potențial se consideră nul, Vj = 0;

pasul 3: se scriu (lEi+ lEic), ecuații de tipul (2.56) pentru potențialele nodurilor adiacente surselor ideale de tensiune;

pasul 4: se aplică relațiile (2.51) în n, = n- 1 -(lEi+ lEic), noduri și suprafețe S, adică pentru i = 1,nr, ținând seama de faptul că în termenii din partea stângă ai relațiilor pot interveni și potențiale pentru care s-au scris ecuațiile de la pasul 3;

pasul 5: sistemului obținut cu relațiile (2.51) și (2.56) i se adaugă ecuațiile de comandă ale surselor comandate, exprimate în funcție de variabilele metodei;

pasul 6: se rezolvă sistemul de la pasul 5 și se obțin valorile celor n-l variabile potențiale ale nodurilor;

pasul 7: cu relația (2.50) se calculează apoi tensiunile la bornele laturilor circuitului;

pasul 8: se determină curenții din laturile circuitului cu ecuația caracteristică a laturii pentru laturile care conțin rezistențe și eventual surse de tensiune înseriate cu acestea, sau cu prima teoremă a lui Kirchhoff pentru cele formate din surse ideale de tensiune;

pasul 9: se verifică bilanțul puterilor.

O altă variantă de aplicare a metodei constă în introducerea ca necunoscute în sistemul de ecuații a curenților prin laturile cu surse ideale de tensiune. Deși are un număr mai mare de ecuații, această metodă cu necunoscute hibride, numită metoda nodală modificată, permite scrierea sistematică a sistemului de ecuații [12].

Pentru formularea matriceală a metodei potențialelor nodurilor în cazul circuitelor reciproce, se adoptă latura standard din figura 2.1. Din cele n noduri ale circuitului se alege un nod de referință, cu potențial zero, iar celorlalte n – 1 noduri li se atribuie potențialele necunoscute V1, V2, …, Vn-1.

Tensiunile laturilor se pot exprima în funcție de potențialele nodurilor cu relația:

(2.57)

Considerând ecuația caracteristică a laturii

(2.58)

unde Gl = R-1l reprezintă matricea pătrată (lxl) a conductanțelor laturilor (se presupune Rl nesingulară), substituind relația (2.57) în (2.58) și ecuația astfel obținută în (2.5), se obține forma matriceală a ecuațiilor circuitului în potențialele nodurilor:

(2.59)

Ecuația (2.59) se poate pune sub forma:

(2.60)

unde

(2.61)

este matricea adrriitanțelor nodale, de ordin (n – 1) x (n – 1) și:

(2.62)

este vectorul curenților injectați în noduri de sursele din laturile incidente în aceste noduri.

Relația (2.60) permite scrierea simplă pe calculator a ecuațiilor circuitului.

Rezolvând sistemul (2.60) se obține vectorul potențialelor necunoscute ale celor n – l noduri, apoi cu relația (2.57) se calculează vectorul tensiunilor la bornele laturilor Ul. Cu relația (2.58) se obține în cele din urmă vectorul curenților laturilor.

Sub forma prezentată mai sus, metoda nu permite rezolvarea matriceală a circuitelor cu laturi alcătuite din surse ideale independente de tensiune, pentru care Rl = 0, deci produsul GlEl, reprezentând curentul de scurtcircuit al sursei, este infinit. Limitele metodei sunt depășite în această situație de metoda nodală modificată [3].

Dacă se alege pentru latura standard schema din figura 2.2, se poate obține formularea matriceală a ecuațiilor circuitului, inclusiv pentru cazul când acesta conține surse ideale independente de tensiune. Partiționând nodurile circuitului în noduri la care laturile cu surse ideale independente de tensiune nu au incidență (a) și noduri adiacente laturilor cu sursă ideală de tensiune (b) și alegând ca nod de referință un nod de tip b, matricea redusă de incidență a laturilor la noduri, A, poate fi pusă sub forma

(2.63)

În aceste conditii ecuatia (2.57) poate fi scrisa ca:

(2.64)

unde UE = -EE, iar UR si UJ sunt necunoscute.

Daca notam l2 = lR + lJ, din ultima ecuatie obtinem vectorul

(2.65)

Prima teorema a lui Kirchhoff devine:

(2.66)

din care putem obtine ecuatiile:

(2.67)

(2.68)

Scriind ecuatiile caracteristice ale laturilor pentru toate laturile de tip lr,

(2.69)

si substituind UR prin relatia corespunzatoare din (2.64), obtinem:

(2.70)

Substituind această relație în ecuația (2.67), cuplând cu primul subsistem din (2.64) și separând variabilele, obținem ecuația matriceală în potențialele nodurilor:

(2.71)

Sistemul (2.71) conține n – 1 variabile independente. Pentru ca sistemul să aibe soluție unică este necesar ca numărul de noduri neadiacente la laturi de tip lE, n , să satisfacă relația n = n – 1 – lE. Această condiție impune ca laturile de tip le să formeze un arbore. Ca urmare, primul subsistem va fi format din n – 1 – lE ecuații, în timp ce al doilea cuprinde lE ecuații simple de tip UE = -EE, exprimate în potențiale la noduri [1]. După rezolvarea sistemului (2.71), se poate calcula vectorul U2, cu relația (2.65), apoi vectorul IR, cu relația (2.69), respectiv vectorul IE cu relația (2.68), matricea AE fiind, în aceste condiții, pătrată și nesingulară. Metoda expusă permite analiza completă a unui circuit reciproc cu surse ideale independente de tensiune, dacă acestea formează un arbore și nodul de referință este adiacent unei laturi cu sursă ideală independentă de tensiune.

Capitolul 3. – Soluționarea numerică a sistemelor de ecuații liniare

3.1. Prezentarea problemei

Așa cum s-a arătat mai sus, multe probleme din domeniul tehnicii, în general, și din cel al aplicațiilor electrotehnice în particular, conduc la modele matematice care implica soluționarea unor sisteme de ecuații liniare de dimensiuni mari: abordarea numerică a funcționării mașinilor electrice, hidraulice sau termice; analiza, optimizarea și estimarea regimurilor de funcționare a sistemelor electroenergetice, modelarea numerică a funcționării aparatelor și echipamentelor electrice. De. altfel, pentru orice problemă se poate elabora, în primă instanță, un model liniar sau, dacă modelul este neliniar, ci poate fi liniamal în primă aproximație, o singură dată sau la fiecare pas al unui proces iterativ de soluționare [8].

Rezolvarea numerică a sistemelor de ecuații liniare constituie una din problemele cele mai frecvent întâlnite în aplicațiile electrotehnice. De exemplu la rezolvarea unui circuit electric (determinarea curenților), scrierea teoremelor lui Kirchhoff conduce la un model matematic care implică în fapt rezolvarea unor sisteme de ecuații liniare de dimensiuni mari. Ca exemplu mai putem da și modelarea numerică și analiza funcționării unor aparate și echipamente electrice.

Se știe că regula lui Cramer nu corespunde cerințelor practice când numărul de ecuații este mai mare ca 3 și când determinantul corespunzător matricii sistemului este zero. Atunci când suntem în acest caz sau atunci când numărul de ecuații nu este egal cu numărul de necunoscute sistemul se rezolvă prin metoda Kroneker – Capelli.

Există 2 categorii de metode de rezolvare a sistemelor de ecuații liniare de forma Ax = b:

metode directe sau “exacte”;

metode indirecte sau “iterative”.

La metodele directe soluția sistemului rezultă printr-o serie de operații care se execută o singură dată, numărul total de operații aritmetice elementare fiind finit și cunoscut de la început. Nu este necesar să se calculeze inversa matricei coeficienților.

La metodele iterative soluția se obține printr-o serie (proces) de aproximații succesive, fiecare secvență de operații aritmetice elementare (mai mic decât la metodele directe) este parcursă de mai multe ori, obtinându-se aproximații din ce în ce mai bune ale soluției până la atingerea unei precizii fixate dinainte (precizie dorită). Aceste metode permit obținerea soluției numerice a unui sistem de ecuații prin generarea unui șir care tinde la soluția exactă.

Un sistem de n ecuații liniare, având coeficienții aij = 1, 2, …, n, cu n neconoscute și termenii liberi bi , i = 1, 2, …. n , cu n necunoscute, xi = , i = 1, 2, …, n, este de forma definită de relația (3.1), care se poate scrie restrâns sub forma (3.2).

(3.1)

, i = 1, 2, …, n (3.2)

Introducând notațiile matriceale și vectoriale pentru matricea A a coeficienților, vectorul coloană b al termenilor liberi și vectorul coloană x al necunoscutelor [11].

A = (3.3)

b = [b1 b2 b3 … bn]t (3.4)

x = [x1 x2 x3 … xn]t (3.5)

rezulta forma matriceală (3.6) a sistemului liniar (sistemului algebric):

A x = b (3.6)

După natura vectorului b se disting următoarele tipuri de sisteme de ecuații liniare: sisteme neomogene, la care b ≠ 0 și sisteme omogene, la care b = 0. Sistemele de ecuații liniare neomogene au o soluție unică dacă și numai dacă matricea A este nesingulară (|A| ≠ 0), acesta fiind, de regulă, cazul de interes pentru aplicațiile din domeniul aplicațiilor electritehnice. Dacă matricea A este singulară (|A| = 0), atunci cel puțin una dintre ecuații nu este liniar independentă, sistemul având o infinitate de soluții (valoarea a cel puțin uneia dintre componentele vectorului x se adoptă arbitrar, celelalte rezultând din soluționarea sistemului care rămâne după renunțarea la ecuațiile care sunt liniar dependente).

Sistemele de ecuații liniare omogene admit soluția banală x = 0 . Ele admit și alte soluții, nebanale, dacă și numai dacă matricea A este singulară. De regulă, aceasta este situația de interes în domeniul aplicațiilor electrotehnice, cu specificația că trebuie determinate și valorile unui parametru care apare în expresiile unor coeficienți ai necunoscutelor, pentru care matricea A este singulară. Soluțiile nebanale se pot calcula în modul indicat anterior.

Din punct de vedere al soluțiilor, pentru un sistem liniar de ecuații pot exista următoarele situații: sistemul nu are nici o soluție, sistemul are soluție unică sau sistemul are o infinitate de soluții. Se utilizează termenul de "sistem consistent" pentru un sistem de ecuații algebrice care are cel puțin o soluție, respectiv cel de "sistem inconsistent" pentru un sistem care nu are nici o soluție [4].

Având în vedere că în multe aplicații din domeniul tehnicii valorile coeficienților și ale termenilor liberi pot fi afectate de erori (au rezultat din determinări experimentale, din alte calcule "aproximative" sau cu ipoteze simplificatoare, etc.), prezintă interes cunoașterea măsurii în care ele influențează soluțiile sistemelor de ecuații liniare. Dacă perturbații sau modificări mici ale elementelor lui A sau b conduc la modificări reduse ale elementelor soluției x, atunci se spune că sistemul este "bine condiționat" în caz contrar, sistemul este "slab condiționat" – modificări relativ reduse ale unor coeficienți sau ale unor termeni liberi au ca efect variații accentuate ale componentelor soluției. Sistemele de ecuații liniare sunt cu atât mai bine condiționate, cu cât elementele diagonale ale matricei A sunt mai "dominante" (mai mari în valoare absolută) în raport cu elementele nediagonale de pe linia respectivă.

În scopul ameliorării gradului de condiționare a sistemelor de ecuații liniare, a stabilității numerice a algoritmelor de soluționare, înainte de începerea efectivă a rezolvării se modifică ordinea ecuațiilor și / sau a necunoscutelor, astfel încât să se asigure dominanța elementelor diagonale ale matricei A. De regulă, în aplicațiile din domeniul energetic, această operație nu este necesară, dominanța elementelor diagonale fiind asigurată prin modul de definire a problemei.

Tot din punct de vedere al asigurării condiționării cât mai bune a sistemului, al stabilității numerice a algoritmilor de soluționare, este important ca variabilele, respectiv termenii liberi, să fie de același ordin de mărime. Dacă aceste cerințe nu sunt îndeplinite, sistemul (3.1) se transformă într-un alt sistem, mai avantajos din punctele de vedere discutate, operația de "scalare" – modificarea corespunzătoare a valorii variabilelor și a termenilor liberi, precum și a elementelor matricei A, prin adoptarea coeficienților de scalare (factorilor de scară) nxj, j = 1, 2, …, n, și nbi, i = 1, 2, …, n:

xj = nxj xj’ (3.7)

bi = nbi bj’ (3.8)

Înlocuind expresiile de mai sus în relația (3.2) rezultă succesiv:

, i = 1, 2, …, n (3.9)

, i = 1, 2, …, n (3.10)

unde coeficienții a’ij au expresiile date de relația (3.11).

a’ij = , i = 1, 2, …, n, j = 1, 2, …, n (3.11)

Ordinea calculelor este următoarea: se adoptă valorile coeficienților de scalare; se determina noile valori ale elementelor lui A, cu relațiile (3.11) și ale lui b, cu relațiile (3.8) rezolvă sistemul (3.10), se calculează valorile necunoscutelor originale cu relațiile (3.7). Se menționează că și calculul în unități relative, foarte mult utilizat în domeniul energetic, constituie de fapt, din punctul de vedere discutat, o operație de scalare [9].

Pentru rezolvarea unei aplicații concrete se alege acea metodă (directă sau iterativă), pentru care timpul total de calcul (numărul total de operații aritmetice elementare) este cel mai redus, în condițiile asigurării preciziei impuse și a posibilității estimării simple a erorii.

3.2. Metode directe de soluționare

3.2.1. Rezolvarea unui sistem n ecuații liniare cu n necunoscute, folosind regula lui Cramer

Considerăm următorul sistem de n ecuații liniare cu n necunoscute:

iar

și

corespunzător, matricea sistemului de ecuații, coloana termenilor liberi și coloana necunoscutelor. Presupunem, că determinantul sistemului de ecuații este:

Dacă acum vom înlocui succesiv în determinantul principal d coloanele coeficienților de pe lângă necunoscutele xj (j = 1, 2, …, n) cu coloana termenilor liberi bi, atunci se vor căpăta următorii determinanți:

;

Teorema lui Cramer. Sistemul de n ecuații liniare cu n necunoscute, determinantul principal al căruia este diferit de zero, are totdeauna soluție unică, determinată astfel:

Formulele de mai sus reprezintă formulele lui Cramer sau regula lui Cramer.

3.2.1.1. Algoritmul și programul rezolvării unui sistem de ecuații liniare după regula lui Cramer

Rezolvarea efectivă a unui astfel de sistem se poate face prin calculul a n+1 determinanți. Dacă cunoaștem algoritmul de calculare a determinantului unei matrice pătratice, atunci n-avem decât să-l transformăm într-o funcție care va întoarce valoarea determinantului unei matrice. Vom apela această funcție pentru calculul tuturor determinanților care intervin în regula lui Cramer [1].

Vom obține programul descris în continuare care rezolvă sisteme de ecuații liniare, comparatile și determinante (determinantul principal este nenul).

Să exemplificăm prin rezolvarea sistemului de ecuații următor:

calculăm determinantul principal:

calculăm determinanții, alcătuiți din coeficienții necunoscutei xi

utilizând regula lui Cramer, determinăm soluțiile:

program rezolvarea_unui_sistem_de_ecuatii_liniare_dupa_Cramer;

uses crt;

const maxn=30;

type matrice=array[1..maxn,1..maxn] of real;

vector=array[1..maxn] of real;

var a:matrice;

liber,x,memo:vector;

n,i,j:integer;

d:real;

procedure citire_matrice(var a:matrice;var n:integer);

var i:integer;

begin

for i:=1 to n do

for j:=1 to n do

begin

write('a[',i,',',j,']='); readln(a[i,j]);

end;

end;

procedure afisare_matrice(var a:matrice; var n:integer);

begin

for i:=1 to n do

begin

for j:=1 to n do write(a[i,j]:5:1,' ');

writeln;

end;

end;

function det(a:matrix;n:integer):real;

{ Determinantul unei matrice nxn }

var i,j,k:integer;d:real;

begin

for i:=1 to pred(n) do

begin

if abs(a[i,i])<1e-15

then begin

det:=0.0;

exit

end;

for j:=succ(i) to n do

begin

d:=a[j,i]/a[i,i];

for k:=i to n do

a[j,k]:=a[j,k]-d*a[i,k];

end;

end;

d:=1.0;

for i:=1 to n do

d:=d*a[i,i];

det:=d;

end;

procedure Cramer(var a:matrice;var n:integer);

var i,j:integer;

d:real;

begin

d:=det(a,n); { determinare determinant }

if d=0 then begin

write('Nu se poate rezolva cu regula lui Cramer!');

exit

end;

for j:=1 to n do

begin

for i:=1 to n do

begin

memo[i]:=a[i,j];

{ memorarea matricei coeficientilor pentru restabilirea ulterioara }

a[i,j]:=liber[i];

{modificarea matricei coeficientilor prin inlocuirea coloanei j prin vectorul

termenilor liberi }

end;

writeln;

afisare_matrice(a,i);

writeln('Determinantul DX','=',det(a,i):5:2);

x[j]:=det(a,n)/d;

write('Solutia x',j,'=',det(a,n):5:2,'/',d:5:2,'=',x[j]:5:2,' ');

readkey;

for i:=1 to n do a[i,j]:=memo[i];

end;

end;

begin { programul de baza }

clrscr;

write('Cate ecuatii liniare ');readln(n);

writeln('Dati matricea coeficientilor: ');

citire_matrice(a,n);

writeln('Dati termenii liberi');

for i:=1 to n do

begin

write('liber[',i,']=');

readln(liber[i]);

end;

clrscr;

afisare_matrice(a,n);

d:=det(a,n);

writeln('Determinantul matricei este d=',d:10:2);

writeln;

Cramer(a,n);

readln;

end.

3.2.2. Metoda inversării matriceale

Se consideră sistemul de ecuații liniare (de ordin n) scris în formă matricială:

Ax = b.

Dacă matricea A este nesingulară (det(A) ≠ 0) – acesta fiind de regulă cazul de interes pentru aplicațiile electrotehnice – relația matricială se poate multiplica cu rezultând:

x = A-1 b (3.12)

Această relație evidențiază cele 2 faze ale metodei:

inversarea matricei A (relații cunoscute);

efectuarea produsului matricei

Metoda inversării matriceale necesită un timp de calcul mare (datorită numărului mare de operații aritmetice elementare), aplicarea ei fiind justificată numai în situațiile în care este necesară soluționarea repetată a sistemului de forma (3.6) pentru diferite valori ale termenilor liberi (inversarea matriceală se face o singura dată, la prima rezolvare, la soluționările următoare fiind necesară numai efectuarea înmulțirii matriceale) [10].

3.2.3. Metoda de eliminare a lui Gauss

Metoda lui Gauss este una din metode tradiționale directe de rezolvare a sistemelor de ecuații liniare. Ideea de bază a metodei constă în aducerea sistemului de ecuații prin transformări elementare la o formă echivalentă, având matrice superior sau inferior triunghiulară, urmată de rezolvarea sistemului rezultat prin procedee recurente specifice, foarte eficiente.

Transformarea sistemului inițial într-un sistem de formă triunghiulară se realizează cu ajutorul a trei operații de elementare sau bază:

interschimabrea a două ecuații între ele;

înmulțirea unei ecuații cu o constantă nenulă;

scăderea unei ecuații din alta și înlocuirea celei de-a doua ecuație cu rezultatul scăderii.

Transformarea sistemului este echivalentă cu eliminarea succesivă a necunoscutelor din ecuații și se numește de aceea faza eliminării.

Rezolvarea sistemului cu matrice triunghiulară constă în determinarea necunoscutelor și substituția lor în ecuațiile sistemului în ordine inversă, fiind denumită din acest motiv faza substituției inverse.

Fie pentru exemplificare este dat un sistem de trei ecuații cu trei necunoscute:

(3.13)

sau sub formă matriceală

Ax=b

cu A = [aij] – matricea sistemului, x = [xi] – matricea coloană a necunoscutelor și b = [bi] – matricea coloană a termenilor liberi [3].

La primul pas al etapei eliminării urmărim eliminarea necunoscutei x1 din toate ecuațiile sistemului, cu excepția primei ecuații. Pentru aceasta, împărțim mai întâi prima linie la elementul pivot a11, presupus nenul, adică a11 0 (dacă nu este îndeplinită această condiție, reordonăm și numerotăm ecuațiile):

Scădem apoi prima ecuație înmulțită cu primul coeficient al celei de-a doua ecuații din această ecuație și, respectiv, înmulțită cu primul coeficient al celei de-a treia ecuații din aceasta din urmă [5]. Obținem astfel sistemul:

cu

Matriceal, primul pas al metodei eliminării lui Gauss conduce la:

(3.14)

La următorul pas, eliminăm necunoscuta x2 din ultima ecuație. Pentru aceasta, împărțim mai întâi a doua ecuație la elementul pivot , diferit de zero (dacă nu este așa, atunci realizăm interschimbarea ecuațiilor a doua și a treia) și apoi scădem linia obținută, înmulțită cu , din ecuația a treia. Matriceal, al doilea pas ne conduce la sistemul:

(3.15)

cu

(3.16)

Faza eliminării se încheie, împărțind cea de a treia ecuație la elementul pivot , care, pentru un sistem cu matrice nesingulară, trebuie să fie diferit de zero. Rezultă după acest pas sistemul:

cu (3.17)

sau matriceal, A(3)x=b(3).

Observăm ca matricea A(3) este superior triunghiulară, iar sistemul căpătat este echivalent cu cel inițial Ax = b, adică are soluția (x1, x2, x3).

Faza substituției inverse (mersul înapoi) presupune parcurgerea în sens invers a ecuațiilor sistemului cu matrice triunghiulară (3.17), rezultat în faza eliminării, și stabilirea soluției sistemului potrivit unui calcul recursiv:

(3.18)

După cum se observă, determinarea componentelor soluției are loc de la indici mari spre indici mici, fiecare nouă componentă depinzând în mod explicit numai de componentele determinate anterior [10].

Observație: metoda de eliminare Gauss permite și calcularea determinantului matricii sistemului. Se observă că, matricea A(3) a sistemului (3.17) fiind triunghiulară, are determinantul egal cu produsul elementelor diagonale, adică det A(3) = 1.

Având în vedere, însă, că prin împărțirea liniilor matricii sistemului la elementele pivot rezultă o matrice având determinantul egal cu determinantul matricii inițiale împărțit la produsul elementelor pivot, rezultă:

adică

Metoda descrisă mai sus poate fi generalizată și aplicată la rezolvarea a m sisteme de n ecuații liniare cu n necunoscute.

(3.19)

sau sub formă matriceală

Atunci putem spune că la fiecare pas k al fazei eliminării se efectuează împărțirea liniei pivot k la elementul diagonal . Deoarece este posibil ca un element pivot să fie egal cu zero (imposibilitate de continuare a fazei eliminării) se procedează în practică astfel:

se determină elementul maxim al matricei A(k–1) situat pe coloana k și liniile l k (pe diagonala principală și sub ea), adică elementul alk;

se aduce elementul maxim alk găsit pe diagonala principală, interschimbând între ele liniile l și k.

Această metodă se numește pivotare parțială pe coloană.

Se descrie în continuare algoritmul eliminării lui Gauss cu pivotare parțială pe coloane:

Date de intrare:

n {ordinul matricei}

m {numărul sistemelor simultane}

aij, i = 1, 2, …, n, j = 1, 2,…, n {matricea sistemului}

bij, i = 1, 2, …, m, j = 1, 2, …, m {matricea termenilor liberi}

Date de ieșire:

bij, i = 1, 2, …, n, j = 1, 2, …, m {soluția sistemului}

Det {determinantul sistemului}

Faza eliminării

Pasul 1. Det = 1;

Pasul 2. Pentru k = 1, …, n execută Pașii 3 – 6

Pasul 3. t = 0 {caută elementul pivot de pe coloana k}

Pentru i=k+1,…, n execută

Dacă t<|aik| Atunci t=|aik|; lin:=i;

Pasul 4. Dacă link { interschimbă liniile lin și k pentru a pune pivotul pe diagonala } Atunci

Det=–Det{ la interschimbarea a două linii Det iși schimbă semnul}

Pentru j=k,…, n execută t=akj; akj=alinj; alinj=t;

Pentru j=1…, m execută t=bkj; bkj=blinj; blinj=t;

Pasul 5. Det=akk*Det { Înmulțește determinantul cu pivotul }

Pentru j=k+1,…, n execută akj=akj/akk { Împarte linia pivot la pivot}

Pentru j=1,…, m execută bkj=bkj/akk;

Pasul 6. Dacă kn Atunci

Pentru i=k+1,…, n execută

Pentru j=k+1,…, n execută aij=aij-aikakj;

Pentru j=1,…, m execută bij=bij-aikbkj

Faza substituției inverse

Pasul 7. Pentru k=n-1,…, 1 execută

Pentru j=1,…, m execută

t=0;

Pentru i=k+1,…, n execută t=t+akibij;

bkj=bkj-t

Programul de realizare a metodei lui Gauss

program metoda_eliminarii_lui_Gasuss;

uses crt;

const nmaxlin=10;

nmaxcol=10;

type vector=array[1..nmaxcol] of real;

matr_coef= array[1..nmaxlin,1..nmaxcol] of real;

var b,x:vector;

a:matr_coef;

n,m:integer;

prima,er:boolean;

procedure schimb(var a,b:real);

var temp:real;

begin

temp:=a;

a:=b;

b:=temp;

end;

procedure citeste(var a:matr_coef;var b:vector;var n,m: integer);

{ citeste valorile lui n si matricele a,b }

var i,j:integer;

begin

writeln;

repeat

write('Dati numarul de ecuatii n='); readln(n);

if not prima then clrscr else prima:=false;

m:=n

until n<nmaxlin;

if n>1

then begin

for i:=1 to n do

begin

writeln('ecuatia',i:3);

for j:=1 to n do

begin

write(j:3,':');

read(a[i,j])

end;

write('termenul liber b:'); read(b[i]);

readln

end;

end;

end;

procedure tipar_matrice(a:matr_coef;b:vector);

var i,j:integer;

begin

writeln('matricea initiala este:');

for i:=1 to n do

begin

for j:=1 to m do write(a[i,j]:8:3);

writeln(':',b[i]:8:3)

end;

writeln

end;

procedure afisare;

{ afisare solutii }

var i:integer;

begin

writeln('solutiile sunt: ');

for i:=1 to m do

write(x[i]:10:3);

writeln

end;

procedure gauss(a:matr_coef;b:vector;var x:vector;ncol:integer);

{ solutia dupa metoda de eliminare Gauss }

var c:matr_coef; { matrice de lucru cu nlin,ncol }

lib:vector; { vector de lucru cu ncol }

i,j,i1,k,l,n:integer;

temp,sum,t,curent,amax:real;

begin

n:=ncol; er:=false;

for i:=1 to n do

begin { copiere in matrice de lucru }

for j:=1 to n do c[i,j]:=a[i,j];

lib[i]:=b[i]

end;

for i:=1 to n-1 do { determina linia pivot avand elementul maxim }

begin { pe coloana i1 }

amax:=abs(c[i,i]); l:=i; i1:=i+1;

for j:=i1 to n do

begin { determina elementul maximal }

curent:=abs(c[j,i]);

if curent>amax then begin

amax:=curent; l:=j

end

end;

if amax=0.0

then er:=true

else begin

if l<>i

then begin

{ interschimba liniile sau ecuatiile pentru a pune

elementul maxim pe diagonala principala }

for j:=1 to n do schimb(c[l,j],c[i,j]);

schimb(lib[l],lib[i]);

end;

for j:=i1 to n do

begin

t:=c[j,i]/c[i,i]; { imparte linia pivot }

for k:=i1 to n do c[j,k]:=c[j,k]-t*c[i,k];

lib[j]:=lib[j]-t*lib[i]

end

end

end;

if c[n,n]=0.0

then er:=true

else begin

x[n]:=lib[n]/c[n,n]; i:=n-1;

{ faza substitutiei inverse }

repeat

sum:=0.0;

for j:=i+1 to n do sum:=sum+c[i,j]*x[j];

x[i]:=(lib[i]-sum)/c[i,i];

i:=i-1

until i=0

end;

if er then writeln('eroare! matrice singulara!');

end;

begin { programul de baza }

prima:=true; clrscr;

repeat

citeste(a,b,n,m); tipar_matrice(a,b);

if n>1 then begin

gauss(a,b,x,n);

if not er then afisare;

end

until n=0;

end.

3.2.3.1. Exemplu numeric. Rezolvarea unui circuit electric

Având circuitul de curent continuu din figura de mai jos cu datele numerice: , să se determine curenții din circuit utilizând teoremele lui Kirchhoff.

Circuit de curent continuu

Se scrie sistemul de ecuații după aplicarea teoremelor lui Kirchhoff:

(3.20)

Notăm curenții necunoscuți cu: și sistemul se rescrie:

(3.21)

Sistemul se va rezolva în continuare cu metoda eliminării versiunea Gauss. Metoda presupune transformarea sistemului inițial , într-un sistem echivalent (soluțiile vor fi similare cu soluțiile sistemului inițial) de forma a cărui rezolvare este foarte simplă, simplitatea fiind asigurată de matricea superior triunghiulară. Deci rezolvarea va avea două etape:

triunghiularizarea matricei A;

soluționarea sistemului echivalent.

pentru cazul particular de rezolvare a circuitului de curent continuu din exemplu matricea A augment (i se adaugă matricii sistemului A o coloană care corespunde cu vectorul b) se scrie:

Prima ecuație a sistemului este numită ecuația pivot iar elementul x1 se numește element pivot și urmărim eliminarea necunoscutei x1 din toate ecuațiile sistemului, cu excepția primei ecuații. Pentru aceasta, împărțim mai întâi prima linie la coeficientul ementului pivot x1 care în acest caz având valoarea 1 nu mai este nevoie de acest pas (dacă este zero, reordonăm și numerotăm ecuațiile).

Scădem prima ecuație înmulțită cu primul coeficient al celei de-a doua ecuații din această ecuație și, respectiv, înmulțită cu primul coeficient al celei de-a treia ecuații din aceasta din urmă [11].

(3.22)

Obținem astfel sistemul:

Ecuația identică se plasează ultima, și se alege ca ecuație pivot în continuare ecuația a treia și x2 element pivot și îl elimin din celelalte ecuații de mai jos:

(3.23)

Deci matricea A s-a transformat în matricea triunghilarizată A*. Se continuă cu faza substituției inverse (mersul înapoi) care presupune parcurgerea în sens invers a ecuațiilor sistemului cu matrice triunghiulară, rezultat în faza eliminării, și stabilirea soluției sistemului potrivit unui calcul recursiv:

(3.24)

După cum se observă, determinarea componentelor soluției are loc de la indici mari spre indici mici, fiecare nouă componentă depinzând în mod explicit numai de componentele determinate anterior [5].

Pașii metodei:

(3.25)

(3.26)

unde:

.

3.2.4. Metode bazate pe factorizarea LR a matricei coeficienților

3.2.4.1. Forma generala a metodei

Această clasă de metode directe se bazează pe factorizarea LR a matricei A a coeficienților sistemului liniar de ecuații (3.6) [8]. Forma generală a factorizării LR, prezentată anterior este definită de relația (1.56), unde se observă că matricea L este inferior triunghiulară iar R este superior triunghiulară. Înlocuind matricea A cu produsul L R, sistemul liniar (3.6) devine:

L R x = b (3.27)

Introducând notația y = R x, unde y este un vector coloană n – dimensional, sistemul (3.27) se poate scrie sub forma a două sisteme liniare de ordinul n:

L y = b (3.28)

R x = y (3.29)

Sistemele (3.28) și (3.29) se rezolvă simplu,.deoarece pentru prima dintre ele matricea coeficienților este inferior triunghiulară, iar pentru a doua este superior, triunghiulară.

(3.30)

y1= (3.31)

yi=, i = 2, 3, …, n (3.32)

(3.33)

(3.34)

xi= , i = n – 1, n – 2, …, 1 (3.35)

În concluzie, algoritmul general al metodei bazate pe factorizarea LR este:

se efectuează factorizarea LR a matricei A;

se rezolvă sistemul (3.28) folosind relațiile (3.31) – (3.32);

se rezolvă sistemul (3.29) folosind relațiile (3.34) – (3.35), soluția obținută fiind chiar cea a sistemului (3.29).

Metodele bazate pe factorizarea LR a matricei coeficienților necesită un timp de calcul de același ordin de mărime cu metoda triunghiularizării [9]. Dacă se utilizează la soluționarea repetată a sistemului de forma (3.6) pentru diferite valori ale termenilor liberi, atunci factorizarea matricei A se face o singură dată, la prima rezolvare. La soluționările următoare se rezolvă numai sistemele (3.28) și (3.29) (punctele b și c ale algoritmului). Evident, este necesară memorarea factorilor L și R care se poate face într-un singur tablou. Metodele particulare de soluționarea a sistemelor de ecuații liniare, abordate în subparagrafele următoare, se diferențiază între ele prin algoritmul de factorizare utilizat: LR Doolittle, LR Crout sau LR Cholesky.

3.2.4.2. Metoda bazată pe factorizarea LR Cholesky (metoda rădăcinii pătrate)

Metoda rădăcinii pătrate se aplică la sistemele care au matricea A a coeficienților simetrică. La punctul a) al algoritmului general din subparagraful 3.2.3.1. se folosește factorizarea LR Cholesky. Ea se caracterizează prin aceea că R = Lt. Calculele decurg conform celor prezentate în algoritmul general.

3.3. Metode iterative de soluționare

3.3.1. Metoda aproximațiilor succesive, versiunea Jacobi

Se consideră sistemul liniar de ecuații de ordinul n definit de relația (3.6), scris sub forma (3.2). Dacă se evidențiază distinct termenul din sumă care corespunde elementului diagonal, rezultă exprimarea (3.36).

aii xi +, i = 1, 2, …, n (3.36)

de unde se poate explicita expresia lui xi;

xi = , i = 1, 2, …, n (3.37)

În aceste condiții algoritmul metodei Jacobi este următorul:

se inițializează x cu x0 (indicele superior se referă la numărul iterației);

x0 = [x10 x20 x30 … xn0]t (3.38)

la un pas oarecare k, k = 1, 2, 3,… , al procesului iterativ de calcul, noile valori ale variabilelor se determină cu relațiile (3.39)

xik =, i = 1, 2, …, n (3.39)

condiția de terminare a procesului de calcul este de forma (3.40).

max{| xik – xik-1 |} , i = 1, 2, …, n (3.40)

Se poate formula condiția suficientă de convergență a procesului de calcul:

(3.41)

ceea ce se traduce prin cerința ca sistemul să fie cât mai bine condiționat. Practic, este suficient ca modulele elementelor diagonale ale matricei A a coeficienților să fie dominante (mai mari) în raport cu modulele elementelor nediagonale, pentru fiecare ecuație în parte [4]. De regulă, această ultimă condiție este îndeplinită în aplicațiile din domeniul aplicațiilor electrotehnice, în caz contrar, sistemul trebuie rearanjat înainte de începerea soluționării, pentru a asigura dominanța elementelor diagonale ale lui A. Rearanjarea se referă la schimbarea ordinii ecuațiilor și / sau a variabilelor.

Dacă convergența este asigurată, numărul de iterații depinde de modul de inițializare a soluției. în absența unor alte informații, inițializarea se poate face sub forma (3.42)

, i = 1, 2, …, n (3.42)

Exemplu numeric

După aplicarea asupra unui circuit funcționând în curent continuu a metodei curenților ciclici se ajunge la un sistem de ecuații liniare cu necunoscutele curenții de buclă notați cu x1, x2, x3:

; ;

;

Aproximația inițială:

Soluția la prima iterație se obține din:

Deci soluția sistemului la prima iterație este:

Soluția la iterația a doua se obține din soluția de la prima iterație:

Metoda converge și soluțiile se stabilizează după a doua iterație și se pot scrie ca fiind:

Se rezolvă sistemul în programul MathCad pentru a vedea toate soluțiile și pentru a stabili câte iterații sunt necesare pentru ca eroarea să fie într-o limită impusă apriori.

Evaluarea erorii în metoda aproximațiilor succesive în funcție de normele lui este:

Știind că procesul iterativ converge, câte iterații trebuie făcute pentru ca eroarea să fie mai mică decât 10-4.

Deci numărul de iterații care trebuie făcute pentru a se atinge precizia impusă este iar soluția la ultima iterație este soluția sistemului:

Metoda se poate aplica și pentru sisteme neliniare.

3.3.2. Metoda aproximațiilor succesive, versiunea Gauss – Seidel

Metoda Gauss – Seidel este foarte asemănătoare cu metoda Jacobi, singura / deosebire dintre ele fiind următoarea: la metoda Jacobi în partea dreaptă a relației (3.40) apar numai valorile "vechi" ale variabilelor (valorile de la iterația k – 1), iar la metoda Gauss – Seidel apar atât unele valori "vechi" (pentru variabilele care nu au fost încă recalculate), cât și unele valori "noi" (valorile de la iterația k, pentru variabilele care au fost recalculate deja).

Pentru metoda Gauss-Seidel relația (3.40) devine:

xik = , i = 1, 2, …, n (3.43)

celelalte elemente ale algoritmului fiind similare cu cele de la metoda Jacobi, prezentată în paragraful anterior. Rămân valabile în totalitate comentariile legate de convergență și de inițializarea soluției. De regulă, versiunea Gauss – Seidel a metodei aproximațiilor succesive aduce un plus de viteză de convergență fată de versiunea Jacobi, ceea ce determină utilizarea ei pe scară largă la soluționarea unor probleme din domeniul aplicațiilor electrotehnice [7].

Se menționează că dacă se elaborează un program de calcul pentru metoda aproximațiilor succesive, unde pentru vectorul coloană x se utilizează un singur tablou, atunci, în mod implicit, se lucrează cu versiunea Gauss – Seidel.

3.3.3. Metoda gradientului (celei mai abrupte coborâri)

Se consideră sistemul liniar de ecuații de ordinul n, scris sub forma (3.36), având matricea A a coeficienților simetrică. Pentru soluția aproximativă xk se definește vectorul rk = [r1k r2k r3k … rnk]t al reziduurilor, la fel ca la metoda Southwell.

În aceste condiții, algoritmul metodei gradientului (numită și metoda celei mai abrupte coborâri, datorita semnificației gradientului) este următorul:

se inițializează x cu x0 (indicele superior se referă la numărul iterației);

x0 = [x10 x20 x30 … xn0]t (3.44)

la un pas oarecare k, k = 1, 2, 3, …, al procesului iterativ de calcul, noile valori ale variabilelor se determină cu relația (3.45)

xk = xk-1 + hk-1 (3.45)

unde vectorul coloană al corecțiilor hk-1 are expresia (3.46)

hk-1 = (3.46)

vectorul coloană al reziduurilor rk-1 se calculează în modul definit la metoda Southwell,

rk-1 = A xk-1 – b (3.47)

iar scalarul are expresia definita de relatia:

(3.48)

condiția de terminare a procesului de calcul este de forma:

max{| hik-1 |} = max{| xik – xik-1 |} , i = 1, 2, …, n (3.49)

Rămân valabile în totalitate comentariile expuse în paragraful 3.3.1 (metoda Jacobi), referitoare la convergență și la inițializarea soluției [9].

Se prezintă si textul sursă al programului de calcul Gradient, cu observația că pentru toți vectorii coloană s-a utilizat câte un singur tablou iar apoi se prezintă o fereastra din programul Delphi cu un exemplu numeric rezolvat.

Program Gradient;

var i,j:integer; {variabile ajutatoare}

n:integer; {ordinul sistemului}

k:integer; {contorul de iteratie}

itermax:integer; {numarul de iteratii maxim admis}

eps:real; {eroarea maxima admisa}

s,u:real; {variabile ajutatoare}

alfa:real; {scalarul alfa}

hmax:real; {corectia maxima in modul}

a:array[1..8,1..8]of real; {matricea coeficientilor}

b:array[1..8]of real; {vectorul coloana al termenilor liberi}

x:array[1..8]of real; {vectorul coloana al necunoscutelor}

r:array[1..8]of real; {vectorul coloana al reziduurilor}

h:array[1..8]of real; {vectorul coloana al corectiilor}

t:array[1..8]of real; {vector coloana ajutator}

begin

writeln('Solutionarea sistemelor de ecuatii liniare prin metoda gadientului');

writeln('***************************************************************');

writeln;

write('Ordinul sistemului (maxim 8): n=');

readln(n);

writeln;

writeln(' Matricea A a coeficientilor');

for i:=1 to n do begin

writeln('Linia ',i);

for j:=1 to n do begin

write('A[',i,',',j,']=');

readln(a[i,j]);

end;

end;

writeln;

writeln(' Vectorul coloana B al termenilor liberi');

for i:=1 to n do begin

write('B[',i,']=');

readln(b[i]);

end;

write('Numarul maxim de iteratii: ');

readln(itermax);

writeln;

write('Eroarea maxima admisa: eps=');

readln(eps);

writeln;

{Afisare A}

writeln(' Matricea A a coeficientilor');

for i:=1 to n do begin

for j:=1 to n do write(a[i,j]:10:5);

if n<8 then writeln;

end;

writeln;

{Afisare B}

writeln(' Vectorul coloana B al termenilor liberi');

for i:=1 to n do writeln(b[i]:10:5);

writeln;

{Inceputul procesului de solutionare a sistemului de ecuatii}

writeln(' Initializarea solutiei');

for i:=1 to n do begin

write('X[',i,']=');

readln(x[i]);

end;

writeln;

{Inceputul procesului iterativ de calcul}

k:=0;

repeat

k:=k+1;

{Calculul reziduurilor}

s:=0;

for i:=1 to n do begin

r[i]:=-b[i];

for j:=1 to n do r[i]:=r[i]+a[i,j]*x[j];

s:=s+r[i]*r[i];

end;

{Calculul numitorului}

u:=0;

for i:=1 to n do begin

t[i]:=0;

for j:=1 to n do t[i]:=t[i]+a[j,i]*r[j];

u:=u+t[i]*r[i];

end;

{Calculul valorii lui alfa}

if abs(u)<1e-10 then begin

writeln('Nu se poate calcula alfa la pasul ',k,' !');

writeln;

writeln('Apasa Enter');

readln;

exit;

end;

alfa:=s/u;

{calculul corectiilor sa al noilor valori ale solutiei}

hmax:=0;

for i:=1 to n do begin

h[i]:=-alfa*r[i];

x[i]:=x[i]+h[i];

if abs(h[i])>hmax then hmax:=abs(h[i]);

end;

{Afisarea rezultatelor de la iteratia k}

writeln;

writeln(' Iteratia ',k);

writeln;

write(' ');

for i:=1 to n do write(i:9);

writeln;

write('X k-1');

for i:=1 to n do write(x[i]-h[i]:9:4);

writeln;

write('R ');

for i:=1 to n do write(r[i]:9:4);

writeln;

write('H ');

for i:=1 to n do write(h[i]:9:4);

writelnl

write('X k ');

for i:=1 to n do write(x[i]:9:4);

writeln;

writeln;

writeln('alfa=',alfa:1:6);

witeln;

writeln('hmax=',hmax:1:6,' eps=',eps:1:6);

writeln;

writeln('Apasa Enter');

readln;

until (hmax<=eps)or(k>itermax);

{Solutia}

if hmax<=eps then begin

writeln('Solutia sistemului, obtinuta in ',k:3,' iteratii');

writeln('––––––––––––––-');

writeln;

for i:=1 to n do writeln('X[',i,']=',x[i]:12:6);

writeln;

end;

else writeln('S-a depasit numarul maxim de iteratii');

writeln;

writeln('Apasa Enter');

readln;

end.

Figura 3.1. – Metoda gradientului de soluționare a unui sistem de ecuații liniare

Aplicarea metodei gradientului este justificata numai în cazul unei convergente mai puternice decat cea a metodei Gauss – Seidel, deoarece acum volumul de calcule pentru o iteratie este mult mai mare.

Considerăm o matrice pătratică A de ordinul n. Se cere să se determine toate cele n valori proprii ale acestei matrice și vectorii proprii corespunzători cu metoda Leverrier – Fadev.

Metoda utilizează noțiunea de urmă a unei matrice, definită ca sumă a elementelor diagonale ale matricei respective, în scopul găsirii valorilor coeficienților polinomului caracteristic .

Program Pol_Car;

type matrice=array[1..8,1..8]of real;

vector=array[1..9]of real;

var i,j,m:integer; {variabile ajutatoare}

n:integer; {ordinul matricei}

k:integer; {contorul de pasi}

A:matrice; {matricea A}

B:matrice; {matricea B}

D:matrice; {matrice ajutatoare}

C:vector; {coeficientii polinomului caracteristic}

{Procedura pentru afisarea matricelor}

procedure VizMat(Mat:matrice);

begin

for i:=1 to n do begin

for j:=1 to n do write(mat[i,j]:10:2);

if n<8 then writeln;

end;

writeln;

end;

begin

writeln(' Determinarea polinomului caracteristic al unei matrice');

writeln(' prin metoda Leverrier-Fadeev');

{citire A}

writeln;

write('Ordinul matricei A (maxim 8): n=');

readln(n);

writeln;

for i:=1 to n do begin

writeln('Linia ',i);

for j:=1 to n do begin

write('A[',i,',',j,']=');

readln(a[i,j]);

end;

end;

writeln;

{afisare A}

writeln('Matricea A');

writeln;

VizMat(a);

{init B}

k:=1;

c[1]:=1;

c[2]:=0;

for i:=1 to n do begin

for j:=1 to n do b[i,j]:=a[i,j];

c[2]:=c[2]-b[i,i];

end;

{tiparirea rezultatelor intermediare pentrul pasul 1}

writeln(' Pasul ',k);

writeln;

{Afisare B}

writeln('Matricea B ',k);

writeln;

VizMat(B);

{afisare coeficienti c1,c2}

writeln('C[',k,']=',c[k]);

writeln('C[',k+1,']=',c[k+1]);

writeln;

{inceputul procesului de calcul al coeficientilor ck, k=3,4..n+1}

for k:=2 to n do begin

{calculul sumei matriceale D=B+ck*I}

for i:=1 to n do begin

for j:=1 to n do d[i,j]:=b[i,j];

d[i,i]:=d[i,i]+c[k];

end;

{calculul matricei B=A*(B+ck*I)=A*D}

c[k+1]:=0;

for i:=1 to n do begin

for j:=1 to n do begin

b[i,j]:=0;

for m:=1 to n do b[i,j]:=b[i,j]+a[i,m]*d[m,j];

end;

c[k+1]:=c[k+1]-b[i,i];

end;

c[k+1]:=c[k+1]/k;

{tiparirea rezultatelor intermediare pentru pasul k}

writeln;

writeln(' Pasul ',k);

writeln;

{Afisare B}

writeln('Matricea B ',k);

writeln;

VizMat(B);

{Afisarea coeficientului c.k+1}

writeln('C[',k+1,']=',c[k+1]);

writeln;

writeln(' Apasa Enter');

readln;

end;

{rezultate finale}

writeln('Coeficientii polinomului caracteristic');

writeln('–––––––––––––');

writeln;

for i:=1 to n+1 do writeln('C[',i,']=',c[i]:1:6);

writeln;

writeln(' Apasa Enter');

readln;

end.

În continuare vom prezenta rezolvarea unei aplicații în MathCAD ale sistemelor de ecuații liniare, metode numerice de rezolvare.

Să se rezolve un sistem de n ecuații liniare cu n necunoscute utilizând funcția lsolve. Rezolvând un circuit electric cu metodele cunoscute ajungem la următorul sistem de ecuații:

unde matricea rezistențelor R respectiv vectorul tensiunilor electromotoare E sunt:

Deci vectorul curenților este:

Rezolvăm același sistem utilizând metoda clasică matricială. Avem relația matricială:

deci:

Valoarea functiei c(R) se numește numărul de condiționare al matricei R și cu cât este mai mare cu atât matricea este mai improprie inversării numerice.

Un alt mod de rezolvare este utilizând blocul de rezolvare Given – Find.

Rezolvând un circuit electric cu metodele cunoscute ajungem la următorul sistem de ecuații. În Mathcad pentru rezolvare ni se pune la dispoziție blocul Given – Find unde pentru aproximarea soluției cautate este necesar să se precizeze un punct inițial (în apropierea soluției căutate) cu care va începe procesul iterativ al funcției find [11].

Pentru sistemul obținut de noi alegem ca aproximații inițiale pentru soluții valorile de mai jos:

Cuvantul given precede ecuațiile sistemului iar semnul egal este egal boolean și se tasteaza CTRL=.

Deci vectorul curenților este:

Sau în loc de find se poate utiliza o altă funcție

.

Dacă aproximația inițială nu este aleasă în domeniul de convergență al procesului iterativ generat de funcțiile find sau minerr procesul nu va converge sau va converge către altă soluție.

Bibliografie

CRET, L., DUMA, D., Teoria circuitelor electrice. Culegere de probleme, UTPres, Cluj-Napoca, 2005.

CZIKER, A., Aplicatii ale metodelor numerice in electrotehnica, Editura Casa Cărții de Știință, , 2002.

*** Mathcad 2001, User's Guide, Mathsoft, , 2002.

Ioan, D. ș.a., Metode numerice în ingineria electrică, Ed. Matrix Rom București, 1998.

BRATCU, A., FILIPESCU, A., Metode numerice utilizate în analiza sistemelor, Ed. Matrix Rom București, 2004.

Comincioli, V., Analisi Numerica. Metodi, modelli, applicazioni, Mc.Grow-Hill Book Co. Milano, 1998.

Rice, J. Numerical methods-Software and Analysis, Mc.Grow-Hill Book , 1982.

TRÂMBITAȘ, R. Analiză Numerică, Ed. , , 2005.

Simion, E., Maghiar, T., Electrotehnică, EDP București, 1981.

Iordache, M., Dumitriu, L., Teoria moderna a circuitelor electrice, Vol. II Editura ALL, Bucuresti, 2000.

BUTA, A., Transportul și distribuția energiei electrice, Litografia Universității Tehnice, Timișoara, 1991.

EREMIA, M., Analiza asistată de calculator a regimurilor sistemelor electroenergetice, Editura Tehnică, București, 1985.

POEATA, A., ARIE, A., CRISAN, O., Transportul și distribuția energiei electrice, Editura Didactică și Pedagogică, București, 1981.

HORTOPAN, G., Aparate și echipamente electrice, Editura Didactica si Pedagogica,Bucuresti, 1999.

Bibliografie

CRET, L., DUMA, D., Teoria circuitelor electrice. Culegere de probleme, UTPres, Cluj-Napoca, 2005.

CZIKER, A., Aplicatii ale metodelor numerice in electrotehnica, Editura Casa Cărții de Știință, , 2002.

*** Mathcad 2001, User's Guide, Mathsoft, , 2002.

Ioan, D. ș.a., Metode numerice în ingineria electrică, Ed. Matrix Rom București, 1998.

BRATCU, A., FILIPESCU, A., Metode numerice utilizate în analiza sistemelor, Ed. Matrix Rom București, 2004.

Comincioli, V., Analisi Numerica. Metodi, modelli, applicazioni, Mc.Grow-Hill Book Co. Milano, 1998.

Rice, J. Numerical methods-Software and Analysis, Mc.Grow-Hill Book , 1982.

TRÂMBITAȘ, R. Analiză Numerică, Ed. , , 2005.

Simion, E., Maghiar, T., Electrotehnică, EDP București, 1981.

Iordache, M., Dumitriu, L., Teoria moderna a circuitelor electrice, Vol. II Editura ALL, Bucuresti, 2000.

BUTA, A., Transportul și distribuția energiei electrice, Litografia Universității Tehnice, Timișoara, 1991.

EREMIA, M., Analiza asistată de calculator a regimurilor sistemelor electroenergetice, Editura Tehnică, București, 1985.

POEATA, A., ARIE, A., CRISAN, O., Transportul și distribuția energiei electrice, Editura Didactică și Pedagogică, București, 1981.

HORTOPAN, G., Aparate și echipamente electrice, Editura Didactica si Pedagogica,Bucuresti, 1999.