Discussion :

[kelyos]

Bonjour à tous, je vais présenter un peu ce que je compte faire.

Il me faut calculer à partir de la matrice (bande), trouver la solution x par la méthode du gradient conjugué sous l'équation Mx=B (M la matrice, B le second membre un vecteur et x aussi).

Par l'algorithme que je trouve sur le net, j'aimerais savoir si parmi vous qui ont codé cet algo me dire si mon code est bon.

(Pour ceux ou celles qui veulent l'intégralité des fichiers et codes je les ai laissé ici http://www.faceweb.fr/Laplace2D.rar)

Code C++ :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
Vecteur MatriceBande::GradientConjugue(const Vecteur &B, double epsilon)
{
    MatriceBande &M = (*this);
    Vecteur r0(N), r1(N), x0(N), p0(N);
 
    double a, b;
    int n = 0, Nmax=100000;
 
    r0 = B - M*x0;
    p0 = r0;
    b = 0;
 
    while (r0.NormeCarree() >= epsilon*epsilon)
    {
        n++;
 
        if (n>Nmax)
        {
            cout << " Le gradient conjugue n'a pas converge en moins de "<< Nmax <<" iterations." << endl;
            cout << " (||r||^2 = " << r0.NormeCarree()<< " )" << endl;
            break;
        }
 
        a = r0.NormeCarree() / (M*p0*p0);
 
        x0 = x0 + a*p0;
        r1 = r0;
        r0 = r0 - a*(M*p0);
        p0 = r0 + b*p0;
 
        b = r0.NormeCarree() / r1.NormeCarree();
    }
 
    return x0;
}

[Pierrot_gre]

Si tu veux résoudre l'équation Mx=B par un algo du gradient, je suppose que tu veux minimiser la fonction:
phi(x) = ||Mx-B||^2
Déjà, si tu as une méthode pour calculer l'inverse d'une matrice, il existe une solution exacte:
x = M^-1B, ou si M est pas inversible:
x = (M^TM)^-1M^TB

Sinon avec le gradient, il te faut calculer le gradient de phi par rapport à x. C'est:
-2M^T(B-Mx)
dans ton code c'est r0 et j'ai l'impression qu'il manque le M^T devant.

Ensuite tu utilises (je crois), un pas "optimal" pour faire converger ton gradient. Je te conseille dans un 1er temps d'utiliser un pas constant et d'accroitre le nombre d'itération max pour voir si ça converge.

Si oui, c'est-à-dire que ton implémentation de la dérivée est bonne, tu peux passer à un pas optimal pour accélérer la convergence.