Metodo di Jacobi

Nel metodo di Jacobi $A=L+D+U$ con L strettamente triangolare inferiore, D diagonale ed U strettamente triangolare superiore

$$\left\{ \begin{array}{lll} M & = & D \\ N & = & -L -U \end{... ...neq 0 \leftrightarrow a_{i\:i} \neq 0 \; \forall i = 1\ldots n$$

visto il precedente schema dei metodi iterativi si ha che

$$D x^{k+1} = -(L+U)x^k +b$$

e procedendo riga per riga

$$a_{i\:i}x_i^{k+1} = b_i - \sum _{j=1, j \neq i} ^{j=n} a_{i\:j} x_j ^ k$$

da cui si ricava

$$x_i^{k+1} = \frac{b_i - \sum _{j=1, j \neq i} ^{j=n} a_{i\:j} x_j ^ k } {a_{i\:i}}$$

x=Jacobi(A,b,$x_0$,epsilon) - Data la matrice A non singolare ed il vettore $\underline{b}$ la funzione risolve il sistema lineare $A\underline{x}=\underline{b}$ mediante il metodo iterativo di Jacobi partendo dal vettore di innesco $\underline{x}_0$ con una tolleranza epsilon; la condizione di arresto è basata sul metodo del residuo.

%JACOBI
% x=Jacobi(A,b,x0,epsilon)
% Pre: A non singolare, A(i,i) non nulli per i=1:n.
%      In realta' perche'il problema sia ben condizionato, 
%      A deve 
% essere  a diagonale dominante.
% La funzione restituisce in x la soluzione del sistema Ax=b, 
% risolto con il metodo semiiterativo di Jacobi.
% Riceve in x0 il vettore di innesco e in epsilon la tolleranza 
% desiderata.
% 
% See also GAUSS_SEIDEL, DD
function x=Jacobi(A,b,x0,epsilon)
n=length(b);
x=x0;
r=b-A*x;

while norm(r)>epsilon*norm(b)
   for i=1:n
      som=0;
      for j=1:n
         if i~=j
            som=som+A(i,j)*x(j);
         end
      end
      som=(b(i)-som)/A(i,i);
      x(i)=som;   
   end
   r=b-A*x;
end
return
   
      

A =

     9     6     2
     3   -12     4
     7    15    -3

» b

b =

     1
     2
     4

» Jacobi(A,b,[1;2;3],1E-9)

ans =

   1.00000000008960
  -0.62500000007290
  -2.12500000015545

» inv(A)*b

ans =

   1.00000000000000
  -0.62500000000000
  -2.12500000000000


» A

A =

3.0000000000 1.0000000000 0.5000000000 0.8000000000
1.5000000000 9.0000000000 0.5000000000 0.6000000000
3.2000000000 4.0000000000 12.7000000000 3.2000000000
0.7000000000 3.2000000000 4.1000000000 11.3000000000

» b

b =

   0.50000000000000
   1.20000000000000
   3.00000000000000
   6.20000000000000

»  inv(A)*b

ans =

  -0.01092762955907
   0.09772671951632
   0.08447073456824
   0.49102600234596

» Jacobi (A,b,[0.01;0.09;0.1;1],1E-9)

ans =

  -0.01092762960712
   0.09772671950406
   0.08447073461875
   0.49102600233408