Intégration numérique (quadrature)

Cours ENPC - Pratique du calcul scientifique

Le problème

Objectif : calculer numériquement

\[ I := \int_{\Omega} u(\mathbf{\boldsymbol{x}}) \, \mathrm d\mathbf{\boldsymbol{x}} \] où $\Omega \subset {\mathbb{{R}}}^d$ et $u\colon \Omega \to {\mathbb{{R}}}$.

Pour simplifier on se concentre sur le cas 1D : $\Omega = [-1, 1] \subset {\mathbb{{R}}}$
Le cas général 1D sur $[a,b]$ exploite la bijection \[ \begin{array}{rrcl} \zeta\colon &[-1, 1] &\to& [a, b] \\ &y &\mapsto& \frac{b+a}{2} + \frac{(b-a)}{2} y \end{array} \] de sorte que \[ \int_{a}^{b} u(x) \, \mathrm dx = \int_{-1}^{1} u\bigl(\zeta(y)\bigr) \, \zeta'(y) \, \mathrm dy = \frac{b-a}{2}\int_{-1}^{1} u \circ \zeta (y) \, \mathrm dy, \]
Question: peut-on approximer $I$ par une somme de la forme

\[ \widehat I = \sum_{i=0}^n w_i \, u(x_i) \quad ? \]

La méthode de Newton-Cotes fermée

Idée principale

Calculer l’interpolation polynomiale $\widehat u$ de $u$ en des nœuds équidistants $-1 = x_0 < x_1 < \dots < x_{n-1} < x_n = 1$.
Puis calculer l’approximation

\[ \widehat I \approx \int_{-1}^{1} \widehat u(x) \, \mathrm dx \]

Expression explicite du polynôme interpolant : \[ \widehat u(x) = \sum_{i=0}^{n} u(x_i) \varphi_i(x), \qquad \text{où} \qquad \varphi_{i}(x) = \prod_{\substack{j = 0\\ j \neq i}}^{n} \frac {x - x_j} {x_i - x_j}. \]
On en déduit \[ \widehat I = \int_{-1}^{1} \widehat u(x) \, \mathrm dx = \sum_{i=1}^n w_i u(x_i), \qquad \text{avec} \qquad w_i := \int_{-1}^{1} \varphi_i(x) \, \mathrm dx. \]

Poids de Newton-Cotes par calcul formel

using SymPy
@vars x

nodes(n) = Sym[-1 + 2i//n for i in 0:n]

lagrange(nodes, i) = prod(x .- nodes[1:end .!= i]) /
                     prod(nodes[i] .- nodes[1:end .!= i]) 

for n ∈ (1,2,4)
  nd = nodes(n) ; println("**********") ; @show n ; println(nd)
  println([integrate(lagrange(nd, i),(x, -1, 1)) for i ∈ eachindex(nd)])
end

**********
n = 1
Sym[-1, 1]
Sym[1, 1]
**********
n = 2
Sym[-1, 0, 1]
Sym[1/3, 4/3, 1/3]
**********
n = 4
Sym[-1, -1/2, 0, 1/2, 1]
Sym[7/45, 32/45, 4/15, 32/45, 7/45]

\[ w_i := \int_{-1}^{1} \varphi_i(x) \, \mathrm dx \]

\[~\]

Méthode	Formule d’intégration
Trapézoïdal	$\widehat I = u(-1) + u(1)$
Simpson	$\widehat I = \frac{1}{3} u(-1) + \frac{4}{3} u(0) + \frac{1}{3} u(1)$
Boole	$\widehat I = \frac{7}{45} u(-1) + \frac{32}{45} u\left(-\frac{1}{2}\right) + \frac{4}{15} u(0) + \frac{32}{45} u\left(\frac{1}{2}\right) + \frac{7}{45} u(1)$

\[~\]

La méthode de Newton-Cotes fermée : exemples

\[~\]

$n$	$d$	Méthode	Formule d’intégration
1	1	Trapézoïdal	$\widehat I = u(-1) + u(1)$
2	3	Simpson	$\widehat I = \frac{1}{3} u(-1) + \frac{4}{3} u(0) + \frac{1}{3} u(1)$
4	5	Boole	$\widehat I = \frac{7}{45} u(-1) + \frac{32}{45} u\left(-\frac{1}{2}\right) + \frac{4}{15} u(0) + \frac{32}{45} u\left(\frac{1}{2}\right) + \frac{7}{45} u(1)$

\[~\]

$n =$ degré polynomial de construction
$d =$ degré de précision $=$ le plus haut degré polynomial tel que l’intégration est exacte

$n=2 ⇒ d=3$ car $\int_{x=-1}^1 x^3 \mathrm dx = \frac{1}{3} (-1)^3 + \frac{4}{3} (0)^3 + \frac{1}{3} (1)^3=0$
En principe, on peut construire des règles d’intégration de degré arbitrairement élevé moyennant d’augmenter le nombre de nœuds d’intégration.

Échec de la méthode de Newton-Cotes

… même avec beaucoup de nœuds car l’interpolation polynomiale ne converge pas toujours.

Code

import Polynomials
import Plots
Plots.default(linewidth=2)
n = 12 ; X = LinRange(-1, 1, n)
f(x) = 1 / (1 + 25 * x^2)
p = Polynomials.fit(X, f.(X))
x = -1:0.005:1
Plots.plot(x, f.(x), label="Fonction de Runge")
Plots.plot!(legend=:topright, size=(900,500))
Plots.plot!(x, p.(x), fillrange = 0 .* x, fillalpha = 0.35, label="Newton-cotes avec " * string(n) * " nœuds")
Plots.scatter!(X, f.(X), label="Nœuds d'interpolation")

Méthode de Newton-Cotes composite

Idée : s’appuyer sur l’interpolation par morceaux

Rappel sur l’erreur pour une interpolation polynomiale, nœuds équidistants \[ \fcolorbox{bleu_C}{bleu_light_C}{$ e_n(x) := u(x) - \widehat u(x) = \frac{u^{(n+1)}(\xi)}{(n+1)!} (x-x_0) \dotsc (x - x_n) $} \]

\[ \fcolorbox{bleu_C}{bleu_light_C}{$ \begin{align} &E_n := \sup_{x ∈ [a, b]} \bigl\lvert e_n(x) \bigr\rvert \leqslant\frac{C_{n+1}}{4(n+1)} \left(\frac{b-a}{n} \right)^{n+1}\\ &\text{avec } C_{n+1} = \sup_{x ∈ [a, b]} \left\lvert u^{(n+1)}(x) \right\rvert \text{ et } h = \frac{b-a}{n} \end{align} $} \]
Rappel sur l’erreur pour une interpolation par morceaux
- intervalle $[a,b]$ discrétisé avec $n+1$ nœuds, $h=\frac{b-a}{n}$
- polynôme de degré $m$ sur $[x_{i},x_{i+1}]$
- erreur \[ \fcolorbox{bleu_C}{bleu_light_C}{$ \sup_{x ∈ [x_{i},x_{i+1}]} \bigl\lvert e_n(x) \bigr\rvert \leqslant\frac{C_{m+1}}{4(m+1)} \left(\frac{h}{m}\right)^{m+1} $} \]
- $n$ a vocation à tendre vers l’infini, pas $m$
- $m$ est fixé petit donc l’erreur est majorée uniformément par $C h^{m+1}$

Règle composite trapézoïdale $\color{red}{m=1}$

$n+1$ nœuds équidistants $a=x_0< x_1< \dotsc < x_n=b$, pas $h=\frac{b-a}{n}$
Intégration trapézoïdale sur $[x_{i},x_{i+1}]$

\[ \int_{x_i}^{x_{i+1}} u(x) \, \mathrm dx = \frac{h}{2} \int_{-1}^{1} u \circ \zeta (y) \, \mathrm dy \approx \frac{h}{2} \bigl( u \circ \zeta(-1) + u \circ \zeta(1) \bigr) = \frac{h}{2} \bigl(u(x_i) + u(x_{i+1})\bigr) \]
application sur chaque sous-intervalle

\[ \begin{align} \int_{a}^{b} u(x) \, \mathrm dx &= \sum_{i=0}^{n-1} \int_{x_i}^{x_{i+1}} u(x) \, \mathrm dx \approx \frac{h}{2}\sum_{i=0}^{n-1} \bigl( u(x_i) + u(x_{i+1}) \bigr) \\ &= \frac{h}{2} \bigl( u(x_0) + 2 u(x_1) + 2 u(x_2) + \dotsb + 2 u(x_{n-2}) + 2 u(x_{n-1}) + u(x_n) \bigr) \end{align} \]

Exercice : démontrer le théorème suivant

On suppose $u \in \mathcal{C}^2([a, b])$ et $\widehat I_h$ désigne l’approximation trapézoïdale de l’intégrale $I$ de $u$. On a alors \[ \bigl\lvert I - \widehat I_h \bigr\rvert \leqslant\frac{b-a}{12} C_2 h^2, \qquad C_2 := \sup_{\xi \in [a, b]} \lvert u''(\xi)\rvert \]

Règle de Simpson composite

Construction de la règle

$n+1$ nœuds équidistants $a=x_0< x_1< \dotsc < x_n=b$, pas $h=\frac{b-a}{n}$
$n$ pair ⇒ nombre pair d’intervalles et nombre impair de nœuds

Découpage \[ \int_{a}^{b} u(x) \, \mathrm dx =\sum_{i=0}^{n/2-1} \int_{x_{2i}}^{x_{2i+2}} u(x) \, \mathrm dx \]
Rappel : Simpson sur $[-1,1]$ \[ \int_{-1}^{1} u(x) \, \mathrm dx \approx \frac{1}{3} \bigl(u(-1) + 4 u(0) + u(1)\bigr) \]
Donc sur $[x_{2i},x_{2i+2}]$ \[ \begin{align} \int_{x_{2i}}^{x_{2i+2}} u(x) \, \mathrm dx & \approx \frac{x_{2i+2}-x_{2i}}{2} × \frac{1}{3} \Bigl(u(x_{2i}) + 4 u(x_{2i+1}) + u(x_{2i+2})\Bigr)\\ & \approx \frac{h}{3} \Bigl(u(x_{2i}) + 4 u(x_{2i+1}) + u(x_{2i+2})\Bigr) \end{align} \]
Règle de Simpson composite \[ \widehat I_h \approx \frac{h}{3} \Bigl( u(x_0) + 4 u(x_1) + 2 u(x_2) + 4 u(x_3) + 2 u(x_4) + \dotsb + 2 u(x_{n-2}) + 4 u(x_{n-1}) + u(x_n)\Bigr) \]

Contrôle de l’erreur

On suppose $u \in \mathcal{C}^4([a, b])$. On a alors \[ \bigl\lvert I - \widehat I_h \bigr\rvert \leqslant\frac{b-a}{180} C_4 h^4, \qquad C_4 := \sup_{\xi \in [a, b]} \lvert u^{(4)}(\xi)\rvert \]

Démonstration dans le polycopié.
La démonstration exploite le fait que le degré de précision de la règle de Simpson est 3 et non 2.

Extrapolation de Richardson

Idée de l’extrapolation de Richardson

Méthode générale pour estimer $J(0)$ sachant que l’on connaît $J(h)$ pour divers $h>0$

Développement de Taylor \[ J(h) = J(0) + J'(0) h + J''(0) \frac{h^2}{2} + J^{(3)}(0) \frac{h^3}{3!} + \dotsb + J^{(n)}(0) \frac{h^n}{n!} + \mathcal O(h^{n+1}) \]
Développement en $\frac{h}{2}$ \[ J(h/2) = J(0) + J'(0) \frac{h}{2} + J''(0) \frac{h^2}{2 × 2^2} + J^{(3)}(0) \frac{h^3}{3! × 2^3} + \dotsb + J^{(n)}(0) \frac{h^n}{n! × 2^n} + \mathcal O(h^{n+1}) \]
On veut éliminer le terme en $h$ \[ J_1(h/2)=2J(h/2)-J(h)=J(0) - J''(0) \frac{h^2}{4} + \mathcal O(h^{3}) \]
Puis le terme en $h^2$ dans $J_1(h/2)$ \[ J_2(h/4) = \frac{4 J_1(h/4) - J_1(h/2)}{4 - 1} = J(0) + \mathcal O(h^3) \]
Généralisation \[ J_k(h/2^k) = \frac{2^k J_{k-1}(h/2^k) - J_{k-1}(h/2^{k-1})}{2^k - 1} = J(0) + \mathcal O(h^{k+1}) \]

Cas d’un développement initial en puissances paires \[ J(h) = J(0) + \sum_{k=0}^n J^{(2k)}(0) \frac{h^{2k}}{(2k)!} + \mathcal O(h^{2n+2}) \]
Généralisation \[ J_k(h/2^k) = \frac{2^{2k} J_{k-1}(h/2^k) - J_{k-1}(h/2^{k-1})}{2^{2k} - 1} \]

Méthode de Romberg

Application de l’extrapolation de Richardson à l’intégration

Méthode trapézoïdale composite avec $h \in \left\{ \frac{b-a}{n}: n \in {\mathbb{{N}}}\right\}$ \[ J(h)= \frac{h}{2} \bigl( u(x_0) + 2 u(x_1) + 2 u(x_2) + \dotsb + 2 u(x_{n-1}) + u(x_n) \bigr) \]

Propriété (non démontrée, cf (Quarteroni et al., 2007)) \[ \forall k \in {\mathbb{{N}}}, \; J(h) = I + \alpha_1 h^2 + \alpha_2 h^4 + \cdots + \alpha_k h^{2k} + \mathcal O(h^{2k+2}) \]

Fonctions $J_k$

J₁(h) = (4J(h) - J(2h))/3
J₂(h) = (16J₁(h) - J₁(2h))/15
J₃(h) = (64J₂(h) - J₂(2h))/63

Exercice

Montrer que si $J(h)$ correspond à la méthode trapézoïdale composite, $J_1(h)$ correspond à la méthode de Simpson composite.

Exemple

\[ I=\int_{0}^{\frac{\pi}{2}} \cos x\, \mathrm dx \]

Code

import Plots
using LaTeXStrings
u(x) = cos(x) ; a, b = 0, π/2 ; I = 1.0
function composite_trapezoidal(u, a, b, h)
    n = Int((b-a)/h)
    x = LinRange(a, b, n + 1)
    ux = u.(x)
    return (h/2) * sum([ux[1]; ux[end]; 2ux[2:end-1]])
end
J(h) = composite_trapezoidal(u, a, b, h)
J₁(h) = (4J(h) - J(2h))/3
J₂(h) = (16J₁(h) - J₁(2h))/15
J₃(h) = (64J₂(h) - J₂(2h))/63
hs = (b-a) * 2.0.^(-8:-3)
colors = Plots.palette(:default)
Plots.plot(hs, 0.1hs.^2, color=colors[1])
Plots.scatter!(hs, abs.(1 .- J.(hs)), color=colors[1], label=L"|I - J(h)|~~~O(h^2)")
Plots.plot!(hs, .005hs.^4, color=colors[2])
Plots.scatter!(hs, abs.(1 .- J₁.(hs)), color=colors[2], label=L"|I - J_1(h)|~~~O(h^4)")
Plots.plot!(hs, .0025hs.^6, color=colors[3])
Plots.scatter!(hs, abs.(1 .- J₂.(hs)), color=colors[3], label=L"|I - J_2(h)|~~~O(h^6)")
Plots.plot!(hs, .003hs.^8, color=colors[4])
Plots.scatter!(hs, abs.(1 .-J₃.(hs)), color=colors[4], label=L"|I - J_3(h)|~~~O(h^8)")
Plots.plot!(xlabel="h", ylabel="erreur relative", xaxis=:log, yaxis=:log, size=(900,450), color=colors[1])
Plots.plot!(bottom_margin=10Plots.mm, left_margin=4Plots.mm, legend=:bottomright)
Plots.xticks!((b-a)*2.0.^(-8:-3), [latexstring("\\frac{b-a}{2^{$(i)}}") for i in 8:-1:3])

Méthode avec nœuds non-équidistants : Gauss-Legendre

Intégration de Newton-Cotes ⇒ nœuds fixés, seuls degrés de liberté = poids
Pour atteindre un degré de précision de $2n+1$ avec Newton-Cotes, il faut $2n+2$ nœuds
Si on optimise aussi les nœuds, on peut espérer un degré de précision de $2n+1$ avec $n+1$ nœuds et $n+1$ poids soit $2n+2$ d.d.l. \[ \sum_{i=0}^{n} w_i x_i^{d} = \int_{-1}^{1} x^d \, \mathrm dx, \qquad d = 0, \dotsc, 2n+1. \]
Produit scalaire \[ \begin{array}{rccl} ⟨•,•⟩ :& \mathcal{C}([-1,1])×\mathcal{C}([-1,1])&→&{\mathbb{{R}}}\\ &(f,g)&↦&⟨f,g⟩=\int_{x=-1}^1 f(x) g(x) \mathrm dx \end{array} \]
Gram-Schmidt sur les polynômes $1$, $X$, $X^2$, … ⇒ polynômes orthogonaux de Legendre $L_n$ (à une constante multiplicative près)
Racines $(r_i)_{(0≤i≤n)}$ de $L_{n+1}$ (scindé à racines simples dans $]-1,1[$) ⇒ les formes linéaires $P↦P(r_i)$ forment une base de ${\mathbb{{R}}}_n[X]^*$
Il existe donc $n+1$ poids tels que $\int_{-1}^{1} P(x) \, \mathrm dx=\sum_{i=0}^{n} w_i P(r_i)$ pour tout $P \in{\mathbb{{R}}}_n[X]$
La formule reste vraie si $P \in{\mathbb{{R}}}_{2n+1}[X]$

Extensions de l’intégration gaussienne

Intégrale avec pondération $w(x)$ \[ \int_{-1}^{1} P(x) {\color{red}w(x)} \, \mathrm dx=\sum_{i=0}^{n} w_i P(r_i) \qquad \forall P \in{\mathbb{{R}}}_{2N+1}[X] \]

Démarche identique à Gauss-Legendre avec le produit scalaire \[ \begin{array}{rccl} ⟨•,•⟩ :& \mathcal{C}([a,b])×\mathcal{C}([a,b])&→&{\mathbb{{R}}}\\ &(f,g)&↦&⟨f,g⟩=\int_{x=-1}^1 f(x) g(x) {\color{red}w(x)} \mathrm dx \end{array} \]

Domaine d’intégration	Pondération $w(x)$	Famille de polynômes orthogonaux
$[–1, 1]$	1	Legendre
$]–1, 1[$	$(1-x)^α (1+x)^β$, $α,β>-1$	Jacobi
$]–1, 1[$	$\frac{1}{\sqrt{1-x^2}}$	Tchebychev (premier type)
$]–1, 1[$	$\sqrt{1-x^2}$	Tchebychev (second type)
${\mathbb{{R}}}^+$	$\mathop{\mathrm{e}}^{-x}$	Laguerre
${\mathbb{{R}}}$	$\mathop{\mathrm{e}}^{-x^2}$	Hermite

References

Quarteroni, A., Sacco, R., Saleri, F., 2007. Numerical mathematics, 2nd ed. ed, Texts in applied mathematics. Springer, Berlin ; New York.

Domaine d’intégration	Pondération \(w(x)\)	Famille de polynômes orthogonaux
\([–1, 1]\)	1	Legendre
\(]–1, 1[\)	\((1-x)^α (1+x)^β\), \(α,β>-1\)	Jacobi
\(]–1, 1[\)	\(\frac{1}{\sqrt{1-x^2}}\)	Tchebychev (premier type)
\(]–1, 1[\)	\(\sqrt{1-x^2}\)	Tchebychev (second type)
\({\mathbb{{R}}}^+\)	\(\mathop{\mathrm{e}}^{-x}\)	Laguerre
\({\mathbb{{R}}}\)	\(\mathop{\mathrm{e}}^{-x^2}\)	Hermite

Méthode	Formule d’intégration
Trapézoïdal	\(\widehat I = u(-1) + u(1)\)
Simpson	\(\widehat I = \frac{1}{3} u(-1) + \frac{4}{3} u(0) + \frac{1}{3} u(1)\)
Boole	\(\widehat I = \frac{7}{45} u(-1) + \frac{32}{45} u\left(-\frac{1}{2}\right) + \frac{4}{15} u(0) + \frac{32}{45} u\left(\frac{1}{2}\right) + \frac{7}{45} u(1)\)