Marc BUFFAT

Marc BUFFAT

Professeur au département de Mécanique, Lyon 1 e-mail

Blog scientifique et pédagogique utilisant des notebooks IPython et Linux

Cours Python scientifique

Table des matières

In [1]:
%matplotlib inline
%autosave 300
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
from matplotlib import rcParams
rcParams['font.family'] = 'serif'
rcParams['font.size'] = 14
from IPython.core.display import HTML
from IPython.display import display,Image
from matplotlib import animation
#from JSAnimation import IPython_display
css_file = 'style.css'
HTML(open(css_file, "r").read())

Autosaving every 300 seconds
Out[1]:

Calcul symbolique avec sympy

Marc BUFFAT, dpt mécanique, Université Lyon 1

Licence Creative Commons
Mise à disposition selon les termes de la Licence Creative Commons
Attribution - Pas d’Utilisation Commerciale - Partage dans les Mêmes Conditions 2.0 France.

Systèmes de calcul symbolique (CAS)

Sympy

utilisation avec une sortie formattée LaTeX

 pprint(expr)      affiche expr
 str(expr)
 latex(expr)

manipulation avec des méthodes

 expr.simplify()
 expr.factor()
 expr.expand()
In [2]:
import sympy as sy
sy.init_printing()

Variables symboliques

In [3]:
# symbole (variable)
x= sy.Symbol('x')
(sy.pi+x)**2
Out[3]:
$$\left(x + \pi\right)^{2}$$
In [4]:
# fonction
f=sy.Function('f')
f(x+2)
Out[4]:
$$f{\left (x + 2 \right )}$$

On peut définir au début du programme une liste de variables symboliques avec leurs types et propriétés (entiere, positive, ..)

In [5]:
x,y,z,t=sy.symbols('x,y,z,t')
i,j,k=sy.symbols("i,j,k",integer=True)
f,g,h=sy.symbols("f,g,h",cls=sy.Function)
In [6]:
yp=sy.Symbol('x',positive=True)
yp<0
Out[6]:
$$\mathrm{False}$$

Liste de variables courantes

In [7]:
import sympy.abc
dir(sympy.abc)
Out[7]:
['A',
 'B',
 'C',
 'D',
 'E',
 'F',
 'G',
 'H',
 'I',
 'J',
 'K',
 'L',
 'M',
 'N',
 'O',
 'P',
 'Q',
 'R',
 'S',
 'T',
 'U',
 'V',
 'W',
 'X',
 'Y',
 'Z',
 '__builtins__',
 '__cached__',
 '__doc__',
 '__file__',
 '__loader__',
 '__name__',
 '__package__',
 '__spec__',
 '_clash',
 '_clash1',
 '_clash2',
 'a',
 'alpha',
 'b',
 'beta',
 'c',
 'chi',
 'd',
 'delta',
 'division',
 'e',
 'epsilon',
 'eta',
 'exec_',
 'f',
 'g',
 'gamma',
 'greeks',
 'h',
 'i',
 'iota',
 'j',
 'k',
 'kappa',
 'l',
 'lamda',
 'm',
 'mu',
 'n',
 'nu',
 'o',
 'omega',
 'omicron',
 'p',
 'phi',
 'pi',
 'print_function',
 'psi',
 'q',
 'r',
 'rho',
 's',
 'sigma',
 'string',
 'symbols',
 't',
 'tau',
 'theta',
 'u',
 'upsilon',
 'v',
 'w',
 'x',
 'xi',
 'y',
 'z',
 'zeta']
In [8]:
from sympy.abc import w,omega
print(type(omega))

<class 'sympy.core.symbol.Symbol'>

nombre complexe

$$ I = \sqrt{-1}$$
In [9]:
from sympy import I
1+1*I
(x+I)**2
Out[9]:
$$\left(x + i\right)^{2}$$

nombre rationnel (fraction)

In [10]:
r1=sy.Rational(3,6)
r2=sy.Rational(4,16)
r1+r2
Out[10]:
$$\frac{3}{4}$$

Evaluation

exp.evalf(n=20)
N(exp,n=20)
In [11]:
from sympy import pi
print(pi.evalf(100))
print(sy.N(pi,100))

3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068
3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068

substitution

exp.subs(x,val)
exp.subs([(x1,val1),(x2,val2),..])
In [12]:
y=(x+pi)**2
y.subs(x,z+pi)
Out[12]:
$$\left(z + 2 \pi\right)^{2}$$

évaluation numérique

