Factorizarea LU

Factorizarea LU

Breviar teoretic

Fie sistemul compatibil determinat

Ax = b.

unde

Factorizarea LU presupune descompunerea matricei A intr-un produs de matrice L U , unde

, ( 2.5 )

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

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

Astfel, rezolvarea sistemului (2.4) se reduce la rezolvarea sistemelor triunghiulare

Ly = b (2.6)

cu solutia

(2.7)

si

Ux = y (2.8)

cu solutia

(2.9)

Problema rezolvata

Exercitiul 1 Sa se determine solutia sistemului urmator, folosind factorizarea LU:

Sistemul se scrie in forma matriceala:

Ax = b

unde

,

Deoarece

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 Crout

Presupunem ca

si ne propunem sa determinam coeficientii l_ij, u_jk. Pentru aceasta, folosim definitia inmultirii matricelor. Astfel, avem:

a = λ₁₁ 1 λ₁₁ = 1

a = λ_{11 12 12} = 1

a = λ_{11 13 13} = −1

a = λ₂₁ 1 λ₂₁ = 2

a = λ_{21 12} + λ₂₂ 1 λ₂₂ = −3

a ₂₃ = −1

a λ₃₁ = 1

a λ₃₂ = 2

a λ₃₃ = 1

sau

B. Rezolvarea sistemelor triunghiulare

Pentru rezolvarea sistemului initial, avem de rezolvat doua sisteme triunghiulare:

a carui solutie este

si respectiv:

a carui solutie este

REZOLVARE FOLOSIND FACTORIZAREA DOOLITTLE

A. Factorizarea Doolittle

Presupunem ca

si ne propunem sa determinam coeficientii l_ij, _jk, la fel ca si in exemplul precedent. Astfel avem:

a = 1 _{11 11} = 1

a = 1 _{12 12} = 1

a = 1 _{13 13} = −1

a = λ_{21 11} λ₂₁ = 2

a = λ_{21 12} + 1 _{22 22} = −3

a ₂₃ = 3

a λ₃₁ = 1

a ₂₂ λ₃₂ =

a ₃₃ = 1

sau

B. Rezolvarea sistemelor triunghiulare

Pentru rezolvarea sistemului initial, avem de rezolvat doua sisteme triunghiulare:

,

a carui solutie este

si respectiv:

a carui solutie este

Implementare

A. Algoritm

Date de intrare: un sistem de ecuatii.

Date 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 =

_ii = 1

pentru i =

pentru j =

pentru j =

b) factorizarea Doolittle

pentru i =

λ_ii = 1

pentru i =

pentru j =

pentru j =

3. Rezolvarea celor doua sisteme triunghiulare

pentru i =

pentru i =

B. Programe MAPLE si rezultate

Deoarece cele doua variante ale descompunerii LU difera doar prin modul de factorizare a matricei sistemului, am implementat separat cele doua variante de factorizare: LUcrout si LUdoolittle, dupa care le-am folosit ca optiuni in 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!');

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

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 (, Doolittle);

Doolittle

lst := n := 3

s := 0 aux ₁ := 2 s := 0 s := −2 aux ₂ := 3

s := 0 s := 4 s := 3 aux ₃ := 2 s := 0 rez ₃ := 2

s := 0 s := 0 rez ₂ := 1

s := 0 s := 1 s := 3 rez ₁ := 1

<-- exit LUsist (now at top level) = z = 1, x = 1, y = 2}

z = 1, x = 1, y = 2

> LUsist(, Crout);

, Crout

lst := n := 3

s := 0 aux ₁ := −2 s := 0 s := −2 aux ₂ := 1

s := 0 s := 4 s := 3 aux ₃ := 2

s := 0 rez ₃ := 2 s := 0 s := 0 rez ₂ := 1

s := 0 s := −1 s := −3 rez ₁ := 1

<-- exit LUsist (now at top level) = z = 1, x = 1, y = 2} z = 1, x = 1, y = 2

Programul pentru factorizarea LU prin limbajul C++

# include<iostream.h>

# include<math.h>

# include<conio.h>

int n, m, i, j, k, p, t;

double a[10][11], y[10], x[10], max, s;

void main (void)

cout<<"Introduceti termenii liberi;"<<endl;

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

k = 1;

t = 0;

do

max = fabs(a[k][k]);

p = k;

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

if (max<fabs(a[i][k]))

if (max = = 0)

t = 1;

else

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

k++;}}

while ((k<=n)&&(t = = 0));

if ((t = =1) | | (a[n][n] = = 0.0))

cout<<"Sistemul nu este compatibil determinat!"<<endl'

else

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

for (i = n; i > 0; i--)

cout.precision(0);

cout<<"Solutia sistemului:"<<end;

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

cout<<"x("<<i<<")="<<setw(12)<<x[i]<<endl;

}

getch ( );

}