Discussion :

[Caroline_1]

Bonjour,

J'ai essayé de programmer la méthode de la puissance itérée inverse comme suit :

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
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
 
      program puissance_iteree_inverse
      implicit double precision (a-h,o-z)
      dimension a(1000,1000),b(1000),p(1000,1000)
      dimension x(1000),y(1000)
 
 
c     ------------------------------------------------------------------
c     SAISIES
c     ------------------------------------------------------------------
 
 
c     Saisie de la taille de la matrice
      write(*,*)"Saisissez la taille de la matrice :"
      read(*,*) m
 
      n=m*m
 
 
c     Saisie de la pr‚cision voulue
      write(*,*)"\nSaisir la pr‚cision:"
      read(*,*) eps
 
 
c     ------------------------------------------------------------------
c     ALGORITHME DE LA PUISSANCE ITEREE INVERSE
c     ------------------------------------------------------------------
      n=3
 
      do j=1,m
      do i=1,m
      a(i,j)=0.000001d0
      enddo
      enddo
 
      do i=1,3
      a(i,i)=i*1.d0
      enddo
 
      do i=1,n
      x(i)=1.d0
      enddo
 
      w=1.d0  ! lambda_anc
      z=0.d0  ! lambda
      cpt=0
      call decomposition_A_LU(a,p,n)
      do j=1,15
c      do while(abs(z-w).GT.eps)
      w=z
      call prod_scal(x,x,r,s,n)
      do i=1,n
      b(i)=(x(i)*1.d0)/s
      print *,b(i)
      enddo
      call descente_triangle(p,b,y,n)
      call remontee_triangle(a,y,x,n)
      call prod_scal(x,b,t,v,n)
      z=(1.d0)/t
      cpt=cpt+1
      enddo
c      enddo
 
 
      write(*,*)"\nLa valeur propre de plus petit module de A est  :"
      print *,z
 
      write(*,*)"\nLe nombre d'it‚rations pour cette pr‚cision est :"
      print *,cpt
 
 
      end
 
 
 
c     ------------------------------------------------------------------
c     DECOMPOSITION A=LU
c     ------------------------------------------------------------------
 
      subroutine decomposition_A_LU(a,p,n)
      implicit double precision (a-h,o-z)
      dimension a(1000,1000),p(1000,1000)
 
c     Algorithme de d‚composition A=LU
c     Au d‚but a=A                A la fin P=L et a=U
 
      do k=1,n-1
      do i=1,k-1
      p(i,k)=0.d0
      enddo
      p(k,k)=1.d0
      do i=k+1,n
      p(i,k)=(a(i,k)*1.d0)/a(k,k)
      enddo
      do i=1,k
      do j=1,n
      a(i,j)=a(i,j)
      enddo
      enddo
      do i=k+1,n
      do j=1,k
      a(i,j)=0.d0
      enddo
      enddo
      do i=k+1,n
      do j=k+1,n
      a(i,j)=a(i,j)-p(i,k)*a(k,j)
      enddo
      enddo
      enddo
      do i=1,n-1
      p(i,n)=0.d0
      enddo
      p(n,n)=1.d0
 
      return
      end
 
c     ------------------------------------------------------------------
c     DESCENTE TRIANGLE
c     ------------------------------------------------------------------
 
      subroutine descente_triangle(p,b,y,n)
      implicit double precision (a-h,o-z)
      dimension p(1000,1000),b(1000),y(1000)
 
      y(1)=b(1);
      do i=2,n
      y(i)=(b(i)-p(i,i-1)*y(i-1))
      enddo
 
      return
      end
 
c     ------------------------------------------------------------------
c     REMONTEE TRIANGLE
c     ------------------------------------------------------------------
 
      subroutine remontee_triangle(u,y,x,n)
      implicit double precision (a-h,o-z)
      dimension u(1000,1000),y(1000),x(1000)
 
      x(n)=(y(n)*1.d0)/u(n,n);
      do i=n-1,1
      x(i)=((y(i)-u(i,i+1)*x(i+1))*(1.d0))/u(i,i)
      enddo
 
      return
      end
 
c     ------------------------------------------------------------------
c     PRODUIT SCALAIRE ET NORME
c     ------------------------------------------------------------------
 
      subroutine prod_scal(x,y,r,s,n)
      implicit double precision (a-h,o-z)
      dimension x(1000),y(1000)
 
      r=0.d0
 
      do i=1,n
        r=r+x(i)*y(i)
      enddo
      s=dsqrt(r)
      return
      end

Pourtant ça ne marche pas et je ne sais pas pourquoi. J'ai utilisé l'algo que l'on peut retrouver notamment à cette adresse :
www-pequan.lip6.fr/~jmc/polycopies/cours5-puissanceiteree.pdf

Pourriez-vous m'aider ??
A bientôt et merci

[Ehouarn]

Bonjour,

ça ne marche pas

Il serait bien plus productif que tu précise ce qui ne marche pas, et en quoi (message d'erreur à la compilation? à l'exécution? Résultats clairement faux?) ...

Attention également aux variables "écrites en dur" du type:

Code :

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

1
2
3
4
5
 
c     ------------------------------------------------------------------
c     ALGORITHME DE LA PUISSANCE ITEREE INVERSE
c     ------------------------------------------------------------------
      n=3

alors que quelques lignes plus haut:

Code :

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

1
2
3
4
5
6
 
c     Saisie de la taille de la matrice
      write(*,*)"Saisissez la taille de la matrice :"
      read(*,*) m
 
      n=m*m

Au delà de ça, je me permets de rappeler quelques méthodes/préceptes pour tester/valider un code:

• Bien tester séparément chaque routine, en affichant chaque étape des calculs, sur des cas simples dont on connait par ailleurs la solution.
• Utiliser des des "implicit none" et d'une façon générale toutes les options de compilations aidant à isoler/repérer les erreurs.

[Caroline_1]

Bonjour,

Merci pour la réponse. Alors ce qui ne marche pas c'est que le programme ne m'indique aucune erreur. Néanmoins il me donne une valeur fausse.
Pour ne n=3 c'était pour un test.
Toutes les subroutines marchent je les ai testées séparemment.

[FR119492]

Salut!

J'ai essayé de programmer la méthode de la puissance itérée inverse comme suit :

Si c'est un exercice dans le cadre de tes études, tu es bien obligée de faire ce qu'on t'a demandé. Sinon, tu aurais meilleur temps d'utiliser des bibliothèques toutes faites, comme on en trouve sur Internet, car elles ont été écrites par des gens beaucoup plus compétents que toi et moi, puis testées par des centaines d'utilisateurs.

Pour revenir à ton code, j'ai quelques remarques:

1. Comme Ehouarn te l'a déjà indiqué, il est prudent d'utiliser l'instruction Implicit None pour éviter l'assignation automatique du type à chaque variable; ainsi, de nombreuses erreurs sont détectées à la compilation et non à l'exécution, et signalées par un message beaucoup plus clair.
2. Tu appliques la méthode de factorisation LU sans pivotage, ce qui peut conduire aux pires choses si un des pivots successifs est nul ou presque.
3. Les méthodes de la puissance itérée ne convergent pas lorsque la valeur propre cherchée est multiple, ce qui est le cas pour ta matrice.
4. Dans le cas que tu as choisi comme exemple, la matrice est symétrique; mais si elle ne l'était pas, tu aurais des valeurs propres complexes et ton programme ne fonctionnerait pas.

Jean-Marc Blanc