Problema metodo gauss seidel

[wecalculus]

Ragà data questa traccia, punto 2, non riesco ad arrivare al risultato, perchè non so richiamare bene la funzione...

http://www.dm.uniba.it/~pugliese/didattica/CalcNumTa1112/prova_lab_2011-06-13.pdf

questo è il file m che uso:


function [xn,k]=gs(A,b,x0,toll,nmax)

% Metodo di Gauss-Seidel
%
% A: matrice del sistema
% b: termine noto
% x0: vettore iniziale
% toll: tolleranza sul residuo normalizzato
% nmax: massimo numero di iterazioni
%
% xn: soluzione ottenuta
% k: numero di iterazioni effettuate

n = length(b);
xn = zeros( n, 1 );
k = 0;

if (( size(A,1)~=n) || (size(A,2)~=n) || (length(x0) ~= n) )
error('dimensioni incompatibili')
end

if (prod(diag(A)) == 0)
error('errore: elementi diagonali nulli')
end

L = tril(A);
xv = x0;
r = b - A * x0;
err = norm(r) / norm(b);

while ( err > toll && k < nmax )
k = k + 1;
z = fwsub(L,r);
xn = xv + z;
r = b - A*xn;
err = norm(r) / norm(b);
xv = xn;
end


e la richiamo cosi:

prima metto A=[2.4 -0.8 -0.7; 0.5 1.5 0.7; -0.1 0.8 2.1]
poi b= [0.9;2.7;2.8]

poi [xn,k]=gs(A,b,x0,10^-5,100)


non so che mettere a x0, pensavo x0=zeros(n,1) dove n è presa dalla b, quindi x0=zeros(3,1)...ma matlab da errore perchè sicuramente sbaglio :S chi mi sa aiutare?
domani avrei l'esame...

[walter891]

io proverei mettendo ones(3,1) perchè il vettore nullo forse non va bene per questo metodo

[wecalculus]

ho risolto con questo nuovo files da me creato:



function [x1,flag] = provagauss (A,b,tol,x0,maxiter)

numiter=0;
err=1+tol; flag=0;
while err>tol && numiter if numiter>0
x0=x1;
end
L=tril(A);
U=L-A;
x1=L\(U*x0+b);
errOld=err;
err= norm(((abs(x1-x0))/(abs(x1))),1);
errNew=err;
numiter=numiter+1;
l=sprintf('%d\t %1.15e\t %1.3e\t %0.3g',...
numiter,x1,err,log(errNew)/log(errOld));
disp(l)
end
if err>tol
disp('Attenzione: raggiunto il massimo numero di iterate')
else
flag=numiter;
end


ho imposto A,b come quelli scritta prima, come tol, quella scritta sulla traccia, come maxiter=100, e con x0 proprio x0=zeros(3,1)
pare funzionare ora...il risultato è pure corretto... secondo voi è corretto??