Cette page montre comment les formules des règles Simpson sont établies. Il n’est pas nécessaire de la lire si on ne cherche qu’à appliquer les règles. Si on souhaite cependant comprendre pourquoi elles sont vraies, on peut approfondir sa connaissance en lisant cette page.
La règle 1/3 et la règle 5/8 reposent toutes les deux sur une approximation quadratique des intervalles de fonction. Elles ont la même assise théorique, mais diffèrent seulement sur le domaine d’intégration. Elles sont donc traitées ensemble ci-dessous.
La règle 3/8 repose quant-à-elle sur l’usage de fonctions cubiques. Elle est traitée dans une section séparée, en dernier.
Représenter une fonction quadratique à partir de trois points
Le point commun de la règle 1/3 et de la règle 5/8 est l’usage d’une fonction quadratique \(y = ax^2+bx+c\) pour approximer la fonction qui nous intéresse. Pour rendre cette fonction utile aux fins d’application, il faut représenter les coefficients \(a, b\) et \(c\) en termes de trois points \(y_{-1}, y_{0}\) et \(y_{1}\), formant deux intervalles.
Le point de départ est donc une fonction quadratique « quelconque », mais dont l’intervalle d’intégration est les trois points \(-\Delta x, 0\) et \(\Delta x\). Une transformation linéaire des axes permet de représenter n’importe quelle fonction sur cet intervalle, si bien qu’il n’y a aucune perte de généralité.
Dans la Figure 1 ci-dessous, l’écart est \(\Delta x = 2\), si bien que l’intervalle de la fonction est [-2, 2].
Voir le calcul.
import numpy as npimport matplotlib.pyplot as pltdef plot_quadratic_area(a, b, c, dx):""" Affiche une fonction quadratique et l'aire sous la courbe sur [-dx, dx]. """# 1. Définition de la fonction et de l'intervalle f =lambda x: a*x**2+ b*x + c x_range = np.linspace(-dx, dx, 500) y_values = f(x_range)# Calcul de l'aire exacte pour l'affichage (via intégration numérique)# On utilise simpson pour la précisionfrom scipy.integrate import simpson area = simpson(y_values, x_range)# 2. Création de la figure fig, ax = plt.subplots(figsize=(5*3/2, 3*3/2)) ax.set_facecolor('#f8f9fa')# 3. Dessin de l'aire sous la courbe# On utilise une couleur douce (lightblue) avec de la transparence ax.fill_between(x_range, 0, y_values, color='skyblue', alpha=0.4, label=f'Aire $\\approx$ {area:.3f}')# 4. Dessin de l'équation quadratique (la courbe)# On utilise une couleur contrastée (navy) ax.plot(x_range, y_values, color='navy', linewidth=2.5, label=f'$f(x) = {a}x^2 {" + "if b >=0else" - "}{abs(b)}x {" + "if c >=0else" - "}{abs(c)}$')# 5. Mise en forme esthétique ax.set_title(f"Fonction quadratique sur $[-\\Delta x, \\Delta x]$ avec $\\Delta x = {dx}$", fontsize=14, fontweight='bold') ax.set_xlabel("$x$", fontsize=12) ax.set_ylabel("$f(x)$", fontsize=12)# Ajout de l'axe X (ligne de base) ax.axhline(0, color='black', linewidth=1) ax.axvline(0, color='black', linewidth=0.5, linestyle='--') # Axe Y pointillé ax.grid(True, linestyle=':', alpha=0.6) ax.legend(loc='upper right', frameon=True, facecolor='white')# Ajustement des limites pour voir les bornes $-\Delta x$ et $\Delta x$ ax.set_xlim(-dx *1.2, dx *1.2)# On ajuste l'axe Y pour qu'il ne soit pas trop écrasé ymin, ymax = np.min(y_values), np.max(y_values) ax.set_ylim(min(0, ymin) -0.5, ymax +0.5) plt.tight_layout() plt.show()# --- PARAMÈTRES À MODIFIER ---A =-1# Coefficient de x^2B =0# Coefficient de xC =5# ConstanteDELTA_X =2# Largeur de l'intervalle de chaque côté de l'origine# -----------------------------plot_quadratic_area(A, B, C, DELTA_X)
Figure 1: Une fonction quadratique.
Pour une fonction quadratique quelconque, soit \(f_{-1}, f_0\) et \(f_1\) les trois valeurs aux \(-\Delta x, 0\) et \(\Delta x\). C’est trois valeurs doivent satisfaire: \[\begin{align}
f_{-1} &= a(\Delta x)^2 - b \Delta x + c & f_0 &= c,\\
f_{1} &= a(\Delta x)^2 + b \Delta x + c.
\end{align}\] Ces trois équations permettent de trouver les paramètres \(a, b,c\) en fonction de \(f_{-1}, f_0\) et \(f_1\): \[\begin{align}
c&= f_0 , & b&=\frac{f_1-f_{-1}}{2\Delta x},\\
a &= \frac{f_{1} + f_{-1} - 2f_0}{2(\Delta x)^2}
\end{align}\]
Ces trois valeurs permettent de calculer \(a, b\) et \(c\) à partir des trois valeurs observées de la fonction.
Démonstration de la règle 1/3
Pour trouver l’aire sous la courbe de la fonction quadratique sur l’intervalle \([-\Delta x, \Delta x]\), on calcule son intégrale: \[\begin{align}
\int_{-\Delta x}^{\Delta x}ax^2+bx+cdx &= \left. \frac{a}{3}x^3 + \frac{b}{2}x^2 + cx\right|_{-\Delta x}^{\Delta x}
\end{align}\]
En subsituant les valeurs de \(a, b\) et \(c\) telles qu’identifiées dans la section précédente, on peut trouver: \[\begin{align}
\int_{-\Delta x}^{\Delta x}ax^2+bx+cdx&= \frac{a}{3}2\Delta x^3 + 2\Delta x f_0,\\
&= \frac{\Delta x}{3}\left[f_n + 4 f_0 + f_p\right].
\end{align}\] On peut remarquer la structure de poids \(1, 4, 1\) propre à la règle 1/3 pour trois points. Cette expression est la valeur exacte de l’intégrale pour une expression quadratique.
Application à des fonctions quelconques
Maintenant, soit une fonction \(g(x)\) quelconque dont on subdivise le domaine en intervalles \(\Delta x\) égaux, à raison de \(2n+1\) points de subdivision.
Soit \(g_0, g_1, g_2, g_3, g_4, g_5, \dots , g_{2n+1}\) les valeurs de \(g\) aux subdivisons. On peut alors décomposer chaque segment en approximation quadratique, dont l’aire est donnée par: \[\begin{align}
\int_a^b g(x)dx &\approx \frac{\Delta x}{3}\left[g_0 + 4g_1 + g_2\right] + \frac{\Delta x}{3}\left[g_2 + 4g_3 + g_4\right] + \frac{\Delta x}{3}\left[g_4 + 4g_5 + g_6\right] + \dots \\
&~~~\dots + \frac{\Delta x}{3}\left[g_{2n-1} + 4g_{2n} + g_{2n+1}\right]\\
&= \frac{\Delta x}{3}\left[g_0 + 4g_1 + 2g_2+ 4g_3 + 2g_4+ 4g_5 + 2g_6+ \dots +2g_{2n-1} + 4g_{2n} + g_{2n+1}\right].
\end{align}\]
On peut remarquer que lorsqu’on regroupe les segments ensemble, la structure de poids \([1, 4, 2, 4, \dots, 4, 1]\) propre à la règle 1/3 apparaît. Bien sûr, on doit pouvoir séparer la fonction en \(2n+1\) points pour être en mesure de l’appliquer.
Démonstration de la règle 5/8
La règle 5/8 emploie la même approximation quadratique, mais ce qui change, c’est l’intervalle d’intégration. On souhaite seulement approximer l’intervalle \([0, \Delta x]\). Ce faisant, l’intégrale qui nous intéresse est: \[\begin{align}
\int_0^{\Delta x} ax^2 + bx + cdx= a \frac{(\Delta x)^3}{3} + b \frac{(\Delta x)^2}{2} + c\Delta x.
\end{align}\] En remplaçant les valeurs pour \(a\), \(b\) et \(c\) identifiées en début de texte et avec un peu d’algèbre, on obtient: \[\begin{align}
\int_0^{\Delta x} ax^2 + bx + cdx = \frac{\Delta x}{12}\left[5f_{1} + 8 f_0 - f_{-1}\right]
\end{align}\]
Pour l’appliquer à une fonction quelconque, on l’approxime par la parabole définie sur un intervalle qui nous intéresse. Si l’intervalle est le dernier d’une suite (comme dans la Figure 2 ci-dessous), il faut ainsi choisir les deux points qui définissent l’intervalle (\(y_{n}, y_{n+1}\)), et le débutant l’intervalle précédent (\(y_{n-1}\)): \[\begin{align}
F \approx \frac{\Delta x}{12}\left[5 y_{n+1} + 8 y_n - y_{n-1}\right].
\end{align}\]
Voir le calcul.
import numpy as npimport matplotlib.pyplot as pltfrom scipy.interpolate import lagrangedef visualize_simpson_5_8_single():# 1. Définition de la fonction réelle f(x)# Une fonction cubique pour que la parabole d'approximation soit visiblement différente f =lambda x: 0.5* x**3-0.5* x**2+2* x +1# 2. Configuration des points (n-2, n-1, n) h =1.0 x_nodes = np.array([0.0, 1.0, 2.0]) # x_{n-2}, x_{n-1}, x_n y_nodes = f(x_nodes)# L'intervalle d'intégration cible pour l'aire est [x_{n-1}, x_n] soit [1.0, 2.0] a, b = x_nodes[1], x_nodes[2]# 3. Génération des courbes# Courbe de la fonction réelle pour l'affichage (haute résolution) x_fine = np.linspace(-0.5, 2.5, 500) y_fine = f(x_fine)# Approximation Quadratique (La parabole de la règle 5/8)# On interpole les 3 points pour obtenir le polynôme de degré 2 poly_quad = lagrange(x_nodes, y_nodes) x_quad_fine = np.linspace(x_nodes[0], x_nodes[2], 200) y_quad_fine = poly_quad(x_quad_fine)# 4. Visualisation plt.figure(figsize=(8, 7*2/3))# Affichage de la fonction réelle plt.plot(x_fine, y_fine, 'k-', label='Fonction réelle $f(x)$', linewidth=2)# Aire sous la courbe (vraie intégrale sur [x_{n-1}, x_n]) x_fill = np.linspace(a, b, 200) plt.fill_between(x_fill, f(x_fill), color='skyblue', alpha=0.3, label='Aire réelle $\int_{x_{n-1}}^{x_n} f(x)dx$')# La courbe quadratique (l'approximation de la règle) plt.plot(x_quad_fine, y_quad_fine, 'r-', label="Courbe de la règle 5/8 (Parabole)", linewidth=2.5)# Marquage des points de contrôle (nodos) plt.scatter(x_nodes, y_nodes, color='black', zorder=5, label='Points de contrôle ($y_{n-2}, y_{n-1}, y_n$)')# Esthétique plt.title(r"Visualisation de la règle $\frac{\Delta x}{12}(5y_{n} + 8y_{n-1} - y_{n-2})$ sur $[x_{n-1}, x_n]$", fontsize=14) plt.xlabel("x", fontsize=12) plt.ylabel("f(x)", fontsize=12) plt.legend(loc='upper left', frameon=True) plt.grid(True, which='both', linestyle='--', alpha=0.5) plt.tight_layout() plt.show()visualize_simpson_5_8_single()
Figure 2: La règle 5/8 approxime un intervalle unique.
Représenter une fonction cubique à partir de quatre points
La règle 3/8 repose sur l’usage d’une expression cubique \(y = ax^3 + bx^2 + cx +d\). Soit trois points \(f_{-1}, f_0, f_1, f_2\) de cette fonction à des intervalles \(\Delta x\) égaux. Sans perte de généralité, supposons que l’intervalle est \([-\Delta x, 2 \Delta x]\). Les valeurs doivent ainsi satisfaire: \[\begin{align}
f_{-1} &= a(-\Delta x)^3 + b(-\Delta x)^2 - c\Delta x +d,\\
f_{0} &= a(0)^3 + b(0)^2 + c\cdot 0 +d,\\
f_{1} &= a(\Delta x)^3 + b(\Delta x)^2 + c\Delta x +d,\\
f_{2} &= a(2\Delta x)^3 + b(2\Delta x)^2 + c2\Delta x +d.
\end{align}\] C’est un système de quatre équations à quatre inconnues. On peut résoudre pour identifier les valeurs \(a, b, c, d\) en fonction de \(f_{-1}, f_0, f_1\) et \(f_2\): \[\begin{align}
a &= \frac{f_2 - 3f_1 + 3f_{0}-f_{-1}}{6(\Delta x)^3},\\
b &= \frac{f_1 + f_{-1}-2f_0}{2(\Delta x)^2},\\
c&=\frac{6f_1 - 2f_{-1}-3f_0-f_2}{6 \Delta x},\\
d &= f_0.
\end{align}\]
Démonstration de la règle 3/8
Avec ces coefficients, on peut maintenant calculer l’intégrale: \[\begin{align}
\int_{-\Delta x}^{2\Delta x} ax^3 + bx^2 + cx + d dx &= \left. \frac{a}{4}x^4 + \frac{b}{3}x^3 + \frac{c}{2}x^2 + dx \right|_{-\Delta x}^{2\Delta x},\\
&=\frac{15}{4}\Delta x^4 + 3b\Delta x^3 + \frac{3}{2}c\Delta x^2 + d \Delta x.
\end{align}\] Avec pas mal d’algèbre (!), on arrive à la solution: \[\begin{align}
\int_{-\Delta x}^{2\Delta x} ax^3 + bx^2 + cx + d dx &= \frac{3}{8}\Delta x\left[f_{-1}+3f_0+3f_1+f_2\right].
\end{align}\]
Pour l’appliquer à une fonction quelconque, on emploie les valeurs aux quatres points d’intervalles égaux et on calcule la formule. Dès que la fonction \(f(x)\) qui nous intéresse n’est pas une fonction cubique, la règle fait une erreur d’approximation.
La Figure 3 illustre l’application de la règle avec \(\Delta x = 1\) et pour une fonction \(f(x)\) quelconque. Les points définissant les intervalles sont \(x=0, x=1, x=2\) et \(x=3\). Les valeurs respectives de la fonction sont \(f_{-1}, f_0, f_1\) et \(f_2\). On note que l’approximation cubique de la fonction génère des écarts à la courbe réelle. Ce faisant, le calcul d’intégrale découlant de la règle de Simpson donnera une valeur légèrement différente de la valeur réelle.
Voir le calcul.
import numpy as npimport matplotlib.pyplot as pltfrom scipy.interpolate import lagrangedef visualize_simpson_3_8():# 1. Définition de la fonction réelle f(x)# On utilise un polynôme de degré 4 pour que l'approximation cubique # de la règle 3/8 montre une différence visible avec la fonction réelle. f =lambda x: x**4-4*x**3+4*x**2+2*x +1# 2. Configuration des points (x_{-1}, x_0, x_1, x_2) h =1.0 x_nodes = np.array([0.0, 1.0, 2.0, 3.0]) # Les 4 points requis pour la règle 3/8 y_nodes = f(x_nodes)# L'intervalle d'intégration cible est [x_{-1}, x_2] soit [0.0, 3.0] a, b = x_nodes[0], x_nodes[3]# 3. Génération des courbes# Courbe de la fonction réelle (haute résolution) x_fine = np.linspace(-0.5, 3.5, 500) y_fine = f(x_fine)# Approximation Cubique (Le polynôme de degré 3 défini par la règle 3/8)# On interpole les 4 points pour obtenir le polynôme de degré 3 poly_cubic = lagrange(x_nodes, y_nodes) x_cubic_fine = np.linspace(x_nodes[0], x_nodes[3], 200) y_cubic_fine = poly_cubic(x_cubic_fine)# 4. Visualisation plt.figure(figsize=(8, 7*2/3))# Affichage de la fonction réelle plt.plot(x_fine, y_fine, 'k-', label='Fonction réelle $f(x)$', linewidth=2)# Aire sous la courbe (vraie intégrale sur [x_{-1}, x_2]) x_fill = np.linspace(a, b, 200) plt.fill_between(x_fill, f(x_fill), color='skyblue', alpha=0.3, label='Aire réelle $\int_{x_{-1}}^{x_2} f(x)dx$')# La courbe cubique (l'approximation de la règle) plt.plot(x_cubic_fine, y_cubic_fine, 'r--', label="Courbe de la règle 3/8 (Cubique)", linewidth=2.5)# Marquage des points de contrôle (nodos) plt.scatter(x_nodes, y_nodes, color='black', zorder=5, label='Points de contrôle ($f_{-1}, f_0, f_1, f_2$)')# Esthétique plt.title(r"Visualisation de la règle $\frac{3\Delta x}{8}(f_{-1} + 3f_0 + 3f_1 + f_2)$ sur $[x_{-1}, x_2]$", fontsize=14) plt.xlabel("x", fontsize=12) plt.ylabel("f(x)", fontsize=12) plt.legend(loc='upper left', frameon=True) plt.grid(True, which='both', linestyle='--', alpha=0.5) plt.tight_layout() plt.show()if__name__=="__main__": visualize_simpson_3_8()
Figure 3: La règle 3/8 approxime une fonction par une cubique.
