USMA1Q — Méthodes Numériques¶

Séance 8 — Intégration numérique I : fondamentaux et quadrature de Gauss¶

Conservatoire National des Arts et Métiers

25 Mars 2026

nicholas-anton.collins-craft@enpc.fr


Un lien vers Google Colab qui peut faire tourner le notebook qui contient les exercices se trouve ici :

https://colab.research.google.com/github/nickcollins-craft/USMA1Q-Methodes-Numeriques/blob/main/

Plan du cours¶

Jour Sujet
03 mars Fondamentaux de la programmation en Python
04 mars (matin) Fonctions, listes approfondies, fichiers et graphiques
04 mars (après-midi) Calcul numérique avec NumPy
11 mars (matin) Tableaux en mémoire, matrices & précision numérique
11 mars (après-midi) Résolution des systèmes linéaires
18 mars (matin) Interpolation & recherche des racines
18 mars (après-midi) Systèmes d'équations non linéaires
25 mars (matin) Intégration numérique I
25 mars (après-midi) Intégration numérique II (Monte Carlo)
02 avril Équations différentielles
08 avril Examen

Plan de la séance (indicatif)¶

Horaire Sujet
0:00 – 0:50 Quadrature de Newton-Cotes : trapèzes et Simpson
0:50 – 1:40 Quadrature de Gauss–Legendre
1:40 – 1:50 ☕ Pause
1:50 – 2:45 Intégration de données discrètes & intégrales cumulatives
2:45 – 3:40 Intégrales impropres et singularités
3:40 – 3:45 Récapitulatif

Motivation — pourquoi intégrer numériquement ?¶

De nombreuses quantités physiques importantes sont des intégrales :

QuantitéIntégrale
Variation d'enthalpie ΔH\(\int_{T_1}^{T_2} C_p(T)\,\mathrm{d}T\)
Ténacité (énergie absorbée)\(\int_0^{\varepsilon_f}\sigma(\varepsilon)\,\mathrm{d}\varepsilon\)
Rigidité d'un élément fini\(\int_0^L EA\left(\frac{\mathrm{d}N}{\mathrm{d}x}\right)^2\mathrm{d}x\)
Chaleur latente (DSC)\(\int_{T_1}^{T_2}\dot{q}(T)\,\mathrm{d}T\)

Problème : $C_p(T)$, $\sigma(\varepsilon)$ ne sont souvent disponibles que sous forme de données tabelées, pas de fonctions analytiques. Il faut donc intégrer numériquement.

1 · Quadrature de Newton-Cotes¶

Idée : remplacer l'intégrande par une fonction simple (polynôme) sur chaque sous-intervalle, puis intégrer ce polynôme exactement.

Avec $N$ points espacés de $h = (b-a)/(N-1)$, on obtient la règle composite.

Règle des trapèzes (méthode des trapèzes)¶

On relie les points par des segments droits :

$$\int_a^b f(x)\,\mathrm{d}x \approx h\left[\frac{f(x_0)}{2} + f(x_1) + f(x_2) + \cdots + f(x_{N-2}) + \frac{f(x_{N-1})}{2}\right]$$

Erreur : $\mathcal{O}(h^2)$ — si on divise $h$ par 2, l'erreur est divisée par 4.

En Python : scipy.integrate.trapezoid(y, x) où trapezoid = trapèze (méthode des trapèzes).

⚠️ Le $y$ est donné en prémier à cette méthode, le $x$ (si on a) est donné deuxième.

In [1]:
import numpy as np
from scipy.integrate import trapezoid

# Exemple : sin(x) sur [0, π] — valeur exacte = 2
N = 10
x = np.linspace(0, np.pi, N)
y = np.sin(x)

resultat = trapezoid(y, x)
print("N =", N, "  résultat =", round(resultat, 8), "  erreur =", round(abs(resultat - 2.0), 4))
# Affichage plus lisible :
print("N = %d   résultat = %.8f   erreur = %.2e" % (N, resultat, abs(resultat - 2.0)))

N = 10   résultat = 1.97965081   erreur = 0.0203
N = 10   résultat = 1.97965081   erreur = 2.03e-02

Règle de Simpson (méthode de Simpson)¶

On relie les points par des paraboles :

$$\int_a^b f(x)\,\mathrm{d}x \approx \frac{h}{3}\bigl[f(x_0) + 4f(x_1) + 2f(x_2) + 4f(x_3) + \cdots + 4f(x_{N-2}) + f(x_{N-1})\bigr]$$

