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...
(pour l'enseigner l'an prochain, par exemple)
jeudi 2 mai 2013
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)
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)
Inscription à :
Commentaires (Atom)