Interpolation et approximation

Cours ENPC - Pratique du calcul scientifique

Cadre général de l’interpolation

Soient

un segment $[a,b]$ de ${\mathbb{{R}}}$
$n+1$ points distincts $a=x_0< x_1< \dotsc < x_n=b$
$n+1$ valeurs $u_0, u_1, \dotsc, u_n$
un sous-espace $\mathop{\mathrm{Vect}}(φ_0,φ_1,\dotsc,φ_n)$ de l’e.v. des fonctions continues sur $[a,b]$

On cherche à identifier un élément de $\mathop{\mathrm{Vect}}(φ_0,φ_1,\dotsc,φ_n)$ soit \[ \fcolorbox{bleu_C}{bleu_light_C}{$ \widehat u(x) = α_0 φ_0(x) + \dotsb + α_n φ_n(x) $} \] tel que \[ \fcolorbox{bleu_C}{bleu_light_C}{$ \forall i \in \{0, \dotsc, n\}, \qquad \widehat u(x_i) = u_i $} \] \[ \fcolorbox{bleu_C}{bleu_light_C}{$ \mathsf A \mathbf{\boldsymbol{\alpha }}:= \begin{pmatrix} \varphi_0(x_0) & \varphi_1(x_0) & \dots & \varphi_n(x_0) \\ \varphi_0(x_1) & \varphi_1(x_1) & \dots & \varphi_n(x_1) \\ \vdots & \vdots & & \vdots \\ \varphi_0(x_n) & \varphi_1(x_n) & \dots & \varphi_n(x_n) \end{pmatrix} \begin{pmatrix} \alpha_0 \\ \alpha_1 \\ \vdots \\ \alpha_n \end{pmatrix} = \begin{pmatrix} u_0 \\ u_1 \\ \vdots \\ u_n \end{pmatrix} := \mathbf{\boldsymbol{b}} $} \]

Nœuds équidistants

\[x_k=a+(b-a)\frac{k}{n}\]

Nœuds de Tchebychev

\[x_k=a + (b-a) \frac{1 - \cos \left(\pi \frac{k + \frac{1}{2}}{n+1} \right)}{2}\]

Familles de polynômes

Base canonique

\[ φ_i(x)=x^i \quad⇒\quad \mathsf A= \begin{pmatrix} 1 & x_0 & \dots & x_0^n \\ 1 & x_1 & \dots & x_1^n \\ \vdots & \vdots & & \vdots \\ 1 & x_n & \dots & x_n^n \end{pmatrix} \]

Code

x_interp = LinRange(0, 1, 200)
pl = Plots.plot(size=(450,300), bottom_margin=4Plots.mm)
for d in 0:10
    p(x) = x^d
    Plots.plot!(p, 0, 1, )
end
Plots.plot!(xlims=(0, 1), xlabel="x",
    legend=false)

Polynômes de Lagrange

\[ φ_{i}(x) = \prod_{\substack{j = 0\\ j ≠ i}}^{n} \frac {x - x_j} {x_i - x_j} \quad⇒\quad \mathsf A=\mathsf I \]

Code

xs = LinRange(0, 1, 11)
x_interp = LinRange(0, 1, 200)
pl = Plots.plot(size=(450,300), bottom_margin=4Plots.mm)
for d in eachindex(xs)
    p(x) = prod(x .- xs[1:end .!= d]) / prod(xs[d] .- xs[1:end .!= d])
    Plots.plot!(p, extrema(xs)...)
end
Plots.scatter!(xs, 0*xs .+ 1, )
Plots.plot!(xlims=(0, 1), xlabel="x",
    title="$(length(xs)) nœuds équidistants",
    legend=false)

Polynômes échelonnés (Gregory-Newton généralisé)

\[ φ_{i}(x) = (x-x_0) (x-x_1) (x-x_2) \dots (x-x_{i-1}) \]

$\quad⇒\quad \mathsf A$ triangulaire inférieure

Code

xs = LinRange(0, 1, 11)
x_interp = LinRange(0, 1, 200)
pl = Plots.plot(size=(450,300), bottom_margin=4Plots.mm)
for d in eachindex(xs)
    p(x) = prod(x .- xs[1:d-1])
    Plots.plot!(p, extrema(xs)...)
end
Plots.scatter!(xs, zero(xs), )
Plots.plot!(xlims=(0, 1), xlabel="x",
    title="$(length(xs)) nœuds équidistants",
    legend=false)

Contrôle de l’erreur

Théorème

