Sisteme tridiagonale

Sisteme tridiagonale

O clasa speciala de sisteme liniare este aceea in care matricea A a sistemului este tridiagonala, adica:

(2.10)

Pentru aceste sisteme se aplica factorizarea LU. Astfel, matricea A se descompune, folosind un caz particular al factorizarii Crout, intr-un produs L U unde:

(2.11)

si

(2.12)

Coeficientii a₂, , a_n sunt cunoscuti din matricea A, iar coeficientii β_i, _j se obtin din definitia inmultirii matricelor:

(2.13)

2 Problema rezolvata

Exercitiul. Sa se rezolve sistemul tridiagonal:

Rezolvare

Matricea sistemului este:

Descompunem aceasta matrice astfel:

Din definitia produsului a doua matrice, obtinem:

.

B. Rezolvarea sistemelor triunghiulare

Pentru a rezolva sistemul initial, avem de rezolvat doua sisteme triunghiulare:

a carui solutie este

si respectiv:

a carui solutie este

Implementare

A. Aloritm

Date de intrare: un sistem de ecuatii tridiagonal

Date de iesire: solutia sistemului

Algoritmul consta in:

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

. n = numarul de linii ale matricei A

2. descompunerea LU aplicata matricei tridiagonale A

pentru i =   

pentru i =

pentru i =

1. rezolvarea sistemelor triunghiulare

pentru i =

pentru i =   

B. Programe MAPLE si rezultate

Observatia 1. Spre deosebire de metodele anterioare, unde ordinea necunoscutelor in sistem nu era esentiala, in cazul sistemelor tridiagonale, daca se schimba ordinea necunoscutelor, atunci matricea sistemului nu va mai fi tridiagonala. De aceea, este necesara construirea unei proceduri, nedeterminate, care sa returneze necunoscutele sistemului in ordinea in care apar ele in 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

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

n := 4

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

aux ₁ := 3, aux ₂ := , aux ₃ :=, aux ₄ := 1,

rez ₄ := 1 , rez ₃ := 1, rez ₂ := 1, rez ₁ := 1

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

x = 1, y = 1, z = 1, t = 1