In [13]:
X=np.linspace(-10,10,51)
Y=np.array([sy.N(y.subs(x,val)) for val in X])
plt.plot(X,Y)
Out[13]:
[<matplotlib.lines.Line2D at 0x7f3f2a65dd68>]

Fonction numérique

lambdify(var,exp,'numpy')

transforme exp en fonction de var

In [14]:
F = sy.lambdify([x],y,'numpy')
Y=F(X)
plt.plot(X,Y)
Out[14]:
[<matplotlib.lines.Line2D at 0x7f3f2a5587f0>]

comparaison de l’éfficacité!!!

utiliser lambdify pour l'évaluation numérique des fonctions

In [15]:
%%timeit
Y=np.array([sy.N(y.subs(x,val)) for val in X])

16.3 ms ± 58.3 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
In [16]:
%%timeit
Y=F(X)

7.78 µs ± 48.6 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

Manipulation algébrique

  • equalité

    Eq(y,x)     equation y=x
y=x         affectation
y==x        test egalité
In [17]:
eq=sy.Eq(y,x**2)
print(eq)
print(sy.latex(eq))
eq

Eq((x + pi)**2, x**2)
\left(x + \pi\right)^{2} = x^{2}
Out[17]:
$$\left(x + \pi\right)^{2} = x^{2}$$
In [18]:
yy=x**2
print(yy)
yy

x**2
Out[18]:
$$x^{2}$$
In [19]:
yy == x**2
Out[19]:
True
  • manipulation d’expression 
    factorisation
expansion
simplify
collect
In [20]:
expr=(1+x)**3;
a=expr.expand()
a
Out[20]:
$$x^{3} + 3 x^{2} + 3 x + 1$$
In [21]:
a.factor()
Out[21]:
$$\left(x + 1\right)^{3}$$
In [22]:
a=(x + x**2)/(x*sy.sin(y)**2 + x*sy.cos(y)**2)
sy.pprint(a)
a.simplify()

                2                  
               x  + x              
───────────────────────────────────
     2⎛       2⎞        2⎛       2⎞
x⋅sin ⎝(x + π) ⎠ + x⋅cos ⎝(x + π) ⎠
Out[22]:
$$x + 1$$
In [23]:
sy.collect(i*x**2 + j*x**2 + i*x - j*x + k, x)
Out[23]:
$$k + x^{2} \left(i + j\right) + x \left(i - j\right)$$
  • dérivation 
       exp.diff(x)    dériv. / x
   exp.diff(x,x)  dériv. 2nd
   exp.diff(x,n)  dériv. ordre n
In [24]:
f(x).diff(x,3)
Out[24]:
$$\frac{d^{3}}{d x^{3}} f{\left (x \right )}$$
  • intégration 
      integrate(exp,x)
    intégrale indéfinie 
      integrate(exp,(x,a,b))
    intégrale définie entre a et b
In [25]:
sy.Integral(sy.sin(x**2),x)
Out[25]:
$$\int \sin{\left (x^{2} \right )}\, dx$$
In [26]:
sy.integrate(sy.sin(x**2),x)
Out[26]:
$$\frac{3 \sqrt{2} \sqrt{\pi} S\left(\frac{\sqrt{2} x}{\sqrt{\pi}}\right)}{8 \Gamma{\left(\frac{7}{4} \right)}} \Gamma{\left(\frac{3}{4} \right)}$$
  • limite 
        limit(exp,x,val)
    limite pour x=val (oo = infinie)
In [27]:
sy.limit(sy.sin(x)/x,x,0)
Out[27]:
$$1$$
  • series de Taylor 
      exp.series(x,val,n)
    développement en serie d’ordre n 
      removeO()
    enleve ordre (polynome)
In [28]:
expr=sy.sin(x)/x
expr.series(x,0,6)
Out[28]:
$$1 - \frac{x^{2}}{6} + \frac{x^{4}}{120} + \mathcal{O}\left(x^{6}\right)$$
In [29]:
expr.series(x,0,6).removeO()
Out[29]:
$$\frac{x^{4}}{120} - \frac{x^{2}}{6} + 1$$
  • approximation par différences finies (version 0.7.6) 
     as_finite_diff(expr)
    calcul différences finies 
     as_finite_diff(expr,[pts])
