Corrigé TD3 ¶
Classe Vecteur¶
In [1]:
class Vecteur:
"""
Une classe de Vecteur
"""
data = []
size = 0
def __init__(self,taille=0,valeur=0.0,liste=[]):
"Constructeur a partir de taille et valeur ou d'une liste"
# tentative de conversion en float plutôt que vérification du type
if int(taille)>0:
self.size = int(taille)
self.data = [float(valeur)]*self.size
elif liste != []:
self.size = len(liste)
self.data = list(liste) # copie la liste
for e in self.data: e = float(e)
def apply(self,f):
"Applique une la fonction f au vecteur: "
from types import FunctionType
if not type(f) is FunctionType:
print "Erreur: l'argument n'est pas une fonction mais",type(f)
return NotImplemented
# par tentative de conversion
for i in range(self.size):
self.data[i] = float(f(self.data[i]))
def __str__(self):
"affichage de l'objet avec le mot-cle print"
s = "[ "
for i in range(self.size):
s += str(self.data[i])+" "
return s+"]"
def __add__(self,other):
"additionne deux vecteurs et retourne le vecteur résultant"
if not other.__class__ is Vecteur:
print "Erreur: l'argument n'est pas un vecteur"
return NotImplemented
if other.size != self.size:
print "Erreur: les vecteurs n'ont pas la même taille"
return NotImplemented
res = Vecteur(taille=self.size,valeur=0.0)
for i in range(self.size):
res.data[i] = self.data[i] + other.data[i]
return res
def __sub__(self,other):
"soustrait deux vecteurs et retourne le vecteur résultant"
if not other.__class__ is Vecteur:
print "Erreur: l'argument n'est pas un vecteur"
return NotImplemented
if other.size != self.size:
print "Erreur: les vecteurs n'ont pas la même taille"
return NotImplemented
res = Vecteur(taille=self.size,valeur=0.0)
for i in range(self.size):
res.data[i] = self.data[i] - other.data[i]
return res
def __iadd__(self,other):
"additionne un réel à tous les élements"
for i in range(self.size):
self.data[i] += float(other)
return self
def __isub__(self,other):
"soustrait un réel à tous les élements"
for i in range(self.size):
self.data[i] -= float(other)
return self
def __mul__(self,other):
"retourne le produit scalaire de deux vecteurs"
if not other.__class__ is Vecteur:
print "Erreur: l'argument n'est pas un vecteur"
return NotImplemented
if other.size != self.size:
print "Erreur: les vecteurs n'ont pas la même taille"
return NotImplemented
res = 0.0
for i in range(self.size):
res += self.data[i] * other.data[i]
return res
def norm(self):
"retourne la norme"
from math import sqrt
return sqrt(self*self)
def __len__(self):
"renvoie la taille"
return self.size
def __getitem__(self,cle):
"accès en lecture aux éléments, marche avec les slices"
return self.data.__getitem__(cle)
def __setitem__(self,cle,valeur):
"accès en écriture aux éléments"
return self.data.__setitem__(cle,valeur)
def norm(v):
return v.norm()
Tests du constructeur et de l'affichage
In [2]:
print Vecteur(), Vecteur(3), Vecteur(taille=3,valeur=1), Vecteur(taille=2,valeur='1.0'), Vecteur(liste=range(3))
a = Vecteur(taille=3,valeur=1)
print type(Vecteur), type(a), a.__class__
Test de l'application d'une fonction
In [3]:
v = Vecteur(liste=range(10))
from math import sin,pi
v.apply(lambda x: sin(pi*x/6.0))
print v
Tests des opérateurs courants
In [4]:
x, y = Vecteur(liste=range(8)), Vecteur(liste=range(0,16,2))
x += 0.5
print x, y, x - y
print x*y, x.norm(), norm(x)
In [5]:
x[1] = 2.1
print x[:3]
Documentation
In [6]:
help(Vecteur)
Algèbre linéaire avec Numpy¶
In [7]:
from numpy import *
set_printoptions(precision = 4, suppress=True)
L = 10
In [8]:
x = arange(1,L+1)
V = sqrt(2/float(L+1.))*sin(pi*x/(L+1))
print V,linalg.norm(V)
In [9]:
W = diag([-2.]*L,k=0) + diag([1.]*(L-1),k=1) + diag([1.]*(L-1),k=-1)
print W
In [10]:
A = dot(W,V)
print A / V
In [11]:
(val,U) = linalg.eigh(W)
print "valeurs propres:", val
print "vecteurs propres:"
for i in range(L): print i,":",U[:,i]
In [12]:
Um1 = linalg.inv(U)
D = diag(val)
print dot(U, dot(D,Um1))
print dot(Um1, dot(W,U))
Application à la dynamique d'une chaîne d'oscillateurs couplés¶
In [13]:
Omega = sqrt(-val); print Omega
Conditions initiales $Y^0 = U^{-1}X^0$
In [14]:
X0 = zeros( (L,) ); X0[0] = 0.2; X0[1] = 0.1; print X0
# en partant d'un état propre
# X0 = U[:,1]
In [15]:
Y0 = dot(Um1,X0); print Y0
Solution des équations: $X(t) = UY(t)$ avec $Y_j(t) = Y^0_j\cos(\omega_j t)$
In [16]:
t = 10
X = dot(U,Y0*cos(Omega*t))
print X
In [19]:
t = linspace(0.0,10.0,101)
allX = zeros( shape=(L,len(t)) )
for i in range(len(t)):
allX[:,i] = dot(U,Y0*cos(Omega*t[i]))
In [20]:
from pylab import *
%matplotlib inline
for x in range(L):
plot(t, allX[x,:] )
xlabel('temps')
ylabel('Position')
title("Dynamique d'oscillateurs couples")
grid(True)
show()
In [ ]: