TP3 Intégration¶
Dans ce TP on va s'intéresser au calcul d'intégrales. On commencera par calculer une intégrale par la méthode des rectangles avant d'utiliser une fonction pré-programmée dans python.
1) La méthode des rectangles¶
Cette méthode consiste à approximer une intégrale en divisant l'intervalle d'intégration en petits sous-intervalles, puis en utilisant des rectangles pour estimer l'aire sous la courbe dans chaque sous-intervalle.
La formule de la méthode des rectangles (à gauche) pour approcher une intégrale avec une grille uniforme sur $[a,b]$ est donnée par :
$$ \int_a^b f(t) \, dt \approx \frac{b-a}{n}\sum_{j=0}^{n-1} f(a_j)\qquad \text{avec}\quad a_j = a + j\frac{b-a}{n}\qquad j\in \{0,\dots,n\}. $$
In [1]:
import numpy as np
import matplotlib.pyplot as plt
a) Ecrire une fonction integrale qui prend en argument le tableau f_grid constituée de $(f(a_0),\dots,f(a_{n-1}))$, et les bornes a et b, puis renvoie la valeur approchée de l'intégrale par la méthode des rectangles.¶
Pour commencer, on pourra prendre $f(x) = \sin(x)$ ou $\cos(x)$ sur l'intervalle $[0,\pi]$ pour vérifier que son code fonctionne correctement (car on connait les valeurs exactes des intégrales).
En effet, on a $$ \int_0^\pi \sin(t) dt = [-\cos(t)]^\pi_0 = 2 \qquad\text{et}\qquad \int_0^\pi \cos(t) dt = [\sin(t)]^\pi_0 = 0. $$
In [2]:
def integrale(f_grid, a, b):
n = len(f_grid)
return sum(f_grid) * (b-a)/n
In [3]:
a, b = 0, np.pi # défini les bornes de l'intervalle
n = 100 # nombre de point de discrétisation de $f$
a_grid = np.arange(a, b, (b-a)/n) # grille de point a_j
# np.arange(point de départ, point d'arrivé, pas entre les points)
In [4]:
f = lambda x: np.sin(x) # définition de la fonction f à intégrer
# équivalent à
# def f(x):
# return np.sin(x)
f_grid = f(a_grid) # évaluation de f sur la grille de point
integrale(f_grid, a, b) # évaluation de l'intégrale
Out[4]:
np.float64(1.9998355038874436)
In [5]:
# comme précédement avec la fonction cos
f = lambda x: np.cos(x)
f_grid = f(a_grid)
integrale(f_grid, a, b)
Out[5]:
np.float64(0.0314159265358977)
b) Illustrer en fonction $n$ la convergence de la méthode des rectangles vers la valeur de l'intégrale.¶
On pourra prendre $f(x) = \sin(x)$ ou $\cos(x)$ sur l'intervalle $[0,\pi]$.
In [6]:
n_grid = range(10, 300, 5) # grille de n pour évaluer l'erreur d'approximation
errors = [] # liste enregistrant les erreurs pour les différents n
f = lambda x: np.sin(x) # on définit la fonction à intégrer
true_int = 2 # vraie valeur de l'intégrale
for n in n_grid: # pour chaque n
a_grid = np.arange(a, b, (b-a)/n) # on crée une grille
f_grid = f(a_grid) # on évalue la fonction sur la grille
err = np.abs(integrale(f_grid, a, b) - true_int) # on calcul l'erreur d'approximation de l'intégrale
errors.append(err) # on enregistre l'erreur dans la liste
In [7]:
plt.plot(n_grid, errors)
Out[7]:
[<matplotlib.lines.Line2D at 0x7849dd4e97c0>]
Comme attendu, l'erreur d'approximation décroit quand $n$ augment, c'est à dire au fur et à mesure qu'on raffine la grille de point sur laquelle on évalue $f$.
2) La fonction quad de scipy.integrate.¶
Cette fonction permet de calculer l'intégrale d'une fonction $f$ sur un intervalle donné $[a,b]$. Elle utilise des techniques d'intégration numérique bien plus avancée que la simple méthode des rectangles.
Pour utiliser cette fonction on procède de la façon suivante pour une fonction $f$ sur l'intervalle $[0,1]$:
In [8]:
from scipy.integrate import quad # on importe la fonction permettant le calcul d'intégrale
In [9]:
# Définir la fonction à intégrer
f = lambda x: x**2
# Calculer l'intégrale de f(x) entre 0 et 1
result, error = quad(f, 0, 1)
result, error
# result contient une valeur approchée de l'intégrale, et
# error une estimation de l'erreur d'approximation
Out[9]:
(0.33333333333333337, 3.700743415417189e-15)
a) Reproduire les résultats de la question 1-a avec la fonction quad.¶
In [10]:
f = lambda x: np.sin(x) # fonction à intégrer
a, b = 0, np.pi # borne de l'intervalle
result, error = quad(f, a, b) # évalution de la fonction quad
result # contient la valeur approchée de l'intégrale
# error contient une estimation de l'erreur d'approximation
Out[10]:
2.0
In [11]:
# même chose pour la fonction cos
f = lambda x: np.cos(x)
a, b = 0, np.pi
result, error = quad(f, a, b)
result
Out[11]:
4.9225526349740854e-17
b) Illustrer numériquement la fonction porte.¶
$$ \Pi(t) = \begin{cases} 1 & \text{ si } |t| \leq T/2 \\ 0 & \text{ sinon.} \end{cases} $$
pour $T=4$.
In [12]:
T = 4
Pi = lambda t: 1.*(np.abs(t) < T/2) # définit la fonction porte
In [13]:
t = np.linspace(-2*T, 2*T, 200) # grille de temps pour illustrer la fonction
plt.plot(t,Pi(t))
Out[13]:
[<matplotlib.lines.Line2D at 0x7849dd3597f0>]
c) Calculer numériquement sa transformée de Fourier $\mathcal{F}(\Pi)$ (on séparera partie réelle et imaginaire pour l'intégration). Comparer avec sa formule théorique.¶
In [14]:
def T_Fourier(f, s):
f_re_tmp = lambda t: f(t) * np.cos(2*np.pi*s*t)
# définit la partie réelle à intégrer dans la transformée de Fourier
tf_re, _ = quad(f_re_tmp, -10, 10)
# évalue la partie réelle de la transformée de Fourier
# on calcul l'intégrale sur un intervalle suffisament grand
# (même si pour cette exemple ca n'a pas vraiment de sens...)
f_im_tmp = lambda t: f(t) * np.sin(-2*np.pi*s*t)
# définit la partie réelle à intégrer dans la transformée de Fourier
tf_im, _ = quad(f_im_tmp, -10, 10)
# évalue la partie réelle de la transformée de Fourier
return np.array([tf_re, tf_im])
# on ranvoie un tableau numpy contenant
# la partie réelle et la partie imaginaire
# de la transformée de Fourier
In [15]:
freq = np.linspace(-4, 4, 100) # grille de fréquence
TF = np.array([T_Fourier(Pi,s) for s in freq])
# créer un tableau numpy contenant la valeur
# de la transformée de Fourier pour chaque $s$
TF_theorie = T*np.sinc(freq*T)
# Formule théorique de la transformée de Fourier pour la comparaison
/tmp/ipykernel_8010/3583195404.py:6: IntegrationWarning: The integral is probably divergent, or slowly convergent. tf_re, _ = quad(f_re_tmp, -10, 10)
In [16]:
TF
Out[16]:
array([[ 7.90117197e-15, 0.00000000e+00],
[-6.90132074e-02, 0.00000000e+00],
[-7.43030631e-02, 0.00000000e+00],
[-8.05234020e-03, 0.00000000e+00],
[ 6.88917232e-02, 0.00000000e+00],
[ 8.26896107e-02, 0.00000000e+00],
[ 1.71373956e-02, 0.00000000e+00],
[-6.80851545e-02, 0.00000000e+00],
[-9.14856876e-02, 0.00000000e+00],
[-2.74016900e-02, 0.00000000e+00],
[ 6.64926497e-02, 0.00000000e+00],
[ 1.00759514e-01, 0.00000000e+00],
[ 3.90402652e-02, 0.00000000e+00],
[-6.39867023e-02, 0.00000000e+00],
[-1.10611142e-01, 0.00000000e+00],
[-5.23186419e-02, 0.00000000e+00],
[ 6.04005620e-02, 0.00000000e+00],
[ 1.21187354e-01, 0.00000000e+00],
[ 6.76073031e-02, 0.00000000e+00],
[-5.55081516e-02, 0.00000000e+00],
[-1.32705560e-01, 0.00000000e+00],
[-8.54379190e-02, 0.00000000e+00],
[ 4.89907768e-02, 0.00000000e+00],
[ 1.45493807e-01, 0.00000000e+00],
[ 1.06599207e-01, 0.00000000e+00],
[-4.03793155e-02, 0.00000000e+00],
[-1.60061317e-01, 0.00000000e+00],
[-1.32309409e-01, 0.00000000e+00],
[ 2.89479490e-02, 0.00000000e+00],
[ 1.77230820e-01, 0.00000000e+00],
[ 1.64547886e-01, 0.00000000e+00],
[-1.35044433e-02, 0.00000000e+00],
[-1.98406833e-01, 0.00000000e+00],
[-2.06748336e-01, 0.00000000e+00],
[-8.06316270e-03, 0.00000000e+00],
[ 2.26176779e-01, 0.00000000e+00],
[ 2.65416119e-01, 0.00000000e+00],
[ 3.98926732e-02, 0.00000000e+00],
[-2.65852763e-01, 0.00000000e+00],
[-3.54517943e-01, 0.00000000e+00],
[-9.13496716e-02, 0.00000000e+00],
[ 3.30277977e-01, 0.00000000e+00],
[ 5.10406429e-01, 0.00000000e+00],
[ 1.89096341e-01, 0.00000000e+00],
[-4.60362714e-01, 0.00000000e+00],
[-8.66442364e-01, 0.00000000e+00],
[-4.51228398e-01, 0.00000000e+00],
[ 8.93478766e-01, 0.00000000e+00],
[ 2.62308213e+00, 0.00000000e+00],
[ 3.83034039e+00, 0.00000000e+00],
[ 3.83034039e+00, 0.00000000e+00],
[ 2.62308213e+00, 0.00000000e+00],
[ 8.93478766e-01, 0.00000000e+00],
[-4.51228398e-01, 0.00000000e+00],
[-8.66442364e-01, 0.00000000e+00],
[-4.60362714e-01, 0.00000000e+00],
[ 1.89096341e-01, 0.00000000e+00],
[ 5.10406429e-01, 0.00000000e+00],
[ 3.30277977e-01, 0.00000000e+00],
[-9.13496716e-02, 0.00000000e+00],
[-3.54517943e-01, 0.00000000e+00],
[-2.65852763e-01, 0.00000000e+00],
[ 3.98926732e-02, 0.00000000e+00],
[ 2.65416119e-01, 0.00000000e+00],
[ 2.26176779e-01, 0.00000000e+00],
[-8.06316270e-03, 0.00000000e+00],
[-2.06748336e-01, 0.00000000e+00],
[-1.98406833e-01, 0.00000000e+00],
[-1.35044433e-02, 0.00000000e+00],
[ 1.64547886e-01, 0.00000000e+00],
[ 1.77230820e-01, 0.00000000e+00],
[ 2.89479490e-02, 0.00000000e+00],
[-1.32309409e-01, 0.00000000e+00],
[-1.60061317e-01, 0.00000000e+00],
[-4.03793155e-02, 0.00000000e+00],
[ 1.06599207e-01, 0.00000000e+00],
[ 1.45493807e-01, 0.00000000e+00],
[ 4.89907768e-02, 0.00000000e+00],
[-8.54379190e-02, 0.00000000e+00],
[-1.32705560e-01, 0.00000000e+00],
[-5.55081516e-02, 0.00000000e+00],
[ 6.76073031e-02, 0.00000000e+00],
[ 1.21187354e-01, 0.00000000e+00],
[ 6.04005620e-02, 0.00000000e+00],
[-5.23186419e-02, 0.00000000e+00],
[-1.10611142e-01, 0.00000000e+00],
[-6.39867023e-02, 0.00000000e+00],
[ 3.90402652e-02, 0.00000000e+00],
[ 1.00759514e-01, 0.00000000e+00],
[ 6.64926497e-02, 0.00000000e+00],
[-2.74016900e-02, 0.00000000e+00],
[-9.14856876e-02, 0.00000000e+00],
[-6.80851545e-02, 0.00000000e+00],
[ 1.71373956e-02, 0.00000000e+00],
[ 8.26896107e-02, 0.00000000e+00],
[ 6.88917232e-02, 0.00000000e+00],
[-8.05234020e-03, 0.00000000e+00],
[-7.43030631e-02, 0.00000000e+00],
[-6.90132074e-02, 0.00000000e+00],
[ 7.90117197e-15, 0.00000000e+00]])
On voit que les parties imaginaires sont nulles, comme attendu.
In [17]:
plt.plot(freq, TF[:,0], label="calcul num")
plt.plot(freq, TF_theorie, label="théorie")
plt.legend()
Out[17]:
<matplotlib.legend.Legend at 0x7849d32a8b60>
On ne voit pas la différence entre la formule théorique
d) Essauyer avec d'autres fonctions (comme celles du TD 3 de math 6 - Analyse)¶
On va faire exactement la même chose pour la fonction
$$ f(t) = \begin{cases} e^{-t}&\text{si }t>0 \\ 0&\text{si }t<=0 \end{cases} $$
illustrée ci-dessous, et dont on sait que la transformée de Fourier est (on l'a calculé en TD)
$$ \mathcal{F}f(s) = \frac{1}{1+2i \pi s}\qquad s\in \mathbb{R}. $$
In [18]:
f = lambda t: np.exp(-t)*(t>0)
t = np.linspace(-2, 10, 200)
plt.plot(t,f(t))
Out[18]:
[<matplotlib.lines.Line2D at 0x7849d31649b0>]
In [19]:
# évaluation de la transformée de Fourier
freq = np.linspace(-4, 4, 200)
TF = np.array([T_Fourier(f,s) for s in freq])
TF_theorie = (1+2j*np.pi*freq)**(-1)
In [20]:
# Illsutration de la partie imaginaire
# dans la 1ère colonne de TF
plt.plot(freq, TF[:,0], label="calcul num")
plt.plot(freq, TF_theorie.real, label="théorie")
plt.legend()
Out[20]:
<matplotlib.legend.Legend at 0x7849d3307350>
In [21]:
# Illsutration de la partie imaginaire
# dans la 1ère colonne de TF
plt.plot(freq, TF[:,1], label="calcul num")
plt.plot(freq, TF_theorie.imag, label="théorie")
plt.legend()
Out[21]:
<matplotlib.legend.Legend at 0x7849d31908f0>
