Livre_silo 30 août 2013 16:32 Page 235
¨
©
¨
©
¨
©
¨
©
C o p y r i g h t E y r o l l e s
235
9 – Résolution numérique d’équations différentielles
Comme on voulait extraire la première composante des tableaux rendus (qui représentent θ), on a écrit :
pl.plot(t, array(theta)[:, 0])
La transformation de theta en tableau array était nécessaire : le slicing ne fonctionne pas
sur les tableaux de tableaux natifs de Python.
On peut, sur ce même exemple, souhaiter représenter le portrait de phase, c’est-à-dire l’ensemble des points des couples (y(t), y
′ (t)) pour t ∈ [0, 10]. C’est quelque chose de facile à
faire, grâce à un double-slicing. Le programme 9 ci-après regroupe les commandes qui ont
permis de réaliser la figure 9.9. On voit en particulier une option qui change l’épaisseur des
traits.
PROGRAMME 9 Un portrait de phase
def fpendule(t, X):
theta, thetap = X
return array([thetap, -theta])
_,theta = euler_vectoriel(fpendule, 0, 30, array([0,1]), 0.1)
pl.plot(array(theta)[:, 0], array(theta)[:, 1])
_,theta = heun(fpendule, 0, 30, array([0,1]), 0.1)
pl.plot(array(theta)[:, 0], array(theta)[:, 1], linewidth = 4)
portrait-pendule.pdf
Figure 9.9
Le portrait de phase du pendule linéarisé avec un pas grossier (h = 0.1) : Euler vs Heun