Série 6 : Interpolation, Splines et Moindres carrés

1. Interpolation

Exercice 1

On considère la fonction \(f(x)=3 x^2 - x^3 - 2\) définie sur l’intervalle [0,2]. Calculer le polynôme de Lagrange \(\Pi_2 f(x)\) de degré 2 interpolant la fonction \(f\) aux nœuds \(x_0=0, x_1=1, x_2=2\) (Sugg. : calculer d’abord les polynômes \(\varphi_1, \varphi_2, \varphi_3\) de la base de Lagrange associée aux nœuds donnés; puis, exprimer \(\Pi_2 f(x)\) comme combinaison linéaire des \(\varphi_i\)).
Calculer maintenant le polynôme de Lagrange \(\Pi_3 f(x)\) de degré 3 interpolant la fonction \(f\) aux nœuds \(x_0=0, x_1=e^{-\sqrt{2}}, x_2=3 e^{-\sqrt{1 / 2}}\) et \(x_3=2\).
On connait les valeurs d’une fonction \(g(x)\) aux points \(x=0, 1, 2, 3\) : \(g(0)=3, g(1)=1, g(2)=0, g(3)=6\). Estimer sa valeur en \(x=1.5\) en interpolant \(g\) par un polynôme de degré 3 aux points \(x=0, 1, 2, 3\).

Solution

a) On peut calculer le polynôme \(\Pi_2\) interpolant la fonction \(f\) en utilisant les polynômes de la base de Lagrange associée aux nœuds \(x_0, x_1, x_2\). Polynômes de la base de Lagrange :

\[\begin{aligned} & \varphi_0(x)=\frac{\left(x-x_1\right)\left(x-x_2\right)}{\left(x_0-x_1\right)\left(x_0-x_2\right)}=\frac{1}{2}\left(x^2-3 x+2\right) \\ & \varphi_1(x)=\frac{\left(x-x_0\right)\left(x-x_2\right)}{\left(x_1-x_0\right)\left(x_1-x_2\right)}=2 x-x^2 \\ & \varphi_2(x)=\frac{\left(x-x_0\right)\left(x-x_1\right)}{\left(x_2-x_0\right)\left(x_2-x_1\right)}=\frac{1}{2}\left(x^2-x\right) \end{aligned}\]

Le polynôme cherché est donc

\[\Pi_2 f(x)=f(0) \varphi_0(x)+f(1) \varphi_1(x)+f(2) \varphi_2(x)=2(x-1)\]

b) Il suffit d’observer que \(f(x)\) est un polynôme de degré trois et donc on a \(\Pi_3 f(x)=f(x)\). c) On peut construire les polynômes de la base de Lagrange et utiliser la même démarche que au point a). On note \(x_0=0, x_1=1, x_2=2\) et \(x_3=3\). On a alors

\[\begin{aligned} \varphi_0(x) & =-\frac{1}{6}(x-1)(x-2)(x-3) & \varphi_1(x) & =\frac{1}{2} x(x-2)(x-3) \\ \varphi_2(x) & =-\frac{1}{2} x(x-1)(x-3) & \varphi_3(x) & =\frac{1}{6} x(x-1)(x-2) \end{aligned}\]

Finalement, le polynôme interpolant est \(\Pi_3 g(x)=x^3-\frac{5}{2} x^2-\frac{1}{2} x+3\) et sa valeur en \(x=1.5\) est \(\Pi_3 g(1.5)=0\).

Notons que, afin de calculer \(\Pi_3 g(x)\), on pourrait imposer directement les conditions d’interpolation. Soit \(\Pi_3 g(x)=a_0+a_1 x+a_2 x^2+a_3 x^3\), avec \(a_j(j=0, \ldots, 3)\) coefficients inconnus. Les conditions d’interpolation \(\Pi_3 g\left(x_i\right)=g\left(x_i\right)(i=0, \ldots, 3)\) nous donnent le système \(4 \times 4\) suivant :

