Programmation de méthode de calcul numérique : Gauss Seidel

[invite1301f9c9]

Bonjour, alors voilà je code actuellement Gauss Seidel, et je ne comprend pas ce qui ne va pas dans mon code
. Mon vecteurX0 devrait aller en convergeant, et ce n'est pas le cas.... Il atteint des valeurs énormes alors que je travaille sur une matrice dont je connais les solutions...
Je sais plus trop quoi faire.

do
    {
        for(i = 0; i<dimension; i++)
        {
            sommeA = 0;
            sommeB = 0;
            for(j = i; j<dimension; j++)
            {
                sommeA += Matrice_A2[i][j]*VecteurX1[j];
            }
            for(j = 0; j<i; j++)
            {
                sommeB += Matrice_A2[i][j]*VecteurX2[j];
            }

            VecteurX2[i] = (Matrice_B[i]-sommeA-sommeB)/Matrice_A2[i][i];
        }

        comp_max = 0;
        for(i = 0; i<dimension; i++)
        {
             VecteurX0[i] = fabs((VecteurX2[i] - VecteurX1[i]));
             VecteurX1[i] = VecteurX2[i];
        }
        comp_max = 0;
        for(k = 0; k < dimension ; k++)
        {

            if(VecteurX0[k] > comp_max)
            {
                  comp_max = VecteurX0[k];
            }
        }


    }while(comp_max > crit_arret);

PS : J'ai publié mon problème partie mathématique également, mais je ne sais trop lequel des deux sujets est le bon.
-----

[invitea8d03b1e]

Bonsoir,

Donnez, peut-etre, un peu plus d'informations
concernant le typage des vecteurs et des variables accessoires,

mais intuitivement je crois que vous tombez içi,
dans [La non représentabilité complète de R],
vu la caractère "discrèt" innérent à n'importe quelle machine de type ~ Von Neumann,...

Par conséquent,
il faut bien avoir à l'esprit que le floating point ne fait que contourner
l'existence de "trous" (éléments) non représentable (...)

Supposition: Problème d'arrondis,... amplifiés à chaque divisions, par ex.

[invite1301f9c9]

Mon code doit être capable de gérer n'importe matrice de dimension N qu'on lui donne

Ma déclarations de variables :
double sommeA, sommeB, comp_max = 0;
double crit_arret = 1e-3;
double* VecteurX0 = NULL;
double* VecteurX1 = NULL;
double* VecteurX2 = NULL;
VecteurX0 = malloc(sizeof(double)*dimensio n);
VecteurX1 = malloc(sizeof(double)*dimensio n);
VecteurX2 = malloc(sizeof(double)*dimensio n);

VecteurX est un tableau 1 dimension de taille 'dimension', tout comme Matrice_B;
et Matrice_A est une matrice carré de taille 'dimension'.

[invitea8d03b1e]

:: comp_max ? --> semble etre globale et modifiée par fabs(...) [/!\ Effets de bord...]
Donc ... Que fait fabs(...) ?

Indépendemment de ce qu'il fait,
je vois que votre critère de sortie est [comp_max > 1e-3] ...

:: a) comp_max est remis à zero dans ton code avant d'etre testé en certain point. C'est une [inconsistence].
:: b) SI comp_max diverge, la boucle (do) principale ne TERMINE PAS.

Que fait donc fabs(...) ?

