CALCUL MATRICIAL(II). Rezolvarea sistemelor de ecuatii neliniare (Metoda iterativa Jacobi, Metoda...

download CALCUL MATRICIAL(II). Rezolvarea sistemelor de ecuatii neliniare (Metoda iterativa Jacobi, Metoda Gauss-Seidel, Metoda lui Newton, Metoda lui Newton-Raphson)

of 5

Transcript of CALCUL MATRICIAL(II). Rezolvarea sistemelor de ecuatii neliniare (Metoda iterativa Jacobi, Metoda...

  • 8/7/2019 CALCUL MATRICIAL(II). Rezolvarea sistemelor de ecuatii neliniare (Metoda iterativa Jacobi, Metoda Gauss-Seidel, Me

    1/5

    Carmen-Sanda Georgescu, Tudor Petrovici, Radu PopaMetode numerice. Fisa de laborator nr. 9:

    CALCUL MATRICIAL(II). Rezolvarea sistemelor de ecuatii neliniare(Metoda iterativa Jacobi, Metoda Gauss-Seidel, Metoda lui Newton,

    Metoda lui Newton-Raphson)

    6. REZOLVAREA SISTEMELOR DE ECUAII NELINIARE

    Fie un sistem de ecuaii neliniare de forma

    ( ) 0xxxf n21i =,,, K , n21i ,,, K= , (9.1)

    care poate fi scris sub forma vectoriala: ( ) 0=XF , F:Rn R, F(X)=[f1(X) f2(X) . . . fn(X) ]si X=[x1 x2 . . . xn] (9.2)

    Sistemul (9.2) este echivalent cu sistemul: ( )XGX = , G:Rn R, G(X)=[g1(X) g2(X) . . . gn(X) ]si X=[x1 x2 . . . xn], unde gi(X) sunt alese astfel incat daca xi= gi(X), pentru i=1,2, . . . ,n atuncifi(X)=0 (9.3)

    6.1. METODA ITERATIV JACOBI PENTRU REZOLVAREA SISTEMELOR NELINIAREAlgoritm: Pornind un proces iterativ cu vectorul de start X(0)=[x(0)1 x

    (0)2 . . . x

    (0)n], vectorul soluie se

    calculeaz iterativ cu relaia:( ) ( )( )k1k XGX =+ ,

    ale carei componente sunt:

    ( ) ( ) ( ) ( )( ) ( )( )kiknk2k1i1ki gxxxgx X==+ ,,, K , n21i ,,, K= . (9.5)In fond, se considera sistemul neliniar sub forma:

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

    ( ) ( ) ( )

    ( )

    =

    =

    =

    +

    +

    +

    0xxxf

    0xxxf

    0xxxf

    1k

    n

    k

    2

    k

    1n

    kn

    1k2

    k12

    kn

    k2

    1k11

    ,,,

    ,,,

    ,,,

    K

    KKK

    K

    K

    ,

    din care se rezolva (explicit sau numeric) fiecare ecuaie pentruif( )1kix+

    in funcie de ,( )kjx n1j ,= si

    .ij Not: Functia iterativa seamana cu metoda substitutiilor succesive pentru si cu metodaiterativa Jacobi pentru sisteme de ecuaii algebrice liniare. Conditia de convergen este:

    ( ) 0xf =

    ( )

    1x

    gn

    1j j

    i

  • 8/7/2019 CALCUL MATRICIAL(II). Rezolvarea sistemelor de ecuatii neliniare (Metoda iterativa Jacobi, Metoda Gauss-Seidel, Me

    2/5

    Dupa ietratia a 4-a solutiile incep sa se departeze usor de solutia exacta (Vezi Problema 1 la APLICATII)si se stabilizeaza dupa iteratia a 19-a la valorile x= 0.50000, y= 0.86603

    6.2. METODA ITERATIV GAUSS-SEIDEL PENTRU REZOLVAREA SISTEMELORNELINIARE

    Algoritm: In scopul imbunatatirii vitezei de convergenta a procesului iterativ Jacobi, metoda luiGauss si Seidel utilizeaza o schema iterativa obtinuta astfel (conditiile de convergenta raman cele de la

    metoda Jacobi): intr-o ecuaie oarecare a sistemului (9.6) toate valorile solutiei deja obtinute laiteratia curenta, adica ecuaia:

    if

    ( ) ( ) ( ) ( ) ( )( ) 0xxxxxf knk1i1ki1k21k1i =++++ ,,,,,, KK (9.8)se rezolva (explicit sau numeric) pentru

    ( )1kix+

    in funcie de( ) ( ) ( ) ( )k

    nk1i

    1k1i

    1k1 xxxx ,,,,, KK +

    +

    +.

    Exemplul 2:Folosind metoda iterativa Gauss-Seidel cu 3 iteratii, sa se rezolve sistemul de ecuatii neliniarx3-3*x*y2+1=03*x2-y3=0folosind solutia de start x0=1, y0=1Pas1: Se expliciteaza x=g1(y) si y=g2(x) din a doua ecuatie si respectiv din prima ecuatieSe obtine sistemul echivalent x=sqrt(y3/3*y), y=sqrt((x3+1)/3*x)Pas2: Daca g1(y)= sqrt(y

    3/3*y) si g2(x)= sqrt((x3+1)/3*x), atunci conditiile de convergenta (9.7) devin

    sqrt(3)/3 , (sqrt(3)/3)*sqrt(x3/(x3+1)) cu

  • 8/7/2019 CALCUL MATRICIAL(II). Rezolvarea sistemelor de ecuatii neliniare (Metoda iterativa Jacobi, Metoda Gauss-Seidel, Me

    3/5

    Pas2: Se construieste schema iterativa Newton generalizata valabila pentru k=0,1,2, . . .x(k+1)=x(k) [ (x(k))2+(y(k))2-1 ]/[ 2*x(k)]y(k+1)=y(k) [(x(k))3-y(k)]/[ -1]Pas3: Se porneste procesul iterativ cu solutia de start x(0)=0.9, y(0)=0.5Rezulta la primul pas de iteratiex(1)=x(o) [(x(0))2+(y(0))2-1]/[ 2*x(0)] x(1)= 0.86667y(1)=y(0) [(x(0))3-y(0)]/[ -1] y(1)= 0.72900Se continua iteratiile pentru k=1,2, . . . , conform schemei iterative de la pasul 2

    6.4. METODA NEWTON-RAPHSON PENTRU REZOLVAREA SISTEMELOR NELINIARE

    Algoritm: Sistemul (9.4) se poate scrie pentru diverse forme ale lui ( )XG . De exemplu, fie:, unde( ) ( ) ( )XFXAXXG = ( ) ( )XXA ija= este o matrice patrata de ordinul n. In metoda

    Newton-Raphson, matricea A se alege a fi inversa matricii Jacobian J , adica: ( ) (XJXA 1= ) , unde

    ( )( )

    =

    j

    i

    x

    f XXJ , functiile fifiind componentele vectorului F.

    Relaia iterativ de calcul a solutiei sistemului devine:( ) ( ) ( )( ) ( )( )kkkk XFXJXX = + 11 . (9.10)

    Daca introducem vectorul de corectie intre doua iteratii succesive:( ) ( ) ( )kkk

    XXX = ++ 11 , se observa ca (9.10) devine:

    ( )( ) ( ) ( )( )kkk XFXXJ = +1 . (9.12)Dupa rezolvarea sistemului liniar (9.12) si aflarea vectorului de corectie X(k+1), se obtine din (9.11) solutiala iteratia curenta:X(k+1).

    Se poate utiliza si un coeficient de pondere plasat in domeniul 10 < :( ) ( ) ( )11 ++ += kkk XXX . (9.13)

    Iteratiile continua pana la realizarea unui criteriu de convergenta impus, cum ar fi:

    ( )

    ( )

    +

    k

    i

    ki

    i x

    x1

    max , sau( )( ) ki

    if Xmax , cu ales suficient de redus.

    Comentariu: O dificultate majora a metodei o constituie necesitatea determinarii Jacobianul functiei F(X)care contine derivatele partiale fi/xj..Algoritmii numerici care implementeaza algoritmul Newton-Raphson au posibilitatea ca in eventualitatea in care utilizatorul nu foloseste Jacobianul functiei F(X) sainlocuiasca il inlocuiasca cu calculul derivatelor partiale prin diferente finite, J(X,h)=i j(X,h).Exemplul 4:Folosind metoda iterativa Newton-Raphson cu 3 iteratii, sa se rezolve sistemul de ecuatiineliniarx3-3*x*y2+1=03*x2-y3=0folosind solutia de start x(0)=1, y(0)=1Pas1: Se formeaza matricea F(X)=[ f1(x,y) f2(x,y)]

    T =[x3-3*x*y2+1 3*x2-y3]T si matricea jacobianJ(X)=[ f1/x(x,y) f1/y(x,y); f2/x(x,y) f2/y(x,y)]=[3*x2-3*y2 -6*x*y; 6*x -3*y2]Pas2: Se calculeaza F(X(0))=[-1 1]T si J-1(X(0))=[-1/12 1/6;-1/6 0] si se parcurge primul pas al procesuluiiterativ cu solutia de start x(0)=1, y(0)=1[x(1) y(1)]T=[x(0)y(0)]T J-1(X(0))* F(X(0)) [x(1) y(1)]T =[1 1]T-[-1/12 1/6;-1/6 0]* [-1 1]T=[11/12 7/6]Pas3: Se calculeaza J-1(X(1)) si se parcurge pasul 2 al procesului iterativ

    [x(2) y(2)]T=[x(1)y(1)]T J-1(X(1))* F(X(1))Pas4: Se continua cautarea solutiei dupa schema iterativa[x(k+1) y(k+1)]T=[x(k)y(k)]T J-1(X(k))* F(X(k)), pentru k=2,3, . . .Functii Octave pentru rezolvarea sistemelor neliniare si/sau transcendente

    [x, info, msg] = fsolve (fcn, x0)Calculeaza solutia ecuatiei transcendente/sitemului transcendent, pornind de la solutia de start x0.Remarcam ca procesul iterativ de cautare a solutiei este sensibil la valoarea de startIntrari: Daca fcn, contine numai numele functiei f(x) si x0 este valoarea de start a solutiei, fsolve ,rezolva ecuatia /sistemul f(x) = 0. Daca fcn este un tablou de doua elemente, primul element estenumele functiei f iar al doilea element este jacobianul j(x)

    df_ijac(i,j) = ----

    dx_jDe remarcat ca se poate folosi si functia fsolve_options care poate contine un set de parametriipentru fsolve

  • 8/7/2019 CALCUL MATRICIAL(II). Rezolvarea sistemelor de ecuatii neliniare (Metoda iterativa Jacobi, Metoda Gauss-Seidel, Me

    4/5

    APLICAII

    Problema 1. Sa se rezolve sistemul de ecuatii neliniare , folosind functia Octave fsolvex3-3*x*y2+1=03*x2-y3=0folosind solutia de start x(0)=1, y(0)=1Sa se compare rezultatul obtinut cu rezultatele obtinute la Exemplele 1 si 2

    Rezolvare: Se editeaza functia y=f(x)octave#1>function y = f (x)> y(1) = x(1).^3 + 3*x(1).*x(2).^2+1;> y(2) = 3*x(1).^2 - x(2).^3 ;>endfunctionapoi se cheama in mediul Octave functia fsolve cu sintaxaoctave#2> [x, info,msg] = fsolve ("f", [1; 1])Rezulta solutia x si info ( info=1 anunta ca procesul prin care s-a gasit solutia, este convergent)x =[ 0.476822160309098 0.880260820452681]T, info = 1, msg = solution converged within specifiedtoleranceVerificare: octave#3>f(x) # Returneaza ans = -3.77475828372553e-15 -1.33226762955019e-15Problema 2. Folosind metodele iterative Jacobi si Gauss-Seidel, cu 4 pasi de iterare, sa se rezolvesistemul de ecuatii neliniarx3-3*x*y2-2*x+2=03*x2*y-y3-2*y=0folosind solutia de start x0=1, y0=1

    Problema 3. Folosind metoda iterativa Newton-Raphson, cu 3 pasi de iterare, sa se rezolve sistemul deecuatii neliniarx2+y2-1=0x3--y=0folosind solutia de start x0=0.9, y0=0.5Sa se compare solutiile obtinute cu solutia care se obtine folosind functia Octave fsolvex= 0.82603 , y=0.56362Problema 4. Folosind metoda iterativa Newton generalizata , cu 2 pasi de iterare, sa se rezolve sistemulde ecuatii neliniare si transcendentex+3*lg(x)-y2=02*x3x*y-5*x+1=0folosind solutia de start x0=3.4, y0=2.2Sa se compare solutiile obtinute cu solutia care se obtine folosind functia Octave fsolvex= 3.4874 , y= 2.2616Problema 5. Folosind metoda iterativa Gauss-Seidel , cu 3 pasi de iterare, sa se rezolve sistemul deecuatii neliniare si transcendente

    sin(x)-y-1.32=0x-cos(y)-0.5=0folosind solutia de start x0=1.3, y0=-0.4Sa se compare solutiile obtinute cu solutia care se obtine folosind functia Octave fsolve x= 1.44679 , y= -0.32768Problema 6. Sa se rezolve sistemul de ecuatii neliniare si transcendente, folosind functia Octave fsolve-2x^2 + 3xy + 4 sin(y) = 63x^2 - 2xy^2 + 3 cos(x) = -4folosind solutia de start x0=1, y0=2Rezolvare: Se editeaza functia y=f(x)octave#1>function y = f (x)> y(1) = -2*x(1). 2 + 3*x(1).*x(2) + 4*sin(x(2)) - 6;> y(2) = 3*x(1).^2 - 2*x(1).*x(2)^2 + 3*cos(x(1)) + 4;>endfunctionapoi se cheama in mediul Octave functia fsolve cu sintaxaoctave#2> [x, info,msg] = fsolve ("f", [1; 2]) Observatiecand nu se specifica Jacobianul lui F

    Rezulta solutia x si info ( info=1 anunta ca procesul prin care s-a gasit solutia, este convergent)x =[ 0.633445666940771 1.482872760937817]T, info = 1, msg = solution converged within specifiedtoleranceVerificare: octave#3>f(x) # Returneaza ans =2.4869e-14 -8.5265e-14Comentariu: Daca se foloseste Jacobianul functiei F(X) (in sensul definit de metoda Newton-Raphson), vatrebui ca inainte de a se chema functia fsolve in mediul Octave, sa se scrie functia f1 care definesteJacobianul sistemului, apoi se cheama in mediul Octave functia fsolve cu sintaxaoctave#3> [x, info] = fsolve (["f" "f1], [1; 2])Rezulta solutia: x=0.633445666940845; y= 1.482872760937504Verificare: octave#3>f(x) # Returneaza ans =2.4869e-14 -8.5265e-14

  • 8/7/2019 CALCUL MATRICIAL(II). Rezolvarea sistemelor de ecuatii neliniare (Metoda iterativa Jacobi, Metoda Gauss-Seidel, Me

    5/5

    Comentariu: Daca se foloseste Jacobianul functiei F(X) (in sensul definit de metoda Newton-Raphson), vatrebui ca inainte de a se chema functia fsolve in mediul Octave, sa se scrie functia care definesteJacobianul si anume:octave#2>function y=f1(x)>y(1,1)=-4*x(1)+3*x(2);>y(1,2)=3*x(1)+4*cos(x(2));>y(2,1)=6*x(1)-2*x(2).^2-3*sin(x(1));>y(2,2)=-4*x(1).*x(2);

    >endfunctionapoi se cheama in mediul Octave functia fsolve cu sintaxaoctave#3> [x, info] = fsolve (["f" ; "f1], [1; 2])Rezulta solutia: x=0.633445666940845; y= 1.482872760937504