$u\colon [a, b] \to {\mathbb{{R}}}$ est une fonction dans $\mathcal{C}^{n+1}([a, b])$
$a=x_0< x_1< \dotsc < x_n=b$ sont $n+1$ nœuds distincts
$\widehat u$ est un polynôme de degré au plus $n$, interpolateur de $u$ aux points $x_0, x_1, \dotsc, x_n$, i.e. $\widehat u(x_i) = u(x_i)$ pour tout $i ∈ \{0,\dotsc,n\}$

Alors on a $\quad∀\, x ∈ [a,b],\quad ∃\, ξ=ξ(x) ∈ [a,b]$ \[ \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) $} \]

Exercice : démontrer le théorème

Examiner le cas où $x$ est l’un des $x_i$.
Posant $ω_n(x) = \prod_{i=0}^n (x - x_i)$ et $g(t) = e_n(t) ω_n(x) - e_n(x) ω_n(t)$ pour un $x$ donné différent des $x_i$, montrer que $g^{(k)}(t)$ avec $0≤k≤n+1$ admet $n+2-k$ racines distinctes dans $[a,b]$.
Conclure.
Que dire du cas où $u$ est un polynôme de degré au plus $n$ ?
En supposant que $u$ est dans $\mathcal{C}^{∞}([a, b])$, fait-on systématiquement tendre l’erreur vers $0$ lorsque le nombre de nœuds tend vers l’infini ?

Corollaire du théorème (mêmes hypothèses)

On pose $C_{n+1} = \sup_{x ∈ [a, b]} \left\lvert u^{(n+1)}(x) \right\rvert$ et $h = \max_{i ∈ \{0, \dotsc, n-1\}} \lvert x_{i+1} - x_i\rvert$

Alors on montre que $E_n := \sup_{x ∈ [a, b]} \bigl\lvert e_n(x) \bigr\rvert \leqslant\frac{C_{n+1}}{4(n+1)} h^{n+1}$

Indications :
- si $x ∈ [x_{i},x_{i+1}]$, $(x-x_{i})(x_{i+1}-x)=\left(\frac{x_{i+1}-x_{i}}{2}\right)^2-\left(x-\frac{x_{i}+x_{i+1}}{2}\right)^2 ≤ \frac{h^2}{4}$
- $\lvert ω_n(x)\rvert ≤ \frac{h^2}{4} × 2h × 3h × 4h × \dotsb × nh$

Remarques

$h$ dépend de $n$ et du choix des nœuds (en $1/n$ pour nœuds équidistants)
$C_{n}$ ne peut a priori être contrôlé
Dans certains cas, $C_{n}$ ne croît pas trop vite avec $n$ de sorte que $E_n\underset{n→∞}{→}0$ (ex : sinus)
Contre-exemple avec la fonction de Runge $u(x)=\frac{1}{1+25x^2}$ pour laquelle le majorant de $E_n$ tend vers l’infini si les nœuds sont équidistants
Optimiser l’interpolation ? → optimisation des nœuds ou interpolation par morceaux

Fonction sinus avec nœuds équidistants

\[ \fcolorbox{orange_C}{orange_light_C}{$ u(x)=\sin{x} $} \] \[ \fcolorbox{pistache_C}{pistache_light_C}{$ x_k=a + (b-a) \frac{k}{n} \quad (0 ≤ k ≤ n) $} \]

Fonction de Runge avec nœuds équidistants

\[ \fcolorbox{orange_C}{orange_light_C}{$ u(x)=\frac{1}{1+25x^2} $} \] \[ \fcolorbox{pistache_C}{pistache_light_C}{$ x_k=a + (b-a) \frac{k}{n} \quad (0 ≤ k ≤ n) $} \]

Optimisation des nœuds : nœuds de Tchebychev

Contrôle du majorant de l’erreur

\[ \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) $} \]

Augmenter le nombre de nœuds $n+1$ change $u^{(n+1)}$ ⇒ difficile à contrôler
Idée : pour $n$ donné, choisir les nœuds pour minimiser $\sup_{x∈[a,b]} \lvert(x-x_0) \dotsc (x - x_n)\rvert$
Pour tout polynôme $p$ unitaire de degré $n$, on montre que $\sup_{x∈[-1,1]} \lvert p(x)\rvert ≥ \frac{1}{2^{n-1}}$
Borne $\frac{1}{2^{n-1}}$ atteinte pour $\frac{T_{n}(x)}{2^{n-1}}$ avec $T_n(x)=\cos(n\arccos x)$ (polynôme de Tchebychev)
Les zéros de $T_n$ sont donc les $x_k=-\cos (\pi \frac{k + \frac{1}{2}}{n} ) \; (0 ≤ k < n)$ (signe “-” pour ordre)
Soit $x_k=-\cos (\pi \frac{k + \frac{1}{2}}{n+1} ) \; (0 ≤ k ≤ n)$ pour $n+1$ points
Dans le cas $[a,b]$, pour $n$ points \[ \fcolorbox{pistache_C}{pistache_light_C}{$ x_k=a + (b-a) \frac{1 - \cos (\pi \frac{k + \frac{1}{2}}{n} )}{2} \quad (0 ≤ k < n) $} \]
Application à la fonction de Runge $u(x)=\frac{1}{1+25x^2}$

