Introduction à l'interpolation numérique

Principales techniques et algorithmes

yahiko

- Titre : Introduction à l'interpolation numérique

- Auteur : yahiko

- Parution : 16 février 2015

- Public visé :

- Licence : Copyright ® 2015 yahiko. Aucune reproduction, même partielle, ne peut être faite de ce site et de l'ensemble de son contenu : textes, documents, images, etc. sans l'autorisation expresse de l'auteur. Sinon vous encourez selon la loi jusqu'à trois ans de prison et jusqu'à 300 000 € de dommages et intérêts.

Sommaire

- 1. Avant-propos

- 2. Introduction

- 3. Interpolation arrondie

- 4. Interpolation linéaire

- 5. Interpolation cosinusoïdale

- 6. Autres changements de variable

- 7. Interpolation polynomiale

- 7-1. Rappels sur les polynômes

- 7-2. Interpolation de Lagrange

- 7-3. Interpolation de Newton

- 8. Splines

- 8-1. Principes

- 8-2. Spline de Hermite

- 8-3. Spline de Catmull-Rom

- 8-4. Spline cardinale

- 8-5. Spline cubique

1 - Avant-propos

Cet article traite de l'interpolation numérique à destination des développeurs en informatique. À cet effet, des algorithmes transcriront les techniques d'interpolation présentées. Le langage d'implémentation sera TypeScript, un surensemble de JavaScript, dont la syntaxe est suffisamment proche du C/C++ ou de Java pour être transposable sans difficulté dans ces langages ou d'autres si nécessaire.

Même si ce n'est pas son objectif premier, cet article peut également intéresser un plus large public souhaitant comprendre l'interpolation puisqu'il tentera de présenter les principales techniques sous un angle pragmatique tout en apportant un éclairage théorique préalable. La plupart des résultats et des algorithmes seront explicités dans leurs principes et leur cheminement.

Néanmoins, afin de ne pas se perdre dans des explications sans fin, cet article présupposera que le lecteur est déjà familiarisé avec les principales techniques de manipulation de fonctions et connaît déjà des notions comme la continuité, la dérivabilité, et dans une moindre mesure les espaces vectoriels et le calcul matriciel.

Il sera fait peu mention des théorèmes associés à ce sujet et encore moins de leur démonstration dans la mesure où tout cela peut être retrouvé dans tout bon ouvrage mathématique sur l'Analyse ou dans les nombreux cours disponibles sur Internet.

2 - Introduction

Le but de l'interpolation est de compléter un ensemble de valeurs données à des positions données avec de nouvelles valeurs plausibles, en cohérence avec les valeurs déjà existantes.

L'ensemble des positions et des valeurs initiales pourraient être représentées par des points sur un repère orthonormé en deux dimensions comme ci-dessous :

Figure 1. Ensemble de points initiaux à interpoler.

0

0.8

1

0.5

2

0.1

3

0.4

4

0.6

5

0.5

6

0.3

En général, les points peuvent ne pas être espacés à intervalle régulier. Cependant, dans la pratique, les systèmes de mesure ou les algorithmes informatiques génèrent très souvent des données réparties à intervalle régulier. C'est pourquoi, pour des raisons de simplicité dans les explications et dans l'implémentation, nous ne considérerons que cette situation.

Une interpolation pourrait consister ici à déterminer pour une position donnée, par exemple , une valeur en .

Intuitivement, compte-tenu des points initiaux de la figure 1, une valeur trop élevée comme , ou bien négative ne serait pas cohérente. Il nous faut donc une règle nous permettant de calculer les différentes valeurs en en fonction de assurant la cohérence de l'ensemble des valeurs.

Une définition un peu plus rigoureuse serait de dire que l'interpolation consiste à déterminer une fonction, simple de préférence, dont la courbe représentative passe par les points initiaux. Cette fonction est appelée fonction interpolatrice ou interpolant. La fonction interpolatrice a pour objet de se rapprocher d'une fonction inconnue dont la courbe représentative passerait elle aussi par les points initiaux.

Si on note la fonction interpolatrice, alors elle doit vérifier où et correspondent respectivement aux positions et aux valeurs initiales.

Les points initiaux indiqués plus haut serviront de jeu d'exemple pour l'ensemble de ce document. Ceci afin de faciliter la compréhension et la comparaison entre les techniques.

L'interpolation est l'un des problèmes fondamentaux de l'analyse numérique avec des applications aussi bien en ingénierie avec par exemple la méthode des éléments finis, qu'en statistique et aussi en informatique, notamment en infographie (Computer Graphics).

À noter que l'interpolation est différente de l'approximation qui elle n'impose pas à la courbe représentative de la fonction approchante de passer par tous les points donnés.

3 - Interpolation arrondie

Figure 2. Courbe représentative d'une interpolation arrondie.

L'interpolation arrondie est également connue comme l'interpolation du plus proche voisin (nearest neighbour). Le principe est très simple puisqu'il s'agit de choisir pour une abscisse donnée, la valeur en du point le plus proche. Cela revient à définir une fonction constante par morceaux, communément appelée fonction en escalier.

Voici une implémentation de cette méthode en TypeScript pour une interpolation entre deux points :

function neighbour(y0: number, y1: number, t: number): number {
return (t < 0.5 ? y0 : y1);
}

Implémentation de l'interpolation arrondie

La fonction neighbour() ainsi définie prend trois paramètres. Les paramètres y0 et y1 et correspondent aux valeurs initiales. Le paramètre t correspond à la position relative de la valeur à interpoler. Si , alors la fonction neighbour() retournera y0. Si , alors la fonction neighbour() retournera y1.

Pour avoir l'interpolation en , nous devons renseigner les valeurs en pour () et pour (). Pour la valeur du paramètre t cela revient et considérer un nouveau repère ayant pour origine ( ; 0). Dans ce nouveau repère, . Le code correspondant serait :

y = neighbour(0.4, 0.6, 0.1);
// y === 0.4

Exemple d'utilisation de la fonction neighbour()

Il est facile de voir qu'en appliquant cette fonction à chaque paire de points successifs, en faisant varier à chaque fois t entre 0 et 1, on obtient une courbe en escalier comme celle de la figure 2.