\[\left\{ \begin{matrix} a_0 & & & & & & & = & 3 \\ a_0 & + & a_1 & + & a_2 & + & a_3 & = & 1 \\ a_0 & + & 2a_1 & + & 4a_2 & + & 8a_3 & = & 0 \\ a_0 & + & 3a_1 & + & 9a_2 & + & 27a_3 & = & 6 \\ \end{matrix} \right.\]

On voit donc que cette méthode pour calculer \(\Pi_3 g(x)\) est beaucoup plus coûteuse que celle basée sur les polynômes \(\varphi_i\).

Exercice 2

On se donne la fonction :

\[f(x)=\sin \left(\frac{x}{3}\right), \quad x \in I=[0,1] .\]

Soit \(\Pi_n f\) le polynôme interpolant la fonction aux nœuds \(x_0, x_1, \ldots, x_n\) équidistribués. Estimer l’erreur d’interpolation \(E_n(f) = \max_{x \in I} |f(x) - \Pi_n f(x)|\) sur l’intervalle \(I\), en fonction du degré \(n\) du polynôme et étudier son comportement quand \(n \rightarrow \infty\).
Trouver le nombre minimal de nœuds équirépartis pour que \(E_n(f) \leq 10^{-4}\). (Sugg. : essayer pour \(n=1, 2, 3, \ldots\))

Solution

a) En général, pour estimer l’erreur d’interpolation d’une fonction continue par un polynôme de degré \(n\) dans le cas de nœuds équidistribués sur l’intervalle \([a, b\)], on peut utiliser l’inégalité

\[E_n(f)=\max _{a \leq x \leq b}\left|f(x)-\Pi_n f(x)\right| \leq \frac{1}{4(n+1)}\left(\frac{b-a}{n}\right)^{(n+1)} \max _{a \leq x \leq b}\left|f^{(n+1)}(x)\right|,\]

où \((b-a)\) est la longueur de l’intervalle d’interpolation. Pour la fonction \(f(x)=\sin \left(\frac{x}{3}\right)\), on a que

\[f^{(1)}=\frac{1}{3} \cos \left(\frac{x}{3}\right), \quad f^{(2)}=-\frac{1}{9} \sin \left(\frac{x}{3}\right), \quad \ldots\]

et donc on obtient immédiatement

\[\max _{0 \leq x \leq 1}\left|f^{(n+1)}(x)\right| \leq \frac{1}{3^{(n+1)}} .\]

Dans ce cas, on a donc l’estimation suivante :

\[E_n(f)=\max _{0 \leq x \leq 1}\left|f(x)-\Pi_n f(x)\right| \leq \frac{1}{4(n+1)(3 n)^{(n+1)}}\]

qui tend vers zéro quand \(n \rightarrow \infty\).

b) On a la condition suivante à satisfaire :

\[\frac{1}{4(n+1)(3 n)^{(n+1)}} \leq \frac{1}{10^4}\]

D’abord, on peut essayer pour \(n=1\) et \(n=2\), mais on voit que l’inégalité (3) n’est pas vérifié. Finalement, pour \(n=3\), on trouve

\[\frac{1}{4 \cdot 4 \cdot 9^4} \sim 9.526 \cdot 10^{-6} < 10^{-4} \text {. }\]

Exercice 5

On veut trouver une formule d’approximation pour l’intégrale

\[I=\int_0^1 f(x) \mathrm{d} x\]

où \(f\) est une fonction quelconque.

Soit \(\Pi_1^H f\) le polynôme linéaire par morceaux interpolant \(f\) aux nœuds \(0=x_0, x_1, \ldots, x_n=\) 1 équirépartis. On a donc \(H=1 / n\) et \(x_i=i H\). Le graphe sur la page suivante montre un exemple avec \(n=4\). On considère l’intégrale de cet interpolant comme approximation de \(I\) :

\[\tilde{I}=\int_0^1 \Pi_1^H f(x) \mathrm{d} x \approx I .\]

Exprimer \(\tilde{I}\) en fonction \(\operatorname{des} f\left(x_i\right), i=0 \ldots n\).
En utilisant l’estimation pour l’erreur de l’interpolant, donner une estimation pour l’erreur de l’integrale

\[E=|I-\tilde{I}| .\]

Suggestion : Utiliser l’inégalité suivante :

