Metode Numerice - Probleme de Seminar Si Lucrari de Laborator

211
Universitatea de Vest din Timi¸ soara Facultatea de Matematicˇa¸ siInformaticˇa METODE NUMERICE PROBLEME DE SEMINAR S ¸I LUCR ˇ ARI DE LABORATOR Simina Mari¸ s LilianaBrˇaescu Timi¸ soara 2007

Transcript of Metode Numerice - Probleme de Seminar Si Lucrari de Laborator

  • Universitatea de Vest din TimisoaraFacultatea de Matematica si Informatica

    METODE NUMERICE

    PROBLEME DE SEMINAR

    SI LUCRARI DE LABORATOR

    Simina Maris Liliana Braescu

    Timisoara

    2007

  • Introducere

    Procesul de restructurare al Invatamantului Superior din Romania si trecereaacestuia pe trei cicluri, a determinat la nivelul ntregii tari elaborarea de noiplanuri de nvatamant si de programe analitice adecvate.

    Metode numerice - Probleme de seminar si lucrari de laborator este unmaterial aditional la cursul de Metode numerice elaborat n acord cu noilecerinte, pe baza programei analitice conceputa la nivelul Departamentului deInformatica si aprobata n Consiliul Profesoral al Facultatii de Matematicasi Informatica de la Universitatea de Vest din Timisoara.

    Problemele si lucrarile de laborator prezentate n aceasta carte se adresea-za n primul rand studentilor de la Facultatea de Matematica si Informatica,fiind abordate toate temele din programa analitica, la nivelul studentilorSectiei de Informatica aflati n semestrul al cincilea de studiu, oferind exemplesi detalii referitoare la metodele numerice prezentate n curs.

    Lucrarea este structurata pe sapte capitole, primul dintre acestea fiindrezervat pentru prezentarea unui set de cunostinte minimale de programaren Maple. Capitolele 2-7 corespund capitolelor din cursul deMetode numericesi sunt organizate dupa cum urmeaza:

    breviar teoretic problema rezolvata probleme propuse implementarePrin aceasta lucrare, autorii pun la dispozitia cititorilor toate cunostintele

    necesare n vederea construirii de algoritmi si proceduri capabile sa ia caargument un obiect matematic si sa returneze un rezultat final.

    Autorii

  • Cuprins

    1 MapleV4 - scurta introducere 91.1 Reguli generale de introducere a comenzilor . . . . . . . . . . 91.2 Pachete de programe . . . . . . . . . . . . . . . . . . . . . . . 111.3 Constante, operatori si functii des utilizate . . . . . . . . . . . 121.4 Structuri de date . . . . . . . . . . . . . . . . . . . . . . . . . 131.5 Calcule cu matrice si vectori. Pachetul linalg . . . . . . . . . 161.6 Grafice . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 201.7 Elemente de programare . . . . . . . . . . . . . . . . . . . . . 22

    2 Rezolvarea sistemelor liniare 292.1 Metoda lui Gauss . . . . . . . . . . . . . . . . . . . . . . . . . 292.2 Factorizarea LU . . . . . . . . . . . . . . . . . . . . . . . . . . 512.3 Sisteme tridiagonale . . . . . . . . . . . . . . . . . . . . . . . 612.4 Factorizarea Cholesky . . . . . . . . . . . . . . . . . . . . . . 682.5 Factorizarea Householder . . . . . . . . . . . . . . . . . . . . . 742.6 Metoda Jacobi . . . . . . . . . . . . . . . . . . . . . . . . . . 852.7 Metoda Gauss-Seidel . . . . . . . . . . . . . . . . . . . . . . . 942.8 Metoda relaxarii succesive . . . . . . . . . . . . . . . . . . . . 100

    3 Ecuatii si sisteme de ecuatii neliniare 1073.1 Metoda punctului fix . . . . . . . . . . . . . . . . . . . . . . . 1073.2 Metoda lui Newton . . . . . . . . . . . . . . . . . . . . . . . . 113

    4 Interpolare polinomiala. Functii spline 1274.1 Polinomul lui Newton cu diferente divizate . . . . . . . . . . . 1274.2 Polinomul de interpolare Lagrange . . . . . . . . . . . . . . . 1414.3 Interpolare spline . . . . . . . . . . . . . . . . . . . . . . . . . 1444.4 Polinoame Bernstein . . . . . . . . . . . . . . . . . . . . . . . 155

    5 Derivare numerica 1595.1 Aproximarea derivatei prin diferente finite . . . . . . . . . . . 159

    7

  • 8 CUPRINS

    5.2 Aproximarea derivatei . . . . . . . . . . . . . . . . . . . . . . 163

    6 Integrare numerica 1656.1 Formule de tip Newton-Cotes . . . . . . . . . . . . . . . . . . 1656.2 Formule de tip Gauss . . . . . . . . . . . . . . . . . . . . . . . 168

    7 Ecuatii diferentiale 1737.1 Metoda diferentelor finite . . . . . . . . . . . . . . . . . . . . 1747.2 Metoda lui Taylor . . . . . . . . . . . . . . . . . . . . . . . . . 1797.3 Metoda Runge-Kutta . . . . . . . . . . . . . . . . . . . . . . . 1857.4 Metoda Adams-Bashforth . . . . . . . . . . . . . . . . . . . . 1917.5 Metoda Adams-Moulton . . . . . . . . . . . . . . . . . . . . . 1977.6 Metoda predictor-corector . . . . . . . . . . . . . . . . . . . . 1987.7 Probleme la limita liniare . . . . . . . . . . . . . . . . . . . . 2027.8 Metoda colocatiei si metoda celor mai mici patrate . . . . . . 209

  • Capitolul 1

    MapleV4 - scurta introducere

    Maple este un sistem de calcul algebric (CAS) dezvoltat de firma Maplesoft(http://www.maplesoft.com), care poate fi utilizat n:

    1. calcule simbolice;

    2. calcule numerice;

    3. programarea unor metode numerice;

    4. reprezentari grafice.

    In cele ce urmeaza, vom prezenta principalele elemente necesare n pro-gramarea unor metode numerice, corespunzatoare softului MapleV4.

    1.1 Reguli generale de introducere a comen-

    zilor

    Un document MapleV4 poate avea patru tipuri de campuri:

    1. comenzi Maple (introduse de catre utilizator);

    2. rezultate Maple (raspunsuri ale CAS-ului la comenzile introduse);

    3. grafice (raspunsuri ale CAS-ului);

    4. texte (introduse de catre utilizator).

    In continuare, vom prezenta cateva reguli de introducere a comenzilor.

    9

  • 10 CAPITOLUL 1. MAPLEV4 - SCURTA INTRODUCERE

    1. Orice comanda se termina cu ; (daca dorim sa afiseze rezultatul) sau :(daca nu dorim ca rezultatul sa fie afisat). De exemplu:

    > sin(Pi);

    0

    > 1+3:

    2. Asignarea se face cu := , iar dezasignarea se face prin asignarea numeluivariabilei. De exemplu, putem avea secventa:

    > x:= 7;

    x := 7

    > x:=x+1:

    > x;

    8

    > x:=x;

    x := x

    > x;

    x

    3. Comentariile sunt precedate de caracterul #. De exemplu:

    > x:=3; # se atribuie lui x valoarea 3

    x := 3

    4. Maple face diferenta ntre litere mici si litere mari:

    > x:=3; X:=5; a:=X-x;

    x := 3

    X := 5

    a := 2

    5. Secventele sunt scrise ntre paranteze rotunde, ( ), listele ntre parantezepatrate, [ ], iar multimile ntre acolade, {}. De exemplu:> secv:=(1,2,3); lista:=[2,1,2,3]; multime:={2,1,2,3};

    secv := 1, 2, 3

    lista := [2, 1, 2, 3]

    multime := {1, 2, 3}

  • 1.2. PACHETE DE PROGRAME 11

    6. Argumentele unei functii se pun ntre paranteze rotunde, () , iar indiciintre paranteze patrate, [] .

    > a:=cos(Pi); b:=lista[2];

    a := 1b := 1

    7. Procentul, % , face referire la ultima comanda executata anterior. Deexemplu:

    > a:=2;

    a := 2

    > b:=%+1;

    b := 3

    8. Daca rezultatul furnizat de Maple este identic cu comanda introdusa(Maple raspunde prin ecou la comanda), atunci aceasta arata ca Maplenu poate interpreta comanda introdusa. Pentru a remedia situatia,verificati daca ati introdus corect comanda sau daca nu cumva functiautilizata face parte dintr-un pachet care trebuie ncarcat n prealabil.

    > arctg(1); # o incercare de a calcula arctangenta de 1

    arctg(1)

    > arctan(1); # apelarea corecta a functiei arctangentapi

    4

    9. Pentru a nu retine eventuale atribuiri anterioare, este util ca pentrurezolvarea unei probleme noi sa ncepem cu instructiunea

    > restart;

    1.2 Pachete de programe

    Pachetele sunt colectii de functii care permit efectuarea de calcule specifice.Apelarea lor se face cu ajutorul comenzii with(nume_pachet). Pentru aapela o anumita functie dintr-un pachet, se foloseste sintaxa:pachet[functie](argumente)

    Printre cele mai utilizate pachete sunt:plots - pentru reprezentari grafice;DEtools - pentru rezolvarea ecuatiilor diferentiale;

  • 12 CAPITOLUL 1. MAPLEV4 - SCURTA INTRODUCERE

    linalg - pentru rezolvarea unor probleme de algebra liniara;

    student - pentru analiza matematica.

    De exemplu, la apelarea pachetului grafic, se obtine lista tuturor functiilorapelabile:

    > with(plots);

    Warning, the name changecoords has been redefined

    [animate, animate3d, changecoords, complexplot, complexplot3d,

    conformal, contourplot, contourplot3d, coordplot,

    coordplot3d, cylinderplot, densityplot, display,

    display3d, fieldplot, fieldplot3d, gradplot, gradplot3d,

    implicitplot, implicitplot3d, inequal, listcontplot,

    listcontplot3d, listdensityplot, listplot, listplot3d,

    loglogplot, logplot, matrixplot, odeplot, pareto,

    pointplot, pointplot3d, polarplot, polygonplot,

    polygonplot3d, polyhedraplot, replot, rootlocus,

    semilogplot, setoptions, setoptions3d, spacecurve,

    sparsematrixplot, sphereplot, surfdata, textplot,

    textplot3d, tubeplot]

    1.3 Constante, operatori si functii des uti-

    lizate

    Constantele folosite de Maple sunt:

    Constanta Semnificatiefalse falstrue adevaratgamma constanta lui Eulerinfinity +Catalan constanta lui CatalanFail valoare de adevar necunoscutaPi piI unitatea imaginaraNULL secventa vida

    Operatorii folositi frecvent sunt:

  • 1.4. STRUCTURI DE DATE 13

    Operator Sintaxa Semnificatie+, - a+b, a-b suma a + b (diferenta a b)* a*b, 2*a produsul a b, sau 2a/ a/b catul

    a

    b^, ** a^b, a**b puterea ab

    ! n! factorialul 1 2 ... nmax, min max(a,b,c), maximul (minimul) dintre a, b, c

    min(a,b,c)

    =,=, operatori booleeni:= f:=expr operatorul de asignare f = expr= a=b ecuatia a = b.. x=a..b a x b

    and, or, xor, operatori logiciimplies, not

    Functii folosite frecvent n Maple:

    Functie Sintaxa Semnificatiesin, cos, tan, cot sin(x) , ... functii trigonometrice

    arcsin, arctan, arccos arctan(x), ...ln, log10 ln(x), log10(x) logaritmi

    exp exp(x), exp(1) functia exponentialasqrt sqrt(x) radicalabs abs(x) modul

    1.4 Structuri de date: secvente, liste, multimi,

    siruri de caractere

    1.4.1 Secvente

    O secventa este o nsiruire de expresii, separate prin virgule. Exista maimulte moduri de a defini o secventa:

    a. direct:> s:=1,2,3,4; t:=(a,b,c);

    s := 1, 2, 3, 4t := a, b, c

    b. cu ajutorul functiei seq:> seq(3*x, x=2..7);

  • 14 CAPITOLUL 1. MAPLEV4 - SCURTA INTRODUCERE

    6, 9, 12, 15, 18, 21

    c. cu ajutorul unui ciclu for (vezi sectiunea 1.7.1)

    Cu ajutorul functiilor min si max se poate calcula minimul, respectiv maxi-mul unei secvente.

    > min(s); max(s,2,15);

    3

    15

    1.4.2 Liste

    O lista este o secventa de expresii, scrisa ntre paranteze patrate, [ ]. Deexemplu, putem avea lista:

    > ll:=[1,2,5,2,4,2,7,2,a,2,c];

    ll := [1, 2, 5, 2, 4, 2, 7, 2, a, 2, c]

    Putem afla numarul de operanzi din lista cu ajutorul functiei nops:

    > nops(ll);

    11

    Al n-lea element din lista poate fi afisat cu una din comenzile op(n,ll) saull[n]:

    > aa:=ll[3]; bb:=op(9,ll);

    aa := 5

    bb := a

    Functia member(elem, ll) returneaza true daca elementul respectiv se aflan lista ll, si false n caz contrar:

    > member(d, ll);

    false

    Lista vida este desemnata prin []:

    > lista:=[];

    lista := [ ]

    Se poate adauga un element nou la lista ll astfel: [op(ll),elem]. Deexemplu:

    > lista:=[op(lista), d,e,f];

    lista := [d, e, f ]

    Se poate sterge al n-lea element din lista ll astfel: subsop(n=NULL,ll). Deexemplu:

    > lista:=subsop(2=NULL, lista);

    lista := [d, f ]

  • 1.4. STRUCTURI DE DATE 15

    1.4.3 Multimi

    O multime este o secventa de expresii, scrisa ntre acolade, {}, n care fiecareelement figureaza o singura data. De exemplu:> ll:={1,2,5,2,4,2,7,2,a,2,c};

    ll := {1, 2, 4, 5, 7, a, c}Adaugarea unui nou element la multime, sau stergerea elementului de pepozitia n se face la fel ca la liste.Multimea vida este desemnata prin {}.Reuniunea a doua multimi se face utilizand operatorul union:> s:={1,2,3} : t:={2,3,4} :

    > s union t;

    {1, 2, 3, 4}Intersectia a doua multimi se realizeaza cu ajutorul operatorului intersect:> s intersect t;

    {2, 3}Diferenta a doua multimi se realizeaza utilizand operatorul minus:> s minus t;

    {1}> s minus s;

    {}

    1.4.4 Siruri de caractere

    Sirurile de caractere sunt delimitate de apostrof invers,, dupa cum urmeaza:> acesta este un sir;

    acesta este un sirSirurile de caractere se pot concatena cu ajutorul comenzii cat. De exemplu,putem avea:> i:=4;

    i := 4> cat( Valoarea lui i este , i);

    V aloarea lui i este 4Atentie! La concatenarea unui sir de cifre, se obtine un sir de caractere, nuun numar:> a:=cat(5,7,9); b:=52;

    a :=579b := 52

    > whattype(a); # afla tipul expresiei a

    symbol> whattype(b); # afla tipul expresiei b

  • 16 CAPITOLUL 1. MAPLEV4 - SCURTA INTRODUCERE

    integer> a+b;

    a :=579+52> whattype(a+b); # afla tipul expresiei a+b

    symbol

    1.5 Calcule cu matrice si vectori. Pachetul

    linalg

    Cu ajutorul cuvantului-cheie array se pot defini vectori si matrice.Un vector se defineste n urmatorul mod:

    > v:=array(1..dim_vect);

    Elementele unui vector se pot defini unul cate unul, sau printr-un ciclufor (vezi sectiunea 1.7.1):> v:=array(1..4);

    v := array(1..4, [ ])> v[1]:=a; v[2]:=b; v[3]:={a,b}; v[4]:=3;

    v1 := av2 := bv3 := {a, b}v4 := 3

    > evalm(v); # evalueaza valoarea lui v

    [a, b, {a, b}, 3]O matrice se defineste astfel:

    > M:=array(1..nr_rand, 1..nr_col);

    Elementele unei matrice se pot defini unul cate unul, sau printr-un ciclufor (vezi sectiunea 1.7.1):> M:=array(1..2,1..2);

    M := array(1..2, 1..2, [ ])> M[1,1]:=1: M[1,2]:=a: M[2,1]:=3: M[2,2]:={}:

    > evalm(M); # evalueaza valoarea lui M[1 a3 { }

    ]

    Un alt mod de a defini matrice si vectori, precum si de a efectua operatiispecifice cu aceste obiecte, este folosirea pachetului linalg. Pachetul linalgse ncarca astfel:

    >with(linalg);

  • 1.5. CALCULE CU MATRICE SI VECTORI. PACHETUL LINALG 17

    Warning, the protected names norm and trace have been

    redefined and unprotected

    [BlockDiagonal, GramSchmidt, JordanBlock, LUdecomp, QRdecomp,

    Wronskian, addcol, addrow, adj, adjoint, angle, augment,

    backsub, band, basis, bezout, blockmatrix, charmat,

    charpoly,cholesky, col, coldim, colspace, colspan,

    companion, concat, cond, copyinto, crossprod, curl,

    definite, delcols, delrows, det, diag, diverge, dotprod,

    eigenvals, eigenvalues, eigenvectors, eigenvects,

    entermatrix, equal, exponential, extend, ffgausselim,

    fibonacci, forwardsub, frobenius, gausselim, gaussjord,

    geneqns, genmatrix, grad, hadamard, hermite, hessian,

    hilbert, htranspose, ihermite, indexfunc, innerprod,

    intbasis, inverse, ismith, issimilar, iszero, jacobian,

    jordan, kernel, laplacian, leastsqrs, linsolve, matadd,

    matrix, minor, minpoly, mulcol, mulrow, multiply, norm,

    normalize, nullspace, orthog, permanent, pivot,

    potential, randmatrix, randvector, rank, ratform, row,

    rowdim, rowspace, rowspan, rref, scalarmul, singularvals,

    smith, stack, submatrix, subvector, sumbasis, swapcol,

    swaprow, sylvester, toeplitz, trace, transpose,

    vandermonde, vecpotent, vectdim, vector, wronskian]

    O matrice se defineste cu comanda matrix:matrix(nr_randuri, nr_coloane, lista_elem sau fc_generatoare)

    Astfel, putem avea:> a:=matrix(3,2,[1,2,3,4,5,6]);

    a :=

    1 23 4

    5 6

    dar si> f:=(i,j)->i+j; # functia generatoare a elem matricei

    f := (i, j) i+ j> b:=matrix(2,3,f); # adica b[i,j]=f(i,j)

    b :=

    [2 3 43 4 5

    ]Elementul aij se scrie a[i,j]. Astfel, pentru matricea a din exemplul ante-rior, putem avea:> a[3,1];

    5

  • 18 CAPITOLUL 1. MAPLEV4 - SCURTA INTRODUCERE

    Un vector se defineste cu ajutorul comenzii vector:vector(nr_elem, lista_elem sau fc_generatoare)

    Astfel, putem avea:v:=vector([2,4,8,2]);

    v := [2, 4, 8, 2]sauf:=x-> 2*x+1;

    f := x 2x+ 1w:=vector(5,f); # adica w[i]=f(i)

    w := [3, 5, 7, 9, 11]Elementul i al vectorului v, vi, se scrie v[i]. Astfel, pentru vectorul v

    din exemplul anterior, putem avea:> v[3];

    8Redam mai jos cele mai utilizate functii din pachetul linalg, mpreuna cu

    descrierea lor. Pentru mai multe detalii referitoare la aceste functii, precumsi la celelalte functii din pachetul linalg, se poate consulta pagina de helpreferitoare la pachetul linalg.

    Functie Descriere

    addcol(A,c1,c2,m) nlocuieste coloana c2 a matricei A cum*c1+c2

    addrow(A,r1,r2,m) nlocuieste linia r2 a matricei A cum*r1+r2

    adj(A), adjoint(A) calculeaza matricea adjuncta a matriceiA

    angle(u,v) calculeaza unghiul vectorilor u si vaugment(A,B) concateneaza (alatura) matricile A si B

    pe orizontalabacksub(U,b) rezolva sistemul Ux=b, prin substitutie

    inversa, unde U este o matrice superiortriunghiulara

    band(b,n) construieste o matrice n x n care arepe diagonala principala elementele vec-torului b, iar celelalte elemente suntnule

    continuare pe pagina urmatoare

  • 1.5. CALCULE CU MATRICE SI VECTORI. PACHETUL LINALG 19

    Pachetul linalg - continuare

    Functie Descrierecholesky(A) efectueaza descompunerea Cholesky a

    matricei Acol(A,i), col(A,i..k) extrage coloana i, respectiv coloanele i

    pana la k, din matricea Acoldim(A) returneaza numarul de coloane ale ma-

    tricei Acrossprod(u,v) returneaza produsul vectorial al vecto-

    rilor u si vdelcols(A,r..s) sterge coloanele de la r la s din ma-

    tricea Adelrows(A,r..s) sterge liniile de la r la s din matricea Adet(A) calculeaza determinantul matricei Adiverge(f) calculeaza divergenta vectorului fdotprod(u,v) calculeaza produsul scalar al vectorilor

    u si vexponential(A) calculeaza eA

    extend(A,m,n,x) adauga m linii si n coloane matricei A,initializate cu x

    forwardsub(L,b) rezolva sistemul Lx=b prin substitutienainte, unde L este o matrice inferiortriunghiulara

    gausselim(A) efectueaza eliminarea gaussiana cusemipivot asupra matricei A

    geneqns(A,x) genereaza un sistem de ecuatii pornindde la matricea A si vectorul necunos-cutelor x

    genmatrix(sist, var) genereaza matricea coeficientilor sis-temului sist, in raport cu multimeavariabilelor var

    grad(expr, vect) calculeaza gradientul expresiei expr, infunctie de variabilele vect

    inverse(A) calculeaza inversa matricei Amatadd(A,B,c1,c2) calculeaza c1*A+c2*Bminor(r,c) calculeaza minorul de ordin (r,c)

    (elimina linia r si coloana c) din ma-tricea A

    continuare pe pagina urmatoare

  • 20 CAPITOLUL 1. MAPLEV4 - SCURTA INTRODUCERE

    Pachetul linalg - continuare

    Functie Descrieremulcol(A,c,expr) multiplica coloana c a matricei A cu ex-

    presia exprmulrow(A,r,expr) multiplica linia r a matricei A cu expre-

    sia exprmultiply(A,B) efectueaza nmultirea matricelor A si Bnorm(A) calculeaza norma matricei Anormalize(v) calculeaza versorul vectorului vrank(A) calculeaza rangul matricei Arow(A,i), row(A,i..j) extrage linia i, respectiv liniile de la i

    la j, ale matricei Arowdim(A) returneaza numarul de linii din ma-

    tricea Ascalarmult(A,s) nmulteste toate elementele matricei A

    cu scalarul sstack(A,B) concateneaza matricele A si B pe verti-

    calasubmatrix(A,r1..r2,c1..c2) extrage o submatrice a matricei A, ntre

    liniile r1, r2, si coloanele c1, c2subvector(A,r1..r2) extrage un subvector al vectorului A, de

    la rangul r1 la rangul r2swapcol(A,c1,c2) interschimba coloanele c1 si c2 ale ma-

    tricei Aswaprow(A,r1,r2) interschimba liniile r1 si r2 ale matricei

    A

    trace(A) calculeaza urma matricei Avectdim(v) returneaza dimensiunea vectorului v

    1.6 Grafice

    Graficul unei functii se realizeaza folosind comanda plot, a carei sintaxa este

    plot(functie, x=x_min..x_max, y_min..y_max)

    unde argumentul y_min..y_max este optional.

    De exemplu, putem avea:

    > plot(sin(x), x=-5..5);

  • 1.6. GRAFICE 21

    1

    0.5

    0

    0.5

    1

    4 2 2 4

    x

    > plot(cos(x)^2, x=-5..5);

    0

    0.2

    0.4

    0.6

    0.8

    1

    4 2 2 4

    x

    > plot([sin(x),cos(x)^2], x=-5..5);

  • 22 CAPITOLUL 1. MAPLEV4 - SCURTA INTRODUCERE

    1

    0.5

    0

    0.5

    1

    4 2 2 4

    x

    Mai multe detalii despre grafice se pot gasi accesand pagina de help refe-ritoare la instructiunea plot, sau la pachetul plots.

    1.7 Elemente de programare

    1.7.1 Conditionarea si ciclarea

    A. ConditionareaSintaxa unei instructiuni conditionale este

    if CONDITIE then EXPRESIE

    [ elif CONDITIE then EXPRESIE ]

    [ else EXPRESIE ]

    fi

    Instructiunile puse ntre paranteze patrate, [ ], sunt optionale.

    De exemplu, putem avea secventa:> if a if x

  • 1.7. ELEMENTE DE PROGRAMARE 23

    [ for CONTOR ] [ from EXPR ] [ by EXPR ] [ to EXPR ]

    [ while EXPR ]

    do INSTRUCTIUNI od;

    sau

    [ for CONTOR ] [ in EXPR ] [ while EXPR ]

    do INSTRUCTIUNI od;

    unde:- from indica punctul de plecare n iteratie (daca este omis, valoarea sa

    implicita este 1);- by indica pasul contorului (daca este omis, se considera implicit ca are

    valoarea 1);- to indica punctul de oprire a iteratiei (daca este omis, se considera

    implicit ca are valoarea + si se obtine o bucla infinita);- while indica o expresie booleana, care trebuie sa poata fi evaluata ca

    adevarata sau falsa;- in indica elementele succesive ale expresiei EXPR.

    De exemplu, pentru a scrie toate numerele pare de la 6 la 100 putemfolosi:> for i from 6 by 2 to 100 do print(i) od;

    Cu ajutorul buclei for se pot defini secvente, liste, multimi, vectori saumatrice.> s:=NULL;

    for i from 1 to 3 do s:=s,2*i+1 od; # definirea unei secvente

    s :=s := 3s := 3, 5s := 3, 5, 7

    > l:=[];

    for i from 1 to 4 do l:=[op(l),i^2] od; # definirea unei liste

    l := [ ]l := [1]l := [1, 4]l := [1, 4, 9]l := [1, 4, 9, 16]

    > v:=vector(3); # definirea vectorului

    for i from 1 to 3 do v[i]:=i^3-i^2+1 od; # definirea elem vect

    evalm(v); # vizualizarea vectorului

    v := array(1..3, [ ])

  • 24 CAPITOLUL 1. MAPLEV4 - SCURTA INTRODUCERE

    v1 := 1v2 := 5v3 := 19[1, 5, 19]

    M:=array(1..3,1..4); # definirea matricei

    M := array(1..3, 1..4, [ ])> for i from 1 to 3 do # definirea elem matricei

    for j from 1 to 4 do

    M[i,j]:=i^j

    od;

    od;

    > evalm(M); 1 1 1 12 4 8 16

    3 9 27 81

    Putem afisa elemetele unei liste (secvente, multimi, matrice, vector) astfel:

    > lista:=[3,2,4,5,1]:

    > for i in lista do print(i) od;

    Mai multe detalii despre instructiunile de conditionare si de ciclare se potgasi accesand pagina de help referitoare la acestea.

    1.7.2 Functii si proceduri

    O functie poate fi definita cu ajutorul operatorului ->. Putem defini functiide o variabila sau functii de mai multe variabile.> f:=x->x^2+1;

    f := x x2 + 1> g:=(x,y)->x^2+y;

    g := (x, y) x2 + y> f(3);

    10> g(3,4);

    13> g(4,3);

    19

    O procedura este un grup de instructiuni, variabile si constante. Sintaxaeste:

  • 1.7. ELEMENTE DE PROGRAMARE 25

    proc (ARGUMENTE)

    local VARIABILE_LOCALE;

    global VARIABILE_GLOBALE;

    options OPTIUNI;

    description SIR_DE_CARACTERE;

    INSTRUCTIUNI;

    end;

    O procedura returneaza ultimul rezultat obtinut. Pentru a forta returnareaunui alt rezultat, se foloseste RETURN. De asemenea, pentru a returna unmesaj de eroare, se foloseste ERROR.De exemplu, putem defini procedura:

    > modul:=proc(a) if a modul(-3);

    3Un alt exemplu de procedura este urmatorul:

    > ec2:=proc(a,b,c)

    local delta,x1,x2;

    description Rezolvarea ecuatiei de gradul 2;

    delta:=b^2-4*a*c;

    if delta>0 then

    x1:=(-b+sqrt(delta))/(2*a);

    x2:=(-b-sqrt(delta))/(2*a);

    RETURN(x1,x2);

    elif delta=0 then RETURN(-b/(2*a));

    else

    RETURN(ecuatia nu ere solutii reale);

    fi;

    end:

    care produce urmatoarele rezultate:> ec2(1,6,9); # ecuatia x^2+6*x+9=0

    3> ec2(1,2,9); # ecuatia x^2+2*x+9=0

    ecuatia nu are solutii reale> ec2(1,2,-3); # ecuatia x^2+2*x-3=0

    1,3

  • 26 CAPITOLUL 1. MAPLEV4 - SCURTA INTRODUCERE

    Pentru a defini tipul unui argument, se foloseste sintaxa argument::tip.De exemplu, sa luam urmatoarea procedura si situatiile care pot aparea:

    > # procedura care returneaza determinantul unei matrice

    > determinant:=proc(A) RETURN(det(A)) end:

    > determinant(2);

    Error, (in linalg:-det) expecting a matrix

    Procedura determinant se poate imbunatati astfel:

    > determinant1:=proc(A)

    if not type(A, matrix)

    then ERROR(argumentul trebuie sa fie matrice!!!)

    fi;

    RETURN(det(A))

    end:

    care produce urmatorul rezultat:

    > determinant1(2);

    Error, (in determinant1) argumentul trebuie sa fie matrice!!!

    Se mai poate defini argumentul A ca fiind de tipul matrice:

    > determinant3:=proc(A::matrix) RETURN(det(A)) end:

    si se obtine urmatorul rezultat:

    > determinant3(2);

    Error, invalid input: determinant3 expects its 1st argument, A,

    to be of type matrix, but received 2

    Mai multe detalii despre tipurile existente se pot gasi accesand pagina dehelp (cuvantul cheie este type).

    Un alt exemplu este procedura rdc, procedura pentru calculul lui1x:

    > rdc:=proc(x)

    if x

  • 1.7. ELEMENTE DE PROGRAMARE 27

    rdc := proc(x) if x < 0 then ERROR(numar negativ!)

    elif x = 0 then RETURN()else simplify(1/(x (1/2))) end if end proc

    > rdc(-1);

    Error, (in rdc) numar negativ!

    > rdc(0);

    > rdc(4);

    1

    2

    Pentru a putea urmari executia unei proceduri, se foloseste debug, iarpentru a stopa urmarirea, se foloseste undebug. De exemplu, putem avea:

    > f:=proc(a,b)

    local y,z;

    y:=a+b/2;

    z:=1/y;

    RETURN(y+z)

    end;

    f := proc(a, b)

    local y, z;

    y := a+ 1/2 b; z := 1/yRETURN(y + z) end proc

    > debug(f);

    f> f(2,4);

    {--> enter f, args = 2, 4

    y := 4

    z :=1

    4 f(0,1);

  • 28 CAPITOLUL 1. MAPLEV4 - SCURTA INTRODUCERE

    {--> enter f, args = 0, 1

    y :=1

    2z := 2

    f(10,20);

    401

    20

    Alte detalii despre functii si proceduri, precum si despre optiunile debugsi undebug, puteti gasi pe paginile de help referitoare la acestea.

  • Capitolul 2

    Rezolvarea sistemelor liniare

    In acest capitol vom prezenta metode de rezolvare a sistemelor liniare detip Cramer (numarul de ecuatii este egal cu numarul de necunoscute, sideterminantul matricei sistemului este nenul):

    a11x1 + a12x2 + . . .+ a1nxn = b1

    a21x1 + a22x2 + . . .+ a2nxn = b2

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

    an1x1 + an2x2 + . . .+ annxn = bn

    (2.1)

    n care aij si bi sunt numere reale date, i = 1 . . . n, j = 1 . . . n, iar x1, x2, . . . , xnsunt numere reale necunoscute.

    Sistemul (2.1) se poate scrie matriceal sub forma:

    Ax = b

    unde: A = (aij)i,j=1,n , b = (b1, b2, . . . , bn)T , x = (x1, x2, . . . , xn)

    T .Daca matricea A este nesingulara, sistemul Ax = b are solutie unica:

    x = A1b.

    Deoarece n cele mai multe cazuri matricea A are numar mare de linii sicoloane, iar calculul matricei A1 este dificil si acumuleaza erori, se impunmetode directe si metode iterative pentru rezolvarea acestor sisteme.

    2.1 Metoda lui Gauss

    2.1.1 Breviar teoretic

    Metoda lui Gauss presupune transformarea sistemului Ax = b ntr-un sistemsuperior triunghiular, si apoi rezolvarea acestuia prin substitutie inversa.

    29

  • 30 CAPITOLUL 2. REZOLVAREA SISTEMELOR LINIARE

    Constructia sistemului superior triunghiular se face astfel: la pasulk se elimina xk din ecuatiile k + 1, ..., n, prin nmultirea ecuatiei k cu

    mik = aikakk

    (elementul akk se numeste pivot) si adunarea acestora la ecuatia i (i > k).In functie de alegerea pivotului, exista urmatoarele variante ale metodei

    lui Gauss:

    1. metoda lui Gauss clasica - n care la fiecare pas, pivotul este ele-mentul akk, k = 1, n;

    2. metoda lui Gauss cu semipivot - n care la fiecare pas, se alegeca pivot elementul aik maxim n valoare absoluta pe coloana, pentrui > k, permutandu-se linia k cu linia i;

    3. metoda lui Gauss cu pivot total - n care la fiecare pas, se alege capivot elementul maxim atat pe linie, cat si pe coloana, pentru i > k,j > k, permutandu-se linia k cu linia i si coloana k cu coloana j;

    In acest fel, sistemul (2.1) se reduce la forma superior triunghiulara

    a11 a12 ... a1,n2 a1,n1 a1,n0 a22 ... a2,n2 a2,n1 a2,n... ... ... ... ... ...0 0 ... 0 an1,n1 an1,n0 0 ... 0 0 ann

    x1x2...xn1xn

    =

    b1b2...

    bn1bn

    (2.2)

    iar rezolvarea sistemului (2.2) se face prin substitutie inversa:

    xn =bnann

    (2.3)

    xk =

    (bk

    nj=k+1

    akj xj) 1akk

    , k = n 1, n 2, . . . , 1

    Observatia 2.1.1. Cu ajutorul eliminarii gaussiene se poate determina siinversa unei matrice. Redam n continuare algoritmul de aflare a inverseiunei matrice A.

    1. generarea matricei B prin concatenarea matricelor A (de dimensiunen) cu matricea In

  • 2.1. METODA LUI GAUSS 31

    2. pentru i = 1, n

    m = Bii

    pentru j = 1, 2n

    Bij =Bijm

    pentru j = 1, n, j 6= im1 = Bjipentru k = 1, 2n

    Bjk = Bjk m1 Bik3. prin stergerea primelor n coloane ale matricei B astfel transformate, se

    obtine inversa matricei A

    2.1.2 Probleme rezolvate

    Exercitiul 2.1.1. Sa se rezolve urmatorul sistem folosind cele trei varianteale eliminarii Gauss:

    x+ y + z = 62x y + 3z = 9x+ 4y + z = 12.

    Matricea sistemului este

    A =

    1 1 12 1 3

    1 4 1

    ,

    iar A este matricea sa extinsa:

    A = (A, b) =

    1 1 1 62 1 3 9

    1 4 1 12

    .

    Deoarece numarul ecuatiilor este egal cu cel al necunoscutelor si

    detA = 3 6= 0,sistemul este compatibil determinat (de tip Cramer), si deci metoda eliminariia lui Gauss este aplicabila.

    In continuare, pentru a efectua operatiile asupra matricei extinse a sis-temului vom nota linia i cu Li, iar coloana j cu Cj .

    Rezolvare utilizand metoda lui Gauss clasicaA. Constructia sistemului superior triunghiular

  • 32 CAPITOLUL 2. REZOLVAREA SISTEMELOR LINIARE

    Pasul 1

    pivot: a11 = 1 m21 = 2

    1= 2

    m31 = 11= 1

    1 1 1 62 1 3 9

    1 4 1 12

    L2L2+m21L1L3L3+m31L1

    1 1 1 60 3 1 3

    0 3 0 6

    Pasul 2

    pivot: a22 = 3 m32 = 33 = 1

    1 1 1 60 3 1 30 3 0 6

    L3L3+m32L2

    1 1 1 60 3 1 3

    0 0 1 3

    In acest moment am ajuns la un sistem de forma Ax = b, echivalent cusistemul initial, n care matricea A este superior triunghiulara, unde:

    A =

    1 1 10 3 1

    0 0 1

    , x =

    xy

    z

    , b =

    63

    3

    .

    B. Rezolvarea sistemului superior triunghiularPrin metoda substitutiei inverse, avem:

    z =3

    1

    y =1

    3(3 1 z)

    x =1

    1(6 1 y 1 z),

    de unde obtinem solutia sistemului: x = 1, y = 2, z = 3.

    Rezolvare cu metoda lui Gauss cu semipivotA. Constructia sistemului superior triunghiular

  • 2.1. METODA LUI GAUSS 33

    Pasul 1

    Ca pivot se ia elementul ai1 de modul maxim de pe coloana 1. Incazul nostru, pivotul este a12, deci se permuta linia 1 cu linia 2, si sefac zerouri pe coloana 1 pentru i > 1:

    1 1 1 62 1 3 91 4 1 12

    L2L1

    2 1 3 91 1 1 6

    1 4 1 12

    2 1 3 91 1 1 6

    1 4 1 12

    L2L2

    12L1

    L3L312L1 2 1 3 90 3

    21

    232

    0 92

    12

    152

    Pasul 2

    Ca pivot se ia elementul ai2 de modul maxim de pe coloana 2, pentrui 2. In cazul nostru, pivotul este a32, deci se permuta linia 2 cu linia3 si se fac zerouri pe coloana 2, pentru i > 2:

    2 1 3 90 32

    12

    32

    0 92

    12

    152

    L3L2

    2 1 3 90 9

    21

    2152

    0 32

    12

    32

    2 1 3 90 9

    21

    2152

    0 32

    12

    32

    L3L3 13L2

    2 1 3 90 9

    21

    2152

    0 0 131

    In acest moment am ajuns la un sistem de forma Ax = b, echivalent cusistemul initial, unde matricea A este superior triunghiulara, iar:

    A =

    2 1 30 9

    21

    2

    0 0 13

    , x =

    xy

    z

    , b =

    915

    2

    1

    .

    B. Rezolvarea sistemului superior triunghiular se face ca si n cazul metodeilui Gauss clasice, si conduce la solutia x = 1, y = 2, z = 3.

    Rezolvare cu metoda lui Gauss cu pivot totalA. Constructia sistemului superior triunghiular

    Pasul 1

  • 34 CAPITOLUL 2. REZOLVAREA SISTEMELOR LINIARE

    ca pivot se alege elementul aij de modul maxim pentru i, j 1.In cazul nostru pivotul este a32, deci se permuta linia 3 cu linia 1, sicoloana 2 cu coloana 1:

    1 1 1 62 1 3 91 4 1 12

    L3L1

    1 4 1 122 1 3 9

    1 1 1 6

    1 4 1 122 1 3 9

    1 1 1 6

    C2C1

    4 1 1 121 2 3 9

    1 1 1 6

    pentru corectitudinea rezultatului final este necesar ca, ori de cateori se permuta coloanele matricei extinse, sa se permute si elementelecorespunzatoare ale vectorului x. Astfel, avem:

    x =

    xy

    z

    x2x1

    yx

    z

    n final, obtinem: 4 1 1 121 2 3 9

    1 1 1 6

    L2L2+

    14L1

    L3L314L1 4 1 1 120 9

    4134

    120 3

    434

    3

    Pasul 2

    ca pivot se alege elementul aij de modul maxim pentru i, j 2.Deoarece pivotul este a23, se permuta coloana 3 cu coloana 2:

    4 1 1 120 9

    4134

    120 3

    434

    3

    C3C2

    4 1 1 120 13

    494

    120 3

    434

    3

    x =

    yx

    z

    x3x2

    yz

    x

    4 1 1 120 13

    494

    120 3

    434

    3

    L3L3 313L2

    4 1 1 120 13

    494

    120 0 3

    13313

  • 2.1. METODA LUI GAUSS 35

    In acest moment am ajuns la un sistem de forma Ax = b, echivalent cusistemul initial, unde matricea A este superior triunghiulara, iar:

    A =

    4 1 10 13

    494

    0 0 313

    , x =

    yz

    x

    , b =

    1212

    313

    .

    B. Rezolvarea sistemului superior triunghiular se face ca si n cazul metodeilui Gauss clasice, si conduce la solutia x = 1, z = 3, y = 2.

    Exercitiul 2.1.2. Sa se gaseasca inversa matricei

    A =

    2 2 32 1 1

    3 1 2

    .

    RezolvareConsideram matricea B obtinuta prin concatenarea matricei A cu ma-

    tricea unitate I3:

    B =

    2 2 3 1 0 02 1 1 0 1 0

    3 1 2 0 0 1

    .

    Folosind metoda eliminarii a lui Gauss, transformam matricea B dupa cumurmeaza:

    2 2 3 1 0 02 1 1 0 1 03 1 2 0 0 1

    L1 1B11 L1

    1 1 32 12 0 02 1 1 0 1 0

    3 1 2 0 0 1

    1 1 32 12 0 02 1 1 0 1 0

    3 1 2 0 0 1

    L2L2B21L1L3L3B31L1

    1 1 32 12 0 00 1 2 1 1 0

    0 2 523

    20 1

    1 1 32 12 0 00 1 2 1 1 0

    0 2 523

    20 1

    L2 1B22 L2

    1 1 32 12 0 00 1 2 1 1 0

    0 2 523

    20 1

    1 1 32 12 0 00 1 2 1 1 0

    0 2 523

    20 1

    L1L1B12L2L3L3B32L2

    1 0 12 12 1 00 1 2 1 1 0

    0 0 32

    12

    2 1

    1 0 12 12 1 00 1 2 1 1 0

    0 0 32

    12

    2 1

    L3 1B33 L3

    1 0 12 12 1 00 1 2 1 1 0

    0 0 1 13

    43

    23

  • 36 CAPITOLUL 2. REZOLVAREA SISTEMELOR LINIARE

    1 0 12 12 1 00 1 2 1 1 0

    0 0 1 13

    43

    23

    L1L1B13L3L2L2B23L3

    1 0 0 13 13 130 1 0 1

    353

    43

    0 0 1 13

    43

    23

    Inversa matricei A va fi matricea C, obtinuta prin stergerea primelor 3coloane ale matricei B:

    C =

    13 13 131

    353

    43

    13

    43

    23

    Intr-adevar, se verifica usor ca

    A C = C A = I3.

    2.1.3 Probleme propuse

    Exercitiul 2.1.3. Sa se rezolve urmatoarele sisteme, folosind cele trei vari-ante ale metodei lui Gauss:

    a)

    x+ 2y + z = 13x y + 5z = 14x+ y z = 2

    b)

    x+ y + z + t = 03x 2y z + t = 8x 2y z + 4t = 1x+ y z + 2t = 1

    Exercitiul 2.1.4. Sa se gaseasca solutia sistemelor anterioare, calculandinversa matricei A a sistemului, si efectuand nmultirea A1b.

    2.1.4 Implementare

    A. AlgoritmAlgoritmii pentru cele 3 metode sunt asemanatori, diferenta dintre ei

    aparand (asa cum se poate vedea si din exemplul rezolvat) n modul derezolvare a eliminarii Gauss.

    Date de intrare: un sistem de ecuatii (scris ca multime de ecuatii)Date de iesire: solutia sistemului

    Algoritmul consta din urmatoarele etape:

    1. generarea matricei extinse a sistemului, A = (aij)i=1,n,j=1,n+1

    n= numarul de ecuatii (numarul de linii ale matricei A);

  • 2.1. METODA LUI GAUSS 37

    2. a) eliminarea Gauss pentru metoda lui Gauss clasica

    - pentru k = 1, n 1- daca akk = 0, atunci se cauta r pentru care akr 6= 0,r = k + 1, n si se schimba linia k cu linia r;

    - daca toti akr = 0, r = k + 1, n atunci se returneaza eroare;

    - pentru i = k + 1, n

    m = aikakk

    , unde akk 6= 0;- pentru j = k, n

    aij = aij +m akj ;b) eliminarea Gauss pentru metoda lui Gauss cu semipivot

    - pentru k = 1, n 1 se cauta elementul de modul maxim pe linie,i.e. daca |akr| > |akk|, r = k + 1, n, se schimba linia k cu linia r- pentru i = k + 1, n

    m = aikakk

    , unde akk 6= 0;- pentru j = k, n

    aij = aij +m akj ;c) eliminarea Gauss pentru metoda lui Gauss cu pivot total

    - pentru k = 1, n 1 se cauta elementul de modul maxim pe linie sicoloana, i.e. daca |apr| > |akk|, p, r = k + 1, n, se schimba coloanap cu coloana k si linia r cu linia k

    - pentru i = k + 1, n

    m = aikakk

    , unde akk 6= 0;- pentru j = k, n

    aij = aij +m akj ;3. rezolvarea sistemului superior triunghiular prin substitutie inversa

    xn =an,n+1ann

    ,

    - pentru i = n 1, 1

    xi =1

    aii

    (ai,n+1

    nj=i+1

    aijxj

    ).

  • 38 CAPITOLUL 2. REZOLVAREA SISTEMELOR LINIARE

    B. Programe MAPLE si rezultateDeoarece diferitele variante ale metodei lui Gauss se deosebesc doar prin

    modul n care se realizeaza eliminarea Gauss, n cele ce urmeaza am imple-mentat separat cele trei variante de eliminare, folosind procedurile cgauss,spgauss, tpgauss. Aceste proceduri vor fi folosite apoi ca optiuni n proce-dura finala gauss.

    restart: with(linalg):

    cgauss:=proc(A::matrix)

    local A1, A2, n, k, r, i, m, j;

    n:=rowdim(A);

    A1:=A;

    A2:=delcols(A1,n+1..n+1);

    if(det(A2)=0) then ERROR(sistemul nu are solutie unica!) fi;

    for k from 1 to n-1 do

    if A1[k,k]=0 then

    for r from k+1 to n

    while A1[k,r]=0 do r=r+1 od;

    if r>n then ERROR(sistemul nu are solutie unica!)

    else A1:=swaprow(A1,k,r);

    fi;

    fi;

    for i from k+1 to n do

    m:=A1[i,k]/A1[k,k];

    for j from k to n+1 do

    A1[i,j]:=A1[i,j]-m*A1[k,j];

    od;

    od;

    od;

    RETURN(evalm(A1));

    end:

    spgauss:=proc(A::matrix)

    local A1, A2, n, k, r, i, m, j, mx;

    n:=rowdim(A);

    A1:=A;

    A2:=delcols(A1,n+1..n+1);

    if(det(A2)=0) then ERROR(sistemul nu are solutie unica!) fi;

    for k from 1 to n-1 do

    mx:=k;

  • 2.1. METODA LUI GAUSS 39

    for r from k to n do

    if (abs(A1[r,k])>abs(A1[k,k])) then mx:=r

    fi;

    od;

    if mxk then A1:=swaprow(A1,k,mx); fi;

    for i from k+1 to n do

    m:=A1[i,k]/A1[k,k];

    for j from k to n+1 do

    A1[i,j]:=A1[i,j]-m*A1[k,j];

    od;

    od;

    od;

    RETURN(evalm(A1));

    end:

    tpgauss:=proc(A::matrix)

    local A1, A2, n, k, r, i, m, j, mx, my, l;

    n:=rowdim(A);

    A1:=A;

    A2:=delcols(A1,n+1..n+1);

    l:=[seq(i), i=1..n];

    if(det(A2)=0) then ERROR(sistemul nu are solutie unica!) fi;

    for k from 1 to n-1 do

    mx:=k; my:=k;

    for i from k to n do

    for j from k to n do

    if (abs(A1[i,j])>abs(A1[k,k])) then mx:=i; my:=j;

    fi;

    od;

    od;

    if mxk then A1:=swaprow(A1,k,mx); fi;

    if myk then

    A1:=swapcol(A1,k,my);

    l:=subsop(k=l[my],my=l[k],l);

    fi;

    for i from k+1 to n do

    m:=A1[i,k]/A1[k,k];

    for j from k to n+1 do

    A1[i,j]:=A1[i,j]-m*A1[k,j];

    od;

    od;

  • 40 CAPITOLUL 2. REZOLVAREA SISTEMELOR LINIARE

    od;

    RETURN(evalm(A1),l);

    end:

    gauss:=proc(eqn::set(equation), opt::symbol)

    local A,A1,l,n,r,k,i,m,j,s,x,rez;

    l:=[op(indets(eqn))];

    n:=nops(l);

    A:=genmatrix(eqn, l, flag);

    if opt=clasic then A1:=cgauss(A);

    elif opt=semipivot then A1:=spgauss(A);

    elif opt=totalpivot then

    rez:=tpgauss(A);

    A1:=rez[1];

    l:=[seq(l[rez[2][i]],i=1..n)];

    else ERROR(optiunile sunt: clasic, semipivot sau totalpivot);

    fi;

    x[n]:=A1[n,n+1]/A1[n,n];

    for i from n-1 by -1 to 1 do

    s:=0;

    for j from i+1 to n do

    s:=s+A1[i,j]*x[j];

    od;

    x[i]:=1/A1[i,i]*(A1[i,n+1]-s);

    od;

    RETURN(seq(l[i]=x[i],i=1..n));

    end:

    Observatia 2.1.2. Instructiunea indets(set_eq) returneaza multimea nede-terminatelor sistemului set_eq. Deoarece ordinea elementelor acestei multiminu este neaparat aceeasi cu ordinea nedeterminatelor din prima ecuatie a sis-temului, pot aparea diferente ntre rezultatele furnizate cu ajutorul coduluiMAPLE si rezultatele calculate pe hartie. Desi matricea sistemului generatacu ajutorul instructiunii indets nu este ntotdeauna aceeasi cu matriceasistemului scrisa pe hartie, rezultatele furnizate de program vor fi aceleasi(eventual ordinea solutiilor va fi schimbata).

    Observatia 2.1.3. Pentru a urmari executia unei proceduri, se folosesteinstructiunea debug. In cazul programelor din exemplele de mai sus, sepoate folosi urmatorul set de instructiuni:

    debug(cgauss):

  • 2.1. METODA LUI GAUSS 41

    debug(spgauss):

    debug(tpgauss):

    debug(gauss):

    Redam mai jos, cu titlu de exemplu, rezultatul urmaririi procedurilorgauss, cgauss, spgauss si tpgauss pentru acelasi sistem de ecuatii.

    > gauss({x+y+z=6,2*x-y+3*z=9,x+4*y+z=12},clasic);{--> enter gauss, args = {x+y+z = 6, 2*x-y+3*z = 9, x+4*y+z = 12},clasic

    l := [x, y, z]

    n := 3

    A :=

    1 1 1 62 1 3 9

    1 4 1 12

    {--> enter cgauss, args = A

    n := 3

    A1 := A

    A2 :=

    1 1 12 1 3

    1 4 1

    m := 2

    A1 2, 1 := 0

    A1 2, 2 := 3A1 2, 3 := 1

    A1 2, 4 := 3m := 1

    A1 3, 1 := 0

    A1 3, 2 := 3

    A1 3, 3 := 0

    A1 3, 4 := 6

    m := 1A1 3, 2 := 0

    A1 3, 3 := 1

    A1 3, 4 := 3

  • 42 CAPITOLUL 2. REZOLVAREA SISTEMELOR LINIARE

    A1 :=

    1 1 1 60 3 1 3

    0 0 1 3

    x3 := 3

    s := 0

    s := 3

    x2 := 2

    s := 0

    s := 2

    s := 5

    x1 := 1

    gauss({x+y+z=6,2*x-y+3*z=9,x+4*y+z=12},semipivot);{--> enter gauss, args = {x+y+z = 6, 2*x-y+3*z = 9, x+4*y+z = 12},semipivot

    l := [x, y, z]

    n := 3

    A :=

    1 1 1 62 1 3 9

    1 4 1 12

    {--> enter spgauss, args = A

    n := 3

    A1 := A

    A2 :=

    1 1 12 1 3

    1 4 1

    mx := 1

    mx := 2

    A1 :=

    2 1 3 91 1 1 6

    1 4 1 12

    m :=1

    2

    A1 2, 1 := 0

    A1 2, 2 :=3

    2

    A1 2, 3 :=12

  • 2.1. METODA LUI GAUSS 43

    A1 2, 4 :=3

    2

    m :=1

    2

    A1 3, 1 := 0

    A1 3, 2 :=9

    2

    A1 3, 3 :=12

    A1 3, 4 :=15

    2mx := 2

    mx := 3

    A1 :=

    2 1 3 90

    9

    2

    12

    15

    2

    03

    2

    12

    3

    2

    m :=1

    3

    A1 3, 2 := 0

    A1 3, 3 :=13

    A1 3, 4 := 1

  • 44 CAPITOLUL 2. REZOLVAREA SISTEMELOR LINIARE

    gauss({x+y+z=6,2*x-y+3*z=9,x+4*y+z=12},totalpivot);{--> enter gauss, args = {x+y+z = 6, 2*x-y+3*z = 9, x+4*y+z = 12},totalpivot

    l := [x, y, z]

    n := 3

    A :=

    1 1 1 62 1 3 9

    1 4 1 12

    {--> enter tpgauss, args = A

    n := 3

    A1 := A

    A2 :=

    1 1 12 1 3

    1 4 1

    l := [1, 2, 3]

    mx := 1

    my := 1

    mx := 2

    my := 1

    mx := 2

    my := 3

    mx := 3

    my := 2

    A1 :=

    1 4 1 122 1 3 9

    1 1 1 6

    A1 :=

    4 1 1 121 2 3 9

    1 1 1 6

    l := [2, 1, 3]

    m :=14

    A1 2, 1 := 0

    A1 2, 2 :=9

    4

  • 2.1. METODA LUI GAUSS 45

    A1 2, 3 :=13

    4

    A1 2, 4 := 12

    m :=1

    4

    A1 3, 1 := 0

    A1 3, 2 :=3

    4

    A1 3, 3 :=3

    4

    A1 3, 4 := 3

    mx := 2

    my := 2

    mx := 2

    my := 3

    A1 :=

    4 1 1 12

    013

    4

    9

    412

    03

    4

    3

    43

    l := [2, 3, 1]

    m :=3

    13

    A1 3, 2 := 0

    A1 3, 3 :=3

    13

    A1 3, 4 :=3

    13

  • 46 CAPITOLUL 2. REZOLVAREA SISTEMELOR LINIARE

    l := [y, z, x]

    x3 := 1

    s := 0

    s :=9

    4x2 := 3

    s := 0

    s := 3

    s := 4

    x1 := 2

    A := matrix([[1, 1, 1, 6], [2, -1, 3, 9], [1, 4, 1, 12]]);

    A :=

    1 1 1 62 1 3 9

    1 4 1 12

    > gausselim(A); 1 1 1 60 3 1 3

    0 0 1 3

    > backsub(%);

    [1, 2, 3]

    Redam programul Maple pentru determinarea inversei unei matrice, pre-cum si un exemplu:

    invmatrix:=proc(A::matrix)

    local n,A1,i,B,j,m,m1,k,C;

    n:=rowdim(A);

    if coldim(A) n then

    ERROR(Matricea trebuie sa fie patratica!)

    fi;

    if det(A)=0 then

    ERROR(Matricea trebuie sa fie nesingulara!)

  • 2.1. METODA LUI GAUSS 47

    fi;

    A1:=matrix(n,n,0);

    for i from 1 to n do

    A1[i,i]:=1

    od;

    B:=concat(A,A1);

    for i from 1 to n do

    m:=B[i,i];

    for j from 1 to 2*n do B[i,j]:=B[i,j]/m od;

    for j from 1 to n do

    if ij then

    m1:=B[j,i];

    for k from 1 to 2*n do

    B[j,k]:=B[j,k]-m1*B[i,k];

    od;

    fi;

    od;

    od;

    evalm(B); # rezultatul eliminarii gaussiene

    C:=delcols(B,1..n); # inversa matricei A

    end:

    > A := matrix([[2, 1, 1], [2, -1, 3], [1, 4, 1]]):

    > debug(invmatrix):

    > C:=invmatrix(A);

    {--> enter invmatrix, args = A

    n := 3

    A1 :=

    0 0 00 0 0

    0 0 0

    A1 1, 1 := 1

    A1 2, 2 := 1

    A1 3, 3 := 1

    B :=

    2 1 1 1 0 02 1 3 0 1 0

    1 4 1 0 0 1

    m := 2

    B1, 1 := 1

    B1, 2 :=1

    2

  • 48 CAPITOLUL 2. REZOLVAREA SISTEMELOR LINIARE

    B1, 3 :=1

    2

    B1, 4 :=1

    2

    B1, 5 := 0

    B1, 6 := 0

    m1 := 2

    B2, 1 := 0

    B2, 2 := 2B2, 3 := 2

    B2, 4 := 1B2, 5 := 1

    B2, 6 := 0

    m1 := 1

    B3, 1 := 0

    B3, 2 :=7

    2

    B3, 3 :=1

    2

    B3, 4 :=12

    B3, 5 := 0

    B3, 6 := 1

    m := 2B2, 1 := 0

    B2, 2 := 1

    B2, 3 := 1

    B2, 4 :=1

    2

    B2, 5 :=12

    B2, 6 := 0

    m1 :=1

    2

    B1, 1 := 1

    B1, 2 := 0

    B1, 3 := 1

  • 2.1. METODA LUI GAUSS 49

    B1, 4 :=1

    4

    B1, 5 :=1

    4

    B1, 6 := 0

    m1 :=7

    2

    B3, 1 := 0

    B3, 2 := 0

    B3, 3 := 4

    B3, 4 :=94

    B3, 5 :=7

    4

    B3, 6 := 1

    m := 4

    B3, 1 := 0

    B3, 2 := 0

    B3, 3 := 1

    B3, 4 :=916

    B3, 5 :=7

    16

    B3, 6 :=1

    4m1 := 1

    B1, 1 := 1

    B1, 2 := 0

    B1, 3 := 0

    B1, 4 :=13

    16

    B1, 5 :=316

    B1, 6 :=14

    m1 := 1B2, 1 := 0

    B2, 2 := 1

    B2, 3 := 0

  • 50 CAPITOLUL 2. REZOLVAREA SISTEMELOR LINIARE

    B2, 4 :=116

    B2, 5 :=116

    B2, 6 :=1

    4

    1 0 013

    16

    316

    14

    0 1 0116

    116

    1

    4

    0 0 1916

    7

    16

    1

    4

    C :=

    13

    16

    316

    14

    116

    116

    1

    4916

    7

    16

    1

    4

    multiply(A,C); 1 0 00 1 0

    0 0 1

    > multiply(C,A); 1 0 00 1 0

    0 0 1

  • 2.2. FACTORIZAREA LU 51

    2.2 Factorizarea LU

    2.2.1 Breviar teoretic

    Fie sistemul compatibil determinat

    Ax = b. (2.4)

    Factorizarea LU presupune descompunerea matricei A ntr-un produs de ma-trice L U , unde

    L =

    11 0 . . . 021 22 . . . 0. . . . . . . . . . . .n1 n2 . . . nn

    U =

    11 12 . . . 1n0 22 . . . 2n. . . . . . . . . . . .0 0 . . . nn

    . (2.5)

    Aceasta descompunere este posibila daca toti determinantii de colt aimatricei A sunt nenuli.

    Pentru a asigura unicitatea descompunerii, trebuie precizate n elementeale matricei L sau U . In mod traditional, se specifica ii sau ii; daca ii = 1atunci factorizarea LU se numeste factorizare Doolittle, iar daca ii = 1se numeste factorizare Crout.

    Astfel, rezolvarea sistemului (2.4) se reduce la rezolvarea sistemelor tri-unghiulare

    Ly = b (2.6)

    cu solutia

    y1 =b111

    yi =

    (bi

    i1j=1

    ijyj

    ) 1ii

    , i = 2, 3, . . . , n(2.7)

    si

    Ux = y (2.8)

    cu solutia

    xn =ynnn

    xi =

    (yi

    nj=i+1

    ijxj

    ) 1ii

    , i = 2, 3, . . . , n.(2.9)

  • 52 CAPITOLUL 2. REZOLVAREA SISTEMELOR LINIARE

    2.2.2 Problema rezolvata

    Exercitiul 2.2.1. Sa se determine solutia sistemului urmator, folosind fac-torizarea LU:

    x+ y z = 22x y + z = 1x+ 3y 2z = 5.

    Sistemul se scrie n forma matriceala:

    Ax = b,

    unde

    A =

    1 1 12 1 1

    1 3 2

    , x =

    xy

    z

    , b =

    21

    5

    .

    Deoarece

    1 6= 0 , 1 12 1

    = 3 6= 0 ,1 1 12 1 11 3 2

    = 3 6= 0 ,rezulta ca matricea A este nesingulara si are toti determinantii de colt nenuli,deci se poate folosi factorizarea LU pentru rezolvarea acestui sistem.Rezolvare folosind factorizarea Crout

    A. Factorizarea CroutPresupunem ca

    A =

    1 1 12 1 1

    1 3 2

    =

    11 0 021 22 0

    31 32 33

    1 12 130 1 23

    0 0 1

    ,

    si ne propunem sa determinam coeficientii lij , ujk. Pentru aceasta, folosimdefinitia nmultirii matricelor. Astfel, avem:

    a11 = 11 1 11 = 1a12 = 11 12 12 = 1a13 = 11 13 13 = 1a21 = 21 1 21 = 2a22 = 21 12 + 22 1 22 = 3a23 = 21 13 + 22 23 23 = 1a31 = 31 1 31 = 1a32 = 31 12 + 32 1 32 = 2a33 = 31 13 + 32 23 + 33 1 33 = 1

  • 2.2. FACTORIZAREA LU 53

    sau

    L =

    1 0 02 3 0

    1 2 1

    , U =

    1 1 10 1 1

    0 0 1

    .

    B. Rezolvarea sistemelor triunghiulare

    Pentru rezolvarea sistemului initial, avem de rezolvat doua sisteme tri-ungiulare:

    1 0 02 3 01 2 1

    y1y2

    y3

    =

    21

    5

    ,

    a carui solutie este

    y =

    y1y2

    y3

    =

    21

    1

    ,

    si respectiv: 1 1 10 1 1

    0 0 1

    xy

    z

    =

    21

    1

    ,

    a carui solutie este

    x =

    xy

    z

    =

    12

    1

    .

    Rezolvare folosind factorizarea Doolittle

    A. Factorizarea Doolittle

    Presupunem ca

    A =

    1 1 12 1 1

    1 3 2

    =

    1 0 021 1 0

    31 32 1

    11 12 130 22 23

    0 0 33

    si ne propunem sa determinam coeficientii lij , jk, la fel ca si n exemplul

  • 54 CAPITOLUL 2. REZOLVAREA SISTEMELOR LINIARE

    precedent. Astfel avem:

    a11 = 1 11 11 = 1a12 = 1 12 12 = 1a13 = 1 13 13 = 1a21 = 21 11 21 = 2a22 = 21 12 + 1 22 22 = 3a23 = 21 13 + 1 23 23 = 3a31 = 31 11 31 = 1a32 = 31 12 + 32 22 32 = 2

    3a33 = 31 13 + 32 23 + 1 33 33 = 1

    sau

    L =

    1 0 02 1 0

    1 23

    1

    , U =

    1 1 10 3 3

    0 0 1

    .

    B. Rezolvarea sistemelor triunghiulare

    Pentru rezolvarea sistemului initial, avem de rezolvat doua sisteme tri-ungiulare:

    1 0 02 1 01 2

    31

    y1y2

    y3

    =

    21

    5

    ,

    a carui solutie este

    y =

    y1y2

    y3

    =

    23

    1

    ,

    si respectiv: 1 1 10 3 3

    0 0 1

    xy

    z

    =

    23

    1

    ,

    a carui solutie este

    x =

    xy

    z

    =

    12

    1

    .

  • 2.2. FACTORIZAREA LU 55

    2.2.3 Probleme propuse

    Exercitiul 2.2.2. Sa se gaseasca solutiile urmatoarelor sisteme, folosind celedoua variante ale factorizarii LU:

    a)

    x+ 2y + z = 13x y + 5z = 14x+ y z = 2

    b)

    3x+ y 2z = 1x+ y + z = 62x y + 4z = 7

    2.2.4 Implementare

    A. AlgoritmDate de intrare: un sistem de ecuatiiDate de iesire: solutia sistemului

    Algoritmul consta din urmatoarele etape:

    1. generarea matricei A a sistemului, si a vectorului coloana b

    n = numarul de linii ale matricei A (numarul de ecuatii ale sistemului)2. a) factorizarea Crout

    pentru i = 1, n

    ii = 1

    pentru i = 1, n

    pentru j = 1, i

    ij = aij j1k=1

    ikkj

    pentru j = i+ 1, n

    ij =1

    ii

    (aij

    i1k=1

    ikkj

    )

    b) factorizarea Doolittle

    pentru i = 1, n

    ii = 1

    pentru i = 1, n

    pentru j = 1, i 1

  • 56 CAPITOLUL 2. REZOLVAREA SISTEMELOR LINIARE

    ij =1

    jj

    (aij

    ik=1

    ikkj

    )

    pentru j = i, n

    ij = aij i1k=1

    ikkj

    3. Rezolvarea celor doua sisteme triunghiulare

    y1 =b111

    pentru i = 2, n

    yi =

    (bi

    i1j=1

    ijyj

    ) 1ii

    xn =ynnn

    pentru i = 2, n

    xi =

    (yi

    nj=i+1

    ijxj

    ) 1ii

    B. Programe MAPLE si rezultateDeoarece cele doua variante ale descompunerii LU difera doar prin modul

    de factorizare a matricei sistemului, am implementat separat cele doua vari-ante de factorizare: LUcrout si LUdoolittle, dupa care le-am folosit caoptiuni n procedura finala LUsist.

    restart: with(linalg):

    LUcrout:=proc(A::matrix)

    local a1,n,l,u,i,s,j,k;

    n:=rowdim(A);

    a1:=A;

    if a1[1,1]=0 then

    ERROR(factorizarea LU nu este aplicabila!);

    fi;

    for i from n by -1 to 2 do

    if det(a1)0 then a1:=delrows(delcols(a1,i..i),i..i);

    else ERROR(factorizarea LU nu este aplicabila!);

  • 2.2. FACTORIZAREA LU 57

    fi;

    od;

    l:=matrix(n,n,0); u:=matrix(n,n,0);

    for i from 1 to n do

    u[i,i]:=1;

    od;

    for i from 1 to n do

    for j from 1 to i do

    s:=0; for k from 1 to j-1 do s:=s+l[i,k]*u[k,j]; od;

    l[i,j]:=A[i,j]-s;

    od;

    for j from i+1 to n do

    s:=0; for k from 1 to i-1 do s:=s+l[i,k]*u[k,j]; od;

    u[i,j]:=1/l[i,i]*(A[i,j]-s);

    od;

    od;

    RETURN(evalm(l), evalm(u));

    end:

    LUdoolittle:=proc(A::matrix)

    local a1,n,l,u,i,s,j,k;

    n:=rowdim(A);

    a1:=A;

    if a1[1,1]=0 then

    ERROR(factorizarea LU nu este aplicabila!);

    fi;

    for i from n by -1 to 2 do

    if det(a1)0 then a1:=delrows(delcols(a1,i..i),i..i);

    else ERROR(factorizarea LU nu este aplicabila!);

    fi;

    od;

    l:=matrix(n,n,0); u:=matrix(n,n,0);

    for i from 1 to n do

    l[i,i]:=1;

    od;

    for i from 1 to n do

    for j from 1 to i-1 do

    s:=0; for k from 1 to i do s:=s+l[i,k]*u[k,i]; od;

    l[i,j]:=1/u[j,j]*(A[i,j]-s);

    od;

    for j from i to n do

  • 58 CAPITOLUL 2. REZOLVAREA SISTEMELOR LINIARE

    s:=0; for k from 1 to i-1 do s:=s+l[i,k]*u[k,i]; od;

    u[i,j]:=A[i,j]-s;

    od;

    od;

    RETURN(evalm(l), evalm(u));

    end:

    LUsist:=proc(l::set(equation), opt::symbol)

    local lst, eqm, A, b, n, lu, L, U,i,s,j,aux, rez, rfin;

    eqm:=genmatrix(l, [op(indets(l))], flag);

    lst:=indets(l);

    n:=nops(lst);

    A:=delcols(eqm,n+1..n+1);

    b:=col(eqm,n+1);

    if opt=Crout then

    lu:=LUcrout(A);

    elif opt=Doolittle then

    lu:=LUdoolittle(A);

    else ERROR(optiunile sunt: Crout sau Doolittle)

    fi;

    L:=lu[1];

    U:=lu[2];

    for i from 1 to n do

    s:=0; for j from 1 to i-1 do s:=s+L[i,j]*aux[j] od;

    aux[i]:=1/L[i,i]*(b[i]-s)

    od;

    for i from n by -1 to 1 do

    s:=0; for j from i+1 to n do s:=s+U[i,j]*rez[j] od;

    rez[i]:=1/U[i,i]*(aux[i]-s)

    od;

    RETURN(seq(lst[i]=rez[i], i=1..n));

    end:

    debug(LUsist);

    LUsist({x+y-z=2,2*x-y+z=1,x+3*y-2*z=5}, Doolittle);

    {--> enter LUsist, args = {x+y-z = 2, 2*x-y+z = 1, x+3*y-2*z = 5},Doolittle

    eqm :=

    1 1 1 21 2 1 12 1 3 5

  • 2.2. FACTORIZAREA LU 59

    lst := {z, x, y}n := 3

    A :=

    1 1 11 2 12 1 3

    b := [2, 1, 5]

    lu :=

    1 0 01 1 0

    213

    1

    , 1 1 10 3 0

    0 0 1

    L :=

    1 0 01 1 0

    213

    1

    U :=

    1 1 10 3 0

    0 0 1

    s := 0

    aux 1 := 2

    s := 0

    s := 2aux 2 := 3

    s := 0

    s := 4

    s := 3

    aux 3 := 2

    s := 0

    rez 3 := 2

    s := 0

    s := 0

    rez 2 := 1

    s := 0

    s := 1

    s := 3

    rez 1 := 1

  • 60 CAPITOLUL 2. REZOLVAREA SISTEMELOR LINIARE

    {--> enter LUsist, args = {x+y-z = 2, 2*x-y+z = 1, x+3*y-2*z = 5},Crout

    eqm :=

    1 1 1 21 2 1 12 1 3 5

    lst := {z, x, y}n := 3

    A :=

    1 1 11 2 12 1 3

    b := [2, 1, 5]

    lu :=

    1 0 01 3 02 1 1

    , 1 1 10 1 0

    0 0 1

    L :=

    1 0 01 3 02 1 1

    U :=

    1 1 10 1 0

    0 0 1

    s := 0

    aux 1 := 2s := 0

    s := 2aux 2 := 1

    s := 0

    s := 4

    s := 3

    aux 3 := 2

    s := 0

    rez 3 := 2

    s := 0

    s := 0

    rez 2 := 1

    s := 0

    s := 1s := 3rez 1 := 1

  • 2.3. SISTEME TRIDIAGONALE 61

  • 62 CAPITOLUL 2. REZOLVAREA SISTEMELOR LINIARE

    2.3.2 Problema rezolvata

    Exercitiul 2.3.1. Sa se rezolve sistemul tridiagonal:

    x +2y = 32x y +z = 2

    3y +2z t = 42z +t = 1.

    RezolvareMatricea sistemului este

    A =

    1 2 0 02 1 1 00 3 2 10 0 2 1

    Descompunem aceasta matrice astfel:

    1 2 0 02 1 1 00 3 2 10 0 2 1

    =

    1 0 0 02 2 0 00 3 3 00 0 2 4

    1 2 0 00 1 3 00 0 1 40 0 0 1

    Din definitia produsului a doua matrice, obtinem:

    b1 = 1 1 1 = 1c2 = 1 2 2 = 2b2 = a2 2 + 2 2 = 5c3 = 23 3 = 1

    5

    b3 = a33 + 3 3 = 135

    c4 = 3 4 4 = 513

    b4 = a4 4 + 4 4 = 313.

    B. Rezolvarea sistemelor triunghiularePentru a rezolva sistemul initial, avem de rezolvat doua sisteme triun-

    ghiulare:

    1 0 0 02 5 0 00 3 13

    50

    0 0 2 313

    y1y2y3y4

    =

    324

    1

    ,

  • 2.3. SISTEME TRIDIAGONALE 63

    a carui solutie este

    y1y2y3y4

    =

    345813

    1

    ,

    si respectiv:

    1 2 0 00 1 1

    50

    0 0 1 513

    0 0 0 1

    xyzt

    =

    345813

    1

    ,

    a cariu solutie este

    xyzt

    =

    1111

    .

    2.3.3 Probleme propuse

    Exercitiul 2.3.2. Sa se rezolve sistemele tridiagonale:

    a)

    x+ y = 32x y + z = 13y z = 5

    b)

    2x+ y = 0x y + 2z = 12y z + t = 5z + 2t = 5.

    2.3.4 Implementare

    A. AlgoritmDate de intrare: un sistem de ecuatii tridiagonalDate de iesire: solutia sistemului

    Algoritmul consta n:

    1. generarea matricei A a sistemului (matrice tridiagonala) si a vectoruluicoloana b

    n = numarul de linii ale matricei A2. descompunerea LU aplicata matricei tridiagonale A

    L = (ij)i,j=1,n, U = (ij)i,j=1,n

  • 64 CAPITOLUL 2. REZOLVAREA SISTEMELOR LINIARE

    pentru i = 1, n

    ii = 1

    pentru i = 2, n

    i,i1 = ai,i1

    11 = a11

    pentru i = 1, n 1i,i+1 =

    ai,i+1ii

    i+1,i+1 = ai+1,i+1 ai+1,i i,i+13. rezolvarea sistemelor triunghiulare

    y1 =b111

    pentru i = 2, n

    yi =

    (bi i,i1yi1

    ) 1ii

    xn = yn

    pentru i = n 1, 1

    xi =

    (yi i,i+1xi+1

    )

    B. Programe MAPLE si rezultate

    Observatia 2.3.1. Spre deosebire de metodele anterioare, unde ordineanecunoscutelor n sistem nu era esentiala, n cazul sistemelor tridiagonale,daca se schimba ordinea necunoscutelor, atunci matricea sistemului nu vamai fi tridiagonala. De aceea, este necesara construirea unei proceduri,nedeterminate, care sa returneze necunoscutele sistemului n ordinea n careapar ele n ecuatii.

    restart: with(linalg):

    tridiagonal:=proc(A::matrix)

    local i,j,n,a1,l,u;

    n:=rowdim(A);

    for i from 1 to n do

    for j from 1 to i-2 do

    if A[i,j]0 then

  • 2.3. SISTEME TRIDIAGONALE 65

    ERROR(matricea nu este tridiagonala!);

    fi;

    od;

    for j from i+2 to n do

    if A[i,j]0 then

    ERROR(matricea nu este tridiagonala!);

    fi;

    od;

    od;

    a1:=A;

    if a1[1,1]=0 then

    ERROR(factorizarea LU nu este aplicabila!);

    fi;

    for i from n by -1 to 2 do

    if det(a1)0 then a1:=delrows(delcols(a1,i..i),i..i);

    else ERROR(factorizarea LU nu este aplicabila!);

    fi;

    od;

    l:=matrix(n,n,0); u:=matrix(n,n,0);

    for i from 1 to n do

    u[i,i]:=1;

    od;

    for i from 2 to n do

    l[i,i-1]:=A[i,i-1];

    od;

    l[1,1]:=A[1,1];

    for i from 1 to n-1 do

    u[i,i+1]:=A[i,i+1]/l[i,i];

    l[i+1,i+1]:=A[i+1,i+1]-A[i+1,i]*u[i,i+1];

    od;

    RETURN(evalm(l), evalm(u));

    end:

    # procedura care returneaza necunoscutele sistemului

    # in ordinea in care apar in ecuatii

    nedeterminate:=proc(l::set(equation))

    local n,i,j,ops,opst;

    n:=nops(l);

    for i from 1 to n do

    ops[i]:=[seq(op(op(l[i])[1])[j] /

  • 66 CAPITOLUL 2. REZOLVAREA SISTEMELOR LINIARE

    coeff(op(l[i])[1],op(indets(op(op(l[i] )[1])[j]))),

    j=1..nops(op(l[i])[1]))];

    od;

    opst:=ops[1];

    for i from 1 to n do

    for j from 1 to nops(ops[i]) do

    if not(ops[i][j] in opst) then

    opst:=[op(opst),ops[i][j]]

    fi;

    od;

    od;

    RETURN(opst);

    end:

    tridiagonalsist:=proc(l::set(equation))

    local eqm, opst, A, b, n, lu, L, U, i, s, j, aux, rez;

    n:=nops(l);

    opst:=nedeterminate(l);

    eqm:=genmatrix(l, opst, flag);

    A:=delcols(eqm,n+1..n+1);

    b:=col(eqm,n+1);

    lu:=tridiagonal(A);

    L:=lu[1];

    U:=lu[2];

    aux[1]:=b[1]/L[1,1];

    for i from 2 to n do

    aux[i]:=1/L[i,i]*(b[i]-L[i,i-1]*aux[i-1])

    od;

    rez[n]:=aux[n];

    for i from n-1 by -1 to 1 do

    rez[i]:=aux[i]-U[i,i+1]*rez[i+1];

    od;

    RETURN(seq(opst[i]=rez[i], i=1..n));

    end:

    debug(tridiagonalsist):

    tridiagonalsist({x+2*y=3,2*x-y+z=2, 3*y+2*z-t=4, -2*z+t=-1});

    {--> enter tridiagonalsist, args = {x+2*y = 3, 2*x-y+z = 2, 3*y+2*z-t= 4, -2*z+t = -1}

  • 2.3. SISTEME TRIDIAGONALE 67

    n := 4

    opst := [x, y, z, t]

    eqm :=

    1 2 0 0 32 1 1 0 20 3 2 1 40 0 2 1 1

    A :=

    1 2 0 02 1 1 00 3 2 10 0 2 1

    b := [3, 2, 4, 1]

    lu :=

    1 0 0 02 5 0 00 3

    13

    50

    0 0 2 313

    ,

    1 2 0 0

    0 115

    0

    0 0 1513

    0 0 0 1

    L :=

    1 0 0 02 5 0 00 3

    13

    50

    0 0 2 313

    U :=

    1 2 0 0

    0 115

    0

    0 0 1513

    0 0 0 1

    aux 1 := 3

    aux 2 :=4

    5

    aux 3 :=8

    13aux 4 := 1

    rez 4 := 1

    rez 3 := 1

    rez 2 := 1

    rez 1 := 1

  • 68 CAPITOLUL 2. REZOLVAREA SISTEMELOR LINIARE

    0 ,

    1 22 5 = 1 > 0 ,

    1 2 12 5 21 2 3

    = 2 > 0 ,

  • 2.4. FACTORIZAREA CHOLESKY 69

    matricea A este pozitiv definita. Aplicand factorizarea Cholesky avem:

    A =

    1 2 12 5 2

    1 2 3

    =

    11 0 021 22 0

    31 32 33

    11 21 310 22 32

    0 0 33

    .

    Folosind definitia produsului a doua matrice, obtinem:

    a11 = 211 11 = 1

    a12 = 11 21 21 = 2a13 = 11 31 31 = 1a22 =

    221 +

    222 22 = 1

    a23 = 21 31 + 22 32 32 = 0a33 =

    231 +

    232 +

    233 33 =

    2.

    Se observa ca pentru gasirea elementelor ij, i = 1, n, j = i, n, este suficientsa calculam dezvoltarile corespunzatoare elementelor aij , i = 1, n, j = i, n.

    B. Rezolvarea sistemelor triunghiulare

    Pentru a determina solutia sistemului initial, avem de rezolvat doua sis-teme triunghiulare:

    1 0 02 1 0

    1 02

    y1y2

    y3

    =

    511

    7

    ,

    cu solutia y1y2

    y3

    =

    51

    2

    ,

    si 1 2 10 1 0

    0 02

    xy

    z

    =

    51

    2

    ,

    cu solutia xy

    z

    =

    21

    1

    .

  • 70 CAPITOLUL 2. REZOLVAREA SISTEMELOR LINIARE

    2.4.3 Probleme propuse

    Exercitiul 2.4.2. Sa se gaseasca solutia sistemului urmator, folosind facto-rizarea Cholesky:

    3x+ y + 3z = 11x+ y + 2z = 63x+ y + 4z = 12.

    2.4.4 Implementare

    A. AlgoritmDate de intrare: un sistem de ecuatiiDate de iesire: solutia sistemului

    Algoritm

    1. generarea matricei A a sistemului (simetrica si pozitiv definita) si avectorului b

    n = numarul de ecuatii ale sistemului (numarul de linii ale matriceiA)

    2. factorizarea Cholesky

    11 = a11

    pentru i = 2, n

    pentru j = 1, i 1

    ij =1

    jj

    (aij

    j1k=1

    ikjk

    )

    ii =

    aii i1k=1

    2ik

    3. rezolvarea sistemelor triunghiulare

    y1 =b111

    pentru i = 2, n

    yi =

    (bi

    i1k=1

    ik yk) 1ii

    xn =ynnn

  • 2.4. FACTORIZAREA CHOLESKY 71

    pentru i = n 1, 1

    xi =

    (yi

    nk=i+1

    ki xk) 1ii

    B. Programe MAPLE si rezultate

    restart: with(linalg):

    desccholesky:=proc(A::matrix)

    local a1,n,l,i,s,j,k;

    n:=rowdim(A);

    a1:=A;

    for i from 1 to n do

    for j from 1 to n do

    if a1[i,j]a1[j,i] then

    ERROR(matricea nu este simetrica!)

    fi;

    od;

    od;

    if a1[1,1]0 then a1:=delrows(delcols(a1,i..i),i..i);

    else ERROR(factorizarea Cholesky nu este aplicabila!);

    fi;

    od;

    l:=matrix(n,n,0);

    l[1,1]:=A[1,1];

    for i from 2 to n do

    for j from 1 to i-1 do

    l[i,j]:=1/l[j,j]*(A[i,j]-sum(l[i,k]*l[j,k],k=1..j-1));

    od;

    l[i,i]:=sqrt(A[i,i]-sum(l[i,k]^2,k=1..i-1));od;

    RETURN(evalm(l));

    end:

    # returnarea necunoscutelor sistemului in ordinea in care

    # apar in ecuatii

  • 72 CAPITOLUL 2. REZOLVAREA SISTEMELOR LINIARE

    nedeterminate:=proc(l::list(equation))

    local n,i,j,ops,opst;

    n:=nops(l);

    for i from 1 to n do

    ops[i]:=[seq(op(op(l[i])[1])[j]/

    coeff(op(l[i])[1],op(indets(op(op(l[i])[1])[j]))),

    j=1..nops(op(l[i])[1]))];

    od;

    opst:=ops[1];

    for i from 1 to n do

    for j from 1 to nops(ops[i]) do

    if not(ops[i][j] in opst) then

    opst:=[op(opst),ops[i][j]]

    fi;

    od;

    od;

    RETURN(opst);

    end:

    choleskysist:=proc(l::list(equation))

    local lst, eqm, A, b, n, lu, L, U,i,s,j,aux, rez, rfin;

    lst:=nedeterminate(l);

    eqm:=genmatrix(l, lst, flag);

    n:=nops(lst);

    A:=delcols(eqm,n+1..n+1);

    b:=col(eqm,n+1);

    lu:=desccholesky(A);

    L:=lu;

    U:=transpose(lu);

    for i from 1 to n do

    s:=0; for j from 1 to i-1 do s:=s+L[i,j]*aux[j] od;

    aux[i]:=1/L[i,i]*(b[i]-s)

    od;

    for i from n by -1 to 1 do

    s:=0; for j from i+1 to n do s:=s+U[i,j]*rez[j] od;

    rez[i]:=1/U[i,i]*(aux[i]-s)

    od;

    RETURN(seq(lst[i]=rez[i], i=1..n));

    end:

  • 2.4. FACTORIZAREA CHOLESKY 73

    debug(choleskysist):

    choleskysist([x+2*y+z=5, 2*x+5*y+2*z=11, x+2*y+3*z=7]);

    {--> enter choleskysist, args = [x+2*y+z = 5, 2*x+5*y+2*z = 11,x+2*y+3*z = 7]

    lst := [x, y, z]

    eqm :=

    1 2 1 52 5 2 11

    1 2 3 7

    n := 3

    A :=

    1 2 12 5 2

    1 2 3

    b := [5, 11, 7]

    lu :=

    1 0 02 1 0

    1 0

    2

    L := lu

    U :=

    1 2 10 1 0

    0 0

    2

    s := 0

    aux 1 := 5

    s := 0

    s := 10

    aux 2 := 1

    s := 0

    s := 5

    s := 5

    aux 3 :=

    2

    s := 0

    rez 3 := 1

    s := 0

    s := 0

    rez 2 := 1

    s := 0

    s := 2

    s := 3

    rez 1 := 2

  • 74 CAPITOLUL 2. REZOLVAREA SISTEMELOR LINIARE

    v2 = e1. (2.17)

    Pentru evitarea ambiguitatilor vom considera vectorul v dat de:

    v = a1 + sign(a11) a1 e1. (2.18)

    Propozitia 2.5.2. Oricare ar fi matricea simetrica A, matricea P definitaprin:

    P = I 2 v vT

    v2 (2.19)

    este simetrica si are proprietatea ca elementele 2, 3, . . . , n de pe prima coloanaa matricei PA sunt nule, unde vectorul v este dat de relatia (2.18).

    Definitia 2.5.1. Se numeste matrice Householder de ordin n1 asociatamatricei A si se noteaza cu Pn1 o matrice de ordin n 1 de forma:

    Pn1 = In1 2 v vT

    v2 (2.20)

  • 2.5. FACTORIZAREA HOUSEHOLDER 75

    unde: v = a1n1 + sign(a21) a1n1 e1 este vectorul format cu componentelevectorului a1 care este prima coloana a matricei A, e1 = ( 1, 0, . . . , 0

    n1

    )T si

    In1 este matricea unitate de ordin n 1.Propozitia 2.5.3. Matricea Hauseholder Pn1 asociata unei matrice sime-trice A este simetrica si are proprietatea ca matricea U1 definita prin:

    U1 =

    1 0 . . . 00

    Pn10

    (2.21)

    este simetrica si verifica relatia:

    A(1) = U1 AU1 =

    a11 1 0 . . . 0

    1 a(1)22 a

    (1)23 . . . a

    (1)2n

    0 a(1)32 a

    (1)33 . . . a

    (1)3n

    0 a

    (1)n2 a

    (1)n3 . . . a

    (1)nn

    (2.22)

    Se considera vectorul coloana a2n2 cu ultimele n2 elemente ale coloaneimatrice A(1). Cu acest vector se construieste o matrice Householder de or-dinul n 2, Pn2. Matricea U2 definita prin:

    U2 =

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

    ... Pn20 0

    (2.23)

    are proprietatea:

    A(2) = U2 A(1) U2 =

    a11 1 0 0 . . . 0

    1 a(1)22 2 0 . . . 0

    0 2 a(2)33 a

    (2)34 . . . a

    (2)3n

    0 0 a

    (2)n3 a

    (2)n4 . . . a

    (2)nn

    (2.24)

    Matricea Pn2 a condus la obtinerea unei noi linii si coloane a matricei tridi-agonale la care vrem sa reducem matricea A.

    Continuand astfel prin n1 transformari, obtinem egalitatea: UAU = Tn care T este matrice tridiagonala.

  • 76 CAPITOLUL 2. REZOLVAREA SISTEMELOR LINIARE

    2.5.2 Problema rezolvata

    Exercitiul 2.5.1. Sa se rezolve urmatorul sistem folosind factorizarea House-holder:

    2x+ 2y + z = 2

    2x y + z = 5x+ y + 2z = 0.

    RezolvareSistemul se poate scrie sub forma Ax = b, unde:

    A :=

    2 2 12 1 1

    1 1 2

    si b :=

    25

    0

    Se observa ca matricea A este simetrica, deci factorizarea Householder esteaplicabila.

    Generarea matricei UCalculam elementele vectorului

    v = a1 + sign(a111) a1 e1unde a1 = (2, 1)T , a1 = 5 si e1 = (1, 0)T . De aici rezulta

    v = (0, 2 +5, 1)T

    si v = 10 + 45. Elementele matricei U sunt date de:

    Uj,k = I3 2vj vkv .

    Dupa efectuarea calculelor, obtinem:

    U =

    1 0 0

    0 2 (2 +5)

    5 + 25

    2 +5

    5 + 25

    0 2 +5

    5 + 25

    2 (2 +5)

    5 + 25

    si

    T = UAU =

    2 5 (2 +5)

    5 + 25

    0

    5 (2 +5)

    5 + 25

    2

    5

    95

    095

    3

    5

    , Ub =

    221

  • 2.5. FACTORIZAREA HOUSEHOLDER 77

    Solutia sistemului tridiagonal

    Ty = Ub

    este

    y =

    1

    20 10521

    5 6515

    ,

    iar solutia sistemului initial este

    x = Uy =

    2 +5

    325

    3

    23

    ,

    adica x =2 +

    5

    3, y =

    253

    z = 23.

    2.5.3 Probleme propuse

    Exercitiul 2.5.2. Sa se gaseasca solutiile urmatoarelor sisteme folosind fac-torizarea Householder:

    a)

    x+ 2y + z = 52x y + 3z = 6x+ 3y 2z = 3

    b)

    x+ y z + t = 3x+ 2y + 3z + t = 9x+ 3y + 2t = 7x+ y + 2z = 5.

    2.5.4 Implementare

    A. AlgoritmDate de intrare: un sistem de ecuatiiDate de iesire: solutia sistemului, obtinuta folosind factorizarea House-

    holder

    Algoritmul consta din urmatoarele etape:

  • 78 CAPITOLUL 2. REZOLVAREA SISTEMELOR LINIARE

    1. generarea matricei tridiagonale pentru i = 1 . . . n 2pentru l = 1 . . . n

    pentru m = 1 . . . n

    daca m = l atunci uml = 1

    daca m 6= l atunci uml = 0//Generam vectorul v

    norm a =

    nj=i+1

    a2ij

    ei+1 = 1

    pentru j = i+ 2 . . . n

    ej = 0

    pentru j = i+ 1 . . . n

    vj = aij + sign(ai,i+1) norm a ej

    norm v=

    nj=i+1

    v2j

    //Generam matricea U

    j = i+ 1 . . . n

    k = i+ 1 . . . n

    ujk = ujk 2 vj vk/ norm v// D=AU

    m = 1 . . . n

    l = 1 . . . n

    dml =

    nk=1

    amk ukl

    //A=UD=UAU

    m = 1 . . . n

    l = 1 . . . n

    aml =

    nk=1

    umk dkl

  • 2.5. FACTORIZAREA HOUSEHOLDER 79

    2. rezolvarea sistemului tridiagonal (vezi paragraful 2.3.4):

    Ty = Ub

    3. gasirea solutiei sistemului initial:

    x = Uy

    B. Programe MAPLE si rezultate

    Observatia 2.5.1. Deoarece o conditie neceesara pentru aplicarea factorizariiHouseholder este ca sistemul sa fie simetric, folosim procedura nedeterminate(prezentata n paragraful 2.4.4). De asemenea, pentru descompunerea matri-cei tridiagonale, am folosit procedura tridiagonal prezentata n paragraful2.3.4.

    householder:=proc(l::list(equation))

    local opst, eqm, A, b, n, i, ii, j, U, UD, e, v, norma, normv,

    k, lu, LL, UU, aux, rez;

    n:=nops(l);

    opst:=nedeterminate(l);

    eqm:=genmatrix(l, opst, flag);

    A:=delcols(eqm,n+1..n+1);

    b:=col(eqm,n+1);

    for i from 1 to n do

    for j from 1 to n do

    if A[i,j]A[j,i] then

    ERROR(sistemul nu este simetric!);

    fi;

    od;

    od;

    rez:=vector(n,0);

    UD:=matrix(n,n,0);

    for ii from 1 to n do

    UD[ii,ii]:=1;

    od;

    for i from 1 to n-2 do

    U:=matrix(n,n,0);

    for ii from 1 to n do

    U[ii,ii]:=1;

    od;

  • 80 CAPITOLUL 2. REZOLVAREA SISTEMELOR LINIARE

    e:=vector(n,0);

    v:=vector(n,0);

    norma:= simplify(sqrt(sum(A[i,k1]^2, k1=i+1..n)));

    e[i+1]:=1;

    for j from i+1 to n do

    v[j]:=A[i,j]+signum(A[i,i+1])*norma*e[j]

    od;

    normv:=simplify(sum(v[k1]^2, k1=i+1..n));

    for j from i+1 to n do

    for k from i+1 to n do

    U[j,k]:=simplify(U[j,k]-2*v[j]*v[k]/normv);

    od;

    od;

    evalm(U);

    A:=simplify( multiply( U, multiply(A,U) ) );

    b:=simplify( multiply(U,b) );

    UD:=multiply(UD,U);

    od;

    lu:=tridiagonal(A);

    LL:=simplify(lu[1]);

    UU:=simplify(lu[2]);

    aux[1]:=simplify(b[1]/LL[1,1]);

    for i from 2 to n do

    aux[i]:=simplify(1/LL[i,i]*(b[i]-LL[i,i-1]*aux[i-1]))

    od;

    rez[n]:=aux[n];

    for i from n-1 by -1 to 1 do

    rez[i]:=simplify(aux[i]-UU[i,i+1]*rez[i+1]);

    od;

    rez:=simplify(multiply(UD,rez));

    RETURN(seq(opst[i]=rez[i], i=1..n));

    end:

    debug(householder):

    householder([x-y+2*z+2*t=0,-x+y+3*z=-1,2*x+3*y-z+2*t=2,2*x+2*z=1]);

    {--> enter householder, args = [x-y+2*z+2*t = 0, -x+y+3*z = -1,2*x+3*y-z+2*t = 2, 2*x+2*z = 1]

    n := 4

    opst := [x, y, z, t]

  • 2.5. FACTORIZAREA HOUSEHOLDER 81

    eqm :=

    1 1 2 2 01 1 3 0 1

    2 3 1 2 22 0 2 0 1

    A :=

    1 1 2 21 1 3 0

    2 3 1 22 0 2 0

    b := [0, 1, 2, 1]rez := [0, 0, 0, 0]

    UD :=

    0 0 0 00 0 0 00 0 0 00 0 0 0

    UD1, 1 := 1

    UD2, 2 := 1

    UD3, 3 := 1

    UD4, 4 := 1

    U :=

    0 0 0 00 0 0 00 0 0 00 0 0 0

    U1, 1 := 1

    U2, 2 := 1

    U3, 3 := 1

    U4, 4 := 1

    e := [0, 0, 0, 0]

    v := [0, 0, 0, 0]

    norma := 3

    e2 := 1

    v2 := 4v3 := 2

    v4 := 2

    normv := 24

    U2, 2 :=13

    U2, 3 :=2

    3

  • 82 CAPITOLUL 2. REZOLVAREA SISTEMELOR LINIARE

    U2, 4 :=2

    3

    U3, 2 :=2

    3

    U3, 3 :=2

    3

    U3, 4 :=13

    U4, 2 :=2

    3

    U4, 3 :=13

    U4, 4 :=2

    3

    1 0 0 0

    013

    2

    3

    2

    3

    02

    3

    2

    3

    13

    02

    3

    13

    2

    3

    A :=

    1 3 0 0

    31

    9

    4

    9

    19

    9

    04

    9

    16

    9

    22

    9

    019

    9

    22

    9

    179

    b :=

    [0,

    7

    3,

    1

    3,23

    ]

    UD :=

    1 0 0 0

    013

    2

    3

    2

    3

    02

    3

    2

    3

    13

    02

    3

    13

    2

    3

    U :=

    0 0 0 00 0 0 00 0 0 00 0 0 0

  • 2.5. FACTORIZAREA HOUSEHOLDER 83

    U1, 1 := 1

    U2, 2 := 1

    U3, 3 := 1

    U4, 4 := 1

    e := [0, 0, 0, 0]

    v := [0, 0, 0, 0]

    norma :=

    377

    9e3 := 1

    v3 :=4

    9+

    377

    9

    v4 :=19

    9

    normv :=754

    81+

    8

    377

    81

    U3, 3 := 4 (4 +

    377)

    377 + 4

    377

    U3, 4 := 19 (4 +

    377)

    377 + 4

    377

    U4, 3 := 19 (4 +

    377)

    377 + 4

    377

    U4, 4 :=4 (4 +

    377)

    377 + 4

    377

    1 0 0 00 1 0 0

    0 0 4 (4 +

    377)

    377 + 4

    37719 (4 +

    377)

    377 + 4

    377

    0 0 19 (4 +

    377)

    377 + 4

    377

    4 (4 +

    377)

    377 + 4

    377

    A :=

    1 3 0 0

    31

    9 377 (4 +

    377)

    9 (377 + 4

    377)0

    0 377 (4 +

    377)

    9 (377 + 4

    377) 2537 (4 +

    377)2

    9 (377 + 4

    377)21122 (4 +

    377)2

    (377 + 4

    377)2

    0 01122 (4 +

    377)2

    (377 + 4

    377)2240

    377

  • 84 CAPITOLUL 2. REZOLVAREA SISTEMELOR LINIARE

    b :=

    [0,

    7

    3,

    34 (4 +

    377)

    3 (377 + 4

    377), 9 (4 +

    377)

    377 + 4

    377

    ]

    UD :=

    1 0 0 0

    013

    46 (4 +

    377)

    3 (377 + 4

    377)10 (4 +

    377)

    377 + 4

    377

    02

    3

    11 (4 +

    377)

    3 (377 + 4

    377)14 (4 +

    377)

    377 + 4

    377

    02

    3 34 (4 +

    377)

    3 (377 + 4

    377)

    9 (4 +

    377)

    377 + 4

    377

    lu :=

    1 0 0 0

    3809

    0 0

    0 377 (4 +

    377)

    9 (377 + 4

    377) 6759 (4 +

    377)2

    80 (377 + 4

    377)20

    0 01122 (4 +

    377)2

    (377 + 4

    377)2240

    377+

    11190080 (4 +

    377)2

    751 (377 + 4

    377)2

    ,

    1 3 0 0

    0 1377 (4 +

    377)

    80 (377 + 4

    377)0

    0 0 1299202253

    0 0 0 1

    LL :=

    1 0 0 0

    3809

    0 0

    0 377 (4 +

    377)

    9 (377 + 4

    377) 6759 (4 +

    377)2

    80 (377 + 4

    377)20

    0 01122 (4 +

    377)2

    (377 + 4

    377)211370320 (393 + 8

    377)

    751 (377 + 4

    377)2

    UU :=

    1 3 0 0

    0 1377 (4 +

    377)

    80 (377 + 4

    377)0

    0 0 1299202253

    0 0 0 1

    aux 1 := 0

    aux 2 :=2180

  • 2.6. METODA JACOBI 85

    aux 3 := 3 (377 + 4

    377)

    751 (4 +

    377)

    aux 4 := 9 (377 + 4

    377) (4 +

    377)

    30160 (393 + 8

    377)

    rez 4 := 9 (377 + 4

    377) (4 +

    377)

    30160 (393 + 8

    377)

    rez 3 := 3 (377 + 4

    377)

    377 (4 +

    377)

    rez 2 :=940

    rez 1 :=27

    40

    rez :=

    [27

    40,

    1

    5,740

    ,116

    ] solve({x-y+2*z+2*t=0,-x+y+3*z=-1,2*x+3*y-z+2*t=2,2*x+2*z=1},> {x,y,z,t});

    {x = 2740

    , y =1

    5, z =

    740

    , t =116}

    2.6 Metoda Jacobi

    2.6.1 Breviar teoretic

    Metoda Jacobi este o metoda iterativa de rezolvare a sistemelor liniare deforma

    Ax = b. (2.25)

    Matricea A se descompune n suma L+D + U , unde

    L =

    0 0 0 . . . 0a21 0 0 . . . 0a31 a32 0 . . . 0 an1 an2 an3 . . . 0

    (2.26)

  • 86 CAPITOLUL 2. REZOLVAREA SISTEMELOR LINIARE

    D =

    a11 0 0 . . . 00 a22 0 . . . 00 0 a33 . . . 0 0 0 0 . . . ann

    (2.27)

    U =

    0 a12 a13 . . . a1n0 0 a23 . . . a2n 0 0 0 . . . an1,n0 0 0 . . . 0

    . (2.28)

    Se defineste traiectoria Jacobi a vectorului x(0) ca fiind vectorul

    x(k+1) = D1[b (L+ U)x(k)] k = 0, 1, 2, . . . (2.29)

    Folosind teorema de convergenta se studiaza daca traiectoria Jacobi con-verge la solutia x() a sistemului (2.25).

    Traiectoria Jacobi converge la solutia x() a sistemului (2.25), daca sinumai daca raza spectrala a matricei

    M = D1(L+ U) (2.30)

    este strict subunitara, adica

    max{|| | det(M In) = 0} < 1. (2.31)

    In caz de convergenta, componentele x(k+1)1 , ..., x

    (k+1)n ale vectorului

    x(k+1), situat pe traiectoria Jacobi a vectorului x(0), sunt date de relatiile:

    x(k+1)i =

    (bi

    nj=1j 6=i

    aij x(k)j) 1aii

    , i = 1, 2, . . . , n; k = 0, 1, . . . (2.32)

    2.6.2 Problema rezolvata

    Exercitiul 2.6.1. Calculati primii trei termeni ai traiectoriei Jacobi asociatevectorului (0, 0, 0) pentru sistemul:

    5x 2y + 3z = 13x+ 9y + z = 22x y 7z = 3.

  • 2.6. METODA JACOBI 87

    RezolvareSistemul se mai poate scrie sub forma Ax = b, unde:

    A =

    5 2 33 9 1

    2 1 7

    , b =

    12

    3

    .

    Matricea A se descompune n suma L+D + U cu

    L =

    0 0 03 0 0

    2 1 0

    , D =

    5 0 00 9 0

    0 0 7

    , U =

    0 2 30 0 1

    0 0 0

    .

    Verificarea conditiei de convergenta a algoritmului presupune calculul valo-rilor proprii ale matricei

    M = D1(L+ U) =

    02

    5

    35

    1

    30

    19

    2

    7

    17

    0

    Calculand maximul n modul al valorilor proprii ale matricei M , obtinem

    (M) = 0.2673998083 < 1,

    si deci algoritmul converge.Aplicam formulele (2.32), plecand de la x(0) = 0, y(0) = 0, z(0) = 0, obtinemsuccesiv:

    x(1) = 0.2000000000y(1) = 0.2222222222

    z(1) = 0.4285714286x(2) = 0.1460317460

    y(2) = 0.2031746032

    z(2) = 0.5174603174x(3) = 0.1917460316

    y(3) = 0.3283950617

    z(3) = 0.4158730159.Pentru comparatie, am rezolvat acest sistem folosind procedura solve

    furnizata de Maple, iar rezultatele sunt:

    x = 0.1861198738, y = 0.3312302839, z = 0.4227129338.

  • 88 CAPITOLUL 2. REZOLVAREA SISTEMELOR LINIARE

    2.6.3 Probleme propuse

    Exercitiul 2.6.2. Sa se verifice daca se poate aplica metoda iterativa a luiJacobi, si n caz afirmativ sa se gaseasca primele 3 elemente ale sirului desolutii partiale. Comparati solutia obtinuta cu solutia exacta:

    a)

    {x+ 3y = 22x+ y = 6

    b)

    x+ 2y + z = 13x y + 5z = 14x+ y z = 2

    2.6.4 Implementare

    A. Algoritm

    Observatia 2.6.1. Deoarece metoda lui Jacobi este o metoda iterativa, tre-buie specificata o conditie de oprire a algoritmului. Algoritmul converge dacasirul (x(k)) este convergent. Convergenta acestui sir poate fi descrisa n modteoretic n diverse moduri. In practica, se foloseste o varianta a criteriuluilui Cauchy, si anume: sirul (x(k)) este convergent, daca

    x(k+1) x(k) <

    unde este o constanta data.In cazul nostru, vom considera ca solutiile sistemului au fost obtinute cu

    eroarea , adica

    x(k+1)i [x()i , x()i + ].

    De aici rezulta o conditie de oprire a algoritmului: ni=1

    (x(k+1)i x(k)i ) <

    n. (2.33)

    Date de intrare: un sistem de ecuatii (o multime de ecuatii), un punctinitial, x(0), o eroare .

    Date de iesire: solutia aproximativa a sistemului, obtinuta n urmaaplicarii traiectoriei Jacobi vectorului x(0) pana cand este ndeplinita conditia(2.33).

    Algoritmul consta n urmatoarele etape:

  • 2.6. METODA JACOBI 89

    1. generarea matricei A a sistemului si a vectorului b

    n - numarul de necunoscute (numarul de linii ale matricei A)2. generarea matricelor L, D, U si verificarea convergentei metodei:

    (D1(L+ U)) < 1

    3. construirea traiectoriei Jacobi

    repeta

    x(k+1)i =

    (bi

    nj=1j 6=i

    aij x(k)j) 1aii

    , i = 1, n

    pana cand

    ni=1

    (x(k+1)i x(k)i ) <

    n

    B. Programe MAPLE si rezultate

    jacobi:=proc(eq::set(equation), init::vector, eps::float)

    local var, n, AA, A, b, l, d, u, i, j, m, lst, xo, test, k, x;

    var:=[op(indets(eq))];

    n:=nops(var);

    if vectdim(init)n then

    ERROR(numarul de necunoscute nu este egal cu

    dimensiunea vectorului initial)

    fi;

    AA:=genmatrix(eq, var, flag);

    A:=delcols(AA,n+1..n+1);

    b:=col(AA,n+1);

    l:=matrix(n,n,0):

    u:=matrix(n,n,0):

    d:=matrix(n,n,0):

    for i from 1 to n do

    for j from 1 to i-1 do

    l[i,j]:=A[i,j];

    od;

    d[i,i]:=A[i,i];

    for j from i+1 to n do

    u[i,j]:=A[i,j];

  • 90 CAPITOLUL 2. REZOLVAREA SISTEMELOR LINIARE

    od;

    od;

    # conditia de convergenta

    m:=multiply(inverse(d),matadd(l,u,-1,-1));

    lst:=[eigenvals(m)];

    if evalf(max(seq(abs(lst[k]),k=1..nops(lst))))>=1 then

    ERROR(Algoritmul nu converge);

    fi;

    # algoritmul propriu-zis

    for i from 1 to n do

    xo[i]:=init[i]

    od;

    test:=1;

    while test>=evalf(eps*sqrt(n)) do

    for i from 1 to n do

    x[i]:=evalf(

    1/A[i,i]*( b[i]-sum(A[i,k]*xo[k],k=1..n)+A[i,i]*xo[i] )

    );

    od;

    test:=evalf(sqrt( sum( (x[k]-xo[k])^2, k=1..n ) ));

    for i from 1 to n do

    xo[i]:=x[i];

    od;

    od;

    RETURN(seq(var[i]=x[i],i=1..n));

    end:

    debug(jacobi):

    jacobi({3*x+y=5, x+2*y=5}, vector(2,[0,0]),0.01);

    {--> enter jacobi, args = {3*x+y = 5, x+2*y = 5}, array(1 ..2,[(1)=0,(2)=0]), .1e-1

    var := [x, y]

    n := 2

    AA :=

    [3 1 51 2 5

    ]

    A :=

    [3 11 2

    ]b := [5, 5]

  • 2.6. METODA JACOBI 91

    l :=

    [0 00 0

    ]

    u :=

    [0 00 0

    ]

    d :=

    [0 00 0

    ]d1, 1 := 3

    u1, 2 := 1

    l2, 1 := 1

    d2, 2 := 2

    m :=

    0

    13

    12

    0

    lst := [

    6

    6,

    6

    6]

    xo1 := 0

    xo2 := 0

    test := 1

    x1 := 1.666666667

    x2 := 2.500000000

    test := 3.004626063

    xo1 := 1.666666667

    xo2 := 2.500000000

    x1 := 0.8333333333

    x2 := 1.666666666

    test := 1.178511303

    xo1 := 0.8333333333

    xo2 := 1.666666666

    x1 := 1.111111111

    x2 := 2.083333334

    test := 0.5007710115

    xo1 := 1.111111111

    xo2 := 2.083333334

    x1 := 0.9722222220

    x2 := 1.944444444

    test := 0.1964185512

    xo1 := 0.9722222220

  • 92 CAPITOLUL 2. REZOLVAREA SISTEMELOR LINIARE

    xo2 := 1.944444444

    x1 := 1.018518519

    x2 := 2.013888889

    test := 0.08346183593

    xo1 := 1.018518519

    xo2 := 2.013888889

    x1 := 0.9953703703

    x2 := 1.990740740

    test := 0.03273642604

    xo1 := 0.9953703703

    xo2 := 1.990740740

    x1 := 1.003086420

    x2 := 2.002314815

    test := 0.01391030679

    xo1 := 1.003086420

    xo2 := 2.002314815

    jacobi({3*x+y=5, x+2*y=5}, vector(2,[0,0]),0.01);x = 1.003086420, y = 2.002314815

    > jacobi({3*x+y=5, x+2*y=5}, vector(2,[0,0]),0.001);x = 0.9998713993, y = 1.999742798

    > jacobi({3*x+y=5, x+2*y=5}, vector(2,[0,0]),0.00001);x = 1.000002381, y = 2.000001786

    De asemenea, se poate compara sirul solutiilor partiale cu solutia exactaobtinuta rezolvand sistemul Ax = b cu ajutorul procedurii linsolve. Pentrusistemul considerat mai sus, a carui solutie exacta este x = 1, y = 2, obtinemurmatoarele grafice:

  • 2.6. METODA JACOBI 93

    comparatie cu solutia exacta

    0

    0.2

    0.4

    0.6

    0.8

    1

    1.2

    1.4

    1.6

    1 2 3 4 5 6 7 8

    comparatie cu solutia exacta

    0

    0.5

    1

    1.5

    2

    2.5

    1 2 3 4 5 6 7 8

  • 94 CAPITOLUL 2. REZOLVAREA SISTEMELOR LINIARE

    2.7 Metoda Gauss-Seidel

    2.7.1 Breviar teoretic

    Metoda Gauss-Seidel este o metoda de rezolvare numerica a si