B. Boutin - benjamin.boutin@univ-rennes1.fr

Introduction à Python

  • Installer au préalable une distribution et un environnement de développement pour Python, par exemple Pyzo.
  • Possibilité d'utiliser un environnement python en ligne (avec notebook Jupyter) : https://jupyter.org/try

Rudiments

Variables et types

  • entiers int
  • flottants float
  • booléens bool
  • chaînes de caractères string
  • tuples tuple (non-modifiables)
  • listes list

Fonctions

Usage compact:

g = lambda x: x*2

Usage standard:

def nom_de_fonction(arg1,arg2)
    """Description de la fonction""""
    INSTRUCTIONS
    return out1,out2,out3

Usage avancé: packing d'un nombre arbitraire de paramètres dans un tuple

def product(*args):
    """Produit d'un nombre arbitraire d'arguments"""
    y = 1.
    for x in args[:]:
        y = y*x
    return y

product(1,2,3,4)
    24.0

Structures de contrôles : for, if, while

  • boucle for
    for element in liste:
          INSTRUCTIONS
    
  • test if (elif) else
    if CONDITION1:
          INSTRUCTIONS1
    elif CONDITION2:
          INSTRUCTIONS2
    elif CONDITION3:
          INSTRUCTIONS3
    else:
          INSTRUCTIONS4
    
  • boucle while
    while CONDITION:
          INSTRUCTIONS
    

Quelques exemples de commandes pour la manipulation des listes

In [1]:
L=[]
L.append(2)
L.append(3)
print("L =",L)
M=[1,2,3,4]
print(M[2:])
N=[i**2 for i in range(11)]
print("carrés =",N)
L = [2, 3]
[3, 4]
carrés = [0, 1, 4, 9, 16, 25, 36, 49, 64, 81, 100]

Exercice : définition de fonction

Définir une fonction qui à $x\in\mathbb{R}$ associe $x^2+1$ si $x\geq 0$ et $\sqrt{-x}$ sinon.

In [2]:
from math import sqrt

def h(x):
    if x>=0:
        return x**2+1
    else:
        return sqrt(-x)

print(h(1))
print(h(1.))
print(h(-1))
2
2.0
1.0

Exercice : calcul des $n$ premiers termes de la suite de Fibonacci

Programme une fonction fibonacci prenant en argument l'entier $n$ et renvoyant la liste des termes $u_0,\ldots,u_n$ de la suite de Fibonacci:

$$ u_0=1,\ u_1=1,\\ \forall n\in\mathbb{N},\ u_{n+2}=u_{n+1}+u_n. $$

In [3]:
# %load Fibonacci.py
def fibonacci(n):
    u = []
    a,b = 1,1
    u.append(a)
    if n > 0:
        u.append(b)    
        for k in range(n-1):
            a,b = b,a+b
            u.append(b)
    return u

fibonacci(20)
Out[3]:
[1,
 1,
 2,
 3,
 5,
 8,
 13,
 21,
 34,
 55,
 89,
 144,
 233,
 377,
 610,
 987,
 1597,
 2584,
 4181,
 6765,
 10946]

Exercice : calcul du pgcd de plusieurs nombres

Programmer la fonction suivante, qui calcule le pgcd de deux entiers $a$ et $b$.

def euclide(a, b):
    """Calcul de pgcd de a et b via l'algorithme d'Euclide"""
    if a < b:
         b,a = a,b
    while b != 0:
         a,b = b,a%b
    return a

euclide(5**3*7**2,9*7*5)

Étant donnés trois entiers $a,b,c$, on peut obtenir leur pgcd en utilisant la propriété suivante:

$${\rm pgcd}(a,b,c)={\rm pgcd}({\rm pgcd}(a,b),c).$$

On peut ainsi calculer le pgcd d'un nombre arbitraire d'entiers par exemple en utilisant le principe de packing présenté précédemment.

In [4]:
# %load pgcdn.py
def euclide(a,b):
    """Calcul du PGCD par l'algorithme d'Euclide"""
    if a<b:
        b,a = a,b
    while b!=0:
        a,b = b,a%b
    return a


