Factorizarea Cholesky

Factorizarea Cholesky

Breviar teoretic

Un caz particular de sisteme liniare este acela in care matricea A a sistemului este simetrica si pozitiv definita (adica toti determinantii de colt sunt strict pozitivi). Pentru astfel de sisteme putem folosi un caz particular al factorizarii LU: descompunem matricea A a sistemului intr-un produs L L^T , unde

(2.14)

Coeficientii λ_ij se obtin din definitia produsului a doua matrice.

2 Problema rezolvata

Exercitiul. Sa se rezolve sistemul:

Rezolvare

1. Factorizarea Cholesky

Matricea sistemului este

Se observa ca a_ij = a_ji, adica matricea A este simetrica. Deoarece

1 > 0, , ,

matricea A este pozitiv definita. Aplicand factorizarea Cholesky avem:

Folosind definitia produsului a doua matrice, obtinem:

Se observa ca pentru gasirea elementelor λ_ij , i = , j = , este suficient sa calculam dezvoltarile corespunzatoare elementelor a_ij , i = , j = .

B. Rezolvarea sistemelor triunghiulare

Pentru a determina solutia sistemului initial, avem de rezolvat doua sisteme triunghiulare:

cu solutia

si

cu solutia

Implementare

A. Algoritm

Date de intrare: un sistem de ecuatii

Date de iesire: solutia sistemului

Algoritm

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

. n = numarul de ecuatii ale sistemului (numarul de linii ale matricei A)

2. factorizarea Cholesky

λ₁₁ = a₁₁

pentru i =

pentru j =

,

3. rezolvarea sistemelor triunghiulare

pentru i =

,

pentru i =

.

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

ERROR('factorizarea Cholesky 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 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

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:

debug(choleskysist):

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

x = 2, y = 1, z = 1