\setcounter{chapter}{2}

\chapterimage{interpol}
\chaptertext{
Premier TP d'analyse numérique (très algébrique...) assisté par Scilab autour des différentes
méthodes d'interpolation polynomiale.
}


\chapter{Interpolations}

\scilab

\section{Schéma de Horner-Ruffini-Holdred-Newton-Al-Tusi-Liu-Hui...}


\subsection{Un peu d'histoire...}



%\vspace{-6.2cm}
\setindent{0cm}
\begin{flushleft}
\begin{minipage}[t]{4.2cm}
\begin{figure}
  \includegraphics[width=2.8cm]{Ruffini}
%\caption{Paolo \textsc{Ruffini}\\(1765 - 1822)}
\end{figure}
\end{minipage}
\end{flushleft}
\setindent{\indentationcours}
\vspace{-5cm}
\begin{minipage}[t]{\linewidth}


\end{minipage}

La  méthode que nous  allons voir  porte le  nom du  britannique William
George  \textsc{Horner} (1786  - 1837)  mais  en fait  elle fut  publiée
presque  10  ans  auparavant  par  un  horloger  londonien,  Theophilius
\textsc{Holdred} et simultanément par l'italien Paolo \textsc{Ruffini} (1765 - 1822)
mais fut déjà utilisée par  \textsc{Newton} 150 ans auparavant et par le
chinois \textsc{Zhu  Shije} cinq siècles  plus tôt (vers 1300)  et avant
lui  par le Persan  \textsc{Sharaf al-Din  al-Muzaffar ibn  Muhammad ibn
al-Muzaffar al-Tusi} vers (1100) et avant lui par le
Chinois \textsc{Liu  Hui} (vers 200)  révisant un des  résultats présent
dans \textit{Les Neuf
Chapitres sur l'art mathématique} publié avant la naissance de JC...



Il faut  cependant noter que  \textsc{Ruffini} l'avait employée  en fait
comme  un moyen  de calculer  rapidement le  quotient et  le  reste d'un
polynôme par $(X-\alpha)$. C'est ce que nous allons (re)découvrir aujourd'hui...


\subsection{Complexité}

Dans toute  la suite, un  polynôme de dégré  $n$ sera représenté  par le
vecteur de ses coefficients. Par exemple, \texttt{[1 2 3]} correspond au
polynôme $1+2x+3x^2$. N'oubliez pas  enfin que Scilab commence à compter
à partir de 1:


\begin{lstlisting}[caption={}]
--> p = [1 2 3]; 
--> p(1)
 ans  =
    1.  
\end{lstlisting}


Nous avons vu le schéma  de \textsc{Horner} au TP précédent. Comparez la
complexité au  pire du calcul de $P(t)$  pour $t\in\bbk$. Vous
prendrez comme \og unité de complexité\fg{} les opérations arithmétiques
de base: \texttt{+ - *}.


\subsection{Dérivées successives}

On note $P_n=a_0+a_1X+\cdots+a_nX^n$ un polynôme de degré $n$.

Déterminez le lien entre  l'algorithme de \textsc{Horner} et la division
euclidienne  de $P_n$  par $(X-t)$.  On utilisera  par la  suite la
notation suivante:

$$P_n=(X-t)(b_nX^{n-1}+b_{n-1}X^{n-2}+\cdots+b_1)+b_0$$


Comment calculer  les dérivées successives $\widetilde{P_n^{(j)}}(t)$
en  utilisant  uniquement  des  schémas  de  \textsc{Horner}?  Quel  est
l'intérêt informatique? Voici un peu d'aide...


Considérons le  développement de $P_n$ selon  les puissances croissantes
de $(X-t)$:

$$P_n=t_n(X-t)^n+t_{n-1}(X-t)^{n-1}+\cdots+t_0$$

Montrez que $\widetilde{P_n^{(j)}}(t)j!t_j$.  Quel est le lien
entre les $t_j$ et nos schémas de \textsc{Horner}?

\subsubsection{Version récursive}

Déterminez une fonction récursive sur Scilab \texttt{der\_j(poly,j,t)} qui calcule
$\widetilde{P_n^{(j)}}(t)$ en utilisant le schéma de \textsc{Horner}.

Pour  cela, on  pourra  commencer par  modifier  légèrement la  fonction
\texttt{horner} pour récupérer le vecteur $\mathtt{[b_0 b_1 ... b_n]}$.

Créez ensuite \texttt{function y = quotient\_horner(poly,t)} qui renvoie
le quotient de \texttt{poly}  par \texttt{(X-t)} puis \texttt{function y
= reste\_horner(poly,t)}.

Déduisez-en    une    fonction    récursive   \texttt{function    y    =
jeme\_poly(poly,j,t)} qui renvoie $\widetilde{P_n^{(j)}}(t)$.


\subsubsection{Version itérative}

Prenons par exemple le  polynôme $P_4=[6,-7,1,-5,2]$. Essayez de trouver
un algorithme simple pour remplir le tableau suivant:


\begin{center}
             \begin{tabularx}{\linewidth}{IZ|ZZZZZ|ZZI}
	     \whline
Restes & 0 & 1 & 2 & 3 & 4 & $j$ & \\
\whline
 / & 6 & $-7$ & 1 & $-5$ & 2 & $a_j$ & $P_4$ \\
$-12$ & $-9$ & $-1$ & $-1$ & 2 & / & $b_j$ & $Q_3$ \\
      & & & & & &$c_j$ & $Q_2$ \\
      & & & & & &$d_j$ & $Q_1$ \\
      & & & & & &$e_j$ & $Q_0$ \\
      & & & & &  & & \\
             \whline
             \end{tabularx}
             \end{center}

En     déduire     l'écriture     de     $P_4$     sous     la     forme
$R_4+R_3(X-2)+R_2(X-2)^2+R_1(X-2)^3+R_0(X-2)^4$  et faire  la  lien avec
les dérivées successives.


Écrivez  ensuite une  fonction  qui renvoie  un  tableau similaire  puis
transformez-la un peu pour obtenir les divisions successives en $t$.


\begin{lstlisting}[caption={}]
-->jeme_poly_mat([6 -7 1 -5 2],5,2)
 ans  = 
    0.     6.  - 7.    1.  - 5.    2.  
  - 12.  - 9.  - 1.  - 1.    2.    0.  
    1.     5.    3.    2.    0.    0.  
    19.    7.    2.    0.    0.    0.  
    11.    2.    0.    0.    0.    0.  
    2.     0.    0.    0.    0.    0.  
\end{lstlisting}

\section{Interpolation de Lagrange}



%\vspace{-6.2cm}
\setindent{0cm}
\begin{flushleft}
\begin{minipage}[t]{4.2cm}
\begin{figure}
  \includegraphics[width=2.7cm]{Waring}
\caption{ Edward \textsc{Waring}\\(1736 - 1798)}
\end{figure}
\end{minipage}
\end{flushleft}
\setindent{\indentationcours}
\vspace{-5cm}
% \begin{minipage}[t]{\linewidth}


% \end{minipage}

Pour la  petite histoire,  même si \textsc{Lagrange}  publia en  1795 la
formule dans \og\textit{Leçons Élémentaires sur les Mathématiques Données a l'École Normale
}\fg{}, elle avait été publiée 16 ans plus tôt par le britannique
\textsc{Waring} dans \og \textit{Problems Concerning Interpolations
}\fg{} mais c'est l'inévitable \textsc{Gauss} qui en fit le
tour  et  la lia  avec  des  travaux  précédents de  \textsc{Newton}  et
\textsc{Euler}, sans nommer personne, bien sûr, dans sa \og{}\textit{Theoria Interpolationis Methodo Nova Tractata}\fg{}.



On travaillera dans cette section avec des polynômes de $\bbr_n[X]$, une
fonction numérique réelle  $f$ définie sur un intervalle  $[a,b]$ et une
famille $\bigl\{x_0,x_1,...,x_n\bigr\}$ de $n+1$ points distincts de $[a,b]$.

\begin{theoreme}
  Il existe un polynôme $P_n$ de $\bbr_n[X]$ et un seul tel que:
$$\forall               i\in               \bigl\{0,1,...,n\bigr\}\qquad
\widetilde{P_n}(x_i)=f(x_i)$$

Ce polynôme s'écrit:

$$P_n(X)=\displaystyle\sum_{i=0}^{n}f(x_i)\LR_i(X)$$

avec:

$$\LR_i(X)=\displaystyle\prod_{j=0,j\neq i}^{n}\frac{X-x_j}{x_i-x_j}$$
\end{theoreme}

La  preuve est  une  formalité  pour de  valeureux  $MP^\star$ (une  des
démonstrations fait le lien avec un déterminant de \textsc{Vandermonde}: comment?)...





D'un point de vue algorithmique,  que se passe-t-il lorsqu'on rajoute un
point pour améliorer l'interpolation?



\section{Interpolation de Newton}


\subsection{L'idée}

Pour remédier  à ce problème,  nous allons procéder  autrement: appelons
$f[x_1,x_2,...,x_n]$   le  coefficient   de  $x^n$   dans  le   polynôme
interpolateur  de \textsc{Lagrange}  (attention  à ne  pas le  confondre
avec le polynôme de \textsc{Lagrange} $\LR$).

On peut  montrer par récurrence  que ces  coefficients sont liés  par la
relation:


$$f[x_1,x_2,...,x_{n+1}]=\frac{f[x_2,...,x_{n+1}]-f[x_1,x_2,...,x_n]}{x_{n+1}-x_1}$$


avec $f[x_i]=f(x_i)$. On visualise la situation à l'aide d'un tableau du
type:





\subsection{Un exemple avec un pas constant}

%\vspace{-6.2cm}
\setindent{0cm}
\begin{flushleft}
\begin{minipage}[t]{4.2cm}
\begin{figure}
  \includegraphics[width=2.8cm]{Isaac}
\caption{ Isaac \textsc{Newton}\\(1642 - 1727)}
\end{figure}
\end{minipage}
\end{flushleft}
\setindent{\indentationcours}
\vspace{-5cm}
\begin{minipage}[t]{\linewidth}


\end{minipage}

Supposons que  nous connaissions quelques valeurs  du logarithme décimal
(par exemple par la méthode de \textsc{Briggs}) et que nous
voulions  connaître des  valeurs intermédiaires.  Nous allons  pour cela
procéder  par interpolation:  étant  données quatre  valeurs, nous  allons
chercher un  polynôme $P$ de degré  3 tel que $P(x)=\log(x)$  pour ces 4
valeurs.

Pour  cela,  nous  allons   utiliser  la  méthode  présentée  par  Isaac
\textsc{Newton} en 1676.

Avec  des notations modernes,  posons $P(x)=a+bx+cx^2+dx^3$.  On connaît
les valeurs du logarithmes pour 1, 2, 3 et 4.

\begin{center}
\begin{tabular}{|l|l|}
  \hline
$x_1=1$ & $p(x_1)=a+b+c+d=f(x_1)$\\
$x_2=2$ & $p(x_2)=a+2b+4c+8d=f(x_2)$\\
$x_3=3$ & $p(x_3)=a+3b+9c+27d=f(x_3)$\\
$x_4=4$ & $p(x_4)=a+4b+16c+64d=f(x_4)$\\
\hline
\end{tabular}
\end{center}

En soustrayant les lignes deux par deux, on fait disparaître $a$:

\begin{center}
\begin{tabular}{|l|}
  \hline
 $f(x_2)-f(x_1)=b+3c+7d=f[x_1,x_2]$\\
$f(x_3)-f(x_2)=b+5c+19d=f[x_2,x_3]$\\
$f(x_4)-f(x_3)=b+7c+37d=f[x_3,x_4]$\\
\hline
\end{tabular}
\end{center}

En soustrayant les lignes deux par deux, on fait disparaître $b$:

\begin{center}
\begin{tabular}{|l|}
  \hline
 $f[x_2,x_3]-f[x_1,x_2]=2c+12d=f[x_1,x_2,x_3]$\\
 $f[x_3,x_4]-f[x_2,x_3]=2c+18d=f[x_2,x_3,x_4]$\\
\hline
\end{tabular}
\end{center}

En soustrayant les lignes deux par deux, on fait disparaître $c$:


$f[x_2,x_3,x_4]-f[x_1,x_2,x_3]=6d=f[x_1,x_2,x_3,x_4]$

ce qui donne:

  $d=\frac{1}{6}f[x_1,x_2,x_3,x_4]$


  $c=\frac{1}{2}f[x_1,x_2,x_3]-f[x_1,x_2,x_3,x_4]$


 $b=f[x_1,x_2]-\frac{3}{2}f[x_1,x_2,x_3]+3f[x_1,x_2,x_3,x_4]-\frac{7}{6}f[x_1,x_2,x_3,x_4]=f[x_1,x_2]-\frac{3}{2}f[x_1,x_2,x_3]+\frac{11}{6}f[x_1,x_2,x_3,x_4]$


$a=f[x_1]-f[x_1,x_2]+\frac{3}{2}f[x_1,x_2,x_3]-\frac{11}{6}f[x_1,x_2,x_3,x_4]-\frac{1}{2}f[x_1,x_2,x_3]+\Delta^3y_1-\frac{1}{6}f[x_1,x_2,x_3,x_4]
=f[x_1]-f[x_1,x_2]+f[x_1,x_2,x_3]-f[x_1,x_2,x_3,x_4]$



Finalement:

\[
P(x)=f[x_1]+(x-1)f[x_1,x_2]+\frac{1}{2}(x-1)(x-2)f[x_1,x_2,x_3]+\frac{1}{6}(x-1)(x-2)(x-3)f[x_1,x_2,x_3,x_4]
\]


\textsc{Newton} avait l'habitude de  réunir les résultats dans un schéma
de cette forme:

\begin{center}
\begin{tabular}{cccc}
$\bm{f[x_1]}$ & $\bm{f[x_1,x_2]}$   & $\bm{f[x_1,x_2,x_3]}$   &  $\bm{f[x_1,x_2,x_3,x_4]}$                       \\
$f[x_2]$          &$f[x_2,x_3]$                      &$f[x_2,x_3,x_4]$ &                         \\
$f[x_3]$          & $f[x_3,x_4]$                     &        &                         \\
$f[x_4]$          &                      &                       &                         \\
\end{tabular}
\end{center}


avec $f[x_1,x_2,...,x_i]=\frac{f[x_2,x_3,...,x_i]-f[x_1,x_2,...,x_{i-1}]}{x_i-x_1}=f[x_2,x_3,...,x_i]-f[x_1,x_2,...,x_{i-1}]$



\subsection{Cas général}


On travaillera dans cette section avec des polynômes de $\bbr[X]$, une
fonction numérique réelle  $f$ définie sur un intervalle  $[a,b]$ et une
famille $\bigl\{x_1,...,x_n\bigr\}$ de $n$ points distincts de $[a,b]$.


On peut montrer que le polynôme interpolateur est défini par:


$$
P_n(t)=f[x_1]+\sum_{k=2}^nf[x_1,...,x_n](x-x_1)(x-x_2)\cdots(x-x_{k+1})
$$




Vous construirez une fonction \verb+function y = newton(X,Y,t)+ 
qui renvoie l'image de \verb+t+  par la fonction polynomiale obtenue par
l'interpolation  de Newton  sur les  points  dont les  abscisses et  les
ordonnées sont les éléments des matrices X et Y.

Le plus pratique est peut-être  de considérer un tableau: déterminez une
relation calculant  \texttt{T(i,j)} en fonction  de \texttt{T(i+1,j-1)},
\texttt{T(i,j-1)}, et des éléments de X...


\begin{lstlisting}[caption={}]
-->p = newton([0 1 2 3 5],[0 1 8 27 125],4)
 p  =
    64.  
\end{lstlisting}


Il peut  s'avérer plus utile  d'obtenir le polynôme lui-même  plutôt que
l'image d'un réel particulier par ce polynôme. 
  
Modifiez  alors
\verb+newton+ pour obtenir une fonction \verb+function y = newton_poly(X,Y)+ 
qui   renvoie  le   polynôme  d'interpolation   lui-même   en  utilisant
\texttt{prod\_mono} et \texttt{add\_poly}:

\begin{lstlisting}[caption={}]
-->newton_poly([0 1 2 3],[0 1 8 27])
 ans  =
    0.    0.    0.    1.  
\end{lstlisting}



Pour construire  une fonction Scilab  associée au problème,  nous aurons
besoin de deux fonctions intermédiaires:
\begin{itemize}
\item \verb+function y  = prod_mono(poly,t) + qui calcule  le produit de
      \verb+poly+ par le monôme de racine \verb+t+;

      \begin{lstlisting}[caption={}]
-->prod_mono([1 3 3 1],-1)
 ans  =
    1.    4.    6.    4.    1.  
      \end{lstlisting}

\item \verb+y = add_poly(p1,p2)+ qui  calcule la somme des polynômes (de
      degrés quelconques...) \verb+p1+ et \verb+p2+.

      \begin{lstlisting}[caption={}]
-->add_poly([1 2 3] , [0 2 4 6 8 10 12])
 ans  =
    1.    4.    7.    6.    8.    10.    12.  
      \end{lstlisting}

 \end{itemize}







Est-ce efficace?


Observons avec trois points de contrôle (approximation parabolique):

\begin{lstlisting}[caption={}]
-->p = newton_poly([42 43 44],[log(42) log(43) log(44)])
 p  =
    2.2608845    0.0465221  - 0.0002705  
 
-->horner(p,43.25)
 ans  =
    3.7669982  
 
-->log(43.25)
 ans  =
    3.7669972  
\end{lstlisting}

Observez en  des lieux \og moins  plats\fg{}. Avec plus  de points? Tout
ceci ne constitue que des  observations.%  Il existe en fait des théorèmes
% dans  le cas où  la fonction  réelle est  suffisamment régulière  et les
% écarts entre les abscisses sont constants.
% Si  $h$  est l'écart  entre  deux  abscisses et  si  $f$  est de  classe
% $\CR^{n}$ sur un $[a,b]$, si $x = x_1 + kh$, alors:

% $$f(x)-P(x) = h^{n}\frac{k(k-1)\cdots(k-n)}{n!}f^{(n)}(\xi_x)$$

% avec $\xi_x$ un élément du  plus petit intervalle contenant $x_1$, $x_n$
% et $x$.

% Des commentaires?



Occupons-nous  de  l'erreur  d'interpolation.  Soit  $P_n$  un  polynôme
interpolateur  en  $x_1$, $x_2$,...,$x_n$  points  distincts de  $[a,b]$
d'une fonction $f$ définie sur $[a,b]$.

Soit    $t\in[a,b]$    distinct    des    $x_i$    (sinon...).    Posons
$e_n(t)=P_n(t)-f(t)$  et  $P_{n+1}$  le  polynôme  qui  interpole  $f$  en
$x_1$,...,$x_n$, $t$.

En pensant à la diagonale descendante de \textsc{Newton}, on obtient:

$$P_{n+1}(t)=f(t)=P_n(t)+f[x_1,x_2,...,x_n,t]\prod_{i=0}^n(t-x_i)$$

donc

$$e_n(t)=f[x_1,x_2,...,x_n,t]\prod_{i=0}^n(t-x_i)$$


Une $n$\eme formule de \textsc{Cauchy} montre que si $f$ est de classe
$\CR^{n+1}$ sur $[a,b]$, il existe $\xi_t\in]a,b[$ tel que:

$$e_n(t)=\frac{f^{(n+1)}(\xi_t)}{(n+1)!}\prod_{i=0}^n(t-x_i)\quad
\text{soit}\quad \left|e_n(t)\right|\leqslant \underset{\xi\in]a,b[}{\text{max}}\left|\frac{f^{(n+1)}(\xi)}{(n+1)!}\prod_{i=0}^n(t-x_i)\right|$$



Le meilleur choix des abscisses pour minimiser l'erreur est donc:


$$\underset{(x_i)}{\text{min}}\left(\underset{t\in[a,b]}{\text{max}}\left|\prod_{i=0}^n(t-x_i)\right|\right)$$


%\vfill\eject

\section{Observation graphique}




Un  peu d'observation  maintenant. Pour  cela nous  aurons  besoin d'une
fonction 

\verb+function y  = poly2x(poly,x)+ qui renvoie  l'image d'un vecteur \verb+x+ par une
fonction polynomiale \verb+poly+.

\begin{lstlisting}[caption={}]
-->poly2x([1 2 1],x)
 ans  =
 
    1.    4.    9.  
\end{lstlisting}


Que fait la boucle suivante:

\begin{lstlisting}[caption={}]
for k = [0.1:0.1:1] do
  x = [0:k:2];
  p = newton_poly(x,1/(1+x^2));
  plot(x,poly2x(p,x))
end
\end{lstlisting}


Commentaires?   La   vie   n'est   pas   si  simple   qu'elle   n'en   a
l'air. Heureusement le cours de $MP^\star$ de M. \textsc{Sauvageot} va vous permettre de résoudre
ce  genre  de problème  de  manière  efficace.  Vous y  croiserez  alors
\textsc{Tchebychev},        \textsc{Bernstein},        \textsc{Fourier},
\textsc{Weierstra\ss}...



\section{Quelques idées pour la frise...}



%\begin{figure}
\begin{center}
  \includegraphics[height=3.5cm]{knuth}
\end{center}
%\end{figure}


Donald \textsc{Knuth} (prix \textsc{Turing 1974}):\og{}\textit{There are ways to amuse yourself while doing things and that's how I look at efficency}\fg{}.



%\begin{figure}
\begin{center}
  \includegraphics[height=3.5cm]{maccarthy}
\end{center}
%\end{figure}

John \textsc{Mac Carthy} (prix \textsc{Turing 1971}):\og\textit{ He who refuses to do arithmetic is doomed to talk nonsense}\fg{}.



%\begin{figure}
\begin{center}
  \includegraphics[height=3.5cm]{dijkstra}
\end{center}
%\end{figure}


Edsger \textsc{Dijkstra} (prix \textsc{Turing 1972}):\og \textit{The question \og{}What is Mathematics?\fg is as anavoidable and as unanswerable as the question \og{}What is Life?\fg. In actual fact I think it's almost the same question }\fg{}.

%%% Localg{} Variables: 
%%% TeX-master: "PolyMapleMPstar_12_13"
%%% End: 