def pgcdn(*n):
    """Calcul du pgcd de plusieurs entiers"""
    p = euclide(n[0], n[1])
    for x in n[2:]:
        p = euclide(p, x)
    return p

print(pgcdn(123,456,789))
3

Exercice : recherche d'un zéro

Par une méthode de dichotomie, déterminer un zéro de la fonction $\cos x - x$ à $10^{-9}$ près.

Programmer la méthode de Newton pour la recherche de ce même zéro.

In [5]:
from math import *

def f(x):
    return(cos(x)-x)

# Mise en place de l'algorithme de dichotomie sur l'exemple
a=0
b=1
while (b-a)>0.001:
    c=(a+b)/2
    if f(a)*f(c)<0:
        a,b = a,c
    else:
        a,b = c,b
print(a,b)

# Encapsulation dans une fonction à paramètres
def dichotomie(a,b,f,epsilon):
    while (b-a)>epsilon:
        c=(a+b)/2
        if f(a)*f(c)<0:
            a,b = a,c
        else:
            a,b = c,b
    return(a,b)

print(dichotomie(0,1,f,10**-16))

def fp(x):
    return(-sin(x)-1)

# Mise en place de l'algorithme de Newton
x=0.73
for i in range(5):
    x = x-f(x)/fp(x)
print(x)
0.73828125 0.7392578125
(0.7390851332151609, 0.7390851332151609)
0.7390851332151607

Utilisation des bibliothèques numpyet matplotlib

Documentation pour les bilbiothèques utilisées:

On pourra parcourir en particulier le tutorial pyplot de matplotlib

Recherche de zéro (suite)

Reprendre l'exercice précédent en utilisant les fonctions scipy.optimize.bisect et scipy.optimize.newton.

In [6]:
from scipy.optimize import *

def f(x):
    return(cos(x)-x)

X1 = bisect(f,0,1)
X2 = newton(f,0)
X3 = root(f,0)

print(X1)
print(X2)
print(X3)
print(X3.x[0])
0.7390851332147577
0.7390851332151607
    fjac: array([[-1.]])
     fun: array([0.])
 message: 'The solution converged.'
    nfev: 9
     qtf: array([-2.66786593e-13])
       r: array([1.67361202])
  status: 1
 success: True
       x: array([0.73908513])
0.7390851332151607

Exercice

Tracer la courbe cardioïde paramétrée suivante

