Gauss, Jacobi et Gauss-Seidel

**JonasB** · 05/03/2009, 12h30

Bonjour,

Pour un projet de math, je dois programmer en matlab les methode de Gauss, Jacobi et Gauss-Seidel.

Je suis débutant sous Matlab et j'ai enormément de mal à rédigé ces méthodes, j'ai essayé pour Gauss mais ça ne donne rien de bon :

Code :

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

1
2
3
4
5
6
7
8
j=0;
for i=1:n
    j=j+1;
    for k=i+1:n
        A(k,i)=A(i,j)*A(k,i)-A(k,j)*A(i,j);
        b(k)=A(i,j)*b(k)-A(k,j)*A(i,j);
    end
end

J'aimerais savoir si vous pouviez m'aider à le faire ou s'il existe les programmes déjà pour que je puisses m'en inspirer.

Merci d'avance

**Jerome Briot** · 05/03/2009, 13h43

Envoyé par JonasB

J'aimerais savoir si vous pouviez m'aider à le faire ou s'il existe les programmes déjà pour que je puisses m'en inspirer.

Quel est le problème exactement ?

Tu ne maîtrise les algorithme ? => Résolution des systèmes linéaires + forum Algorithme > Mathématiques

Dans le cas contraire, le code ne produit pas les bons résultats ?
Il y a des messages d'erreurs ?

**JonasB** · 05/03/2009, 13h57

J'ai déjà les algorithmes mais je n'arrive pas à les programmer !

Pour mon code, il donne bien les 0 sur les lignes qu'il faut mais il ne modifie pas les autre chiffre.

Code :

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

A=[3 5 3 ; 7 9 19 ; -4 5 11]

je trouve

Code :

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

A=[3 5 3 ; 0 9 19 ; 0 0 11]

au lieu de

Code :

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

A=[3 5 3 ; 0 -2 9 ; 0 0 81]

il donne une triangulaire inférieur alors que c'est pas le but !

**FR119492** · 05/03/2009, 15h00

Salut!
J'ai encore quelques routines que j'avais écrites en Fortran. Il te suffira de les traduire.

Code :

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
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
C **********************************************************************
C
      Subroutine L101(A,LDA,NLA,Kerr)
C
C **********************************************************************
C
C     Système linéaire
C     Factorisation
C
C     Matrice générale
C     Stockage dans un tableau à 2 indices
C
C     Jean-Marc Blanc, janvier 2005
C **********************************************************************
C
C     A l'entrée:
C
C     A      Matrice à factoriser
C
C     LDA    Nombre de lignes du tableau A tel que déclaré dans le
C            programme appelant.
C
C     NLA    Taille effective de la matrice.
C
C
C     A la sortie:
C
C     A      Matrices triangulaires résultant de la factorisation,
C            l'une triangulaire inférieure (termes diagonaux égaux
C            à 1) et l'autre triangulaire supérieure.
C
C     Kerr   Paramètre d'erreur: Kerr vaut 0 si la factorisation
C            s'est déroulée normalement. Dans le cas contraire, sa
C            valeur indique quel pivot a été trouvé nul.
C
C     Pour résoudre un système linéaire dont la matrice est stockée
C     dans le tableau A et les colonnes de seconds membres dans le
C     tableau V, on utilisera successivement les sous-programmes
C     L001, L102 et L004:
C
C        Call L101(A,LDA,NLA,Kerr)
C        If (Kerr.Gt.0) Stop
C        Call L002(A,LDA,V,LDV,NLA,NCV)
C        Call L004(A,LDA,V,LDV,NLA,NCV,Kerr)
C	
C **********************************************************************
C
      Implicit None
C
      Integer LDA,NLA,KErr
      Real*8 A(LDA,1)
C
      Integer IL,IC,II
	Real*8 S
C
C **********************************************************************
C
      Kerr=0
      If (A(1,1).Eq.0.d0) Then
        Write (*,'(1X,A)') 'L101: Pivot 1 nul'
        Kerr=1
        Return
      End If
C
      Do IL=2,NLA
        A(IL,1)=A(IL,1)/A(1,1)
      End Do
C
      Do IC=2,NLA
        Do II=IC,NLA
          S=0.d0
          Do IL=1,IC-1
            S=S+A(IL,II)*A(IC,IL)
          End Do
          A(IC,II)=A(IC,II)-S
        End Do
