Jacobi matlab

[andretop00]

Salve, vi sembra corretto questo algoritmo per il metodo di jacobi in matlab:

function x = jac(A,b,x0,tol,it)
x=x0;
xold=x0;
n=size(A);
res=norm(b-A*x);
k=1;

while res>tol && k<it
for i=1:n
s=0;
for j=1:n
if j~=i
s=s+A(i,j)*xold(j);
end 
end
x(i)=(b(i)-s)/A(i,i);
end
res=norm(b-A*x);
k=k+1;
xold=x;
end
end

Dove it è il numero massimo di iterazioni e tol il minimo residuo che posso accettare.

Potreste dirmi se è corretto
Grazie

[feddy]

Ciao, se usi i tag per i codici questo diventa molto piu leggibile sia a te che a chi deve leggere. Hai provato a testarlo su un sistema lineare di cui conosci la soluzione e vedere se converge e come decade l'errore ?

[andretop00]

Ciao, potresti gentilmente spiegarmi come utilizzare i tag?
Ho provato, ma escono numeri elevatissimi, non so il motivo.
Ti sembra corretto il codice?

[feddy]

Devi inserire il tuo codice tra i tag che trovi nell'editor. Per vedere come si fa, premi "Cita" su questo messaggio e guarda il sorgente. Il tuo codice formattato e'

function x = jac(A,b,x0,tol,it)
x=x0;
xold=x0;
n=size(A);
res=norm(b-A*x);
k=1;

while res>tol && k<it
    for i=1:n
        s=0;
        for j=1:n
            if j~=i
                s=s+A(i,j)*xold(j);
            end
        end
        x(i)=(b(i)-s)/A(i,i);
    end
    res=norm(b-A*x);
    k=k+1;
    xold=x;
end
end

L'ho testato su due esempi per i quali si puo' testare e funziona. Invece di scrivere il tutto per componenti, puoi provare a scriverne la versione matriciale, che e' ancora piu' semplice da scrivere, e verificare che trovi li risultati identici. Probabilmente non hai utilizzato una matrice che soddisfa le condizioni necessarie di convergenza.

clear all
close all

%% primo test con matrice presa da wiki
A = [10,-1,2,0;
    -1,11,-1,3;
    2,-1,10,-1;
    0,3,-1,8];
b = A*ones(4,1);
x0 = zeros(4,1);
tol = 1e-6;
it = 2000;

x = jac(A,b,x0,tol,it)





%% altro test con matrice che discretizza l'operatore di laplace


h=0.1;
n = 10;
B = (1/h^2)*toeplitz(sparse(1,[1,2],[-2,1],1,n),sparse(1,[1,2],[-2,1],1,n));
bb = B*ones(n,1);
x0 = zeros(n,1);

[xx,res] = jac(B,bb,x0,tol,it);
iter = [1:length(res)];
semilogy(iter,res,'-o')

function [x,res] = jac(A,b,x0,tol,it)
x=x0;
xold=x0;
n=size(A);
k=1;
res(k)=norm(b-A*x);

while res(k)>tol && k<it
    for i=1:n
        s=0;
        for j=1:n
            if j~=i
                s=s+A(i,j)*xold(j);
            end
        end
        x(i)=(b(i)-s)/A(i,i);
    end
   
    k=k+1;
    res(k)=norm(b-A*x);
    xold=x;
end
end

[andretop00]

Grazie mille

[feddy]

Prego!

[andretop00]

Nel caso mi chiedesse come metodo di arresto il residuo relativo, va bene se aggiungo una variabile res0=norm(b-A*x0) e nel while scrivo res/res0>tol?