View file src/j/deflation.ijs - Download

NB. Calcul des valeurs et vecteurs propres d'une matrice par la méthode de la puissance itérée et de la déflation
NB. Référence : 
NB.  http://www.bibmath.net/dico/index.php?action=affiche&quoi=./m/methodepuissance.html
NB.  http://log.chez.com/text/math/methodepuissance.pdf
NB.  http://www.normalesup.org/~pastre/meth-num/MN/9-val-pro/cours-valeurspropres.pdf
NB.  http://log.chez.com/text/math/cours-valeurspropres.pdf

and =: *.
or =: +.
not =: -.

transpose =: |:
extprod =: */
matprod =: +/ . *

avg =: 3 : 0
 (+/ y) % #y
)

pavg =: 4 : 0
 +/ y * x % +/ x
)

puissance =: 3 : 0
 A =. y
 n =. #A
 NB. u =. 1 , (n-1) $ 0
 u =. (?. n $ 1000) % 1000
 u =. u % (u matprod u)^0.5
 i =. 0
 while. (i < 10000)
 do.
  i =. i + 1
  u1 =. u
  u =. A matprod u
  u =. u % (u matprod u)^0.5
  d =. u - u1
  if. (d matprod d) < 1e_30
  do. break.
  end.
 end.
 NB. l =. avg (A matprod u) % u
 l =. (u * u) pavg (A matprod u) % u
 l; u; i
)

deflation =: 3 : 0
 A =. y
 NB. B = transpose A
 n =. #A
 l =. 0 $ 0
 u =. (n,0) $ 0
 for. A
 do.
  NB. echo 'A='; A
  lu1 =. puissance A
  NB. echo lu1
  l1 =. > 0 { lu1
  u1 =. > 1 { lu1
  l =. l , l1
  u =. u ,. u1
  lu2 =. puissance (transpose A)
  l2 =. > 0 { lu2
  u2 =. > 1 { lu2
  A =. A - l1 * u1 extprod u2 % u1 matprod u2
  NB. B =. B - l1 * u2 extprod u1 % u2 matprod u1
 end.
 l ; u
)