C
        If (IC.Lt.NLA) Then
          If (A(IC,IC).Eq.0.d0) Then
            Write (*,'(1X,A,I3,A)') 'L101: Pivot ',IC,' nul'
            Kerr=IC
            Return
          End If
          Do IL=IC+1,NLA
            S=0.d0
            Do II=1,IC-1
              S=S+A(IL,II)*A(II,IC)
            End Do
            A(IL,IC)=(A(IL,IC)-S)/A(IC,IC)
          End Do
        End If
      End Do
      Return
C
      End

Code :

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
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
C
C **********************************************************************
C
      Subroutine L002(A,LDA,V,LDV,NLA,NCV)
C
C **********************************************************************
C
C     Système linéaire
C     Résolution par substitution
C
C     Matrice triangulaire gauche
C     Tous les termes diagonaux égaux à 1
C     Stockage dans un tableau à 2 indices
C
C     Jean-Marc Blanc, janvier 2005
C
C **********************************************************************
C
C     A l'entrée:
C
C     A      Matrice du système. Aucune référence n'est faite aux
C            termes situés au-dessus de la diagonale principale et
C            sur celle-ci.
C
C     LDA    Nombre de lignes du tableau A tel que déclaré dans le
C            programme appelant.
C
C     V      Les vecteurs (colonnes) de seconds membres.
C
C     LDV    Nombre de lignes du tableau V tel que déclaré dans le
C            programme appelant.
C
C     NLA    Ordre du système.
C
C     NCV    Nombre de vecteurs de seconds membres.
C
C
C     A la sortie:
C
C     V      Les vecteurs (colonnes) de la solution.
C
C **********************************************************************
C
      Implicit None
C
      Integer LDA,LDV,NLA,NCV,Kerr
      Real*8 A(LDA,1),V(LDV,1)
C
      Integer IL,IS,IC
	Real*8 S
C
      If (NLA.Eq.1) Return
C
      Do IL=2,NLA
        Do IS=1,NCV
          S=0.d0
          Do IC=1,IL-1
            S=S+A(IL,IC)*V(IC,IS)
          End Do
          V(IL,IS)=V(IL,IS)-S
        End Do
      End Do
      Return
C
      End

Code :

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
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
C
C **********************************************************************
C
      Subroutine L004(A,LDA,V,LDV,NLA,NCV,Kerr)
C
C **********************************************************************
C
C     Système linéaire
C     Résolution par substitution
C
C     Matrice triangulaire droite
C     Stockage dans un tableau à 2 indices
C
C     Jean-Marc Blanc, janvier 2005
C
C **********************************************************************
C
C     A l'entrée:
C
C     A      Matrice du système. Aucune référence n'est faite aux
C            termes situés au-dessous de la diagonale principale.
C
C     LDA    Nombre de lignes du tableau A tel que déclaré dans le
C            programme appelant.
C
C     V      Les vecteurs (colonnes) de seconds membres.
C
C     LDV    Nombre de lignes du tableau V tel que déclaré dans le
C            programme appelant.
C
C     NLA    Ordre du système.
C
C     NCV    Nombre de vecteurs de seconds membres.
C
C
C     A la sortie:
C
C     V      Les vecteurs (colonnes) de la solution.
C
C     Kerr   Paramètre d'erreur: Kerr vaut 0 si la résolution s'est
C            déroulée normalement. Dans le cas contraire, sa valeur
C            indique quel pivot a été trouvé nul.
C
C **********************************************************************
C
      Implicit None
C
      Integer LDA,LDV,NLA,NCV,Kerr
      Real*8 A(LDA,1),V(LDV,1)
C
      Integer IL,IS,IC
	Real*8 P,S
C
      Kerr=0
      If (A(NLA,NLA).Eq.0.d0) Then
        Write (6,'(1X,A,I3,A)') 'L013: Pivot ',NLA,' nul'
        Kerr=NLA
        Return
      End If
C
      Do IS=1,NCV
        V(NLA,IS)=V(NLA,IS)/A(NLA,NLA)
      End Do
      If (NLA.Eq.1) Return
C
      Do IL=NLA-1,1,-1
        P=A(IL,IL)
        If (P.Eq.0.d0) Then
          Write (6,'(1X,A,I3,A)') 'L013: Pivot ',IL,' nul'
          Kerr=IL
          Return
        End If
C
        Do IS=1,NCV
          S=0.d0
          Do IC=IL+1,NLA
            S=S+A(IL,IC)*V(IC,IS)
          End Do
          V(IL,IS)=(V(IL,IS)-S)/P
        End Do
      End Do
      Return
C
      End

Jean-Marc Blanc

**JonasB** · 05/03/2009, 15h28

