B. Boutin - benjamin.boutin@univ-rennes1.fr
int
float
bool
string
tuple
(non-modifiables)list
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
for
, if
, while
¶for
for element in liste:
INSTRUCTIONS
if
(elif
) else
if CONDITION1:
INSTRUCTIONS1
elif CONDITION2:
INSTRUCTIONS2
elif CONDITION3:
INSTRUCTIONS3
else:
INSTRUCTIONS4
while
while CONDITION:
INSTRUCTIONS
Quelques exemples de commandes pour la manipulation des listes
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)
Définir une fonction qui à $x\in\mathbb{R}$ associe $x^2+1$ si $x\geq 0$ et $\sqrt{-x}$ sinon.
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))
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. $$
# %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)
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.
# %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))
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.
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)
numpy
et matplotlib
¶Documentation pour les bilbiothèques utilisées:
On pourra parcourir en particulier le tutorial pyplot de matplotlib
Reprendre l'exercice précédent en utilisant les fonctions scipy.optimize.bisect
et scipy.optimize.newton
.
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])
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$.
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()
# %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()
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))
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()
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))
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()
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()
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()
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}.$$
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()
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()
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()
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()
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))$.
[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()
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()
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()
[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()
# %load EquationDifferentielle.py
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\}$.
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)$.
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()
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()
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.
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()
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)
# %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()