In [30]:
h=sy.Symbol('h')
sy.Derivative.as_finite_difference(f(x).diff(x,2),h)
Out[30]:
$$- \frac{2}{h^{2}} f{\left (x \right )} + \frac{1}{h^{2}} f{\left (- h + x \right )} + \frac{1}{h^{2}} f{\left (h + x \right )}$$
In [31]:
sy.Derivative.as_finite_difference(f(x).diff(x),h)
Out[31]:
$$- \frac{1}{h} f{\left (- \frac{h}{2} + x \right )} + \frac{1}{h} f{\left (\frac{h}{2} + x \right )}$$
In [32]:
sy.Derivative.as_finite_difference(f(x).diff(x,4),h)
Out[32]:
$$\frac{6}{h^{4}} f{\left (x \right )} + \frac{1}{h^{4}} f{\left (- 2 h + x \right )} - \frac{4}{h^{4}} f{\left (- h + x \right )} - \frac{4}{h^{4}} f{\left (h + x \right )} + \frac{1}{h^{4}} f{\left (2 h + x \right )}$$
  • résolution d’equations algébriques 
    solve(eq,var)
    résoud l’équation eq par rapport à var 
    solve([eq],[var])
    résoud un système d’équations
In [33]:
sy.solve(x**2-y,x)
Out[33]:
$$\left [ - \frac{\pi}{2}\right ]$$
In [34]:
y=sy.Rational(1,4)
sy.solve(sy.cos(x)-y,x)
Out[34]:
$$\left [ - \operatorname{acos}{\left (\frac{1}{4} \right )} + 2 \pi, \quad \operatorname{acos}{\left (\frac{1}{4} \right )}\right ]$$

  • Equations de Lagrange

     euler_equations(L,fonct.,var.)
In [35]:
m,g=sy.symbols('m,g')
L=m*sy.diff(x(t),t)**2/2 - g*x(t)
sy.euler_equations(L,x(t),[t])
Out[35]:
$$\left [ - g - m \frac{d^{2}}{d t^{2}} x{\left (t \right )} = 0\right ]$$

Formulation de Lagrange pour le pendule double

In [36]:
m,g,l=sy.symbols('m,g,l')
theta1,theta2=sy.symbols('theta1,theta2')
dtheta1=theta1(t).diff(t)
dtheta2=theta2(t).diff(t)
Ec=m*l**2/6*(4*dtheta1**2+dtheta2**2+3*dtheta1*dtheta2*sy.cos(theta1(t)-theta2(t)))
print("énergie cinétique E:")
Ec

énergie cinétique E:
Out[36]:
$$\frac{l^{2} m}{6} \left(3 \cos{\left (\theta_{1}{\left (t \right )} - \theta_{2}{\left (t \right )} \right )} \frac{d}{d t} \theta_{1}{\left (t \right )} \frac{d}{d t} \theta_{2}{\left (t \right )} + 4 \left(\frac{d}{d t} \theta_{1}{\left (t \right )}\right)^{2} + \left(\frac{d}{d t} \theta_{2}{\left (t \right )}\right)^{2}\right)$$
In [37]:
Up=-m*g*l/2*(3*sy.cos(theta1(t))+sy.cos(theta2(t)))
print("énergie potentielle U:")
Up

énergie potentielle U:
Out[37]:
$$- \frac{g l}{2} m \left(3 \cos{\left (\theta_{1}{\left (t \right )} \right )} + \cos{\left (\theta_{2}{\left (t \right )} \right )}\right)$$
In [38]:
L=Ec-Up
print("Lagrangien L:")
L

Lagrangien L:
Out[38]:
$$\frac{g l}{2} m \left(3 \cos{\left (\theta_{1}{\left (t \right )} \right )} + \cos{\left (\theta_{2}{\left (t \right )} \right )}\right) + \frac{l^{2} m}{6} \left(3 \cos{\left (\theta_{1}{\left (t \right )} - \theta_{2}{\left (t \right )} \right )} \frac{d}{d t} \theta_{1}{\left (t \right )} \frac{d}{d t} \theta_{2}{\left (t \right )} + 4 \left(\frac{d}{d t} \theta_{1}{\left (t \right )}\right)^{2} + \left(\frac{d}{d t} \theta_{2}{\left (t \right )}\right)^{2}\right)$$
In [39]:
print("Equations de lagrange")
eq=sy.euler_equations(L,{theta1(t),theta2(t)},[t])
eq