Le coût en calcul de la fonction neighbour() est très faible puisqu'elle ne comporte en tout et pour tout qu'un seul test logique. C'est probablement l'interpolation la plus simple et la plus rapide. Si l'intervalle entre deux valeurs initiales est faible, réduisant d'autant l'erreur potentielle d'interpolation, cette méthode peut être une option tout à fait pertinente. Elle est d'ailleurs régulièrement utilisée en animation 3D au niveau du rendu des textures (texture filtering).

Cependant, il apparaît qu'une telle courbe en escalier n'est pas des plus « esthétique » et peut ne pas être en cohérence avec les valeurs attendues. L'interpolation linéaire est une étape vers plus de vraisemblance.

4 - Interpolation linéaire

Figure 3. Courbe représentative d'une interpolation linéaire

L'interpolation linéaire remplace les fonctions constantes par morceaux par des fonctions linéaires (en toute rigueur, affines) dont la courbe représentative est composée de segments de droites (dans le cas d'une courbe représentative en deux dimensions).

Pour rappel, une droite a pour équation tels que et sont des nombres réels où est le coefficient directeur (ou encore la pente) et est l'ordonnée à l'origine.

Pour obtenir une valeur de entre deux points en fonction de , il faut donc déterminer les valeurs de et .

Cela revient à résoudre un système linéaire de deux équations à deux inconnues.

Sur l'intervalle [0 ; 1] on a :

function lerp(y0: number, y1: number, t: number): number {
return (y1 - y0) * t + y0;
}

Implémentation non stable de l'interpolation linéaire

Tout comme la fonction neighbour() vue précédemment dans le cadre de l'interpolation arrondie, la fonction lerp() prend les trois mêmes paramètres.

Pour avoir l'interpolation en , le code correspondant serait :

var y = lerp(0.4, 0.6, 0.1);
// y === 0.42000000000000004

Exemple d'utilisation de la fonction lerp()

Bien que théoriquement juste, le code de l'interpolation linéaire fourni précédemment ne garanti pas les relations lerp(y0, y1, 0) === y0 et lerp(y0, y1, 1) === y1 qui définissent le principe de base de l'interpolation : où est un point initialement donné.

À cause de la représentation même des nombres décimaux au sein de la machine, l'addition ou la soustraction entre deux nombres ayant un écart important d'ordre de grandeur, peut entraîner des aberrations de calcul. Ce n'est donc pas un bug lié au langage TypeScript. Pour s'en convaincre, avec l'implémentation précédente, l'appel lerp(10e20, 0.6, 1) renvoie 0 au lieu de la valeur attendue qui est 0.6.

Pour remédier à ce problème d'implémentation, il est possible de reformuler la fonction lerp() de la sorte :

function lerp(y0: number, y1: number, t: number): number {
return y0 * (1.0 - t) + y1 * t;
}

Implémentation stable de l'interpolation linéaire

L'interpolation linéaire est utilisée dans d'innombrables applications de par sa simplicité et sa rapidité. De plus, par rapport à l'interpolation arrondie, elle améliore significativement la vraisemblance des valeurs interpolées.

Si on considère sa forme stable, cette méthode ne met en jeu que deux multiplications et deux additions. Ceci peut s'avérer plus performant même que le test logique de l'interpolation arrondie sur certaines implémentations et plate-formes.

Néanmoins, l'interpolation linéaire a ses défauts, notamment les « discontinuités » présentes aux points initiaux, peu naturelles. Pour être plus exact, la dérivée première de l'interpolation linéaire par morceau n'est pas continue.

Pour modéliser plus précisément de nombreux phénomènes du monde réel, il peut être souhaitable de rendre le résultat plus lisse. Pour cela, il existe plusieurs techniques.

5 - Interpolation cosinusoïdale

Figure 4. Courbe représentative d'une interpolation cosinusoïdale

Le but de l'interpolation cosinusoïdale est de lisser les angles formés par les segments de droite de l'interpolation linéaire en recourant à un changement de variable.

Plutôt que d'injecter directement la variable dans la fonction lerp(), l'idée est de « ralentir » la variation de au voisinage des points initiaux, adoucissant ainsi la courbe représentative résultante. L'essentiel étant que pour tout dans l'intervalle [0 ; 1], la nouvelle variable soit également dans l'intervalle [0 ; 1].

Cela peut être obtenu en remaniant légèrement la fonction cosinus. En ne considérant que sa demi-période croissante et en la mettant à l'échelle pour que l'image soit comprise dans l'intervalle [0 ; 1], on peut obtenir une nouvelle variable exprimée en fonction de dont les variations ralentissent à proximité des points initiaux.

Figure 5. Variable u en fonction de t sur l'intervalle [0 ; 1]

L'implémentation de l'interpolation cosinusoïdale fait donc appel à l'interpolation linéaire via la fonction lerp() après avoir transformé la variable fournie en paramètre avec la fonction cosinus que nous venons de voir.

function cosine(y0: number, y1: number, t: number): number {
var u = (1.0 - Math.cos(t * Math.PI)) * 0.5;
return lerp(y0, y1, u);
}

Implémentation de l'interpolation cosinusoïdale

Cette nouvelle fonction cosine() prend les mêmes paramètres que les fonctions précédentes. Le lecteur pourra donc se reporter aux précédents chapitres pour avoir plus d'information sur ces paramètres.

La courbe représentative de cette nouvelle fonction est assez satisfaisante dans la mesure où elle est relativement lisse même si elle présente quelques artefacts au niveau de certaines jonctions.

Cependant, si la fonction cosine() est appelée intensivement, cela peut poser des soucis de performance, à cause du coût relativement élevé de la fonction cosinus.

C'est pourquoi d'autres changements de variable que celui que nous venons de présenter ont été mis au point.

6 - Autres changements de variable

Pour améliorer les performances ou pour appliquer un autre type de transition entre deux points initiaux à interpoler, il est possible de faire appel à d'autres changements de variable que celui de la fonction cosinus. L'essentiel étant que pour tout compris dans l'intervalle [0 ; 1], la nouvelle variable soit elle aussi comprise dans l'intervalle [0 ; 1].

Les interpolations les plus courantes se basant sur un changement de variable de l'interpolation linéaire sont les suivantes :

- Smoothstep :

- Smootherstep :

- La parabole concave (accélération) :

- La parabole convexe (décélération) :

Figure 6. Variable u en fonction de t sur l'intervalle [0 ; 1]

