Interpolation.

Soit une fonction y = f(x) connue seulement par quelques couples (x,y).
Le problème qui se pose est le suivant: comment déterminer la valeur de y pour d'autres x que ceux donnés?

La technique permettant de répondre à cette question s'appelle "interpolation".

On peut considérer deux approches différentes:

la détermination d'un polynôme vérifiant les couples (x,y) de degré le plus petit (méthodes de Lagrange, Newton, ...). Le polynôme est choisi à cause de sa facilité de calcul (seulement des additions et des multiplications) et le degré minimum assure l'unicité de la solution;

la détermination d'une succession de polynômes de degré donné (en général, de degré 3) passant par deux points et tels que les raccords entre ceux-ci soient souples. Cette méthode dite des fonctions splines a un avantage énorme. Elle assure la tension minimale dans une pièce qui serait réalisée en respectant sa solution. Elle est donc utilisée pour déterminer des fuselages d'avion ou des pare-chocs d'automobiles.

Note: Il existe un autre problème fondamentalement différent du précédent mais très proche dans sa formulation: déterminer une fonction ne vérifiant pas les couples (x,y) donnés mais passant "au mieux" à travers les points donnés. La méthode de résolution la plus connue est la méthode des moindres carrés dont le principal inconvénient est de devoir déterminer au préalable le type de fonction qui répond le mieux au problème (exponentielle, logarithmique, ...).

Méthode des splines:

Soient n+1 points donnés (x_i;y_i) i=0,..,n. Nous recherchons une fonction y=f(x) qui passe par ces points et dont les dérivées f'(x) et f"(x) sont continues sur [x₀;x_n]. Les points x_i sont appelés noeuds et h_i=x_i+1-x_i le pas (pas nécessairement constant mais positif: x₀<x₁<...<x_n).

Il faut donc rechercher n fonctions S_i(X) i=0,1,..,n-1 telles que:

les couples (x_i;y_i) doivent vérifier les 2 fonctions dont ils font partie sauf les points (x₀;y₀) et (x_n;y_n) qui ne font partie que de la première et de la dernière fonction.
```
S₀(x₀)=f(x₀)=y₀                (1)

S_i(x_i)=S_i-1(x_i)=f(x_i)=y_i	(1)   pour tout i=1,...,n-1

S_n-1(x_n)=f(x_n)=y_n               (1)
```
les dérivées premières de 2 fonctions consécutives doivent être égales aux points communs (elles partagent la même tangente).
```
S'_i(x_i)=S'_i-1(x_i)		(2)    pour tout i=1,...,n-1
```
les dérivées deuxièmes de 2 fonctions consécutives doivent être égales aux points communs (elles doivent avoir la même concavité).
```
S"_i(x_i)=S"_i-1(x_i)		(3)    pour tout i=1,...,n-1
```

Comme les S_i(x) sont des polynômes de degré <= 3, nous devons déterminer, pour chaque cubique, les 4 coefficients des puissances de x. Il y a donc au total 4n nombres à déterminer grâce à:

2 n conditions (1)
(n-1) conditions (2)
(n-1) conditions (3)

Nous disposons donc de 4n-2 conditions linéaires en les coefficients des S_i(x).
Comme il nous manque 2 conditions, nous pourrons imposer 2 conditions supplémentaires: par exemple:

S'₀(x₀)=f'(x₀) et S'_n-1(x_n)=f'(x_n)

S'₀(x₀)=f'(x₀) et S"₀(x₀)=f"(x₀)

S"₀(x₀)=S"_n-1(x_n)=0  (splines naturelles)

........

Comme S_i(x) est de degré <=3, S"_i(x) est linéaire et donc combinaison linéaire de f"_i=f"(x_i) et f"_i+1=f"(x_i+1):

S"_i(x)=(x_i+1-x)/h_i*f"_i+(x-x_i)/h_i*f"_i+1

Deux intégrations successives permettent de faire apparaître:

S_i=y_i*(x_i+1-x)/h_i

  +y_i+1*(x-x_i)/h_i

  -h²_i/6*f"_i[(x_i+1-x)/h_i-((x_i+1-x)/h_i)³]

  -h²_i/6*f"_i+1[(x-x_i)/h_i-((x-x_i)/h_i)³]

En dérivant, on a:


S'_i(x_i)

=(y_i+1-y_i)/h_i-h_i/6*(f"_i+1+2*f"_i)

=S'_i-1(x_i)  même tangente aux points de contact 

=(y_i-y_i-1)/h_i-1+h_i-1/6*(f"_i-1+2*f"_i)

Ceci fournit n-1 équations. Comme ce nombre est insuffisant pour déterminer les n+1 f"_i, nous ajouterons f"₀=f"_N=0.
Ce système sera résolu par la méthode de Gauss-Seidel qui dit que si un système d'équations linéaires peut être écrit sous la forme
x=f(x, y, ...)
y=g(x, y, ...)
...
avec la matrice des coefficients des inconnues des seconds membres ayant une norme <1, la résolution se fait par calculs successifs des x, y, ... à partir de n'importe quelles valeurs initiales.

Méthode de Newton

Déterminons le polynôme d'interpolation de degré <= n admettant les valeurs y_i pour x=x_i i=0,..,n. Il n'est pas exigé que x_i+1-x_i=cste.
Ecrivons le polynome sous la forme:


P_n(X)=a₀

     +a₁*(x-x₀)

     +a₂*(x-x₀)*(x-x₁)

     +...

     +a_n*(x-x₀)*(x-x₁)*...*(x-x_N-1)

Si x=x₀:

P_n(x₀)=y₀=a₀

dont on déduit a₀=y₀
Si x=x₁:

P_n(x₁)=y₁=a₀+a₁*(x₁-x₀)=y₀+a₁*(x₁-x₀)

dont on déduit a₁=(y₁-y₀)/(x₁-x₀)=D¹y₀ différence divisée.
En continuant avec x=x₂, x=x₃, ..., on obtient


a_k=(D^k-1y_i+1-D^k-1y_i)/(x_i+k-x_i)=D^ky_i

La formule finale est donc


P_n(X)=y₀

     +Dy₀*(x-x₀)

     +D²y₀*(x-x₀)*(x-x₁)

     +...

     +D^Ny₀*(x-x₀)*(x-x₁)*...*(x-x_N-1)

En privilégiant x_N, on aurait obtenu:


P_n(X)=y_N

     +Dy₀*(x-x_N)

     +D²y_N-2*(x-x_N)*(x-x_N-1)

     +...

     +D^Ny₀*(x-x_N)*(x-x_N-1)*...*(x-x₁)

J'ai écrit un programme en Turbo Pascal qui interpole par les 2 méthodes et surtout dessine les graphes de chaque solution ainsi que la fonction initiale. Je sais que cela n'a pas de sens d'interpoler une fonction qu'on connaît mais cela permet de visualiser le degré de précision des 2 méthodes et entre autres, que quand le nombre de points donnés est grand et que les points sont équidistants, il apparaît des variations très grandes aux bords du domaine.
Pour charger une archive contenant le source, cliquer ici

Page d'accueil.