Nécessite un nombre impair de points (nombre pair d'intervalles).

Erreur : $\mathcal{O}(h^4)$ — deux ordres meilleurs que les trapèzes pour une même régularité !

En Python : scipy.integrate.simpson(y, x) où simpson = Simpson

⚠️ Le $y$ est donné en prémier à cette méthode, le $x$ (si on a) est donné deuxième.

In [11]:
from scipy.integrate import trapezoid, simpson

exact = 2.0    # valeur exacte de sin(x) sur [0, π]

print("N    | erreur trapèzes | erreur Simpson")
print("-----|-----------------|---------------")
for N in [5, 11, 21, 51, 101, 201]:
    x = np.linspace(0, np.pi, N)
    y = np.sin(x)
    err_trap = abs(trapezoid(y, x) - exact)
    err_simp = abs(simpson(y, x) - exact)
    print("%4d | %.3e       | %.3e" % (N, err_trap, err_simp))

N    | erreur trapèzes | erreur Simpson
-----|-----------------|---------------
   5 | 1.039e-01       | 4.560e-03
  11 | 1.648e-02       | 1.095e-04
  21 | 4.114e-03       | 6.784e-06
  51 | 6.580e-04       | 1.733e-07
 101 | 1.645e-04       | 1.082e-08
 201 | 4.112e-05       | 6.765e-10

Convergence : erreur en fonction de $h$¶

Sur un graphique log-log, l'erreur décroit en droite avec la pente :

  • Trapèzes : pente $-2$ (erreur $\propto h^2$)
  • Simpson : pente $-4$ (erreur $\propto h^4$)

À vous de faire les exercices 1.1 et 1.2 !¶

2 · Quadrature de Gauss–Legendre¶

Idée : au lieu d'imposer des points équirépartis, on choisit à la fois les nœuds $x_i$ et les poids $w_i$ de façon optimale.

$$\int_{-1}^{1} f(t)\,\mathrm{d}t \approx \sum_{i=1}^{n} w_i f(t_i)$$

Propriété remarquable : une quadrature de Gauss à $n$ points intègre exactement tout polynôme de degré $\leq 2n-1$.

La méthode des trapèzes composite à $n$ points n'intègre exactement que les polynômes de degré 1. Simpson n'intègre exactement que de degré 3.

Nœuds et poids de Gauss–Legendre¶

Les nœuds $t_i \in [-1, 1]$ sont les racines des polynômes de Legendre. En Python : np.polynomial.legendre.leggauss(n).

t, w = np.polynomial.legendre.leggauss(n)

Changement de variable pour passer de $[-1, 1]$ à $[a, b]$ :

$$x = \frac{b-a}{2}\,t + \frac{a+b}{2}, \qquad \mathrm{d}x = \frac{b-a}{2}\,\mathrm{d}t$$

donc :

$$\int_a^b f(x)\,\mathrm{d}x = \frac{b-a}{2} \int_{-1}^{1} f\!\left(\frac{b-a}{2}t + \frac{a+b}{2}\right) \mathrm{d}t \approx \frac{b-a}{2}\sum_{i=1}^{n} w_i\, f(x_i)$$

In [4]:
import numpy as np

# Nœuds et poids sur [-1, 1]
t3, w3 = np.polynomial.legendre.leggauss(3)
print("n=3 — poids :", w3)
print("Somme des poids :", sum(w3), " (doit être 2.0 — longueur de [-1,1])")

n=3 — poids : [0.55555556 0.88888889 0.55555556]
Somme des poids : 2.0000000000000004  (doit être 2.0 — longueur de [-1,1])
In [5]:
def gauss_legendre(f, a, b, n):
    # Intégrale de f sur [a, b] par quadrature de Gauss-Legendre à n points.
    t, w = np.polynomial.legendre.leggauss(n)
    # Changement de variable : [-1,1] → [a,b]
    x = (b - a) / 2.0 * t + (a + b) / 2.0
    return (b - a) / 2.0 * np.sum(w * f(x))

# Démonstration : intégrer sin(x) sur [0, π]
print("Gauss-Legendre (n=3) :", round(gauss_legendre(np.sin, 0, np.pi, 3), 10))
print("Gauss-Legendre (n=5) :", round(gauss_legendre(np.sin, 0, np.pi, 5), 10))
print("Valeur exacte         : 2.0")

Gauss-Legendre (n=3) : 2.0013889136
Gauss-Legendre (n=5) : 2.0000001103
Valeur exacte         : 2.0

scipy.integrate.quad — intégration adaptative¶

quad (abréviation de quadrature) utilise la quadrature de Gauss-Kronrod adaptative : il raffine automatiquement là où la fonction varie rapidement.

from scipy.integrate import quad

valeur, erreur_estimee = quad(f, a, b)
  • Renvoie deux valeurs : le résultat et une estimation de l'erreur
  • Toujours vérifier que erreur_estimee est acceptable !
  • Gère les limites infinies : quad(f, 0, np.inf)
In [6]:
from scipy.integrate import quad

# Intégrer sin(x) sur [0, π]
valeur, erreur = quad(np.sin, 0, np.pi)
print("Résultat       :", valeur)
print("Erreur estimée :", erreur)   # doit être très petite
print("Erreur réelle  :", abs(valeur - 2.0))

Résultat       : 2.0
Erreur estimée : 2.220446049250313e-14
Erreur réelle  : 0.0

Gauss en pratique : éléments finis¶

Dans la méthode des éléments finis (MEF), les matrices de rigidité sont calculées (en pratique) par quadrature de Gauss sur chaque élément.

Pour un élément de barre 1D de longueur $L$, avec module $E$ et section $A$ :

$$k_{11} = \int_0^L EA \left(\frac{\mathrm{d}N_1}{\mathrm{d}x}\right)^2 \mathrm{d}x$$

Avec des fonctions de forme linéaires $N_1(x) = 1 - x/L$, on a $\mathrm{d}N_1/\mathrm{d}x = -1/L$.

Donc $k_{11} = EA/L$ pour un élément fini linéaire.

Gauss et sous-intégration¶

Le nombre de points de Gauss requis dépend du degré des fonctions de forme : l'intégrande $k_{ij} = EA\,(\mathrm{d}N_i/\mathrm{d}x)(\mathrm{d}N_j/\mathrm{d}x)$ est de degré $2(p-1)$ si les fonctions de forme $N_i$ sont de degré $p$.

Élément Degré des $N_i$ Degré de l'intégrande Pts de Gauss min.
Linéaire (2 nœuds) 1 0 (constant) 1
Quadratique (3 nœuds) 2 2 2
Cubique (4 nœuds) 3 4 3

Sous-intégration : utiliser trop peu de points introduit des erreurs dans la matrice de rigidité, la rendant trop souple (parfois très utile dans le contexte d'un intégration dit "reduit").

In [ ]:
import numpy as np
from numpy.polynomial.legendre import leggauss

# Exemple MEF 1D : terme k11 d'un élément quadratique (3 nœuds)
# N1(xi) = 0.5*xi*(xi-1)  => dN1/dx = (2/L)*(xi - 0.5), 
# k11 = ∫_{-1}^{1} EA*(dN1/dx)^2 * (L/2) dxi

E, A, L = 1.0, 1.0, 1.0

def integrande_k11(xi, E=1.0, A=1.0, L=1.0):
    dN1_dx = (2.0 / L) * (xi - 0.5)
    jac = L / 2.0
    return E * A * dN1_dx**2 * jac

def gauss_integrate_1d(func, n, E=1.0, A=1.0, L=1.0):
    xi, w = leggauss(n)
    return np.sum(w * func(xi, E, A, L))

k11_exact = 7.0 * E * A / (3.0 * L)

for n in [1, 2, 3]:
    k11_n = gauss_integrate_1d(integrande_k11, n, E, A, L)
    rel_err = abs(k11_n - k11_exact) / abs(k11_exact)
    print("n = %d  ->  k11 = %.10f   erreur relative = %.3e" % (n, k11_n, rel_err))

n = 1  ->  k11 = 1.0000000000   erreur relative = 5.714e-01
n = 2  ->  k11 = 2.3333333333   erreur relative = 1.903e-16
n = 3  ->  k11 = 2.3333333333   erreur relative = 0.000e+00

À vous de faire l'exercice 2 !¶

☕ Pause — 10 minutes¶

3 · Intégration de données discrètes¶

Situation courante : votre capteur enregistre des mesures à intervalles réguliers. Vous ne pouvez pas repositionner les points — vous êtes limités à Newton-Cotes, ou à interpoler d'abord, puis intégrer.

Trois stratégies :

Stratégie Outil Quand l'utiliser
Directe trapezoid(y, x) Données brutes, précision modérée
Spline puis intégration cs.integrate(a, b) Données propres, besoin de précision
Spline puis quad quad(cs, a, b) Idem, vérification indépendante

Intégrale cumulative¶

scipy.integrate.cumulative_trapezoid(y, x) (= intégrale cumulative par les trapèzes) permet de calculer non pas la valeur totale de l'intégrale, mais la valeur en fonction de la borne supérieure :

$$F(x) = \int_a^x f(t)\,\mathrm{d}t$$

Très utile pour :

  • Énergie absorbée en fonction de la déformation
  • Déplacement à partir de la vitesse
  • Enthalpy cumulée en fonction de la température

⚠️ Le $y$ est donné en prémier à cette méthode, le $x$ (si on a) est donné deuxième.

Comparaison des stratégies d'intégration¶

Pour des données bien résolues (pas $h$ petit), les trois méthodes donnent des résultats très proches. Les différences apparaissent pour des données grossières ou irrégulièrement espacées.

Règle pratique :

  • Si les données sont denses et régulières → trapezoid suffit
  • Si les données sont éparses ou irrégulières → interpoler d'abord, puis intégrer

À vous de faire les exercice 3.1 et 3.2 !¶

4 · Intégrales impropres et singularités¶

Une intégrale est dite impropre si :

  • Un (ou les deux) limite(s) d'intégration est infinie
  • L'intégrande a une singularité (divergence) dans l'intervalle

scipy.integrate.quad gère ces cas :

from scipy.integrate import quad

# Limite infinie
quad(f, 0, np.inf)        # ∫₀^∞ f(x) dx

# Limite doublement infinie
quad(f, -np.inf, np.inf)  # ∫₋∞^∞ f(x) dx

La syntaxe est quad(nom_de_function, borne_inferieur, borne_superieur).

In [9]:
from scipy.integrate import quad
import numpy as np

# Exemple 1 : ∫₀^∞ e^(-x) dx = 1
def f_exp(x):
    return np.exp(-x)

val, err = quad(f_exp, 0, np.inf)
print("∫₀^∞ e^(-x) dx =", val, "  (exact : 1.0)")

# Exemple 2 : ∫₋∞^∞ e^(-x²) dx = √π ≈ 1.7725...
def gaussienne(x):
    return np.exp(-x**2)

val2, err2 = quad(gaussienne, -np.inf, np.inf)
print("∫₋∞^∞ e^(-x²) dx =", round(val2, 6), "  (exact : √π =", round(np.sqrt(np.pi), 6), ")")

∫₀^∞ e^(-x) dx = 1.0000000000000002   (exact : 1.0)
∫₋∞^∞ e^(-x²) dx = 1.772454   (exact : √π = 1.772454 )

Singularités intégrables¶

Certaines fonctions divergent en un point, mais leur intégrale est finie.

Exemple : $\int_0^1 x^{-1/2}\,\mathrm{d}x = [2\sqrt{x}]_0^1 = 2$

quad le gère automatiquement dans les cas modérés. Pour une singularité forte, on peut indiquer son emplacement :

quad(f, 0, 1, points=[0.5])   # Signale une singularité en x=0.5
In [10]:
from scipy.integrate import quad

# Singularité en x=0 : f(x) = 1/sqrt(x)
def f_sing(x):
    return 1.0 / np.sqrt(x)

val, err = quad(f_sing, 0, 1)
print("∫₀¹ x^(-1/2) dx =", round(val, 8), "  (exact : 2.0)")
print("Erreur estimée  =", err)

∫₀¹ x^(-1/2) dx = 2.0   (exact : 2.0)
Erreur estimée  = 5.773159728050814e-15

À vous de faire l'exercice 4 !¶

Points clés de la séance¶

Concept Ce qu'il faut retenir
Trapèzes Erreur $\mathcal{O}(h^2)$ ; scipy.integrate.trapezoid(y, x)
Simpson Erreur $\mathcal{O}(h^4)$ — bien meilleur pour les fonctions lisses ; scipy.integrate.simpson(y, x)
Gauss–Legendre $n$ points intègrent exactement les polynômes de degré $2n-1$ ; acceder aux points et poids avec numpy.polynomial.legendre.leggauss(n)
scipy.integrate.quad Gauss-Kronrod adaptatif ; renvoie (valeur, erreur_estimée) ; toujours vérifier l'erreur
Cumulative scipy.integrate.cumulative_trapezoid(y, x, initial=0) — intégrale en fonction de la borne supérieure
Données discrètes Interpoler d'abord (spline) puis intégrer si les données sont éparses
Impropre scipy.integrate.quad(f, 0, np.inf) — quad gère les limites infinies et les singularités modérées
Vérification Comparer au moins deux méthodes ; vérifier erreur_estimée de quad