1> % risolvere un sistema lineare mal condzionato

1> A = hilb(13);

2> % A e' simmetrica e definita positiva

2> % partiamo dalla soluzione precisa e ci calcoliamo i termini

2> % noti b

2> x = ones(13,1);

3> b = A * x;

4> cond(A)
ans =  7.4103e+18

5> cond(A,1)
warning: inverse: matrix singular to machine precision, rcond = 1.95902e-19
ans =  5.1046e+18

6> cond(A,'fro')
warning: inverse: matrix singular to machine precision, rcond = 1.95902e-19
ans =  2.2152e+18

7> [L, U, P] = lu(A);

8> function y = avanti(L,b)
> y = b(1);
> for k = 2 : length(b)
> y = [y; b(k) - L(k, 1: k-1) * y];
> end
> end

9> function x = indietro(U,y)
> n = length(y);
> x = y(n)/U(n,n);
> for k = n-1 : -1 : 1
> x = [(y(k)-U(k,k+1 : n)*x)/U(k,k); x];
> end
> end

10> y = avanti(L, P*b);

11> x_lu = indietro(U, y)
x_lu =

   1.00000
   1.00001
   0.99948
   1.00842
   0.92615
   1.39309
  -0.35066
   4.09304
  -3.76710
   5.88476
  -2.19004
   2.20141
   0.80144


12> eps_x_lu = norm(x - x_lu)/norm(x)
eps_x_lu =  2.3171

13> norm(x - x_lu)
ans =  8.3544

14> norm(x)
ans =  3.6056

15> x_div = A\b
warning: matrix singular to machine precision, rcond = 1.24226e-18
warning: matrix singular to machine precision, rcond = 1.34211e-18
x_div =

   1.00000
   1.00000
   1.00002
   0.99972
   1.00194
   0.99217
   1.01866
   0.97528
   1.01194
   1.01292
   0.97618
   1.01439
   0.99679


16> eps_x_div = norm(x - x_div)/norm(x)
eps_x_div =  0.012767

17> help itermeth
warning: /home/achilles/programmi_matlab/Programmi_CS4a/itermeth.m: possible Matlab-style short-circuit operator at line 37, column 17
`itermeth' is a function from the file /home/achilles/programmi_matlab/Programmi_CS4a/itermeth.m

ITERMETH    Un metodo iterativo generale
  [X, ITER] = ITERMETH(A,B,X0,NMAX,TOL,P) cerca di
  risovere iterativamente il sistema di equazioni
  lineari A*X=B su X. La matrice A di N-per-N coef-
  ficienti deve essere non singolare ed il termine
  noto B deve avere lunghezza N. Se P='J' viene usato
  il metodo di Jacobi, se P='G' viene invece selezio-
  nato il metodo di Gauss-Seidel. Altrimenti, P e'
  una matrice N-per-N non singolare che gioca il
  ruolo di precondizionatore in un metodo di
  Richardson a parametro dinamico.
  Il metodo si arresta quando il rapporto
  fra la norma del residuo corrente e quella del
  residuo iniziale e' minore di TOL e ITER e' il
  numero di iterazioni effettuate. NMAX prescrive
  il numero massimo di iterazioni consentite. Se P
  non viene precisata, viene usato il metodo di
  Richardson non precondizionato con parametro di
  rilassamento uguale a 1.

Additional help for built-in functions and operators is
available in the on-line version of the manual.  Use the command
`doc <topic>' to search the manual index.

Help and information about Octave is also available on the WWW
at http://www.octave.org and via the help@octave.org
mailing list.

18> [x_gs, iter] = itermeth(A,b,zeros(13,1),100000,1.e-16,'G')
x_gs =

   0.99995
   1.00106
   0.99468
   1.00705
   1.00271
   0.99604
   0.99447
   0.99737
   1.00175
   1.00488
   1.00502
   1.00133
   0.99364

iter =  100000

19> [x_gs, iter] = itermeth(A,b,zeros(13,1),200000,1.e-16,'G')
x_gs =

   0.99998
   1.00032
   0.99872
   1.00040
   1.00312
   0.99960
   0.99693
   0.99745
   1.00000
   1.00254
   1.00337
   1.00140
   0.99614

iter =  200000

20> eps_x_gs = norm(x - x_gs)/norm(x)
eps_x_gs =  0.0021923

21> quit