$$ \left\{\begin{aligned} x(t) &= - \cos(t) + \cos(t)^2\\ y(t) &= - \sin(t) + \cos(t)\sin(t) \end{aligned}\right. $$

et sa tangente au point de paramètre $t=\pi/2$, puis l'ensemble des tangentes aux points de paramètres $s=2k\pi/10$ pour $k=0,\ldots,9$.

In [8]:
from numpy import *
import matplotlib.pyplot as plt

t = arange(0, 2*pi, 0.05)
x = -cos(t)+cos(t)**2
y = -sin(t)+cos(t)*sin(t)
plt.axis('equal')
plt.plot(x, y)
plt.show()
In [9]:
# %load Cardioide.py
from numpy import *
import matplotlib.pyplot as plt

t = linspace(0,2*pi,100)
x = - cos(t) + cos(t)**2
y = - sin(t) + cos(t)*sin(t)
plt.plot(x,y)

n = 12
p = linspace(-1,1,100)
for s in linspace(0.,2.*pi,n+1):
    xs = - cos(s) + cos(s)**2
    ys = - sin(s) + cos(s)*sin(s)
    dxs = sin(s) - 2*cos(s)*sin(s)
    dys = - cos(s) - sin(s)*sin(s) + cos(s)*cos(s)
    if sqrt(dxs**2+dys**2)!=0:
        Lx = [xs+k*dxs/sqrt(dxs**2+dys**2) for k in p]
        Ly = [ys+k*dys/sqrt(dxs**2+dys**2) for k in p]
        plt.plot(Lx,Ly,'r')
        plt.plot(xs,ys,'r',marker='.')

plt.axis('equal')
plt.show()

Intégration numérique

Afin d'introduire quelques exemples élémentaires de syntaxe Python concernant les boucles et les définitions de fonctions, on souhaite calculer une approximation numérique de l'intégrale suivante:

$$I = \int _0^3 \sin(10 x + e^x) dx.$$

Pour ce faire, on emploie la méthode des trapèzes décrite ci-après. Étant donné un entier naturel non-nul $n$, on doit évaluer la quantité:

$$I_n = \dfrac{3}{n} \sum_{i=0}^{n-1} \dfrac{f(x_i)+f(x_{i+1})}{2},$$

ayant posé $x_i=3i/n$ pour $0\leq i \leq n$.

Compléter le code suivant pour mener ce calcul à bien :

def f(x):
    y = ...
    return y

def Trapeze(f,n):
    S = 0
    h = 3/n
    x,y = 0,h
    for i in range(0,n):
        S = S + ...
        x,y = y,y+h
    return S

print(Trapeze(f,100))
In [10]:
from numpy import *
import matplotlib.pyplot as plt

def f(x):
    y=sin(10*x+exp(x))
    return y

x = arange(0,3,0.005)
y = f(x)
plt.plot(x,y)
plt.title('La fonction à intégrer')
plt.show()
In [11]:
def Trapeze(f,n):
    """ Méthode des trapèzes pour intégrer f sur [0,3] avec n sous-intervalles """
    S = 0
    h = 3/n
    x,y = 0,h
    for i in range(0,n):
        S = S + h/2*(f(x)+f(y))
        x,y = y,y+h
    return S

print(Trapeze(f,100))
print(Trapeze(f,500))
print(Trapeze(f,1000))
0.01903430191679251
0.017301610805672312
0.017248340424185184

Étude de convergence

Estimer l'erreur commise dans l'intégration de la fonction précédente par la méthode des trapèzes pour différentes valeurs de l'entier $n$. En guise de valeur de référence pour l'intégrale $I$, on utilisera l'approximation obtenue par la commande intégrée quad de la bibliothèque scipy.integrate.

from scipy import integrate
I = integrate.quad(f,0,3)[0]
print(I)

On cherchera ensuite à illustrer graphiquement (par un tracé en échelle logarithmique) la majoration suivante:

$$ \exists C>0,\ \forall n\in\mathbb{N}^\star,\ |I_n-I|\leq \dfrac{C}{n^2}.$$

On pourra utiliser le code suivant qui met en œuvre le calcul de l'erreur $|I_n-I|$ pour différentes valeurs de $n$ et représente l'évolution de cette erreur en échelle logarithmique.

Err = []
ListN = np.arange(1,200)
for n in ListN:
    In = Trapeze(f,n)
    Err.append(np.abs(In-I))

plt.plot(np.log(ListN),np.log(Err))
plt.plot(np.log(ListN),-2*np.log(ListN))
plt.title('Courbe de convergence')
plt.xlabel('log n')
plt.ylabel('log erreur')
plt.show()
In [12]:
from numpy import *
from scipy import integrate
I = integrate.quad(f,0,3)[0]
print(I)

Err=[]
ListN = arange(1,200)
for n in ListN:
    In = Trapeze(f,n)
    Err.append(abs(In-I))

plt.clf
plt.plot(np.log(ListN),np.log(Err))
plt.plot(np.log(ListN),-2*np.log(ListN))
plt.title('Courbe de convergence - Méthode des trapèzes')
plt.xlabel('log n')
plt.ylabel('log erreur')
plt.show()
0.017230595206954503
In [13]:
plt.plot(ListN,Err,ListN,[n**(-2.) for n in ListN])
plt.title('Courbe de convergence - Méthode des trapèzes')
plt.xlabel('n')
plt.ylabel('erreur')
plt.xscale('log')
plt.yscale('log')
plt.show()

Méthode de Simpson

Reprendre le code précédent pour utiliser maintenant la méthode de quadrature de Simpson :

$$J_n = \dfrac{3}{n}\sum_{i=0}^{n-1} \dfrac{1}{6}\left(f(x_i)+f(\tfrac{x_{i}+x_{i+1}}{2}) +f(x_{i+1})\right).$$

Estimer numériquement l'ordre de convergence effectif, c'est-à-dire la plus grande constante $\alpha>0$ telle que

$$ \exists C>0,\ \forall n\in\mathbb{N}^\star,\ |J_n-I|\leq \dfrac{C}{n^\alpha}.$$

In [14]:
def Simpson(f,n):
    """ Méthode de Simpson pour intégrer f sur [0,3] avec n sous-intervalles """
    S = 0
    h = 3/n
    x,y = 0,h
    for i in range(0,n):
        S = S + h/6*(f(x)+4*f(0.5*(x+y))+f(y))
        x,y = y,y+h
    return S

ErrS=[]
ListN = arange(1,200)
for n in ListN:
    In = Simpson(f,n)
    ErrS.append(abs(In-I))
    
plt.plot(ListN,Err,ListN,ErrS,ListN,[n**(-2.) for n in ListN],ListN,[n**(-4.) for n in ListN])
plt.title('Courbe de convergence')
plt.xlabel('n'), plt.ylabel('erreur')
plt.xscale('log'), plt.yscale('log')
plt.legend(['Trapeze','Simpson','ordre 2', 'ordre 4'])

plt.show()

Équations différentielles

Afin de résoudre des équations différentielles ordinaires, on peut utiliser la fonction odeint de la bibliothèque Scipy.

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

Voici un exemple simple pour la résolution de l'équation logistique

$$y'(t) = ry(t)\left(1-\dfrac{y(t)}{K}\right).$$

Dans les simulations, on pourra par exemple utiliser les paramètres $r=1.5$ et $K=6$ ainsi que la donnée initiale $y_0=1$.

def f(y,t):
    return 1.5*y*(1.-y/6.)

y0 = 1.0
t = np.linspace(0,5,201)
sol = odeint(f,y0,t)

plt.plot(t,sol)
plt.axis('equal')
plt.show()
In [15]:
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

def f(y,t):
    return 1.5*y*(1.-y/6.)

y0 = 1.0
t = np.linspace(0,5,201)
sol = odeint(f,y0,t)

plt.plot(t,sol)
plt.plot(t[0],sol[0],'.b')
plt.axis('equal')
plt.show()
In [16]:
t = np.linspace(0,5,201)
for y0 in linspace(0,12,20):
    sol = odeint(f,y0,t)
    plt.plot(t,sol)
    plt.plot(t[0],sol[0],'.b')
plt.show()

Modèle proie-prédateur

Adapter l'exemple précédent de sorte à résoudre le modèle proie-prédateur de Lotka-Volterra:

$$ \left\{ \begin{aligned} x'(t) &= +x(t)(a-by(t))\\ y'(t) &= -y(t)(c-dx(t)) \end{aligned} \right. $$

On choisira pour les simulations les paramètres $a=2$, $b=1$, $c=1$ et $d=0.3$ et pour ce qui concerne la donnée initiale : $x(0)=1$, et $y(0)=1$. La solution sera évaluée sur l'intervalle de temps $[0,30]$.

À l'aide de la fonction quiver de matplotlib décrite dans la documentation, tracer le portrait de phase du système différentiel, ainsi que la trajectoire solution vue comme la courbe paramétrée $t\mapsto (x(t),y(t))$.

In [17]:
[a,b,c,d]=[2,1,1,0.3]
def f(y,t):
    return [y[0]*(a-b*y[1]), -y[1]*(c-d*y[0])]

y0 = [1.0,1.0]
t = np.linspace(0,30,1000)
sol = odeint(f,y0,t)

plt.subplot(211)
plt.plot(t,sol)
plt.axis('equal')
plt.show()

plt.subplot(212)
plt.plot(sol[:,0],sol[:,1])
plt.axis('equal')
plt.show()
In [18]:
X = np.linspace(0, 10, 15)
Y = np.linspace(0, 6, 20)
U, V = np.meshgrid(X, Y)

Q1 = 0*U
Q2 = 0*U

for j in range(len(X)):
    for i in range(len(Y)):
        Q = f((U[i][j],V[i][j]),0)
        Q1[i][j] = Q[0]
        Q2[i][j] = Q[1]
        
fig, ax = plt.subplots()
q = ax.quiver(X, Y, Q1, Q2)
plt.plot(sol[:,0],sol[:,1])
plt.axis('equal')
plt.show()

Attracteur de Lorenz

Résoudre numériquement quelques trajectoires pour le système dynamique de l'attracteur de Lorenz en dimension 3 donné par

$$ \left\{ \begin{aligned} x'(t) &= \sigma(y(t)-x(t))\\ y'(t) &= \rho x(t) - y(t) -x(t)z(t)\\ z'(t) &= x(t)y(t) - bz(t). \end{aligned} \right. $$

On retiendra les paramètres $\sigma=10$, $\rho=30$ et $b=2$. Les trajectoires pourront être calculées sur l'intervalle de temps $t\in[0,10]$ pour les données initiales $(20,0,0)$, $(10,0,0)$ et $(0.1,0,0)$.

Pour le tracé des trajectoires en dimension 3, la commande suivante permet un tracé de courbe paramétrée en dimension 3 d'espace. Il faut préalablement avoir stocké le vecteur sol dans un ndarray à 3 colonnes.

from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure()
ax = pl.axes(projectio='3d')
ax.plot3D(sol[:,0],sol[:,1],sol[:,2])
plt.show()
In [19]:
[s,r,b]=[10,30,2]

def Lorenz(y,t):
    return [s*(y[1]-y[0]), r*y[0]-y[1]-y[0]*y[2],y[0]*y[1]-b*y[2]]

y0 = [2,0,0]
t = np.linspace(0,30,2000)
sol = odeint(Lorenz,y0,t)


from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure()
ax = plt.axes(projection='3d')
ax.plot3D(sol[:,0],sol[:,1],sol[:,2])
plt.show()
In [ ]:
# %load EquationDifferentielle.py

Interpolation de Lagrange

Theorème. _Soient $f$ une fonction à valeurs réelles définie sur un intervalle $I=[a,b]$ et $n+1$ points deux à deux distincts $a_0,a_1,\ldots,a_n$ de l'intervalle $I$. Il existe un unique polynôme $P_n$ de degré au plus $n$ vérifiant_

$$ \forall i\in\{0,\ldots,n\},\ \ p_n(a_i)=f(a_i). $$

_Ce polynôme d'interpolation de Lagrange s'obtient dans la base des polynômes de Lagrange associée aux points $a_0,\ldots,a_n$_

$$ P_n(x) = \sum_{i=0}^n f(a_i)L_i(x),\ \ \ L_i(x) = \prod_{j\neq i}\dfrac{x-a_j}{a_i-a_j}. $$

Pour le choix des points $a_i$, on peut retenir entre autres possibilités:

  • les points équidistants de l'intervalle $I$: $a_i = a+i(b-a)/(n+1)$, $i\in\{0,\ldots,n\}$ ;

  • les points de Chebyshev : $a_i = \tfrac{a+b}{2} + \tfrac{b-a}{2}\cos\left(\tfrac{(2i+1)\pi}{2n+2}\right)$, $i\in\{0,\ldots,n\}$.

Exercice :

Définir une fonction Lagrange(a,i,x) qui prend en argument le tuple (ou la liste) des points a et renvoie la valeur $L_i(x)$.

Définir une fonction Interpolation(f,a,x) qui renvoie la valeur $P_n(x)$.

In [20]:
def Lagrange(a,i,x):
    s = [(x-a[j])/(a[i]-a[j]) for j in range(len(a)) if j != i]
    return prod(s)

n = 3
a = linspace(0,1,n+1)
t = linspace(0,1,200)
for i in range(n+1):
    plt.plot(t,[Lagrange(a,i,s) for s in t],a[i],1,'.')
plt.title("Base de Lagrange pour n="+str(n))
plt.show()
In [21]:
def Interpolation(f,a,x):
    y = 0
    for i in range(a.size):
        y = y + f(a[i])*Lagrange(a,i,x)
    return y

def f(x):
    return cos(10*x)

n = 4
a = linspace(0,1,n)
t = linspace(0,1,200)
plt.plot(t,[Interpolation(f,a,x) for x in t],t,f(t))
plt.title("Interpolation de Lagrange pour n="+str(n))
plt.show()

Exercice : le phénomène de Runge

Soient $I$ l'intervalle $I=[-1,1]$ et $f$ la fonction définie sur $I$ par

$$f(x) = \dfrac{1}{1+25x^2}.$$

À l'aide des fonctions précédemment définies, construire et visualiser le polynôme d'interpolation de Lagrange pour $n=6$ obtenu dans le cas des points équidistants ainsi que pour les points de Chebyshev.

In [22]:
def f(x):
    return 1./(1+25*x**2)

n = 6
a = linspace(-1,1,n)
c = array([cos((2*i+1)*pi/(2*n+2)) for i in range(n+1)])
t = linspace(-1,1,200)
plt.plot(t,[Interpolation(f,a,x) for x in t],t,[Interpolation(f,c,x) for x in t],t,f(t))
plt.title("Interpolation de Lagrange pour n="+str(n))
plt.legend(["Equidistant","Chebyshev","Fonction"])
plt.show()

n = 12
a = linspace(-1,1,n)
c = array([cos((2*i+1)*pi/(2*n+2)) for i in range(n+1)])
t = linspace(-1,1,200)
plt.plot(t,[Interpolation(f,a,x) for x in t],t,[Interpolation(f,c,x) for x in t],t,f(t))
plt.title("Interpolation de Lagrange pour n="+str(n))
plt.legend(["Equidistant","Chebyshev","Fonction"])
plt.show()

Séries de Fourier

Ayant noté $[x]$ la partie entière d'un réel $x$, on considère la fonction $1$-périodique

$$f(x)=e^{x-[x]}.$$ Après quelques longs calculs, on pourrait vérifier que les sommes partielles de la série de Fourier de $f$ sont données par

$$ S_N(x) = (e-1) + \sum_{n=1}^N \dfrac{2(e-1)}{1+4\pi^2 n^2}\cos(2\pi n x) - \sum_{n=1}^N \dfrac{4\pi n (e-1)}{1+4\pi^2 n^2}\sin(2\pi n x). $$

Sur une même figure, comparer les graphes de $f$, $S_3$, $S_5$ et $S_9$.

Pour la série de Fourier de $g(x) = |\sin(\pi x)|$ sur l'intervalle $[0,1]$:

$$ R_N(x) = \dfrac{2}{\pi} - \sum_{n=1}^N \dfrac{4}{\pi(4n^2-1)} \cos(2 \pi n x), $$

examiner numériquement la vitesse de convergence de l'erreur $\sup_{[0,1]} |g-R_N|$ en fonction de $N$.

Remarque : quelques possibilités de calcul symbolique sont proposées par la bibliothèque sympy:

from sympy import fourier_series, pi
from sympy.abc import x
s = fourier_series(x**2,(x,-pi,pi))
s.truncate(n=4)
In [27]:
# %load Fourier.py
from numpy import *
import matplotlib.pyplot as plt

def f(x):
    y = exp(x-floor(x))
    return y

def SF(n,x):
    m = exp(1.)-1.
    A = [m] + [2.*m/(1.+4.*pi**2*k**2) for k in range(1,n+1)]
    B = [0] + [-4.*pi*k*m/(1.+4.*pi**2*k**2) for k in range(1,n+1)]
    return A[0] + sum(A[k]*cos((2.*k*pi)*x) for k in range(1,n+1)) + sum(B[k]*sin((2.*k*pi)*x) for k in range(1,n+1))


x = linspace(-0.5,1.5,2000)
plt.plot(x,SF(11,x))
plt.plot(x,f(x),'r')
plt.show()


############################################

def g(x):
    y = np.abs(np.sin(np.pi*x))
    return y

    
def SG(n,x):
    A = [2./np.pi] + [-4./(np.pi*(4.*k**2-1)) for k in range(1,n+1)]
    return A[0] + sum(A[k]*np.cos((2.*k*np.pi)*x) for k in range(1,n+1))


x = np.linspace(-0.5,1.5,2000)
plt.plot(x,SG(11,x))
plt.plot(x,g(x),'r')
plt.show()
/Users/boutin/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:13: DeprecationWarning: Calling np.sum(generator) is deprecated, and in the future will give a different result. Use np.sum(np.from_iter(generator)) or the python sum builtin instead.
  del sys.path[0]
/Users/boutin/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:31: DeprecationWarning: Calling np.sum(generator) is deprecated, and in the future will give a different result. Use np.sum(np.from_iter(generator)) or the python sum builtin instead.