Méthode de Gauss Seidel - scilab

Chaka

Bonsoir,
Je me permets de vous consulter car je bloque sur une question de mon TP d'analyse numérique.
Dans ce TP, nous devons implémenter la méthode de Jacobi et de Gauss Seidel.
J'ai donc réussi pour la méthode de Jacobi :

function[x,iter] = jacobi(A,b,x0,nmax,tol)
    n = size(A,1);
    iter = 0;
    r = b-A*x0; normb = norm(b); err = norm(r)/normb; x=x0;
    while err > tol & iter < nmax
        iter = iter + 1;
        xnew = zeros(n,1);
        for i = 1:n
            xnew(i)=(b(i)-A(i,[1:i-1 i+1:n])*x([1:i-1 i+1:n]))/A(i,i)
        end 
    x = xnew; r=b-A*x; err = norm(r)/normb;
    end
endfunction

À présent, j'ai besoin de modifier ce script pour obtenir la méthode de Jacobi, voici une de mes tentatives :

function[x,iter] = gaussSeidel2(A,b,x0,nmax,tol)
    n = size(A,1)
    x0=zeros(n,1)
    iter = 0;
    r = b-A*x0; normb = norm(b); err = norm(r)/normb; x=x0;
    while err > 10^(-8) & iter < 1000
        iter = iter + 1;
        xnew = zeros(n,1);
        for i = 1:n
            xnew(i)=(1/A(i,i))*(b(i)-A(i,1:i-1)*xnew(1:i-1)-A(i,i+1:n)*x(i+1:n))
        end 
    x = x; r=b-A*x; err = norm(r)/normb;
    end
endfunction

Cependant Scilab me renvoie une erreur...
Pouvez-vous me donner une piste de réflexion et/ou corriger mon précédent script ?
Merci par avance !
Chaka

tatof

Je commente juste ton deuxième code pour commencer:

while err > 10^(-8) & iter < 1000

Il faut que tu remplaces $10^{-8}$ et $1000$ par les variables de ta fonction. Deuxième chose ce n'est pas x=x mais x=x_new si je compare à ton premier script. De même pour $r=b-Ax$. De plus, tu peux mettre ces $2$ codes ensemble et choisir une variable de sélection de la méthode.

En regardant la page Wikipédia concernée, j'ai fais un code Octave (ne maitrisant pas scilab) qui te montre l'implémentation de ces méthodes. Tu peux essayer de le convertir en Scilab.

function [sol,iter]=meth_jacobi_gauss(A,b,x0,iter_max,seuil_err,choix_meth)

%% Paramètres
% Nous partons du système Ax=b. La matrice A est décomposée: A=M-N.
% Pour la méthode de Jacobi, M=D (matrice diagonale), N=E+F (E matrice
% triangulaire inférieure à diagonale nulle et F matrice triangulaire supérieure à diagonale nulle). 
% Pour la méthode de Gauss-Seidel, M=D-E, N=F.

% iter_max: nombre d'itérations maximales (nombre entier positif)
% seuil_err: seuil d'erreur maximale (nombre réel positif)
% choix_meth=0 (Jacobi), 1 (Gauss-Seidel)

%%
D=diag(diag(A));
E=tril(A)-D;
F=triu(A)-D;

%% Vérification taille
if (size(A,2)~=size(b,1))
    disp('Revoir taille de b');
end

if (size(x0,1)~=size(b,1))
    disp('Revoir la taille de x0 et b');
end

%% Vérification que A est une matrice à diagonale strictement dominante pour avoir convergence
S=abs(E)+abs(F);
S2=sum(S,2);
D2=diag(D);
for i=1:size(S2,1)
    if (D2(i)<S2(i))
        disp('A pas à diagonale strictement dominante');
    end
end

%% Choix méthode
if (choix_meth==0)
    M=D;
    N=-(E+F);
elseif (choix_meth==1)
    M=D+E;
    N=-F;
else
    disp('Méthode non disponible');
end

%% Algorithme
B=inv(M)*N;
invM=inv(M);
sol_pre=x0;

for iter=1:iter_max
    sol_ac=B*sol_pre+invM*b;
    err=abs(sol_ac-sol_pre);
 
    %% Erreur
    r=D*err; %vecteur résidu
    r2=r/norm(b);
    if (r2<repmat(seuil_err,size(r2,1),1))
        break; %arrêt de l'algorithme car condition sur l'erreur remplie
    else
        sol_pre=sol_ac; %continuer l'itération
    end
end
sol=sol_ac;
end