Je connais pas du tout le fortran, j'en avais jamais entendu parler, ça m'aide pas beaucoup mais merci quand meme !

**FR119492** · 05/03/2009, 23h11

Salut!

Je connais pas du tout le fortran, j'en avais jamais entendu parler

Pour ton information, c'est le premier langage conçu spécifiquement pour le calcul numérique (1956), donc pour appliquer des algorithmes comme Gauss, Jacobi, Gauss-Seidel, etc., et il n'a cessé d'évoluer depuis lors, ce qui fait que c'est le plus efficace; la preuve en est qu'Intel diffuse gratuitement un compilateur Fortran, afin de mettre en évidence les performances de ses processeurs.

De toute manière, quand tu connais 42 langages, le 43ème est facile à apprendre.

ça m'aide pas beaucoup

Au contraire, ça peut t'aider beaucoup: il te suffit de le traduire ligne par ligne dans le langage de programmation de Matlab. Tu obtiendras un programme qui ne sera peut-être pas optimal, mais qui devrait fonctionner.

Alors, pour commencer, toutes lignes qui ont un C en première position sont des commentaires. Ensuite, chaque fois que tu rencontreras une ligne que tu ne comprends pas, on sera là pour t'expliquer.

Au travail!
Jean-Marc Blanc

**Jerome Briot** · 06/03/2009, 08h23

Envoyé par FR119492

on sera là pour t'expliquer.

Dans le forum Fortran bien entendu

Sinon, si tu veux continuer à corriger ton code MATLAB, tu devrais regarder l'évolution des indices et des variables dans les boucles et les comparer simplement aux valeurs que tu trouverais à la main :

Code :

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

1
2
3
4
5
6
7
8
9
10
11
12
j=0;
for i=1:n
    i
    j=j+1
    for k=i+1:n
        k
        A(i,j)
        A(k,i)
        A(k,j)
        A(k,i)=A(i,j)*A(k,i)-A(k,j)*A(i,j)
    end
end

Cette méthode brutale de recherche d'erreur est possible ici car ton problème est de taille réduite (3x3)

**JonasB** · 07/03/2009, 00h10

J'ai essayé de regarder le Fortran, j'ai un peu de mal à comprendre !
Je me suis dit que j'allais d'abord essayé la méthode de Dut donc après quelques bugs en tout genre, j'ai réussi à produire le code pour Gauss et il semble fonctionner !
Pareil pour Jacobi et Gauss-seidel !
Merci pour votre aide !
J'ai besoin de tracer pour jacobi le vecteur erreur en fonction des itération sur une échelle logarytmique ! J'ai les x, les itérations !
Pour le vecteur erreur j'ai une matrice comme celle-ci

Code :

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

1
2
3
4
5
6
7
YX=Y-SX
 
YX =
 
    4.0000   -1.0222    0.8000   -0.0277    0.1388    0.0377    0.0324    0.0169    0.0100    0.0059
    6.4000   -2.1333    0.0978   -0.1659   -0.0946   -0.0351   -0.0323   -0.0138   -0.0099   -0.0052
    3.3333    0.2667    0.0148    0.2504    0.0184    0.0620    0.0184    0.0162    0.0079    0.0050

J'aimerais calculer la norme de chaque colonne et la mettre dans une matrice donc avec une ligne et d'autant de colonne ! Un truc comme ça :

[norme(1er colonne) norme(2emcolonne)

Je n'arrive pas à isolé chaque colonne et à calculer sa norme !

**Jerome Briot** · 07/03/2009, 05h53

Jette un oeil ici : http://briot-jerome.developpez.com/m...tion/#LIII-D-2

**JonasB** · 07/03/2009, 18h17

Oui, c'est bon merci, j'ai fait une boucle en prenant chaque fois la norme de la colone indice i pour la mettre dans une autre !

Par contre, je dois donc tracer une courbe sur une échelle logarithmique de l'erreur en fonction des itérations ! La c'est plus une question de compréhension ! Il faut mettr log(erreur)??

Et sinon, il n'y a pas un beau tutoriel pour la mise en forme des courbes sur matlab parce que je ne connais que plot !

**Jerome Briot** · 07/03/2009, 18h30

Tu peux utiliser la fonction SEMILOGY à la place de PLOT ici

Ensuite pour modifier le graphique, deux possibilités

Soit tu traces le graphique et ensuite tu utilises le menu Tools > Edit Plot (ou l'icone en forme de flèche blanche dans la barre d'outil). une fois que tu as la main, tu double clique sur la ligne.

Soit tu utilises SEMILOGY comme ceci :

Code :

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