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

samedi 16 mars 2013

développements limités, graphiques avec sympy


Des développements limités avec sympy (le code python est ici).

Les fonctions rencontrées dans la suite:
series
limit
subs (plus un peu de programmation objet)
factor
solve
plot

1. Calculs

Allons-y (dans Python(x,y) par exemple)

from sympy import *
x, y, z, t = symbols('x y z t')

# la fonction à utiliser pour calculer un développement limité en sympy
# est series:
series(exp(x),x,0,10)

# comme sympy utilise $O(x^n)$ au lieu de $o(x^n)$ dans les développements
# limités, ajouter 1 à l'ordre usuel
# pour un DL de $e^x$ en 0 à l'ordre 3:
series(exp(x),x,0,4)

# mais python est un langage de programmation, on peut écrire une
# fonction convenable:
def DL(f,v,a,n):
   return (series(f,v,a,n+1))

# quelques exemples
series((sin(x)-x)/x**3, x, 0, 3)
DL((sin(x)-x)/x**3, x, 0, 2)
# on en déduit la limite:
limit((sin(x)-x)/x**3, x, 0)

# DL de $\sqrt{1+x^3}$ à l'ordre 6 en 0
DL(sqrt(1+x**3), x, 0, 6)

# DL de $\tan$ en 0 à l'ordre 5
dl=DL(tan(x), x, 0, 5)
print(dl)

Transformons le développement limité en polynôme en éliminant le $O(x^6)$.
Pour cela on utilise la programmation objet de python: l'objet dl possède une fonction propre à sa classe (comme une classe sociale ou une classe de lycée) qui s'appelle subs: on l'utilise sans donner l'objet dl en paramètre puisqu'elle le connaît déjà. Ce genre de fonction de classe s'appelle méthode.

p=dl.subs(O(x**6),0)
print(p)

# factorisons p
p1=factor(p)
print(p1)

#cherchons les racines de p1
sols=solve(p1,x)
print('polynome de degre '+str(p1.as_poly().degree())
      + ', '+str(len(sols))+ ' racines trouvees',
      sols)

# un DL nul en 0+
DL(exp(-1/x), x, 0, 3)
# en effet la limite suivante est nulle:
limit(exp(-1/x), x, 0, dir="+")

2. Graphiques

On va utiliser le module de tracé de courbes de sympy (documentation ici), qui permet de tracer des fonctions données par des expressions formelles, ce que matplotlib ne sait pas faire:


from sympy import *
x, y, z, t = symbols('x y z t')

La fonction de tracé de sympy s'appelle simplement plot,
elle rend en argument la fonction, puis un triplet qui indique la variable de la fonction et ses bornes de variation:

plot(sin(x), (x,-4*pi,4*pi), title='$\$$\sin(x)$\$$')

donne

on peut faire des tracés multiples, il suffit de mettre à la suite les fonctions. Ici $\sin$ et son développement limité à l'ordre 3:

plot(sin(x), series(sin(x),x,0,4).subs(O(x**4),0), (x,-pi,pi),
     title='$\$$\sin(x)$\$$ et son DL en 0 a l\'ordre 3')

donne

même chose à l'ordre 5:

plot(sin(x), series(sin(x),x,0,6).subs(O(x**6),0), (x,-pi,pi),
     title='$\$$\sin(x)$\$$ et son DL en 0 a l\'ordre 5')

donne

enfin avec plus de DL, où on voit bien l'approximation:

plot(sin(x),
     *[series(sin(x),x,0,2*n).subs(O(x**(2*n)),0) for n in
     range(1,5)]
     +[(x,-pi,pi)],
     title='$\$$\sin(x)$\$$ et ses developpements limites en 0
     aux ordres $1,3,5,7,9$')

donne

Je ne sais pas encore mettre des couleurs différentes aux courbes...

3 commentaires:

  1. Pour les couleurs sur tes courbes. Un quick hack :

    p = plot(sin(x),
    *[series(sin(x),x,0,2*n).subs(O(x**(2*n)),0) for n in range(1,4)]
    +[(x,-pi,pi),],
    show=False)
    p[0].line_color=(1.0,0.0,0.0)
    for n in range(1,4):
    p[n].line_color=(0.5-(n-1)/6.0,1.0-(n-1)/3.0,(n-1)/6.0+0.5)
    p.show()

    Mais si tu veux faire comme Maple voila :
    red = (1.0,0.0,0.0)
    green = (0.0,1.0,0.0)
    blue = (0.0,0.0,1.0)
    yellow = (1.0, 1.0, 0.0)
    magenta = (1.0, 0.0, 1.0)
    cyan = (0.0, 1.0, 1.0)
    base_colors = [ red, green, blue, yellow, magenta, cyan,
    (1.0, 0.5, 0.5), (0.5, 1.0, 0.5), (0.5, 0.5, 1.0),
    (1.0, 0.5, 0.0), (1.0, 0.0, 0.5), (0.5, 1.0, 0.0),
    (0.0, 1.0, 0.5), (0.5, 0.0, 1.0), (0.0, 0.5, 1.0) ]


    plage = (x,-pi,pi)
    p = plot(sin(x), plage, show=False)
    for n in range(1,4):
    p.extend( plot(series(sin(x),x,0,2*n).subs(O(x**(2*n)),0),
    plage, show=False) )
    for i, courbe in enumerate(p):
    courbe.line_color = base_colors[i]
    p.show()

    bien sur en commentaire c'est moche, alors voila sur pastebin : http://pastebin.com/WCNYg625

    RépondreSupprimer
    Réponses
    1. Sinon, il me reste le problème des accents qui ne passent pas de Python(x,y) à plot (en fait il y a le même problème avec matplotlib). C'est un problème qui risque de nécessiter des heures d'investigations ingrates, alors je retarde! De toutes façons ce n'est pas très grave.

      Supprimer