Interpolation par morceaux

Interpolation continue par morceaux

\[ \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)} h^{n+1}\\ &\text{avec } C_{n+1} = \sup_{x ∈ [a, b]} \left\lvert u^{(n+1)}(x) \right\rvert \text{ et } h = \max_{i ∈ \{0, \dotsc, n-1\}} \lvert x_{i+1} - x_i\rvert \end{align} $} \]

Idée pour contrôler le majorant :

découper l’intervalle $[a,b]$ avec $n+1$ nœuds, par exemple $h=\frac{b-a}{n}$
interpoler avec un polynôme de degré $m$ sur $[x_{i},x_{i+1}]$

\[ \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} $} \]
$m$ est fixé petit donc l’erreur est majorée uniformément par $C h^{m+1}$

Régularité supérieure

L’interpolation précédente n’assure que la continuité $\widehat u(x_i^-)=\widehat u(x_i^+)$
Pour améliorer la régularité, on peut imposer des conditions sur les dérivées supérieures

→ splines cubiques :
- polynôme de degré 3 sur chaque sous-intervalle : $4n$ inconnues
- interpolation aux nœuds $x_i$ ($0 ≤ i ≤ n$) : $2n$ équations
- raccord des dérivées aux nœuds $x_i$ ($0 < i < n$) : $n-1$ équations
- raccord des dérivées secondes aux nœuds $x_i$ ($0 < i < n$) : $n-1$ équations
- nombre total d’équations : $4n-2$ ⇒ il en faut 2 de plus (en général dérivées secondes nulles aux extrémités)

Approximation par moindres carrés

Position du problème

$n+1$ nœuds distincts $a=x_0< x_1< \dotsc < x_n=b$
$n+1$ valeurs $u_0, u_1, \dotsc, u_n$
un sous-espace $\mathop{\mathrm{Vect}}(φ_0,φ_1,\dotsc,φ_m)$ de l’e.v. des fonctions continues sur $[a,b]$ (en général des polynômes) mais ici $\color{red}{m<n}$

On cherche à identifier un élément de $\mathop{\mathrm{Vect}}(φ_0,φ_1,\dotsc,φ_m)$ soit \[ \fcolorbox{bleu_C}{bleu_light_C}{$ \widehat u(x) = α_0 φ_0(x) + \dotsb + α_m φ_m(x) $} \] tel que \[ \fcolorbox{bleu_C}{bleu_light_C}{$ \forall i \in \{0, \dotsc, n\}, \qquad \widehat u(x_i) \, {\color{red}{\approx}} \, u_i $} \] \[ \fcolorbox{bleu_C}{bleu_light_C}{$ \mathsf A \mathbf{\boldsymbol{\alpha }}:= \begin{pmatrix} \varphi_0(x_0) & \varphi_1(x_0) & \dots & \varphi_m(x_0) \\ \varphi_0(x_1) & \varphi_1(x_1) & \dots & \varphi_m(x_1) \\ \varphi_0(x_2) & \varphi_1(x_2) & \dots & \varphi_m(x_2) \\ \vdots & \vdots & & \vdots \\ \varphi_0(x_{n-2}) & \varphi_1(x_{n-2}) & \dots & \varphi_m(x_{n-2}) \\ \varphi_0(x_{n-1}) & \varphi_1(x_{n-1}) & \dots & \varphi_m(x_{n-1}) \\ \varphi_0(x_n) & \varphi_1(x_n) & \dots & \varphi_m(x_n) \end{pmatrix} \begin{pmatrix} \alpha_0 \\ \alpha_1 \\ \vdots \\ \alpha_m \end{pmatrix} {\color{red}{\approx}} \begin{pmatrix} u_0 \\ u_1 \\ u_2 \\ \vdots \\ u_{n-2} \\ u_{n-1} \\ u_n \end{pmatrix} =: \mathbf{\boldsymbol{b}} $} \]

