Méthode de Gauss Seidel - scilab — Les-mathematiques.net The most powerful custom community solution in the world

Méthode de Gauss Seidel - scilab

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

Réponses

  • 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
    
    
Connectez-vous ou Inscrivez-vous pour répondre.
Success message!