(pour l'enseigner l'an prochain, par exemple)

jeudi 2 mai 2013

diagonalisation d'une matrice

Avec numpy:

>>> from numpy import *

on donne une matrice $3\times 3$ ligne par ligne :


>>> A = matrix([[1,2,3],[4,5,6],[7,8,1]])

>>> A
matrix([[1, 2, 3],
        [4, 5, 6],
        [7, 8, 1]])


le tableau de ses valeurs propres :
>>> eigvals(A)
array([ 12.45416474,  -0.37976219,  -5.07440255])

$3$ racines réelles simples, $A$ est donc diagonalisable, voici ses valeurs propres suivies de la matrice de passage dans la base des vecteurs propres :
>>> eig(A)
(array([ 12.45416474,  -0.37976219,  -5.07440255]), matrix([[-0.29373774, -0.73967882, -0.29720654],
        [-0.69005397,  0.66500848, -0.39870229],
        [-0.66147083, -0.10314536,  0.86758559]]))
>>> P=eig(A)[1]

on construit une matrice diagonale à partir du tableau des valeurs propres :
>>> D=diag(eig(A)[0])

et on vérifie que $A = P D P^{-1}$
>>> A-P*D*inv(P)
matrix([[ -2.22044605e-15,  -1.33226763e-15,  -4.44089210e-16],
        [ -1.77635684e-15,   0.00000000e+00,   1.77635684e-15],
        [ -4.44089210e-15,  -5.32907052e-15,  -2.66453526e-15]])

comme mumpy fait des calculs approchés, le résultat n'est pas exactement $0$, mais pas loin...

lundi 29 avril 2013

calcul approché de $\pi$ avec la moyenne arithmético-géométrique

Suivant un vieux sujet du CAPES (1995), voici programmé en python (avec sympy) un algorithme inspiré de [Salamin 1976]:


from sympy import *
import sympy as sp

x = sp.symbols('x')

def fu(n,x):
    if n==0: return x
    else: return (fu(n-1,x)+fv(n-1,x))/2

def fv(n,x):
    if n==0: return 1
    else: return sqrt(fu(n-1,x)*fv(n-1,x))

def fy(n,y):
    return (fu(n,x)/fv(n,x)).subs(x,y)
    
def fz(n,y):
    return (diff(fv(n,x),x)/diff(fu(n,x),x)).subs(x,y)
    
def fpi(n):
    if n==0: return (2+sqrt(2))
    else: return fpi(n-1) * (1 + fy(n,1/sqrt(2))) / (1 + fz(n,1/sqrt(2)))
    
for n in range(10): 
    p = 1000
    print('n='+ str(n) 
          + ': erreur 10**('
          + str(N(ln(N(fpi(n)-pi,p))/ln(10),5)) 
          + ')')


en faisant les calculs approchés avec $1000$ chiffres significatifs, l'exécution de ce code donne $\pi$ avec 700 décimales au bout de seulement 8 étapes (la précision double à chaque étape, en fait):

n=0: erreur 10**(-0.56444)
n=1: erreur 10**(-2.9939)
n=2: erreur 10**(-8.1322)
n=3: erreur 10**(-18.737)
n=4: erreur 10**(-40.262)
n=5: erreur 10**(-83.619)
n=6: erreur 10**(-170.64)
n=7: erreur 10**(-344.98)
n=8: erreur 10**(-693.95)