Equations de lagrange
Out[39]:
$$\left [ - \frac{g l}{2} m \sin{\left (\theta_{2}{\left (t \right )} \right )} - \frac{l^{2} m}{6} \left(- 3 \left(\frac{d}{d t} \theta_{1}{\left (t \right )} - \frac{d}{d t} \theta_{2}{\left (t \right )}\right) \sin{\left (\theta_{1}{\left (t \right )} - \theta_{2}{\left (t \right )} \right )} \frac{d}{d t} \theta_{1}{\left (t \right )} + 3 \cos{\left (\theta_{1}{\left (t \right )} - \theta_{2}{\left (t \right )} \right )} \frac{d^{2}}{d t^{2}} \theta_{1}{\left (t \right )} + 2 \frac{d^{2}}{d t^{2}} \theta_{2}{\left (t \right )}\right) + \frac{l^{2} m}{2} \sin{\left (\theta_{1}{\left (t \right )} - \theta_{2}{\left (t \right )} \right )} \frac{d}{d t} \theta_{1}{\left (t \right )} \frac{d}{d t} \theta_{2}{\left (t \right )} = 0, \quad - \frac{3 g}{2} l m \sin{\left (\theta_{1}{\left (t \right )} \right )} - \frac{l^{2} m}{6} \left(- 3 \left(\frac{d}{d t} \theta_{1}{\left (t \right )} - \frac{d}{d t} \theta_{2}{\left (t \right )}\right) \sin{\left (\theta_{1}{\left (t \right )} - \theta_{2}{\left (t \right )} \right )} \frac{d}{d t} \theta_{2}{\left (t \right )} + 3 \cos{\left (\theta_{1}{\left (t \right )} - \theta_{2}{\left (t \right )} \right )} \frac{d^{2}}{d t^{2}} \theta_{2}{\left (t \right )} + 8 \frac{d^{2}}{d t^{2}} \theta_{1}{\left (t \right )}\right) - \frac{l^{2} m}{2} \sin{\left (\theta_{1}{\left (t \right )} - \theta_{2}{\left (t \right )} \right )} \frac{d}{d t} \theta_{1}{\left (t \right )} \frac{d}{d t} \theta_{2}{\left (t \right )} = 0\right ]$$
In [40]:
print("forme numérique ")
print(eq)

forme numérique 
[Eq(-g*l*m*sin(theta2(t))/2 - l**2*m*(-3*(Derivative(theta1(t), t) - Derivative(theta2(t), t))*sin(theta1(t) - theta2(t))*Derivative(theta1(t), t) + 3*cos(theta1(t) - theta2(t))*Derivative(theta1(t), t, t) + 2*Derivative(theta2(t), t, t))/6 + l**2*m*sin(theta1(t) - theta2(t))*Derivative(theta1(t), t)*Derivative(theta2(t), t)/2, 0), Eq(-3*g*l*m*sin(theta1(t))/2 - l**2*m*(-3*(Derivative(theta1(t), t) - Derivative(theta2(t), t))*sin(theta1(t) - theta2(t))*Derivative(theta2(t), t) + 3*cos(theta1(t) - theta2(t))*Derivative(theta2(t), t, t) + 8*Derivative(theta1(t), t, t))/6 - l**2*m*sin(theta1(t) - theta2(t))*Derivative(theta1(t), t)*Derivative(theta2(t), t)/2, 0)]

Résolution d’EDO

soit l’EDO $$ \frac{d^2y}{d x^2} - \omega^2 y = x^2 $$ avec les conditions $$y(0)=1, \frac{d y}{dx}(1)=0 $$

tracer de la solution analytique pour $\omega=2\pi$

In [41]:
x,y,C1,C2 = sy.symbols('x,y,C1,C2')
omega = sy.Symbol('omega',positive=True)
eq=y(x).diff(x,x)+omega**2*y(x)-x**2
eq
Out[41]:
$$\omega^{2} y{\left (x \right )} - x^{2} + \frac{d^{2}}{d x^{2}} y{\left (x \right )}$$
In [42]:
sol=sy.dsolve(eq,y(x))
sol
Out[42]:
$$y{\left (x \right )} = C_{1} \sin{\left (\omega x \right )} + C_{2} \cos{\left (\omega x \right )} + \frac{x^{2}}{\omega^{2}} - \frac{2}{\omega^{4}}$$
In [43]:
print("extraction terme de gauche et droite")
sol.lhs, sol.rhs