\[\left|\int_a^b g(x) \mathrm{d} x\right| \leq \int_a^b|g(x)| \mathrm{d} x\]
Soit \(\Pi_0^H f\) le polynôme constant par morceaux interpolant \(f\) aux nœuds \(0=x_0, x_1, \ldots, x_n=\) 1 équirépartis. On a donc \(H=1 / n\) et \(x_i=i H\). Le polynôme constant par morceaux est défini comme suit :

\[\Pi_0^H f(x)=f\left(x_i\right) \quad \text { si }\left|x-x_i\right| < H / 2\]

On considère l’intégrale de cet interpolant comme approximation de \(I\) :

\[\hat{I}=\int_0^1 \Pi_0^H f(x) \mathrm{d} x \approx I\]

Exprimer \(\hat{I}\) en fonction des \(f\left(x_i\right), i=0 \ldots n\). Comparer avec l’expression trouvée à la première question.

Solution

Exercice 6

On considère la fonction \(\sin (\mathrm{x})\) définie sur l’intervalle \([0,3 \pi\)].

En utilisant les commandes polyfit et polyval, calculer le polynôme d’interpolation \(\Pi_n \sin (x)\) pour une distribution de nœuds uniforme dans \([0,3 \pi\)], dans les cas \(n=1, \ldots, 5\) ( \(n\) étant le degré du polynôme). Comparez graphiquement le résultat avec la fonction donnée.
Evaluer l’erreur \(E_n(\sin )=\max _{x \in[0,3 \pi\)}\left|\sin (x)-\Pi_n \sin (x)\right|] et visualiser le graphe de \(E_n\) en fonction de \(n\).
En observant que \(\left|\sin ^{(n)}(x)\right| \leq 1, \forall n\), comparer l’erreur obtenue au point b) avec l’estimation théorique donnée au cours :

\[E_n(f) \leq \frac{1}{4(n+1)} \left(\frac{b-a}{n}\right)^{n+1} \max _{x \in[a, b]}\left|f^{(n+1)}(x)\right| .\]
calculer le polynôme linéaire par intervalles \(\Pi_1^H \sin (x)\), sur \(N\) sous-intervalles de longueur \(H=\frac{3 \pi}{N}\). Considérer \(N=2,4,6,8,10\) et comparer graphiquement les résultats obtenus avec la fonction donnée.
Evaluer l’erreur \(E_1^H(\sin )=\max _{x \in[0,3 \pi]}\left|\sin (x)-\Pi_1^H \sin (x)\right|\) et visualiser le graphe de \(E_1^H\) en fonction du nombre d’intervalles \(N\).

Solution

2. Moindres carrés

Exercice 3

On considère la fonction

\[f(x)=\frac{1}{1+x}, \quad x \in I=[0,1]\]

Déterminer le nombre minimal \(N\) d’intervalles uniformes (de longueur \(h=1 / N\)) pour que le polynôme linéaire par morceaux qui interpole la fonction par intervalles donne une erreur \(\leq 10^{-5}\).
On considère les points \(x_0=0, x_1=1 / 2, x_2=1\). Trouver les coefficients \(a_0\) et \(a_1\) du polynôme \(p(x)=a_0 + a_1 x\) de degré 1 interpolant \(f\) aux points \(x_i (i=0, 1, 2)\) au sens des moindres carrés (droite de régression), en minimisant la fonction "distance polynôme-données" suivante :

\[\Phi\left(a_0, a_1\right) = \sum_{i=0}^2 \left[f\left(x_i\right) - p\left(x_i\right)\right]^2 = \sum_{i=0}^2 \left[f\left(x_i\right) - a_0 - a_1 x_i\right]^2 .\]

Solution

a) On utilise l’estimation de l’erreur qui suit :

\[E_n(f)=\max _{a \leq x \leq b}\left|f(x)-\Pi_n f(x)\right| \leq \frac{1}{4(n+1)}\left(\frac{b-a}{n}\right)^{(n+1)} \max _{a \leq x \leq b}\left|f^{(n+1)}(x)\right|,\]

Soit \(N\) le nombre de sous-intervalles de \(I=[0,1\)] de longueur \(h=1 / N\). On rappelle que pour l’erreur d’interpolation linéaire par intervalles, on peut fournir l’estimation suivante :