Sur la figure 6, on peut constater que la fonction smoothstep (courbe rouge) est très proche de la fonction cosinus (courbe noire) mais comme la fonction smoothstep est un polynôme cubique ne mettant en jeu que quelques additions et multiplications, ses performances seront meilleures. Ce type d'interpolation est d'ailleurs implémentée nativement dans les langages de shading pour le rendu d'images 3D.

La fonction smootherstep a été proposée par Ken Perlin (l'inventeur du bruit de Perlin) et a pour avantage d'avoir non seulement une dérivée première qui s'annule pour et comme pour smoothstep, mais également une dérivée seconde qui s'annule aussi en 0 et en 1 accentuant ainsi la courbure en 0 et en 1.

À la formule du changement de variable près, l'implémentation de ces quatre nouvelles formes d'interpolations est identique à celle de l'interpolation cosinusoïdale.

function smoothstep(y0: number, y1: number, t: number): number {
var u = t*t*(3 - 2*t);
return lerp(y0, y1, u);
}

Implémentation de l'interpolation Smoothstep

Figure 7. Courbe représentative d'une interpolation Smoothstep

function smootherstep(y0: number, y1: number, t: number): number {
var u = t*t*t*(t*(t*6 - 15) + 10);
return lerp(y0, y1, u);
}

Implémentation de l'interpolation Smootherstep

Figure 8. Courbe représentative d'une interpolation Smootherstep

function acceleration(y0: number, y1: number, t: number): number {
var u = t*t;
return lerp(y0, y1, u);
}

Implémentation de l'interpolation Acceleration

Figure 9. Courbe représentative d'une interpolation Acceleration

function deceleration(y0: number, y1: number, t: number): number {
var u = t * (2 - t);
return lerp(y0, y1, u);
}

Implémentation de l'interpolation Deceleration

Figure 10. Courbe représentative d'une interpolation Deceleration

Bien évidemment, il est possible d'imaginer bien d'autres changements de variable selon les besoins, néanmoins, les artefacts aux jonctions persisteront. Seules d'autres approches peuvent gommer ce qui peut être considéré comme un défaut.

7 - Interpolation polynomiale

Ce chapitre plutôt théorique et nécessitant un certain niveau en mathématiques peut être ignoré en première lecture même s'il peut faciliter la compréhension du chapitre sur les splines.

7-1 - Rappels sur les polynômes

Tout ce que nous avons vu jusqu'à présent consistait à relier des paires de points à l'aide de fonctions indépendantes pour obtenir au final une courbe plus ou moins « lisse ». C'est ce qu'on appelle l'interpolation par morceau.

Dans ce chapitre, nous montrerons qu'il existe une fonction interpolatrice capable de relier un nombre arbitraire de points en un seul tenant.

L'idée est de se dire qu'il doit peut-être exister un polynôme dont la courbe représentative passe par tous les points initiaux.

Pour rappel, un polynôme est une fonction usuellement présentée sous la forme d'une somme de monômes tel que sont les coefficients de ce polynôme et où la plus grande puissance de correspond au degré de ce polynôme. À noter qu'ici les coefficients seront classés par ordre de puissance croissante par souci de cohérence avec les algorithmes.

Par exemple, une équation de droite (comme pour l'interpolation linéaire) est un polynôme de degré 1 et de coefficients .

La fonction smoothstep vue précédemment est un polynôme de degré 3 et de coefficients .

Supposons que pour points initiaux, il existe un polynôme de degré dont la courbe représentative passe par tous ces points. En se basant sur la forme usuelle d'un polynôme, cela revient à déterminer les coefficients de ce polynôme à partir d'un système de Vandermonde :

Dans le cas général, un tel système est mal conditionné, c'est-à-dire, que les erreurs d'arrondi deviennent rapidement importantes et faussent significativement le calcul.

Pour remédier à cet écueil, on préférera utiliser une autre formulation des polynômes.

7-2 - Interpolation de Lagrange

Figure 11. Courbe représentative d'une interpolation polynomiale de Lagrange

Pour contourner le système de Vandermonde, l'idée est de réécrire le polynôme de degré comme une combinaison linéaire de polynômes de degré possédant les propriétés suivantes :

1. Les polynômes forment une base dans l'espace vectoriel des polynômes de degré inférieur ou égal à .

Ces polynômes sont communément nommés les polynômes d'interpolation de Lagrange.

La propriété 2 doit s'interpréter comme le fait que chaque polynôme a pour racines les points initiaux lorsque (le polynôme s'annule). Lorsque , le polynôme ne s'annule pas et voit sa valeur normalisée à 1.

Chaque polynôme de degré doit donc avoir racines parmi les points .

Ces polynômes peuvent donc s'écrire sous la forme d'un produit à termes faisant apparaître les racines. Par exemple pour , un tel produit peut s'écrire , pour , le produit s'écrit , et ainsi de suite.

Plus généralement, pour un donné, ce produit peut s'écrire ainsi :

Ce produit s'annule lorsqu'un de ses termes vaut zéro et donc lorsque pour un donné, tel que . Par exemple, pour , si le produit s'annule puisque le facteur .

Par contre, pour un donné, lorsque , le produit ne s'annule pas. Autrement dit, pour chaque indice , il existe un point initial parmi les points pour lequel le produit ne s'annule pas. Par exemple, pour , si le produit ne s'annule pas et ce produit vaut .

Or, puisque nous souhaitons avoir une base normalisée et obtenir 1 pour pour un donné, il est nécessaire de diviser l'expression du produit par la valeur de ce produit quand .

Les polynômes ont donc la forme suivante :

Ceci peut également s'écrire de façon plus compacte avec la notation du produit indexé :

Il ne reste plus qu'à faire passer chaque polynôme par l'ordonnée du point initial correspondant et à en faire la somme pour obtenir finalement le polynôme interpolateur de degré  :

Ce polynôme vérifie par construction .

L'interpolation de Lagrange démontre donc l'existence d'un polynôme de degré passant par points, assertion uniquement présupposée à la section précédente.

Il est possible de transcrire littéralement en langage informatique l'expression du polynôme interpolateur de Lagrange à l'aide de deux boucles imbriquées.

function lagrange(X: number[], Y: number[], x: number): number {
var sum = 0;
for (var i = 0; i < X.length; i++) {
var product = 1;
for (var k = 0; k < X.length; k++) {
if (k != i) {
product *= (x - X[k]) / (X[i] - X[k]);
}
} // for k
sum += Y[i] * product;
} // for i
return sum;
}

Implémentation de l'interpolation de Lagrange

Le paramètre X est la liste des abscisses des points initiaux. Le paramètre Y est la liste des ordonnées des points initiaux tel que Y[0] correspond à l'ordonnée en X[0], Y[1] est la valeur en X[1] et ainsi de suite jusqu'à Y[n] où est le nombre de valeurs initiales. Le paramètre x est l'abscisse du point à interpoler sur l'intervalle [X[0] ; X[n]]. La fonction retourne la valeur de l'interpolation de Lagrange en x.

var y = lagrange([0, 1, 2, 3, 4, 5, 6], [0.8, 0.5, 0.1, 0.4, 0.6, 0.5, 0.3], 2.4);
// y === 0.17979647999999995

Exemple d'utilisation de la fonction lagrange()

On peut noter que la complexité de cet algorithme est élevé. Il y a additions et multiplications ce qui peut se révéler très coûteux en temps de calcul à mesure que le degré et donc que le nombre de points à interpoler augmente. Même s'il est possible de faire appel à des algorithmes plus efficaces comme celui de Neville pour des petits degrés, il est souvent préférable de déterminer au préalable les coefficients du polynôme et ensuite obtenir l'interpolation en .

Pour calculer les coefficients de , une possibilité est de calculer les coefficients de chaque polynôme , de les multiplier par et ensuite de faire la somme pour chaque puissance de .

Plus formellement, chaque polynôme peut s'écrire sous la forme usuelle :

tel que sont ses coefficients.

Par conséquent, chaque coefficient du polynôme peut s'exprimer ainsi :

Pour calculer les coefficients d'un polynôme , posons :

À une racine près, ce nouveau polynôme a les même racines que celles de chaque polynôme , de telle sorte qu'il est possible d'exprimer en fonction de  :

Cette relation permet de déduire les coefficients de en fonction des coefficients de . qui eux peuvent s'obtenir par un simple développement successif des termes du produit.

Pour obtenir les coefficients de , il suffit pour chacune des puissances de , de sommer les coefficients des polynômes correspondants multipliés par un facteur .

Armé de toutes ces relations, nous pouvons donc écrire un algorithme permettant de calculer les coefficients du polynôme interpolateur de Lagrange .

function lagrangeCoefficients(X: number[], Y: number[]): number[] {
var n = X.length - 1; // degrée
// coeffs[] sont les coefficients du polynôme interpolateur f(x)
var coeffs = new Array(n + 1);
// beta[] sont les coefficients des polynômes L[i](x) (x-X[0])...(x-X[i-1])(x-X[i+1])...(x-X[n])
var beta = Array(n + 1);

// initialisation des coefficients
for (var i = 0; i < n + 1; i++) {
coeffs[i] = 0.0;
} // for i

// gamma[] sont les coefficients de (x-X[0])(x-X[1])...(x-X[n])
var gamma = new Array(n + 2);
gamma[0] = 1.0;
for (var i = 0; i < n + 1; i++) {
for (var j = i; j > 0; j--) {
gamma[j] = gamma[j-1] - gamma[j] * X[i];
} // for j
gamma[0] *= -X[i];
gamma[i+1] = 1;
} // for i

// calculs des coefficients
for (var i = 0; i < n + 1; i++) {
// dénominateur du polynôme d'interpolation de Lagrange Li(x)
// denom = (X[i]-X[0])...(X[i]-X[i-1])(X[i]-X[i+1])...(X[i]-X[n])
var denom = 1;
for (var j = 0; j < n + 1; j++) {
if (i != j) {
denom *= X[i] - X[j];
}
} // for j

var lambda = Y[i] / denom;
beta[n] = gamma[n + 1]; // gamma[n1] = 1
coeffs[n] += lambda * beta[n];
for (var j = n - 1; j >= 0; j--) {
beta[j] = gamma[j+1] + beta[j+1] * X[i];
coeffs[j] += lambda * beta[j];
} // for j
} // for i

return coeffs;
}

Implémentation du calcul des coefficients du polynôme de Lagrange

Le paramètre X est la liste des abscisses des points initiaux. Le paramètre Y est la liste des ordonnées des points initiaux tel que Y[0] correspond à l'ordonnée en X[0], Y[1] est la valeur en X[1] et ainsi de suite jusqu'à Y[n] où est le nombre de valeurs initiales. La fonction lagrangeCoefficients() retourne la liste des coefficients du polynôme interpolateur de Lagrange par ordre de puissance croissante.

var fCoeffs = lagrangeCoefficients([0, 1, 2, 3, 4, 5, 6], [0.8, 0.5, 0.1, 0.4, 0.6, 0.5, 0.3]);
// fCoeffs === [
//   0.8,
//   1.256666666666667,
//   -3.013333333333334,
//   1.9250000000000007,
//   -0.5333333333333337,
//   0.06833333333333334,
//   -0.003333333333333334]

Exemple d'utilisation de la fonction lagrangeCoefficients()

Une fois que nous avons les coefficients d'un polynôme, il est relativement simple et efficace de calculer sa valeur en un point à l'aide par exemple de la méthode de Horner.

function horner(A: number[], x: number): number {
var n = A.length - 1; // degré
var result = 0;

for (var i = n; i >= 0; i--) {
result = result * x + A[i];
} // for i

return result;
}

Implémentation de la méthode de Horner

Le paramètre A représente la liste des coefficients par degré croissant du polynôme à évaluer. Le paramètre x est l'abscisse du point à évaluer par le polynôme. La fonction horner() renvoie donc la valeur en x du polynôme de coefficients A.

Connaissant les coefficients, la complexité est en . Ainsi, il est donc possible de calculer efficacement un grand nombre de valeurs d'interpolation.

Cependant, bien qu'il soit plus efficace dans l'absolu de pré-calculer les coefficients du polynôme interpolateur que d'utiliser naïvement la formule de Lagrange en chaque point à interpoler, une rapide évaluation de la complexité de la fonction lagrangeCoefficients() indique tout de même une complexité en . Autrement dit, à partir d'un certain nombre de points, l'interpolation de Lagrange devient coûteuse.

D'autant plus qu'ajouter un nouveau point nous oblige recalculer chaque polynôme et chaque coefficient sans pouvoir réutiliser des résultats antérieurs à l'ajout de ce nouveau point.

Enfin, bien que préférable à la résolution d'un système de Vandermonde, cette méthode reste sensible aux erreurs d'arrondi.

7-3 - Interpolation de Newton

L'interpolation de Newton respecte les même conditions que l'interpolation de Lagrange, à savoir où pour tout point initial . Le polynôme sous-jacent est le même que celui obtenu par la forme de Lagrange. C'est pourquoi il est plus correct de parler de forme de Newton car seule la formulation diffère. La forme de Newton est donc en quelque sorte un autre algorithme pour déterminer le polynôme interpolateur avec pour but avoué de calculer plus efficacement les valeurs du polynôme interpolateur.

Nous avons vu avec l'interpolation de Lagrange qu'il était possible d'exprimer le polynôme interpolateur à partir d'une certaine base de polynômes. Avec la forme de Newton, l'idée est de se donner une autre base de telle sorte que l'ajout d'un nouveau point initial ne nous obligera pas de recalculer les polynômes de cette base, mais simplement d'adjoindre un nouveau polynôme à la base déjà existante.

Pour cela, commençons par ne considérer que le premier point initial .

Cherchons un polynôme de degré 0 tel que .

C'est très simple. C'est le polynôme constant de coefficient . On a donc .

Ajoutons ensuite un nouveau point .

Cherchons un polynôme de degré 1 vérifiant les conditions d'interpolation et et s'exprimant en fonction de .

Ce polynôme peut s'écrire ainsi : .

Comme , avec coefficient réel, la première condition d'interpolation est donc vérifiée.

Déterminons le coefficient de façon à vérifier la seconde condition d'interpolation :

Typiquement, est la pente de la droite d'équation (et l'ordonnée à l'origine).

Ajoutons ensuite un nouveau point et cherchons un polynôme en poursuivant le processus.

On a bien par construction . Reste à déterminer le coefficient pour que .

Après un rapide calcul on obtient :

En généralisant, ajoutons un ème point à un ensemble de points déjà existants.

Pour un ensemble de points et son polynôme interpolateur , il est donc possible d'ajouter un nouveau point et d'exprimer le nouveau polynôme interpolateur en fonction . Les produits des termes faisant apparaître les racines sont appelés polynômes de Newton et forment une base.

Le polynôme interpolateur peut ainsi s'écrire comme une combinaison linéaire de ces polynômes de Newton :

Il ne reste à présent qu'à expliciter les coefficients que nous nommerons coefficients de Newton. Ils peuvent s'exprimer sous la forme de différences divisées, usuellement notées […], tel que .

Les différences divisées sont une succession de divisions définies de manière récursive :

D'où .

Pour implémenter l'interpolation de Newton, nous avons donc besoin d'une fonction calculant tout d'abord les différences divisées puis d'une autre calculant la valeur interpolée proprement dite via la formule de et des polynômes .

function divDiffsRec(X: number[], Y: number[]): number {
if (X.length === 1) {
return Y[0];
}
else {
return (
divDiffsRec(X.slice(1, X.length), Y.slice(1, Y.length)) -
divDiffsRec(X.slice(0, X.length - 1), Y.slice(0, Y.length - 1))
) / (X[X.length - 1] - X[0]);
}
}

Implémentation récursive des différences divisées

Le paramètre X est la liste des abscisses des points initiaux. Le paramètre Y est la liste des ordonnées des points initiaux tel que Y[0] correspond à l'ordonnée en X[0], Y[1] est la valeur en X[1] et ainsi de suite jusqu'à Y[n]. La fonction divDiffsRec() renvoie la valeur des différences divisées selon les vecteurs X et Y.

var alphaN = divDiffsRec([0, 1, 2, 3, 4, 5, 6], [0.8, 0.5, 0.1, 0.4, 0.6, 0.5, 0.3]);
// alphaN === -0.0033333333333333353

Exemple d'utilisation de la fonction divDiffsRec()

Nous pouvons donc définir une fonction permettant de calculer l'interpolation en un point donné à l'aide des différences divisées.

function newtonRec(X: number[], Y: number[], x: number): number {
var sum = 0;
for (var i = 0; i < Y.length; i++) {
var product = 1;
for (var k = 0; k < i; k++) {
if (k != i) {
product *= x - k;
}
} // for k
sum += divDiffsRec(X.slice(0, i+1), Y.slice(0, i+1)) * product;
} // for i
return sum;
}

Implémentation de l'interpolation de Newton à l'aide des différences divisées récursives

La fonction newton() appelée avec les même paramètres que ceux utilisés plus haut pour la fonction lagrange() renvoie sans surprise le même résultat.

var y = newton([0, 1, 2, 3, 4, 5, 6], [0.8, 0.5, 0.1, 0.4, 0.6, 0.5, 0.3], 2.4);
// y === 0.17979647999999995

Exemple d'utilisation de la fonction lagrange()

Cette implémentation de l'interpolation de Newton utilisant les différences divisées récursives, bien qu'absolument correcte, a un défaut qui est de réaliser des calculs redondants sur les différences divisées intermédiaires. En effet dans la boucle principale, à l'itération , dépend de qui a été calculé à l'itération .

Pour optimiser, nous pourrions simplement mémoriser les différences divisées de l'itération précédente. Cependant, il est possible de faire mieux en reconsidérant le calcul des différences divisées non plus de façon récursive, transcription directe de sa définition, mais de façon itérative, en construisant au fur et à mesure une liste des coefficients de Newton.

Coefficients

Initialisation

Itération 0

Itération 1

Itération 2

Ces derniers pourront ainsi être utilisés efficacement soit par l'implémentation précédente de l'interpolation de Newton, soit par un algorithme de type méthode de Horner comme vue précédemment qui s'applique aisément à la forme de Newton.

function divDiffsIter(X: number[], Y: number[]): number[] {
var n = X.length - 1; // degré
var alphas = Y; // initialisation des différences divisées

for (var i = 0; i < n; i++) {
for (var j = n; j > i; j--) {
alphas[j] = (alphas[j] - alphas[j - 1]) / (X[j] - X[j - i - 1]);
} // for j
} // for i

return alphas;
}

Implémentation itérative des différences divisées

Le paramètre X est la liste des abscisses des points initiaux. Le paramètre Y est la liste des ordonnées des points initiaux tel que Y[0] correspond à l'ordonnée en X[0], Y[1] est la valeur en X[1] et ainsi de suite jusqu'à Y[n].

var alphas = divDiffsIter([0, 1, 2, 3, 4, 5, 6], [0.8, 0.5, 0.1, 0.4, 0.6, 0.5, 0.3]);
// alphas === [
//   0.8,
//   -0.30000000000000004,
//   -0.04999999999999999,
//   0.13333333333333333,
//   -0.06666666666666668,
//   0.01833333333333334,
//   -0.0033333333333333353]

Exemple d'utilisation de la fonction divDiffsIter()

À l'aide de cette nouvelle implémentation des différences divisées et en utilisant la méthode de Horner, nous pouvons réécrire l'interpolation de Newton pour un point donné.

function newtonHorner(X: number[], alphas: number[], x: number): number {
var n = X.length - 1; // degré
var result = 0;

for (var i = n; i >= 0; i--) {
result = result * (x - X[i]) + alphas[i];
} // for i

return result;
}

Implémentation de l'interpolation de Newton à l'aide de la méthode de Horner

Le paramètre X est la liste des abscisses des points initiaux. Le paramètre alphas correspond aux coefficients de Newton calculés préalablement, par exemple à l'aide des différences divisées. Le paramètre x est l'abscisse du point à interpoler sur l'intervalle [X[0] ; X[n]]. La fonction newtonHorner() renvoie la valeur de l'interpolation en x.

var alphas = divDiffsIter([0, 1, 2, 3, 4, 5, 6], [0.8, 0.5, 0.1, 0.4, 0.6, 0.5, 0.3]);
var y = newtonHorner([0, 1, 2, 3, 4, 5, 6], alphas, 2.4);
// y === 0.17979648000000004

Exemple d'utilisation de la fonction newtonHorner()

C'est cette méthode qui est la plus couramment utilisée pour l'interpolation polynômiale.

Bien que nous n'ayons pas réellement besoin de calculer les coefficients du polynôme de , cela reste tout à fait possible via les coefficients de Newton . Cela est tout de même coûteux puisque la complexité de l'algorithme correspondant reste en comme son homologue pour la méthode de Lagrange.

function lagrangeCoefficientsNewton(X: number[], alphas: number[]): number[] {
var n = X.length - 1; // degré
// coeffs[] sont les coefficients du polynôme interpolateur f(x)
var coeffs = new Array(n + 1);

// initialisation des coefficients
for (var i = 0; i < n + 1; i++) {
coeffs[i] = 0.0;
} // for i

coeffs[0] = alphas[n];
for (var i = n - 1; i >= 0; i--) {
for (var j = n - i; j > 0; j--) {
coeffs[j] = coeffs[j-1] - X[i] * coeffs[j];
} // for j
coeffs[0] = alphas[i] - X[i] * coeffs[0];
} // for i

return coeffs;
}

Implémentation du calcul des coefficients du polynôme de Lagrange à l'aide des coefficients de Newton.

Qu'elle utilise la méthode de Lagrange ou de Newton, l'interpolation polynomiale produit une courbe représentative plus régulière, plus lisse, que les interpolations par morceaux vues précédemment. Néanmoins, malgré les méthodes et algorithmes pour calculer efficacement les valeurs, elle nécessite de manipuler des polynômes de degré potentiellement élevé, directement fonction du nombre de points initiaux. De plus, dans le cas particulier mais pourtant fréquent où les points sont équirépartis comme c'est l'hypothèse dans cet article, les valeurs tendent à osciller à mesure que le degré du polynôme augmente. C'est le phénomène de Runge qui rend les valeurs de moins en moins cohérentes. Ce qui est un peu paradoxal.

Par conséquence, dans les situations où le nombre de points initiaux est élevé, par exemple 100, 1 000, 10 000 voire davantage, une interpolation polynomiale de Lagrange ou de Newton est à exclure.

8 - Splines

8-1 - Principes

Les techniques que nous avons rencontrées jusqu'à présent comportent certes des qualités, mais également des inconvénients majeurs.

Les fonctions simples par morceaux présentent des irrégularités aux jonctions, elles ne sont pas dérivables sur l'ensemble de l'intervalle [ ; ]. Cela est dû à leur caractère local.

L'interpolation polynomiale de Lagrange (ou de Newton), bien que dérivable sur [ ; ], a quant à elle une complexité croissante à mesure que le nombre de points à interpoler augmente. Cela est dû à son caractère global.

L'idéal serait de conjuguer la simplicité et la flexibilité de l'interpolation par morceaux, avec la régularité et la dérivabilité de l'interpolation polynomiale.

Pour cela, il existe une technique au caractère intermédiaire entre le niveau local d'une interpolation simple par morceau et le niveau global d'une interpolation polynomiale de Lagrange. Il s'agit du concept de spline, fonction polynomiale par morceaux, qui se veut « lisse » sur l'ensemble de l'intervalle, et en particulier aux points de jonction qu'on appelle aussi des nœuds.

Les splines sont basées sur l'interpolation de Hermite. De façon schématique, ce type d'interpolation étend non seulement les conditions sur les valeurs aux nœuds de la fonction d'interpolation comme vu au chapitre précédent, mais aussi sur ses dérivées.

C'est pourquoi l'interpolation par splines peut être considérée comme une technique semi-locale ou semi-globale, c'est selon.

Figure 12. Notations sur les splines

Plus formellement, une spline est une fonction définie par morceaux sur un intervalle [ ; ]. Au ème sous-intervalle [ ; ] est défini un morceau de la spline que nous noterons .

Chaque fonction est un polynôme devant vérifier les conditions d'interpolations et .

Pour assurer continuité et la dérivabilité sur l'ensemble de l'intervalle (et notamment aux jonctions), une spline impose également les conditions suivantes :

Pour assurer la dérivabilité d'ordre 1, il est nécessaire que le polynôme soit de degré supérieur ou égal à deux. En généralisant, on montre que pour assurer une dérivabilité jusqu'à un ordre , il est nécessaire que soit un polynôme de degré supérieur ou égal à .

Plus le degré des morceaux sera élevé, plus la spline aura un aspect lisse, mais en contrepartie, plus les calculs seront coûteux. C'est pourquoi en pratique on se limite à des splines cubiques constituées de polynômes de degré 3.

Il est donc possible d'imposer à de telles splines des conditions de dérivabilité aux nœuds jusqu'à l'ordre 2. Au-delà, les dérivées s'annulent.

8-2 - Spline de Hermite

Les splines de Hermite sont une famille de splines cubiques se limitant à être une fois continûment dérivable (classe ) sur l'intervalle [ ; ].

Cela rend cette famille de splines plus simple à manipuler par rapport au cas général d'une continuité en .

Posons le changement de variable . On peut obtenir ainsi une fonction dans l'intervalle [0 ; 1] équivalente à la fonction dans le ème sous-intervalle [ ; ].

L'expression générale d'un morceau de spline, qui est ici un polynôme de degré 3, et de sa dérivée première peut donc s'écrire :

En et , on obtient les valeurs de et de sa dérivée :

Ceci nous permet d'exprimer les quatre coefficients , , et de notre morceau de spline :

En utilisant les conditions d'interpolation on obtient :

avec et les pentes des tangentes et (cf. figure 13) aux nœuds et . Il est facile de deviner que ces pentes influencent l'aspect de la spline en chaque nœud.

Figure 13. Tangentes d'une spline de Hermite

Il est courant d'utiliser la représentation matricielle pour caractériser les polynômes d'une spline où chaque ligne de la matrice correspond au développement de chaque coefficient , , et  :

La matrice principale est appelée matrice de base. À partir de ses colonnes, il est possible d'extraire les fonctions de base de Hermite permettant de réécrire les polynômes de la spline :

avec

Cette forme présente les polynômes indépendamment des valeurs et des pentes des tangentes aux nœuds. Elle est intéressante puisqu'elle permet d'optimiser les calculs de tracé de splines. Si on se donne un pas de tracé constant ce qui correspond à une subdivision fixe de l'intervalle [0 ; 1] de la variable , les valeurs des fonctions de base de Hermite pour chacun des pas ne seront à calculer qu'une seule fois pour un sous-intervalle , et pourront ensuite être réutilisées pour tous les autres sous-intervalles. La valeur d'interpolation n'étant ensuite qu'une somme de quatre multiplications.

L'implémentation ci-dessous de l'interpolation par spline cubique de Hermite n'intègre pas cette optimisation. Il s'agit simplement d'une transcription directe des formules.

function hermiteSpline(y0: number, y1: number, m0: number, m1: number, t: number): number {
var t2 = t * t;
var t3 = t2 * t;

// Calcul des fonctions de base d'une spline de Hermite
var b0 = 1 - 3 * t2 + 2 * t3;
var b1 = 3 * t2 - 2 * t3;
var b2 = t - 2 * t2 + t3;
var b3 = - t2+ t3;

return b0 * y0 + b1 * y1 + b2 * m0 + b3 * m1;
}

Implémentation de l'interpolation par spline cubique de Hermite

Les paramètres y0 et y1 correspondent aux valeurs des nœuds en et . Les paramètres m0 et m1 correspondent aux pentes respectives. Et le paramètre t compris entre 0 et 1 correspond à la position relative sur l'axe l'abscisse entre les deux nœuds concernés. La fonction hermiteSpline() renvoie l'interpolation par spline cubique de Hermite en t.

Les pentes et des deux tangentes peuvent ne pas être connues au préalable.

Il existe diverses façons de les définir en fonction du rendu final souhaité. Les principales sont la forme de Catmull-Rom et celle plus générale dite cardinale ou canonique.

8-3 - Spline de Catmull-Rom

Figure 14. Courbe représentative d'une spline de Catmull-Rom

Dans l'hypothèse où les nœuds sont répartis de façon plus ou moins homogène, une spline cubique de Catmull-Rom définie chaque tangente à un nœud comme parallèle à la droite passant par les nœuds et . Autrement dit, la pente de la tangente en chaque nœud s'exprime ainsi :

Figure 15. Tangentes d'une spline cubique de Catmull-Rom

Les paramètres étant déduits des nœuds, cela signifie que l'inclinaison locale de la spline suit la direction des nœuds avoisinants.

Dans la situation où les nœuds sont équirépartis, l'expression des pentes se simplifie et donne :

Ce qui permet d'obtenir la forme courante des polynômes d'une spline cubique de Catmull-Rom avec dans l'intervalle [0 ; 1] :

Ces polynômes peuvent s'écrire sous leur forme matricielle :

D'où l'implémentation suivante :

function catmullromSpline(y0: number, y1: number, y2: number, y3: number, t: number): number {
var t2 = t * t;
var t3 = t2 * t;

// Calcul des fonctions de base d'une spline de Catmull-Rom
var b0 = -t + 2 * t2 - t3;
var b1 = 2 - 5 * t2 + 3 * t3;
var b2 = t + 4 * t2 - 3 * t3;
var b3 = -t2 + t3;

return (b0 * y0) + (b1 * y1) + (b2 * y2) + (b3 * y3);
}

Implémentation de l'interpolation par spline cubique de Catmull-Rom

Les paramètres y0, y1, y2 et y3 correspondent respectivement aux valeurs des nœuds , , et . Le paramètre t correspond à l'abscisse relative entre les nœuds correspondant à y1 et y2. La fonction catmullromSpline() renvoie l'interpolation selon la spline de Catmull-Rom.

Comme précédemment avec les splines de Hermite, la spline de Catmull-Rom n'étant qu'un cas particulier, il est possible d'optimiser en pré-calculant les différentes valeurs des fonctions de base sur l'intervalle [0 ; 1].

Même si nous avons à présent l'algorithme de la spline de Catmull-Rom, il nous reste encore à traiter le cas des bords.

Pour le premier sous-intervalle [ ; ] (resp. le dernier sous-intervalle [ ; ]), la valeur de la pente (resp. ) ne peut pas se déduire de la formule de la pente car le nœud (resp ) n'existe pas. Pour se doter de nœuds imaginaires complémentaires, plusieurs conventions possibles peuvent alors être adoptées dont les suivantes :

- L'une des plus simples d'entre elles consiste à se donner un nœud imaginaire ayant la même valeur que le nœud voisin. Ce qui revient à définir :

- Une autre convention consiste à placer un nœud imaginaire (resp. ) sur la droite passant par les nœuds et (resp. et ), de telle sorte que le nœud (resp. ) soit le milieu du segment ayant pour extrémités les nœuds et (resp. et ). Ce qui revient à définir :

8-4 - Spline cardinale

Figure 16. Courbe représentative d'une spline cardinale de tension 0.5

Une spline cardinale, nommée également spline canonique, est une forme plus générale d'une spline de Catmull-Rom où est introduit un paramètre supplémentaire, la tension.

La pente d'une tangente en plus de dépendre des points voisins, dépend ici d'une tension notée  :

De façon imagée, la tension peut être assimilée à la longueur de la tangente. Plus la tension d'une tangente est élevée, et plus cette tangente se rapprochera de la perpendiculaire à la droite parallèle passant par les deux nœuds immédiatement voisins.

On constate aussi qu'une spline de Catmull-Rom n'est qu'une spline cardinale particulière où la tension est nulle.

En se plaçant dans le cas de nœuds équirépartis et en posant , on a :

Ce qui nous permet de déduire la matrice de base d'une spline cardinale à partir d'une spline de Hermite (cf. § 8.2) :

On déduit aisément de la matrice de base l'implémentation suivante en TypeScript.

function cardinalSpline(y0: number, y1: number, y2: number, y3: number, t: number, w: number): number {
var t2 = t * t;
var t3 = t2 * t;
var u = (1-w) * 0.5;

// Calcul des fonctions de base d'une spline cardinale
var b0 = -u * t + 2 * u * t2 - u * t3;
var b1 = 1 + (u - 3) * t2 + (2 - u) * t3;
var b2 = u * t  + (3 - 2 * u) * t2 + (u - 2) * t3;
var b3 = -u * t2 + u * t3;

return (b0 * y0) + (b1 * y1) + (b2 * y2) + (b3 * y3);
}

Implémentation de l'interpolation par spline cardinale

Les paramètres ici sont identiques à la fonction catmullromSpline() vu précédemment, à l'exception du paramètre w qui correspond à la tension aux tangentes. La fonction cardinalSpline() renvoie l'interpolation entre les nœuds correspondant à y1 et y2 suivant t.

8-5 - Spline cubique

Les splines cubiques présentées jusqu'à présent étaient des splines continûment dérivable une fois (). Cependant, le terme de spline cubique sous-entend en général une fonction de classe , ayant une dérivée seconde continue et non nulle.

Une continuité de classe peut être utile car non seulement plus esthétique, mais pouvant également satisfaire par exemple des contraintes physiques dans le domaine de la conception industrielle :

- la propriété d'une surface non plane à réfléchir la lumière de façon homogène (comme la carrosserie d'une voiture) ou bien,

