Méthode de la puissance pour le calcul des valeurs propres -- code phyton

methode_puissance.py

# coding: utf-8

## Méthode de la puissance pour le calcul de valeurs propres

# In[1]:

from numpy import *
import matplotlib as plt
import pylab as pl


# Une matrice 2x2 à coefficients positifs

# In[2]:

A = array([[ 0.00356101,  0.53584758],
       [ 0.33824661,  0.02571842]])


# On génère 200 points le long du cercle unité $\{ {}^t(\cos(\theta), \sin(\theta) ) | \theta \in [0, 2\pi] \}$

# In[3]:

theta = linspace(0,2*pi,200)
v = array([cos(theta), sin(theta)])


# On itère les puissances successives normalisées $$\frac{A^i v}{|| A^i v ||},$$ pour <var>i</var> allant de 1 à 50. 

# In[4]:

for i in range(0,50): # on itere Av/||v||
    v = dot(A,v) # 1 produit Av
    mv = max(linalg.norm(v,axis=0))
    v = v/mv
    # plot un r sur 10
    if i % 10 == 0:
        pl.plot(v[0,:],v[1,:])

pl.show()


# Les vecteurs de <var>v</var> doivent converger vers un vecteur <var>x</var> qui satisfait $Ax = \lambda x$. Le quotient de Rayleigh en $v$ est $r(v) = {}^tv A v/{}^t v v$. On prend par exemple comme approximation de <var>x</var>, le vecteur <code>v[:,0]</code>, mais on aurait pu en prendre un autre vecteur, sachant que certains sont plus près du vrai vecteur propre. 

# In[5]:

vk = v[:,0]
rv = dot(vk.T,dot(A,vk))/dot(vk.T,vk)


# In[6]:

print "vecteur propre {v}\nvaleur propre {w}\n".format(v = vk, w = rv)
ws, vs = linalg.eig(A)
print "linalg.eig(A):\nvecteurs propres {v}\nvaleurs propres {w}".format(v = vs, w = ws)