Méthode de Gauss Seidel

[Chaka]

Bonsoir, je bloque sur une question de mon TP d'analyse numérique...
Je me permets de mettre mon script ci-dessous pour que vous puissiez me le corriger (le langage utilisé est du Scilab):
 

function [x, iter]=gaussSeidel(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)=(A(i,1:i-1)*x0(1:i-1))
            xnew(i)=b(i)-xnew(i)-(A(i,i+1:n)*x(i+1:n))
            xnew(i)=xnew(i)/A(i,i)
        end 
    x = xnew; r=b-A*x; err = norm(r)/normb;
    end
endfunction

Le problème viendrai apparemment de la ligne 11 :  xnew(i)=xnew(i)/A(i,i) Indice invalide et Somme avec une matrice vide ! 

Ci-joint les scripts associés à mon TP.
Merci par avance pour vos retours.

Chaka

jacob1i.sce jacob31i.sce jacobi.sce jacobi0.sce main.sce

[julesx]

Bonsoir,

Tu auras peut-être du mal à avoir une réponse détaillée car j'ai bien l'impression que personne, sur ce site, ne maitrise Scilab. C'est également mon cas, j'essaie cependant de te donner un semblant d'explication.

Déjà, dans cette ligne

xnew(i)=(A(i,1:i-1)*x0(1:i-1)

ce doit être x, pas x0, mais c'est peut être simplement une faute de frappe.

Un essai de ta démarche avec Python conduit également à un échec car xnew est un réel et, pour ce logiciel, le produit partiel A(i,1:i-1)*x(1:i-1) retourne un autre type de donnée.

J'ai trouvé sur la toile un script Python où cette partie est traitée ainsi :

for i in range(len(A)):
            s=0
            xprim=x
            for j in range(len(A)):
                  if i!=j:
                        s=s+A(i,j)*(j)
            x=w*((b-s)/A)+(1-w)*x w est ce qu'ils appellent le facteur de relaxation

En d'autres termes, on calcule la somme d'un seul coup au lieu de la décomposer, ce qui retourne bien un réel. A noter que, pour Python, les indices sont entre crochets, mais l'éditeur de e-bahut interprète i entre crochets comme la mise en italique de ce qui suit !

Bon, peut-être que quelqu'un d'autre t'aidera davantage.

 

[Chaka]

Merci pour ton aide
J'ai finalement réussi et vous mets ci-joint le script final.

FINAL GAUSS SEIDEL.sce

[julesx]

Merci pour le retour, par contre, ne disposant pas de Scilab, je ne peux pas afficher ton script. Je sais que c'est un logiciel libre, mais je n'ai pas le courage de m'attaquer à ce dernier.

[Chaka]

J'ai fait le calcul des sommes séparément.

function [S1]=somme1(A, i, x)
    S1=0
    n=size(A,1);
    for j=1:n
        if j>=i+1
            S1=S1+A(i,j)*x(j)
        end
    end
endfunction

function [S2]=somme2(A, i, xnew)
    S2=0
    n=size(A,1);
    for j=1:n
        if j<i
            S2=S2+A(i,j)*xnew(j);
        end
    end
endfunction

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

 

[julesx]

Bonjour,

Comme dit, je ne connais pas Scilab, donc je ne peux pas apprécier. Mais l'essentiel est que ton script te donne le bon résultat.

Bonne continuation.