[Matlab] Metodo iterativo di Jacobi

[Dalfi1]

Salve ragazzi, ho creato il seguente algoritmo per implementare il metodo iterativo di Jacobi per la risoluzione di sistemi lineari

function [ sol ] = Jacobiter( A,b,x0,tol,imax )
%UNTITLED Summary of this function goes here
%   Detailed explanation goes here

[row,column]=size(A);
n=length(b);
iter=0;
r=b-A*x';
err=norm(r);
while err>tol & iter<imax
    for i=1:n
        s=sum(A(i,1:i-1).*x0(1:i-1));
        s=s+sum(A(i,i+1:end).*x0(i+1:end));
        x(i)=(b(i)-s)/A(i,i);
        iter=iter+1;
    end
    r=b-A*x';
    err=norm(r);
    x0=x;
end
sol=x';

    
end

Ho però un problema alquanto strano: se faccio agire l'algoritmo su matrici e vettori random funziona (ho verificato con A\b), se invece lo uso per sistemi presi da libri di testo funziona solo se come vettore iniziale do il vettore soluzione. Riuscireste a farmi capire quale parte del codice è sbagliata gentilmente?

PS: L'algoritmo è ancora in stato di bozza: mancano ancora i vari casi da escludere (matrici non singolari, con elementi diagonali nulli, ecc.) ma per i test ho preso solo i casi possibili, quindi dovrebbe funzionare.

PPS: Anche se il controllo dell'errore lo faccio su norm(x-x0) non cambia niente...stessi errori

[Lory314]

Non ho esaminato a fondo il codice, ma l'altro giorno scrivere una cosa del genere in MatLab mi aveva dato dei problemi.

while err>tol & iter<imax

Prova con

while(and(err>tol, iter<imax))

Non ti assciuro che sia questo l'errore!

[Dalfi1]

Grazie mille per la risposta ma non credo sia quello l'errore. Se uso l'algoritmo su A=rand(n)+q*eye(n), b=rand(n,1) va tutto bene, altrimenti no
Idem per l'algoritmo di Gauss Seidel appena creato e presumo per il metodo di rilassamento che sto per andare a scrivere.
Ecco Gauss Seidel:

function [ sol,iter ] = GausSeidel( A,b,x0,tol,imax )
%UNTITLED2 Summary of this function goes here
%   Detailed explanation goes here

[row,column]=size(A);
n=length(b);
if column~=n
    error('Errore: la matrice A e il vettore b non hanno dimensioni compatibili');
end
D=diag(A);
if any(D)==0
    error('La matrice ha elemento/i diagonale/i nullo/i')
end
iter=0;
err=inf;
x0=x0';
x=x0;
while (err>tol & iter<imax)
    for i=1:n
        s1=A(i,1:i-1)*x(1:i-1);
        s2=A(i,i+1:n)*x0(i+1:n);
        x(i) = (b(i)-s1-s2)/A(i,i);
        iter=iter+1;
    end
    err=norm(x-x0);
    x0=x;
end
sol=x;


end

[Raptorista1]

[xdom="Raptorista"]Sposto nella sezione più appropriata! [/xdom]

[Cerrone1]

Ho un problema analogo ... nel senso che il programma in matlab è strano mi restituisc e ogni giorno un numero di iterazioni diverso e a volte nn mi da la soluzione esatta pure se la matrice è strett a diag dominante...molto strano... posto il codice magari qlk è csì gentile da aiutarmi?...

function [ x, itnum, res, ratio, err, rap, flag ] = jacobi( A, b, w, x0, tol, maxit, xsol )

% Metodo iterativo di Jacobi rilassato per sistemi lineari
%
%Input:
% A: Matrice del sistema
% b: Termine noto (vettore colonna)
% w: Parametro di rilassamento ( se è uguale a 1 stiamo lavorando con
% Jacobi non rilassato)
% x0: Vettone iniziale
% tol: Tolleranza sul test d'arresto
% maxit: Numero massimo di iterazioni
% xsol: Soluzione del sistema

%Output:
% x: Vettore soluzione del metodo
% itnum: Numero di iterazioni effettuate
% res: Residuo al passo k
% ratio: Rapporto tra i residui ai passi consecutivi
% err: Errore al passo k
% rap: Rapporto ta gli errori ai passi consecutivi
% flag: Se flag = 0 converge nel numero massimo di iterazioni
% fissato altrimenti flag=1



if ( det(A) == 0 )
disp('** La matrice è singolare');
else
[n,n] = size (A);
I = eye(n);
L = -tril(A);
U = -triu(A);
D = diag(diag(A));
L = L+D;
U = U+D;
x = x0;
flag = 0;
itnum = 1;
T = (1-w)*I+w*D^(-1)*(L+U);
C = w*D^(-1)*b;

res(1) = norm(A*x-b);
err(1) = norm(xsol-x0);

while ( itnum <= maxit && res(itnum) >= tol )
itnum = itnum+1;
x = T*x+C;
res(itnum )= norm (A*x-b);
ratio(itnum-1)= res(itnum)/res(itnum - 1);
err(itnum)= norm(xsol-x);
rap(itnum-1)= (err(itnum))/(err(itnum-1));
end
if (res(itnum)>= tol)
flag=1;
else
disp('** Jacobi non converge nel numero di iterazioni fissato');
end
end


% Grafico rappresentate il rapporto tra gli errori ai passi
% consecutivi-numero di iterazioni in scala logaritmica sull'asse y
figure(1)
y=[2:itnum];
semilogy(y,rap,'r-o');
% Nomino gli assi
xlabel('num iterazioni');
ylabel('rapporto errori');

end

[Raptorista1]

Stai infrangendo in un colpo solo una marea di punti del regolamento, ed inoltre ti ammonisco con il mio cartellino preferito