Jacobi estrapolato

Implementare in una M–function il metodo di Jacobi estrapolato e applicarlo a un problema per cui si ha convergenza.
File: esercizio.m

%Esercizio 8 file "esercizi_iterativi"
%8. Implementare in una M–function il metodo di Iacobi estrapolato e applicarlo a un problema per cui si ha convergenza.
%Bonafè Antonio 058223
%Mesin Alberto 066711

clear all;
close all;
clc;
%
%definizione A, b, x0
%
s=[5 -3 -2 10 -7 -3 -2 -6 9 -1 -9 10 -1 -3 -4 7 -1 4];
i=[1 1 1 2 2 2 3 3 3 3 4 4 4 5 5 5 6 6];
j=[1 4 5 2 5 6 1 2 3 6 3 4 5 1 2 5 5 6];
A=sparse(i,j,s);

b=[1;6;2;5;3;4];

x0=zeros(6,1);
%
%ricerca di w -> minimo raggio spettrale
%
D1=inv(diag(diag(A)));
I=speye(size(A));

ws=warning('off'); %escludo visualizzazione Warning
i=1;
ro=zeros(201,2);
for wg=0:0.01:2
    Gw=(I-wg*D1*A);
    ro(i,2)=wg;
    ro(i,1)=max(abs(eigs(Gw)));
    i=i+1;
end;
warning(ws); %reimposto il precedente stato di visualizzazione Warning

ro_index=find(ro(:,1)<1); %Quando converge
%plot(ro(:,2),ro(:,1));
subplot(2,1,1); %%%%%%%%%%%%%%%%%%%%%%%%%
plot(ro(ro_index,2),ro(ro_index,1));
title('raggio spettrale in funzione di w');
xlabel('w');
ylabel('raggio spettrale');
%ro_index=find(ro(:,1)==min(ro(:,1)));

%definizione w, tol, maxit
%
%w=1.17;
w=ro((find(ro(:,1)==min(ro(:,1)))),2);
tol=1e-10;
maxit=500;
%
%Soluzione sistema con JACOBI
%
%[xk,k]=myjacobi(A,b,x0,maxit,tol);
%display('JACOBI');
%display(xk);
%display(k);
[xk,k,r]=myjacobiestrapolato(A,b,x0,1,maxit,tol);
display('JACOBI');
fprintf('\t%.4f \n',full(xk));
fprintf('Iterazioni=%g \n',k);

fprintf('\n');
%
%Soluzione sistema con JACOBI ESTRAPOLATO
%
[xke,ke,re]=myjacobiestrapolato(A,b,x0,w,maxit,tol);
fprintf('JACOBI ESTRAPOLATO (parametro=%g)\n',w);
fprintf('\t%.4f \n',full(xke));
fprintf('Iterazioni=%g \n',ke);
fprintf('\n');
fprintf('Risparmiate %g iterazioni\n',k-ke);
%
%Costruzione grafico
%
mi=10+max(max(find(r>0)),max(find(re>0)));
%x=(1:maxit);
x=(1:mi);
%figure; %nuova finestra grafica
%hold on;
subplot(2,1,2);
plot(x,r(1:mi),'r',x,re(1:mi),'g');
axis([0,mi,-0.25,max(max(r),max(re))]);
legend('Jacobi (normale)','JACOBI ESTRAPOLATO');
title('Residuo');
xlabel('iterazioni');
ylabel('||r^k||');

File: myjacobiestrapolato.m

function [x,k,g]=myjacobiestrapolato(A,b,x,w,maxit,tol);
% Metodo di JACOBI ESTRAPOLATO
% [x,k]=myjacobiestrapolato(A,b,x,w,maxit,tol)
% A -> Matrice dei coefficenti
% b -> Vettore termini noti
% x -> Prima approssimazione delle soluzioni
% w -> Parametro di rilassamento
% maxit -> Massimo numero di iterazioni
% tol -> Tolleranza
% ------------------------------------------------
% x -> Vettore approssimazione soluzioni
% k -> Numero iterazioni effettuate
% g -> Vettore per il grafico (norma del residuo)

%
%Preparazione di G(w) e c(w)
%
D=diag(diag(A));
D1=inv(D);
I=speye(size(A));

Gw=(I-w*D1*A);
cw=w*D1*b;

%
%Preparo g
%
%g=zeros(maxit,1);
g=-0.5*ones(maxit,1);

%
% Iterazioni
%
for k=1:maxit
    %calcolo iterazione successiva
    xtemp=x;
    x=Gw*xtemp+cw;
    %grafico
    g(k)=norm(b-A*x);
    %controllo sulla convergenza
    %g(k)=norm(xtemp-x,inf);
    %if g(k)<tol*norm(x,inf)
    if norm(xtemp-x,inf)<tol*norm(x,inf)
        break
    end;
end;
% controllo sul supereamento delle massime iterazioni
if k==maxit
    error('Convergenza non raggiunta')
end;