Autre chose: le 'double' implique la meme problèmatique que j'ai mentionnée pour 'float'.
Ces représentations sont binaires, autant pour les mantisses, exposant et évidemments les signes.
Incompatible donc avec des représentation décimales complète. (ce que j'entends par "trous" ds |R)

do
    {
        for(i = 0; i<dimension; i++)
        {
            sommeA = 0;
            sommeB = 0;
            for(j = i; j<dimension; j++)
            {
                sommeA += Matrice_A2[i][j]*VecteurX1[j];
            }
            for(j = 0; j<i; j++)
            {
                sommeB += Matrice_A2[i][j]*VecteurX2[j];
            }

            VecteurX2[i] = (Matrice_B[i]-sommeA-sommeB)/Matrice_A2[i][i];
        }

        comp_max = 0;
        for(i = 0; i<dimension; i++)
        {
             VecteurX0[i] = fabs((VecteurX2[i] - VecteurX1[i])); //-----> J'espère qu'il n'y a pas de "jump/out" ds fabs(...) ... [sinon faut tout revoir :( ]
             VecteurX1[i] = VecteurX2[i];
        }
//---------> Que faire de comp_max ici ? (a moins d'un bad-jump -a éviter-, dans fabs(...) \\supposition! )
        comp_max = 0; // ---------------> La valeur de comp_max n'est pas exploitée dans ce "scope"...
        for(k = 0; k < dimension ; k++)
        {

            if(VecteurX0[k] > comp_max)
            {
                  comp_max = VecteurX0[k];
//--------------> Ici on MAXIMISE (constat) la variable comp_max, OK
            }
        }


    } while (comp_max > crit_arret);

[A voir en vidéo sur Futura]

[invitea8d03b1e]

J'oubliais...

La "MAXIMISATION" de comp_max ne porte strictement que sur la dernière phase d'une itération globale.
C'est voulu ?


REFS / Normes sur les représentation binaires des 'floating point'...
http://fr.wikipedia.org/wiki/IEEE_754

[invite1301f9c9]

comp_max permet d'enregistrer la composante maximale recensé en faisant fabs(X2 - X1)
Fabs permet de prendre la valeur absolu.
Donc je vérifie que la différence entre X2 et X1 est inférieur à mon critère de précision, et si c'est pas le cas, je refait un tour, pour être plus précis, donc comp_max redevient 0, et si c'est assez précis je sors de la boucle, et dans ce cas je pourrais récupérer X2 comme vecteur casi solution.

[invitea8d03b1e]

D'accord,

mais, si com_max n'est pas référencée dans fabs(...),
la "remise-à-zéro" court-circuite comp_max qui semble n'etre pas utilisée après sa première [re]init à 0...
soit comp_max est un "alias" ou quelque chose d'opaque, soit fabs(...) l'utilise [c'est normale de se le demander dans ce context partiel], soit une intention n'était pas clair quand à sa "première utilisation locale" et peux notamment sous-entendre un "oubli" de traitement avant sa réinitialisation (avant le troisième 'for' ...) ?

Mais là,

je stagne sur la logique que je vois (celle que tu montres), et je ne peux qu'analyser sur la ~sémantique, assumant que le reste, y compris le modèle algorithmique "serait" assumé.

Alors que tu (vous?) semblez avoir un problème de résultat.

D'une manière générale je "rationnalise" toujours mon code, avant de remettre en question son "intention", en cas de ce genre de problème.

Une fois rationnalisé (TErm-graph rewriting cycles, travaux de refactoring --> simplification ou modulation du paradigme représentatif du modèle algorithmique, etc.), puis seulement, SI et seulement SI le fontionnement demeure "hors-domaine", je remets en question le typage ++profondément, puis s'il y a persistence, je révise l'"axiome" original (ou "intention"/algo-structure)... Est-il vraiment tangible, correcte et adapté (domaines mis en jeu, typage, etc.), logique, on retombe en générale quasi toujours sur des raisonnements diagrammatiques, et/ou des raisonnements "par récurrence",...

SSI (si et seulement si), l'on est certain/certaine d'etre parti(e) sur des base de programmation parfaitement rigoureuses: initialisation des variables, éviter l'effet spaguetti (structure de control + structure de pointeurs, etc.), problème de l'arrêt (sortie des boucles... Sauf en RT, mais ça c'est une autre histoire, etc.), modularité, layout du code (intelligibilité visuelle, effet de blocs imbriqués, procédure et fonctions indépendantes de global-scope!!!), transparence référentielle (pas d'effet de bord: redondant.), justement: éviter toute redondance (fonctionnelle, structure de données, alias, etc.

Ne pas hésiter à "Fonctoriser" les fonctions (c'est toujours possible en C++ et en Java par exemple...). La modélisation objet devient aussi très pratique pour des forme de Fonctorisation ++complexes au service d'une banalisation adéquatte, etc. etc. etc.

- c.f. règles de bases pour le "génie logiciel" -

:: Sorry d'etre un peu pompeux

En dehors de ça, je ne comprends pas forcément l'intention mathématique dans sont contexte et dans son "entier", donc je ne peux que "machinalement"/mécaniquement raisonner sur le code "visible" (celui que je vois içi).

[invite1301f9c9]

int GaussSeidel(int dimension, float Matrice_A3[dimension][dimension], float Matrice_B3[dimension][1])
{

    int i, j, k, condition = 0;
    double sommeA, sommeB, comp_max = 0;
    double crit_arret = 1e-3;
    double* VecteurX0 = NULL;
    double* VecteurX1 = NULL;
    double* VecteurX2 = NULL;
    VecteurX0 = malloc(sizeof(double)*dimension);
    VecteurX1 = malloc(sizeof(double)*dimension);
    VecteurX2 = malloc(sizeof(double)*dimension);
    if(VecteurX0 == NULL || VecteurX1 == NULL || VecteurX2 == NULL)
    {
        printf("Erreur dans l'allocation dynamique \n");
        return false;
    }

    for( k = 0; k<dimension; k++)
    {
        VecteurX1[k] = 1; // Cette boucle nous permet d'initialiser le premier tour du vecteur x
    }
 do
    {
        condition++;
        for(i = 0; i<dimension; i++)
        {
            sommeA = 0;
            sommeB = 0;
            for(j = i; j<dimension; j++)
            {
                sommeA += Matrice_A2[i][j]*VecteurX1[j];
            }
            for(j = 0; j<i; j++)
            {
                sommeB += Matrice_A2[i][j]*VecteurX2[j];
            }

            VecteurX2[i] = (Matrice_B[i]-sommeA-sommeB)/Matrice_A2[i][i];
        }

        comp_max = 0;
        for(i = 0; i<dimension; i++)
        {
             VecteurX0[i] = fabs((VecteurX2[i] - VecteurX1[i]));
             VecteurX1[i] = VecteurX2[i];
        }
        comp_max = 0;
        for(k = 0; k < dimension ; k++)
        {

            if(VecteurX0[k] > comp_max)
            {
                  comp_max = VecteurX0[k];
            }
        }


    }while(comp_max > crit_arret);
    printf("%d BOUCLES ! \n", condition);

    for(i=0;i<dimension;i++) {printf("V0 : %f \t  V1 : %f \t  V2: %f \n", VecteurX0[i], VecteurX1[i], VecteurX2[i]);}
return true;
}

C'est mon code entier, je suis casiment sur que mon problème vient du calcul de norme...

[invitea8d03b1e]

Je me base pour le moment sur ce papier (wiki)...

http://fr.wikipedia.org/wiki/M%C3%A9...e_Gauss-Seidel

Qui parle bien d'un algo qui évite les ....

je cite:

...mais de résoudre « en cascade » le système triangulaire inférieur (D-L)x^{(k)} = [f], ce qui est plus rapide, plus simple et numériquement moins sujet aux arrondis.

Dans l'ensemble: Est-ce que votre implémentation prends vraiment en compte les prédicâts de cette méthode.

Je ne la connais pas personnellement, et mes études dans le domaine des "méthodes numériques" sont dépassée... J'ai lu quelque chose comme l'existence d'une fonction udp(...) - ma mémoire fourche! - liée à ce genre d'itérations sur des matrices/vecteurs abordant les limites machines dans leur représentation des nombres "réels". Je ne peux pas vraiment aider plus

Si je retombe sur ce papier, j'y reviendrais.
(J'avais horreur des calculs numériques portant sur des matrices, à l'époque "forte" de FORTRAN ~ 1977-80...)
+tard je me suis "récusé", à la vue des méthodes par "L'élimination de Gauss-Jordan", pivot de Gauss, etc.