extraction terme de gauche et droite
Out[43]:
$$\left ( y{\left (x \right )}, \quad C_{1} \sin{\left (\omega x \right )} + C_{2} \cos{\left (\omega x \right )} + \frac{x^{2}}{\omega^{2}} - \frac{2}{\omega^{4}}\right )$$
In [44]:
# conditions aux limites en 0 et 1
ys=sol.rhs
eq1=ys.subs(x,0)-1
eq2=sy.diff(ys,x).subs(x,1)
eq1,eq2
Out[44]:
$$\left ( C_{2} - 1 - \frac{2}{\omega^{4}}, \quad C_{1} \omega \cos{\left (\omega \right )} - C_{2} \omega \sin{\left (\omega \right )} + \frac{2}{\omega^{2}}\right )$$
In [45]:
# valeurs des coefficients
sol=sy.solve({eq1,eq2},{C1,C2})
sol
Out[45]:
$$\left \{ C_{1} : \frac{1}{\omega^{4} \cos{\left (\omega \right )}} \left(- 2 \omega + \left(\omega^{4} + 2\right) \sin{\left (\omega \right )}\right), \quad C_{2} : 1 + \frac{2}{\omega^{4}}\right \}$$
In [46]:
# d'où la solution
ysol=ys.subs([(C1,sol[C1]),(C2,sol[C2])])
ysol
Out[46]:
$$\left(1 + \frac{2}{\omega^{4}}\right) \cos{\left (\omega x \right )} + \frac{x^{2}}{\omega^{2}} + \frac{\sin{\left (\omega x \right )}}{\omega^{4} \cos{\left (\omega \right )}} \left(- 2 \omega + \left(\omega^{4} + 2\right) \sin{\left (\omega \right )}\right) - \frac{2}{\omega^{4}}$$
In [47]:
# conversion en fonction numérique
Ys=sy.lambdify(x,ysol.subs(omega,2*pi),'numpy')
X=np.linspace(0,1,50)
Y=Ys(X)
plt.plot(X,Y)
Out[47]:
[<matplotlib.lines.Line2D at 0x7f3f2a1d4fd0>]

Algébre linéaire

Matrix matrice ou vecteur
résolution (solve)
inversion (inv)
décomposition QR
In [48]:
M=sy.Matrix([[1,x],[x,2]])
In [49]:
M.inv()
Out[49]:
$$\left[\begin{matrix}\frac{2}{- x^{2} + 2} & - \frac{x}{- x^{2} + 2}\- \frac{x}{- x^{2} + 2} & \frac{1}{- x^{2} + 2}\end{matrix}\right]$$
In [50]:
b=sy.Matrix([1,2])
M.solve(b)
Out[50]:
$$\left[\begin{matrix}- \frac{2 x}{- x^{2} + 2} + \frac{2}{- x^{2} + 2}\- \frac{x}{- x^{2} + 2} + \frac{2}{- x^{2} + 2}\end{matrix}\right]$$
In [51]:
A=sy.Matrix([[1,1,3],[4,5,6],[7,8,9]])
Q,R = A.QRdecomposition()
sy.pprint(Q)
Q*R

⎡ √66   -√6   3⋅√11⎤
⎢ ───   ────  ─────⎥
⎢  66    6      11 ⎥
⎢                  ⎥
⎢2⋅√66   √6    √11 ⎥
⎢─────   ──    ─── ⎥
⎢  33    3      11 ⎥
⎢                  ⎥
⎢7⋅√66  -√6   -√11 ⎥
⎢─────  ────  ─────⎥
⎣  66    6      11 ⎦
Out[51]:
$$\left[\begin{matrix}1 & 1 & 3\4 & 5 & 6\7 & 8 & 9\end{matrix}\right]$$

Pattern matching

calcul de paramêtres par correspondance de modèle

In [52]:
expr=(1+5*x+x**2)**2
expr
Out[52]:
$$\left(x^{2} + 5 x + 1\right)^{2}$$
In [53]:
p = sy.Wild('p', exclude=[x])
expr.match((p*x+1+x**2)**2)
Out[53]:
$$\left \{ p : 5\right \}$$

Documentation

Bibliothéques utilisées

In [54]:
print("\t\tSystème utilisé")
import sys
print("Système :\t\t",sys.platform)
import platform
print(platform.platform())
print("Ordinateur:\t\t",platform.machine())
print("Version de Python:\t",sys.version)
import IPython
print("Version de IPython:\t",IPython.__version__)
import numpy
print("Version de numpy:\t",numpy.version.version)
import scipy
print("Version de scipy:\t",scipy.version.version)
import matplotlib
print("Version de matplotlib:\t",matplotlib.__version__)
import sympy
print("Version de sympy:\t",sympy.__version__)

		Système utilisé
Système :		 linux
Linux-4.4.0-112-generic-x86_64-with-Ubuntu-16.04-xenial
Ordinateur:		 x86_64
Version de Python:	 3.5.2 (default, Nov 23 2017, 16:37:01) 
[GCC 5.4.0 20160609]
Version de IPython:	 6.2.1
Version de numpy:	 1.11.0
Version de scipy:	 1.0.0
Version de matplotlib:	 1.5.1
Version de sympy:	 1.1.1