Minimisation de \[ \fcolorbox{bleu_C}{bleu_light_C}{$ J(\mathbf{\boldsymbol{\alpha}})= \frac{1}{2}\,\sum_{i=0}^{n} {\lvert {\widehat u(x_i) - u_i} \rvert}^2 = \frac{1}{2}\,\sum_{i=0}^{n} \left( \sum_{j=0}^{m} \alpha_j \varphi_j(x_i) - u_i\right)^2 =\frac{1}{2}\,{\lVert {\mathsf A \mathbf{\boldsymbol{\alpha}}-\mathbf{\boldsymbol{b}}} \rVert}^2 $} \]

On cherche donc $\alpha$ tel que \[ \nabla J(\mathbf{\boldsymbol{\alpha}})=\frac{1}{2}\,\nabla \Bigl( (\mathsf A \mathbf{\boldsymbol{\alpha }}- \mathbf{\boldsymbol{b}})^T(\mathsf A \mathbf{\boldsymbol{\alpha }}- \mathbf{\boldsymbol{b}}) \Bigr)=\mathbf{\boldsymbol{0}} \] soit \[ \mathrm dJ=(\mathsf A \mathrm d\mathbf{\boldsymbol{\alpha}})^T(\mathsf A \mathbf{\boldsymbol{\alpha }}- \mathbf{\boldsymbol{b}}) =\mathrm d\mathbf{\boldsymbol{\alpha}}^T\mathsf A^T(\mathsf A \mathbf{\boldsymbol{\alpha }}- \mathbf{\boldsymbol{b}}) \Rightarrow \nabla J(\mathbf{\boldsymbol{\alpha}})=\mathsf A^T(\mathsf A \mathbf{\boldsymbol{\alpha }}- \mathbf{\boldsymbol{b}}) =\mathbf{\boldsymbol{0}} \] ce qui donne le système à résoudre \[ \fcolorbox{bleu_C}{bleu_light_C}{$ \mathsf A^T\mathsf A \mathbf{\boldsymbol{\alpha }}= \mathsf A^T\mathbf{\boldsymbol{b}} $} \] avec $\mathsf A^T\mathsf A$ matrice $m×m$ inversible si $\mathsf A$ est de rang $m$ (colonnes indépendantes)

Remarque

Sous Julia le résultat du calcul $\alpha=(\mathsf A^T\mathsf A)^{-1} \mathsf A^T\mathbf{\boldsymbol{b}}$ peut simplement être obtenu par

α = A\b

Bibliothèque `Polynomials.jl`

Dépôt GitHub https://github.com/JuliaMath/Polynomials.jl
Doc https://juliamath.github.io/Polynomials.jl/stable/

Exemple d’utilisation de fit

using Plots, Polynomials
xs = range(0, 10, length=10)
ys = @. exp(-xs)
f = fit(xs, ys) # degree = length(xs) - 1
f2 = fit(xs, ys, 2) # degree = 2

scatter(xs, ys, label="Data")
plot!(f, extrema(xs)..., label="Fit")
plot!(f2, extrema(xs)..., label="Quadratic Fit")
width, height = 750, 500 ; plot!(size = (width, height), legend = :topright)

Bibliothèque `Interpolations.jl`

Dépôt GitHub https://github.com/JuliaMath/Interpolations.jl
Doc http://juliamath.github.io/Interpolations.jl/stable/

Exemple d’utilisation ici (section Convenience Constructors)

using Interpolations, Plots
a = 1.0 ; b = 10.0 ; xᵢ = a:1.0:b # bounds and knots
F(x) = cos(x^2/9)
yᵢ = F.(xᵢ) # function application by broadcasting
itp_linear = linear_interpolation(xᵢ,yᵢ) ; itp_cubic = cubic_spline_interpolation(xᵢ,yᵢ)
F_linear(x) = itp_linear(x) ; F_cubic(x) = itp_cubic(x)
x_new = a:0.1:b # smoother interval, necessary for cubic spline
scatter(xᵢ, yᵢ, markersize=10,label="Data points")
plot!(F_linear, x_new, w=4,label="Linear interpolation")
plot!(F_cubic, x_new, linestyle=:dash, w=4, label="Cubic Spline interpolation")
plot!(F, x_new, linestyle=:dot, w=4, label="Initial function")
width, height = 750, 500 ; plot!(size = (width, height), legend = :bottomleft)

Interpolation et approximation

Cadre général de l’interpolation

Familles de polynômes

Contrôle de l’erreur

Fonction sinus avec nœuds équidistants

Fonction de Runge avec nœuds équidistants

Optimisation des nœuds : nœuds de Tchebychev

Interpolation par morceaux

Approximation par moindres carrés

Bibliothèque Polynomials.jl

Bibliothèque Interpolations.jl

Bibliothèque `Polynomials.jl`

Bibliothèque `Interpolations.jl`