- la discrétion et les qualités mécaniques d'une jonction entre deux objets assemblés.

Les splines cubiques de classe qu'on nommera simplement splines cubiques dans cette section, doivent donc vérifier les propriétés de continuité et de dérivabilité déjà évoquées à la section §8.1.

Étant donné nœuds, la spline cubique ayant pour morceaux les polynômes définis pour les sous-intervalles doit être deux fois continûment dérivable sur . Autrement dit les propriétés suivantes doivent être vérifiées :

P1. pour compris entre et et .

P2. pour compris entre et .

P3. pour compris entre et .

P4. pour compris entre et .

Chaque morceau de la spline est un polynôme de degré 3 et a donc 4 coefficients constants inconnus . Il y a donc inconnues à déterminer.

La propriété P1 correspond aux contraintes d'interpolations des nœuds. Elle fournit équations.

Les propriétés P2, P3 et P4 correspondent aux différents niveaux de continuité souhaités. Elles fournissent pour chacune d'elles équations et non pas équation à cause des deux nœuds aux extrémités en et qui ne sont pas concernés par la continuité.

Ces équations sous une forme développée donnent ceci :

Nous avons ainsi équations pour inconnues. Pour avoir une solution unique, il nous manque 2 équations concernant chacune des deux extrémités en et . Ces équations peuvent être définies selon l'aspect souhaité aux extrémités. On distingue généralement deux grandes familles de splines cubiques, même si de nombreuses variantes ont été imaginées :

