Metoda Jacobi

Metoda Jacobi

Breviar teoretic

Metoda Jacobi este o metoda iterativa de rezolvare a sistemelor liniare de forma

Ax = b. (2.25)

Matricea A se descompune in suma L + D + U , unde

(2.26)

(2.27)

(2.28)

Se defineste traiectoria Jacobi a vectorului x⁽⁰⁾ ca fiind vectorul

x^(k+1) = D⁻¹[b − (L + U )x^(k)] k = 0, 1, 2, (2.29)

Folosind teorema de convergenta se studiaza daca traiectoria Jacobi converge la solutia x^(*) a sistemului (2.25).

Traiectoria Jacobi converge la solutia x^(*) a sistemului (2.25), daca si numai daca raza spectrala ρ a matricei

M = −D⁻¹(L + U) (2.30)

este strict subunitara, adica

max < 1. (2.31)

In caz de convergenta, componentele ale vectorului , situat pe traiectoria Jacobi a vectorului x⁽⁰⁾ , sunt date de relatiile:

(2.32)

Problema rezolvata

Calculati primii trei termeni ai traiectoriei Jacobi asociate vectorului (0, 0, 0) pentru sistemul:

Rezolvare

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

Matricea A se descompune in suma L + D + U cu

Verificarea conditiei de convergenta a algoritmului presupune calculul valorilor proprii ale matricei

M = −D⁻¹(L + U) =

Calculand maximul in modul al valorilor proprii ale matricei M , obtinem

(M ) = 0.2673998083 < 1,

si deci algoritmul converge.

Aplicam formulele (2.32), plecand de la x⁽⁰⁾ = 0, y⁽⁰⁾ = 0, z⁽⁰⁾ = 0, obtinem succesiv:

x y⁽¹⁾ = 0.2222222222, z⁽¹⁾ = −0.4285714286

x y⁽²⁾ = 0.2031746032 z⁽²⁾ = −0.5174603174

x = 0.1917460316 y⁽³⁾ = 0.3283950617 z⁽³⁾ = −0.4158730159

Pentru comparatie, am rezolvat acest sistem folosind procedura solve furnizata de Maple, iar rezultatele sunt:

x = 0.1861198738, y = 0.3312302839, z = −0.4227129338.

Implementare

A. Algoritm

Observatia. Deoarece metoda lui Jacobi este o metoda iterativa, trebuie specificata o conditie de oprire a algoritmului. Algoritmul converge daca sirul (x^(k)) este convergent. Convergenta acestui sir poate fi descrisa in mod teoretic in diverse moduri. In practica, se foloseste o varianta a criteriului lui Cauchy, si anume: sirul (x^(k)) este convergent, daca

unde ε este o constanta data.

In cazul nostru, vom considera ca solutiile sistemului au fost obtinute cu eroarea ε, adica

De aici rezulta o conditie de oprire a algoritmului:

(2.33)

Date de intrare: un sistem de ecuatii (o multime de ecuatii), un punct initial, x⁽⁰⁾, o eroare ε.

Date de iesire: solutia aproximativa a sistemului, obtinuta in urma aplicarii traiectoriei Jacobi vectorului x⁽⁰⁾ pana cand este indeplinita conditia (2.33).

Algoritmul consta in urmatoarele etape:

1. generarea matricei A a sistemului si a vectorului b

. n - numarul de necunoscute (numarul de linii ale matricei A)

2. generarea matricelor L, D, U si verificarea convergentei metodei:

(−D⁻¹(L + U )) < 1

3. construirea traiectoriei Jacobi

repeta

pana cand

B. Programe MAPLE si rezultate

jacobi := proc (eq::set (equation), init::vector, eps::float)

local var, n, AA, A, b, l, d, u, i, j, m, lst, xo, test, k, x;

var := [op (indets (eq))];

n := nops (var);

if vectdim (init) <> n then

ERROR ('numarul de necunoscute nu este egal cu dimensiunea vectorului initial')

fi;

AA := genmatrix (eq, var, flag);

A := delcols (AA, n+1 .. n+1);

b := col (AA, n+1);

l := matrix (n, n, 0):

u := matrix (n, n, 0):

d := matrix (n, n, 0):

for i from 1 to n do

for j from 1 to i-1 do

l[i, j] := A[i, j];

od;

d[i, i] := A[i, i];

for j from i+1 to n do

u[i, j] := A[i, j];

od;

od;

# conditia de convergenta

m := multiply (inverse (d), matadd (l, u, -1, -1));

lst := [eigenvals (m)];

if evalf (max (seq (abs (lst[k]), k=1 .. nops (lst))))>=1 then

ERROR('Algoritmul nu converge');

fi;

# algoritmul propriu-zis

for i from 1 to n do

xo[i] := init[i]

od;

test := 1;

while test >= evalf (eps_*sqrt(n)) do

for i from 1 to n do

x[i] := evalf (1/A[i, i]*( b[i]-sum(A[i, k]*xo[k], k=1 .. n)+A[i,i]*xo[i]));

od;

test := evalf (sqrt (sum ((x[k]-xo[k])^2, k=1 .. n )));

for i from 1 to n do

xo[i] := x[i];

od;

od;

RETURN (seq (var[i] = x[i], i=1 .. n));

end:

debug(jacobi):

jacobi(, vector(2, [0,0]), 0.01);

, array(1 .. 2, [(1)=0, (2)=0]), .1e-1

var := [x, y]

n

, b

d u_{1, 2} := 1 l_{2, 1} := 1 d_{2, 2} := 2

xo₁ := 0 xo₂ := 0 test := 1

x x₂ := 2.500000000 test := 3.004626063

xo₁ := 1.666666667 xo₂ := 2.500000000

x := 0.8333333333 x₂ := 1.666666666 test := 1.178511303

xo₁ := 0.8333333333 xo₂ := 1.666666666

x x₂ := 2.083333334 test := 0.5007710115

xo₁ := 1.111111111 xo₂ := 2.083333334

x := 0.9722222220 x₂ := 1.944444444 test := 0.1964185512

xo₁ := 0.9722222220 xo₂ := 1.944444444

x x₂ := 2.013888889 test := 0.08346183593

xo₁ := 1.018518519 xo₂ := 2.013888889

x := 0.9953703703 x₂ := 1.990740740 test := 0.03273642604

xo₁ := 0.9953703703 xo₂ := 1.990740740

x := 1.003086420 x₂ := 2.002314815 test := 0.01391030679

xo₁ := 1.003086420 xo₂ := 2.002314815

<-- exit jacobi (now at top level) = x = 1.003086420, y = 2.002314815}

x = 1.003086420, y = 2.002314815

Se pot compara rezultatele obtinute daca eroarea ε se modifica:

> jacobi (, vector (2, [0,0]), 0.01);

x = 1.003086420, y = 2.002314815

> jacobi (, vector (2, [0,0]), 0.001);

x = 0.9998713993, y = 1.999742798

> jacobi (, vector (2, [0,0]), 0.00001);

x = 1.000002381, y = 2.000001786

De asemenea, se poate compara sirul solutiilor partiale cu solutia exacta obtinuta rezolvand sistemul Ax = b cu ajutorul procedurii linsolve. Pentru sistemul considerat mai sus, a carui solutie exacta este x = 1, y = 2, obtinem urmatoarele grafice: