Interpolation

Interpolation de Lagrange

Pour illustrer

On considère une suite de points \((x_0,\dots x_n)\) et une suite de valeurs \((y_0,\dots y_n)\).

1. Donner une valeur à n, et utiliser seq et rand pour créer deux listes lx et ly de \(n+1\) éléments aléatoires dans \([-5,5]\).

   solution ⇲ ⇱ solution

   n := 10;
   lx := [ seq(rand(-10,10),k=0..n) ];
   ly := [ seq(rand(-10,10),k=0..n) ];
   

Soient les fonctions \(f:x\mapsto \frac1{1+x^2}\) et \(g:x\mapsto \sin(x)\) définies sur l’intervalle \(I=[-5,5]\).

1. Ecrire des fonctions pf(n) et pg(n) qui renvoient des listes de \(n\) évaluations de \(f\) ou \(g\) en des points équirépartis de \(I\).

On rappelle que l’interpolateur de Lagrange de \((y)\) aux points \((x)\) est l’unique polynôme \(P(X)\) de degré inférieur ou égal à \(n\) tel que \(P(x_i)=y_i\).

On a la formule explicite

\[ \begin{align}\begin{aligned}\begin{split}P(X) &= \sum_{i=0}^n y_i L_i(X) \\\end{split}\\\text{où } L_i(X) &= \prod_{j\neq i} \frac{X-x_j}{x_i-x_j}\end{aligned}\end{align} \]

Calcul direct

1. Écrire une fonction LagrangeNaif(lx,ly) qui calcule \(P(X)\) à l’aide de la formule explicite.

   solution ⇲ ⇱ solution

   LagrangeNaif(lx,ly):={
     local n,P,k,j;
     n := size(lx);
     for(k:=0;k<n;k++) {
       local(Lk:=1);
       for(j:=0;j<n;j++) {
         if(j==k) continue;
         Lk *= (x-lx[j])/(lx[k]-lx[j]);
       }
       P += ly[k]*Lk;
     }
     return(P);
   }
   

2. Quelle est la complexité de l’algorithme en fonction de \(n\) ?

   solution ⇲ ⇱ solution

   L’algorithme effectue une double boucle contenant la mutiplication d’un polynôme de taille \(n\) par un monôme, soit \(O(n^3)\) opérations.

3. Calculer les interpolateurs de Lagrange \(F_n\) et \(G_n\) des fonctions \(f\) et \(g\) en les points équirépartis pf(n) et pg(n).

4. Quels sont les comportements de ces interpolateurs ? Est-ce prévisible ?

   solution ⇲ ⇱ solution

   L’interpolation de \(\sin(x)\) est excellente, tandis que celle de \(f(x)\) diverge violemment en norme infinie. En effet \(\sin(x)\) est donné par une série entière rapidement convergente, donc s’approxime très bien par des polynômes, tandis que \(f\) est une fraction rationnelle qui n’est pas proche des polynômes.

Algèbre linéaire

Notons \(V(x_0,\dots x_n)\) la matrice de passage de la base \((L_0(X),\dots L_n(X))\) des interpolateurs de Lagrange aux points \((x_i)\), à la base canonique \((1,X,\dots X^n)\).

(la matrice de passage de \(\mathcal B\) à \(\mathcal B'\) est la matrice de l’application identité \((E,\mathcal B') \to (E,\mathcal B)\)).

1. Définir cette matrice, et déterminer l’instruction xcas permettant de la créer rapidement.

2. L’utiliser pour calculer les coefficients de l’interpolateur de Lagrange \(P(X)\) dans la base canonique.

3. En déduire une fonction LagrangeVdm(lx,ly). Quelle est sa complexité ?

   solution xcas ⇲ ⇱ solution xcas

   On sait que \(P(X)\) a pour coordonnées \((y_0,\dots y_n)\) dans la base des interpolateurs \(L_i\). Pour repasser dans la base canonique, il suffit de multiplier ce vecteur par l’inverse de \(V(x_0,\dots x_n)\).

   D’où la fonction suivante, dans laquelle on a transposé l’équation et on renverse la liste à la fin pour avoir l’écriture usuelle en partant du coefficient de \(X^n\):

   LagrangeVdm(lx,ly) := {
     return revlist(ly*(vandermonde(lx)^-1));
   }
   

   La complexité est à nouveau de \(O(n^3)\) pour l’inversion matricielle.

Différences divisées

On essaie de résoudre le problème par induction en introduisant la suite de polynômes \(P_k\) tels que \(P_k\) soit l’interpolateur de Lagrange des \((y_0,\dots y_k)\) aux points \((x_0,\dots x_k)\).

On peut écrire

\[P_k(X) = P_{k-1}(X)+\lambda_{x_0,\dots x_k} \prod_{j<k}(X-x_j)\]

où \(\lambda_{x_0,\dots x_k}\) est le coefficient dominant de \(P_k\) et peut se calculer par récurrence par la formule

\[\lambda_{x_0,\dots x_k} = \frac{\lambda_{x_1,\dots x_k}-\lambda_{x_0,\dots x_{k-1}}}{x_k-x_0}\]

en partant des \(\lambda_{x_i}=y_i\).

1. Écrire une fonction lambda(ly) qui renvoie la suite des \(\lambda_{x_0,\dots x_k}\) pour \(0\leq k<n\).

2. En déduire une fonction LagrangeDivise de calcul de l’interpolant de Lagrange \(P(X)\).

   solution xcas ⇲ ⇱ solution xcas

   lx := [-3,-1,0,2,3];
   ly := [2,3,1,-1,2];
   
   difdiv(lx,ly) := {
     n := size(lx)-1;
     //fig :=
     seq(point(lx[k],ly[k]),k=0..n);
     L := ly;
     for(d:=0;d<=n;d++) {
       seq(plot(L[k],x,lx[k]-0.5,lx[k+d]+0.5),k,0,n);
       //fig
       for(k:=n-d;k>0;k--) {
         // L[k] interpolates x_k,...x_{k+d}
         // L[k-1] interpolates x_{k-1},... x_{k+d-1}
         L[k-1] := (L[k]*(x-lx[k-1])-L[k-1]*(x-lx[k+d]))/(lx[k+d]-lx[k-1]);
       }
     }
   }:;
   

Application au calcul de polynôme caractéristique

1. Écrire une matrice \(M\) aléatoire de taille 20, et calculer son polynôme caractéristique à l’aide d’une interpolation de Lagrange en les \(\det(k*I-M)\) pour \(0\leq k < 20\).

Moindres carrés

Par projection

Soient \(n\) points \((x_i)\) dans \(I\), on considère la forme quadratique \(P,Q\mapsto \sum_i P(x_i)Q(x_i)\) définie sur \(R_n[X]\).

1. Montrer que c’est un produit scalaire pour lequel les interpolateurs de Lagrange forment une base orthonormale de \(R_n[x]\).

On considère la famille de points donnée par pf(20).

1. Orthonormaliser la base canonique de \(R_n[X]\) pour ce produit scalaire. On pourra utiliser la fonction gramschmidt de xcas.

   solution xcas ⇲ ⇱ solution xcas

   n := 20;
   X := pf(20);
   pscal(P,Q) := sum(P(X[k])*Q(X[k]),k,0,20);
   can := seq(x^j,j,0,n);
   gramschmidt(can, pscal);
   

2. En déduire, pour \(m\leq n\), le polynôme \(P_m(X)\) qui minimise \(\sum_i \abs{f(x_i)-P_m(x_i)}\), et le représenter graphiquement.

3. Que se passe-t-il quand \(m\) augmente ?

Méthode directe

On veut trouver un polynôme \(P_m(X)=\sum_{j=0}^m a_jX^j\) minimisant

\[\sum_{i=1}^n \abs{y_i-P_m(x_i)}^2\]

Soit \(X\) la matrice rectangulaire \(n\times (m+1)\) de coefficients \(x_i^j\) pour \(1\leq i\leq n\) et \(0\leq j\leq m\), et \(y\) le vecteur de composantes \(y_i=f(x_i)\).

1. Montrer que le problème consiste à minimiser \(\norm{Xa-y}^2\) pour \(a\in\R^{m+1}\).

   solution ⇲ ⇱ solution

   En effet, en notant \(a=(a_j)\), \(Xa\) est le vecteur des \(P_m(x_i)\).

2. Montrer que le minimum existe et est atteint en un vecteur \(a\) tel que

   \[^tXXa = ^tXy\]
   solution ⇲ ⇱ solution

   En effet, \(Xa\to\infty\) pour \(a\to\infty\), donc par compacité l’inf est atteint. Par théorème de Rolle, la différentielle s’annule en le minimum, donc soit en \(a\) tel que \(h\mapsto 2\langle Xh, Ax-y\rangle=0\), soit l’équation de l’énoncé.

3. En déduire une fonction \(LS(m,lx,ly)\) qui renvoie une apprimation des moindres carrés de degré \(\leq m\) pour des points \((x_i,y_i)\).

Évaluation de Hörner

Note

Soit

\[P(X) = \sum_{k=0}^n a_k X^{n-k}\]

un polynôme que l’on souhaite évaluer en \(\alpha\).

Le schéma d’évaluation de Hörner consiste à écrire

\[P(\alpha)= \bigl(\dots (a_n \alpha + a_{n-1})\alpha+\dots +a_1\bigr)\alpha+a_0\]

ce qui permet de calculer la valeur \(P(\alpha)\) à l’aide de \(n\) multiplications et additions.

Remarquer la similitude avec l’exponentiation rapide.

1. Écrire une fonction Horner(P,a) qui évalue le polynôme \(P\) en \(a\). La tester sur le polynôme formel P:=sum(a[k]*X^(8-k),k,0,8).

Division par un monôme

On définit précisément la suite de Hörner

\[ \begin{align}\begin{aligned}\begin{split}u_0(\alpha) &= a_0 \\\end{split}\\u_k(\alpha) &= \alpha u_{k-1}(\alpha)+a_k \text{ pour } 1\leq k\leq n\end{aligned}\end{align} \]

de sorte que \(u_n(\alpha)=P(\alpha)\).

1. Montrer que l’on a

   \[P(X) = (X-\alpha) \sum_{k=0}^{n-1} u_k(\alpha) X^{n-1-k} + P(\alpha)\]
   solution ⇲ ⇱ solution

   Si l’on aligne les puissances de \(X\) dans le membre de droite (en mettant \(P(\alpha)=u_n(\alpha)\) comme terme d’indice \(k=n\) dans la première somme), on obtient

   \[u_0(\alpha)X^n + \sum_{k=1}^n ( u_k(\alpha) - \alpha u_{k-1}(\alpha) ) X^{n-k}\]

   ce qui vaut bien \(\sum_{k=0}^n a_k X^{n-k}=P(X)\) d’après la récurrence des \(u_k\).

2. En déduire une fonction horner_divrem(P,a) qui renvoie le quotient et le reste de la division euclidienne de \(P(X)\) par \((X-a)\).

   solution ⇲ ⇱ solution

   D’après la question précédente, le schéma de Hörner calcule les coefficients du quotient. Il suffit de les conserver.

   horner_divrem(P,a):={
     local n,k,u,Q;
     n := degree(P);
     Q := 0;
     u := coeff(P,n);
     for(k:=n-1; k>=0; k--) {
       Q := x*Q + u;
       u := a*u + coeff(P,k);
       }
     return( Q,u );
     }
   }
   

3. Vérifier les résultats en utilisant la fonction xcas quorem.

   solution ⇲ ⇱ solution

   horner_divrem(X^8-3,5)
   quorem(X^8-3,X-5)
   

Application au calcul d’une interpolation

Pour économiser les calculs, on propose la stratégie suivante :

• calculer d’abord le produit complet

  \[A(X) = \prod_{j=0}^n (X-x_i)\]

• puis calculer chacun des \(L_i\) à l’aide d’une division selon le schéma de Hörner.

• Au fur et à mesure de l’étape précédente, on peut construire \(P(X)\).

1. Écrire une fonction LagrangeHorner(X,Y) selon ce principe.

   solution ⇲ ⇱ solution

   En partant de \(A(X)\), \(L_i(X)\) vaut

   \[L_i(X) = \frac{A_i(X)}{A_i(x_i)} \text{ où } A_i(X) = \frac{A(X)}{X-x_i}\]

   On commence par modifier la fonction horner_divrem pour qu’elle calcule à la fois \(A_i(X)\) et \(A_i(x_i)\) (qui vaut simplement \(u_{n-1}\)).

   horner_divrem_A(P,a):={
     local n,k,u,Q;
     n := degree(P);
     Q := 0;
     u := 0;
     for(k:=n; k>0; k--) {
       u := a*u + coeff(P,k);
       Q := x*Q + u;
       }
     return( Q,u );
     }
   

   LagrangeHorner(lx,ly) := {
     local n,A,P,k;
     n := size(lx);
     A := product(x-lx[k],k,0,n-1);
     P := 0;
     for(k:=0;k<n;k++) {
       local Ak,ak;
       Ak,ak := horner_divrem_A(A,lx[k]);
       P += ly[k]*Ak/ak;
     }
     return(P);
   }
   

2. Quelle est la complexité obtenue ?

   solution ⇲ ⇱ solution

   quadratique