- La spline naturelle (natural spline). A défaut de connaître la première ou la seconde dérivée aux extrémités, on peut poser et .

- La spline encastrée (clamped spline) lorsque la valeur de la dérivée première aux extrémités est non nulle.

Une fois qu'on a bien équations pour inconnues, on obtient un système linéaire pouvant s'exprimer sous la forme d'une matrice inversible.

Cependant, inverser naïvement une telle matrice, en utilisant par exemple la méthode du pivot de Gauss dont la complexité est en , n'est pas plus enviable que de résoudre un système de Vandermonde dans le cas de l'interpolation polynomiale globale.

En choisissant utilisant une telle matrice nous perdrions l'avantage semi-locale des splines puisque pour déterminer coefficients d'un seul morceau de la spline, nous devrions prendre en compte la totalité des nœuds.

Heureusement, il est possible de réarranger et simplifier ce système pour aboutir à un système tridiagonal qui peut être résolu en une complexité .

Les calculs qui vont suivre ont pour but d'expliciter l'obtention de ce système tridiagonal.

Chaque morceau étant un polynôme de degré 3, par conséquent la dérivé seconde est un polynôme linéaire de degré 1 qui peut s'écrire comme une interpolation linéaire lagrangienne entre et  :

Cette formule peut se réécrire en posant et , où désigne le « pas » qui est ici constant et non nul comme c'est le cas depuis le début de ce document où on se donne des points initiaux équirépartis. On obtient alors :

En intégrant deux fois la dérivée seconde pour revenir au polynôme interpolateur on obtient l'expression suivante :

Les coefficients et ne sont qu'une simplification de la notation, tandis que et sont les constantes d'intégration résolues en utilisant les conditions d'interpolation et (propriété P1).

Il ne nous manque que les valeurs pour connaître . Pour obtenir ces valeurs, nous pouvons dériver la nouvelle expression de  :

L'évaluation de et en donne :

En faisant appel à la propriété P3, on en déduit :

avec pour compris entre et .

De là, nous pouvons en déduire le système tridiagonal associé :