\[\max _{0 \leq x \leq 1}\left|\Pi_h^1 f(x)-f(x)\right| \leq \frac{\max _{0 \leq x \leq 1}\left|f^{(2)}\right|}{8} h^2 .\]

On a

\[\max _{0 \leq x \leq 1}\left|f^{(2)}\right|=\max _{0 \leq x \leq 1}\left|\frac{2}{(1+x)^3}\right| \leq 2\]

Pour satisfaire la condition demandée par l’exercice on doit donc prendre

\[h^2 \leq 4 \cdot 10^{-5}\]

et donc \(N \geq \sqrt{0.25 \cdot 10^5} \sim 159\).

b) On calcule les dérivées de la fonction \(\Phi\left(a_0, a_1\right)\) :

\[\begin{aligned} \frac{\partial \Phi}{\partial a_0} & =-2 \sum_{i=0}^2\left[f\left(x_i\right)-\left(a_0+a_1 x_i\right)\right]=-2\left(\sum_{i=0}^2 f\left(x_i\right)-3 a_0-a_1 \sum_{i=0}^2 x_i\right) \\ & =-2\left(\frac{13}{6}-3 a_0-\frac{3}{2} a_1\right) ; \\ \frac{\partial \Phi}{\partial a_1} & =-2 \sum_{i=0}^2 x_i\left[f\left(x_i\right)-\left(a_0+a_1 x_i\right)\right]=-2\left(\sum_{i=0}^2 x_i f\left(x_i\right)-a_0 \sum_{i=0}^2 x_i-a_1 \sum_{i=0}^2 x_i^2\right) \\ & =-2\left(\frac{5}{6}-\frac{3}{2} a_0-\frac{5}{4} a_1\right) . \end{aligned}\]

Donc les coefficients \(a_0\) et \(a_1\) du polynôme sont les solutions du système linéaire que l’on obtient en imposant les conditions de minimisation de la "distance polynôme-données", c’est-à-dire \(\frac{\partial \Phi}{\partial a_0}=0\) et \(\frac{\partial \Phi}{\partial a_1}=0\) :

\[\left\{\begin{array}{l} 3 a_0+\frac{3}{2} a_1=\frac{13}{6} \\ \frac{3}{2} a_0+\frac{5}{4} a_1=\frac{5}{6} \end{array} .\right.\]

La solution est

\[a_0=-\frac{1}{2}, \quad a_1=\frac{35}{36} .\]

Le polynôme cherché est donc \(p(x)=-\frac{35}{36}-\frac{1}{2} x\).

On peut remarquer que la matrice \(A=\left(\begin{array}{cc}3 & 3 / 2 \\ 3 / 2 & 5 / 4\end{array}\right)\) du système (4) est bien la matrice 3 \(B^{\top} B\) des équations normales (voir notes du cours), \(B\) étant donnée dans notre cas par

\[B=\left[\begin{array}{cc} x_0^0 & x_0^1 \\ x_1^0 & x_1^1 \\ x_2^0 & x_2^1 \end{array}\right]=\left[\begin{array}{cc} 1 & 0 \\ 1 & 1 / 2 \\ 1 & 1 \end{array}\right]\]

Exercice 4

On considère la sortie numérique d’un instrument dont la mesure est affectée par un bruit de fond. On indique cette donnée par \(f\) et on la définit par les commandes :

import numpy as np
x = np.linspace(0, 1, 10)
f = 10 * x + 2 * np.random.rand(x.size)

où np.random.rand est une fonction Python à valeurs aléatoires dans \(\left(0, 1\right)\). On veut effectuer un "lissage" de la donnée \(f\), en utilisant la méthode d’approximation au sens des moindres carrés.

A l’aide des commandes polyfit et polyval, calculer le polynôme de degré 9 interpolant les données \(\mathrm{f}\), et les polynômes approximant de degré 2,3,4 au sens des moindres carrés.
Comparer graphiquement les polynômes calculés avec la donnée initiale.
Utiliser le polynôme interpolant et les polynômes aux moindres carrés pour extrapoler la valeur de \(f\) au point \(x=1.2\). Comparer les valeurs obtenues avec la valeur "exacte" (sans bruit) de f, notamment 10*(1.2)^2. Quelles conclusions peut-on tirer?

Si l’on répète l’exercice, le vecteur rand et la donnée \(f\) changent. Observer que les changements des polynômes d’approximation au sens des moindres carrés sont petits par rapport aux changements des données (propriété de stabilité).

Solution

a) En procédant comme dans les exercices précédents, on trouve le code Python suivant pour calculer les polynômes interpolants de différents degrés :

