Factorizarea Householder

Factorizarea Householder

Breviar teoretic

Factorizarea Householder este o metoda de rezolvare numerica a sistemelor de tip Cramer simetrice, si consta in determinarea unei matrice simetrice nesingulare U , astfel incat UAU = T sa fie o matrice tridiagonala. Atunci solutia sistemului Ax = b este data de

x = Uy (2.15)

unde y este solutia sistemului

T y = Ub (2.16)

Factorizarea Householder se bazeaza pe urmatoarele rezultate.

Propozitia 1. Oricare ar fi A o matrice patratica de ordinul n si si-metrica, exista un vector v = (v₁, v₂, . . . , v_n)^T astfel incat vectorul coloana a¹ = Ae₁, e₁ = (1, 0, . . . , 0)^T (a¹ este prima coloana a matricei A) are proprietatea

(2.17)

Pentru evitarea ambiguitatilor vom considera vectorul v dat de:

v = a¹ + sign (a₁₁) ||a¹|| e₁. (2.18)

Propozitia 2. Oricare ar fi matricea simetrica A, matricea P definita prin:

(2.19)

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

Definitia 1. Se numeste matrice Householder de ordin n−1 asociata matricei A si se noteaza cu P_n−1 o matrice de ordin n − 1 de forma:

unde: v = a¹_n−1 + sign(a₂₁) ||a¹_n−1|| e₁ este vectorul format cu componentele vectorului a¹ care este prima coloana a matricei A, e₁= ()^T si I _{n −1} este matricea unitate de ordin n − 1.

Propozitia 3. Matricea Hauseholder P_{n − 1} asociata unei matrice simetrice A este simetrica si are proprietatea ca matricea U₁ definita prin:

(2.21)

este simetrica si verifica relatia:

(2.22)

Se considera vectorul coloana a²_n−2 cu ultimele n − 2 elemente ale coloanei matrice A⁽¹⁾. Cu acest vector se construieste o matrice Householder de ordinul n − 2, P_n−2. Matricea U₂ definita prin:

(2.23)

are proprietatea:

(2.24)

Matricea P_n−2 a condus la obtinerea unei noi linii si coloane a matricei tridiagonale la care vrem sa reducem matricea A.

Continuand astfel prin n −1 transformari, obtinem egalitatea: UAU = T in care T este matrice tridiagonala.

2 Problema rezolvata

Exercitiul. Sa se rezolve urmatorul sistem folosind factorizarea Householder:

Rezolvare

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

si

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

Generarea matricei U

Calculam elementele vectorului

v = a¹ + sign (a¹₁₁) ||a¹|| e₁.

unde a¹ = (2, 1)^T , si e₁ = (1, 0)^T. De aici rezulta v = (0, 2 + )^T si . Elementele matricei U sunt date de:

Dupa efectuarea calculelor, obtinem:

si ,

Solutia sistemului tridiagonal

T y = Ub

este

iar solutia sistemului initial este

, adica

Implementare

A. Algoritm

Date de intrare: un sistem de ecuatii

Date de iesire: solutia sistemului, obtinuta folosind factorizarea Householder

Algoritmul consta din urmatoarele etape:

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

pentru l = 1 . . . n

pentru m = 1 . . . n

daca m = l atunci u_ml = 1

daca m ≠ l atunci u_ml = 0

//Generam vectorul v

, e_i

pentru j = i + 2 . . . n

e_j

pentru j = i + 1 . . . n

v_j = a_ij + sign(a_i,i+1) norm a e_j

//Generam matricea U

j = i +1 n

k = i + 1 n

u_jk = u_jk − 2 v_j v_k / norm v // D=AU

m = 1 n

l = 1 n

//A=UD=UAU

m = 1 n

l = 1 n

2. rezolvarea sistemului tridiagonal: Ty = Ub

3. gasirea solutiei sistemului initial: x = Uy

B. Programe MAPLE si rezultate

Observatie. Deoarece o conditie necesara pentru aplicarea factorizarii Householder este ca sistemul sa fie simetric, folosim procedura nedeterminate. De asemenea, pentru descompunerea matricei tridiagonale, am folosit procedura tridiagonal

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;

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]);

Pentru comparatie, am rezolvat acelasi sistem folosind procedura solve:

> solve (, );