J'ai essayé de programmer la méthode de la puissance itérée inverse comme suit :
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
