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 |
Puissance itérée inverse
Répondre avec citation
Partager