import numpy as np
from numpy.polynomial.polynomial import Polynomial
import plotly.graph_objects as go

# Données
x = np.linspace(0, 1, 10)
f = 10 * x**2 + 2 * (np.random.rand(x.size) - 0.5)

# Calcul des polynômes interpolants
c2 = np.polyfit(x, f, 2)
c3 = np.polyfit(x, f, 3)
c4 = np.polyfit(x, f, 4)
c9 = np.polyfit(x, f, 9)

# Noms des polynômes pour la légende
names = ['n=2', 'n=3', 'n=4', 'n=9', 'data']
colors = ['blue', 'green', 'black', 'red', 'black']

Results

b) Tracé des Graphes

Voici le code Python utilisant Plotly pour tracer les graphes des polynômes interpolants :

# Points pour l'évaluation des polynômes
x100 = np.linspace(0, 1, 100)

# Tracé avec Plotly
fig = go.Figure()
fig.add_trace(go.Scatter(x=x100, y=np.polyval(c2, x100), mode='lines', name='n=2', line=dict(color='blue')))
fig.add_trace(go.Scatter(x=x100, y=np.polyval(c3, x100), mode='lines', name='n=3', line=dict(color='green')))
fig.add_trace(go.Scatter(x=x100, y=np.polyval(c4, x100), mode='lines', name='n=4', line=dict(color='black')))
fig.add_trace(go.Scatter(x=x100, y=np.polyval(c9, x100), mode='lines', name='n=9', line=dict(color='red')))
fig.add_trace(go.Scatter(x=x, y=f, mode='markers', name='data', marker=dict(color='black')))
fig.update_layout(title='Polynômes interpolants de degré 2 et 9 pour les données f')
fig.show()

Results

b) Extrapolation

On calcule la valeur de f au point x=1.2 :

x = 1.2
f_exact = 10 * x**2
estimations = [np.polyval(c, x) for c in [c2, c3, c4, c9]]
print(f"Valeur exacte : {f_exact}")
print(f"Estimations : {estimations}")

Results

Valeur exacte : 14.399999999999999
Estimations : [np.float64(11.88640889014812), np.float64(12.20174762673811), np.float64(13.239912549778104), np.float64(-2034.7378534462014)]

On voit donc que seules les approximations aux moindres carrés nous donnent des extrapolations fiables (respectivement 15.310, 15.69, 16.54, la valeur "exacte" sans bruit étant 14.4); le polynôme d’interpolation de Lagrange, par contre, lorsque on l’utilise pour extrapoler la valeur correspondante à x=1.2, donne une prévision complètement erronée.

# Points pour l'évaluation des polynômes
x100 = np.linspace(0, 1.3, 100)

# Tracé avec Plotly
fig = go.Figure()
fig.add_trace(go.Scatter(x=x100, y=np.polyval(c2, x100), mode='lines', name='n=2', line=dict(color='blue')))
fig.add_trace(go.Scatter(x=x100, y=np.polyval(c3, x100), mode='lines', name='n=3', line=dict(color='green')))
fig.add_trace(go.Scatter(x=x100, y=np.polyval(c4, x100), mode='lines', name='n=4', line=dict(color='black')))
fig.add_trace(go.Scatter(x=x100, y=np.polyval(c9, x100), mode='lines', name='n=9', line=dict(color='red')))
fig.add_trace(go.Scatter(x=x, y=f, mode='markers', name='data', marker=dict(color='black')))
fig.update_layout(title='Polynômes interpolants de degré 2 et 9 pour les données f')
fig.show()

Results