                                                    < M A T L A B >
                                        Copyright 1984-2005 The MathWorks, Inc.
                                        Version 7.1.0.183 (R14) Service Pack 3
                                                    August 02, 2005

  To get started, select MATLAB Help or Demos from the Help menu.

>> % controlliamo di avere nel path le funzioni che ci servono, e cioe'
>> % faia, jacobistep, gershcircles, e proviamole
>> n=5
n =
     5
>> gamma=1
gamma =
     1
>> a=faia(n,gamma)
a =

     3    -1     0     0     0
    -1     3    -1     0     0
     0    -1     3    -1     0
     0     0    -1     3    -1
     0     0     0    -1     3
>> % sembra a posto
>> % facciamo un passo di Jacobi
>> b=ones(n,1)
b =
     1
     1
     1
     1
     1
>> y=b
y =
     1
     1
     1
     1
     1
>> y=jacobistep(n,gamma,b,y)
y =
    0.6667
    1.0000
    1.0000
    1.0000
    0.6667
>> % un altro passo
>> y=jacobistep(n,gamma,b,y)
y =
    0.6667
    0.8889
    1.0000
    0.8889
    0.6667
>> % la soluzione di a*x=b viene esattamente
>> x=inv(a)*b
x =
    0.6111
    0.8333
    0.8889
    0.8333
    0.6111
>> % facciamo una ventina di passi di Jacobi
>> for k=1:20
y=jacobistep(n,gamma,b,y);
end
>> y
y =
    0.6111
    0.8333
    0.8889
    0.8333
    0.6111
>> % vediamo la differenza tra la soluzione esatta e quella approssimata
>> norm(x-y)
ans =
   2.6611e-06
>> % sembra soddisfacente .
>> % infine vediamo i cerchi di gershgoring
>> gershcircles(a)
>> % in questo caso si calcolano prima a mano che non con il matlab  !
>> % ora che abbiamo fatto i vari test delle routines, andiamo al punto 4 dell'esercizio
>> n=100
n =
   100
>> % ora sara' opportuno non dimenticarsi il ; alla fine delle righe di comando !!
>> b=ones(n,1);
>> y0=zeros(n,1);
>> e=b;
>> % iniziamo col caso in cui gamma = 1.
>> gamma=1;
>> b(1)=1+gamma;
>> b(n)=1+gamma;
>> % notiamo en passant che il vettore e risolve il sistema:
>> a=faia(n,gamma);
>> norm(a*e-b)
ans =
     0
>> % 20 passi dell'algoritmo di Jacobi salvando i risultati in gg
>> y=y0;
>> gg(1:n,1)=y;
>> for k=2:20
y=jacobistep(n,gamma,b,y);
gg(1:n,k)=y;
end
>> size(gg)
ans =
   100    20
>> % ripetiamo la procedura con gamma=0.01 salvando i risultati in ff
>> gamma=0.01;
>> b(1)=1+gamma;
>> b(n)=1+gamma;
>> % notiamo en passant che il vettore e risolve il sistema:
>> a=faia(n,gamma);
>>  norm(a*e-b)
ans =
     0
>> y=y0;
>> % 20 passi dell'algoritmo di Jacobi salvando i risultati in ff
>> ff(1:n,1)=y;
>> for k=2:20
y=jacobistep(n,gamma,b,y);
ff(1:n,k)=y;
end
>> size(ff)
ans =
   100    20
>> % ora per confrontare la velocita' di convergenza facciamo il grafico degli errori relativi
>> for k=1:20
zz(k)=norm(gg(1:n,k)-e,inf)/norm(e,inf);
ww(k)=norm(ff(1:n,k)-e,inf)/norm(e,inf);
end
>> % in realta' siccome norm(e,inf) vale 1 non occorreva scriverlo !
>> s=1:20;
>> figure(1);plot(s,zz);figure(2);plot(s,ww);
>> % nel secondo caso la convergenza e' piu' veloce.