Interpolations polynomiales et cubic splines

Sylvain Barthélémy, Avril 1998

Vous pouvez tlcharger la librairie qui vous permettra d'exprimenter les mthodes prsentes sur cette page sur GitHub, l'adresse suivante https://github.com/sylbarth/interp


Interpolations de Lagrange (1735-1796)

Selon le thorme de Weiestrass' (et la preuve de Benstein) :

Si f(x) est continue sur l'intervalle [a,b] alors il quelque soit e>0 rel positif, il existe un entier n et un polynome pn(x) tels que |f(x)-pn(x)| < e

Supposons que vous disposez de N+1 points initiaux et que vous dsirez effectuer une interpolation polynomiale sur(entre) ces points. La mthode de Lagrange propose de choisir N points parmis les N+1 et de rsoudre le systme de N quations N inconnues form par les quation polynomiales passant par ces points. Si vous disposez de 4 points, vous rechercherez donc un polynome de degr 3 passant par ces points.

double lagrange(const double& u, const double* x, const double* y, const int& n)
{
  double interp = 0.0;

  for(int i=0; i < n; i++) {
    double p = 1.0;
    for(int j=0; j < n; j++) {
      if(i!=j) p *= (u-x[j])/(x[i]-x[j]);
    }
    interp += p * y[i];
  }

  return interp;
}

Le problme pos par cette mthode est qu'elle ne permet pas de choisir le degr du polynome. Un degr trop lev engendrera des oscillations excessives et inversement. Les mthodes travaillant sur les diffrences et permettent de corriger ce dfaut en choisissant un seuil de convergence "acceptable".


Mthodes des diffrences

La mthode standard des diffrences peut donner le mme rsultat que la mthode de Lagrange, s'il n'est pas possible de trouver un degr polynomiale infrieur N sans perte de prcision excessive. Dans certains cas elle remplacera cependant avantageusement l'interpolation de Lagrange en ajustant un polynome de degr infrieur.

void divdiff_coef(double* dd, const double* x, const double* y, const int& n)
{
  double temp1, temp2;
  for ( int i=0; i < n; i++ )
    dd[i] = y[i];
  for ( int j=1; j < n; j++ ) { 
    temp1 = dd[j-1]; 
    for ( int k=j; k < n; k++ ) {
      temp2 = dd[k];
      dd[k] = (dd[k] - temp1)/(x[k] - x[k-j]);
      temp1 = temp2;   
    }
  }
}

double divdiff(const double& u, const double* dd, const double* x, const int& n)
{
  double interp = 0.0;
   for ( int i=(n-1); i>=1; i-- )
      interp = (interp + dd[i]) * (u - x[i-1]);
   interp += dd[0];
  return interp;
}

Cubic Splines

A la diffrence de la mthode prcdente, le cubic spline n'est pas une mthode d'interpolation polynomiale proprement parler. Nous cherchons plutot lisser une srie de points (un peu comme on le ferai avec une moyenne mobile) et non pas trouver quel polynome reprsente le mieux ces points. Effectuer une interpolation par un cubic spline revient rsoudre un systme tridiagonal. Comme pour la mthode des diffrences, on recherche tout d'abord les valeurs des coefficients avant de calculer les points interpols.

L'algorithme utilis pour rechercher un cubic spline sur une srie de points est un peu plus compliqu que les deux prcdents. Si vous tre intrss par une implmentation en C++, vous pouvez tlcharger enl_interp.zip. Ce code utilise la rsolution d'une matrice tridiagonale propose par les fameux Numerical Recipies in C.


Comparaison

Pour dmontrer les diffrences qui peuvent exister entre une interpolation polynomiale et par les splines, nous avons gnr 9 points non quidistants de facon (pseudo) alatoire. Le cubic spline donne des rsultats assez differents de l'interpolation de Lagrange, avec un aspect bien plus liss (courbe bleue).


Comme vous pourez le constater sur le tableau suivant, la mthode de Lagrange donne bien les mmes rsultats que la mthode des diffrences.

uLagrangeDivDiffCSpline
0.000.01790.01790.0179
0.031.57181.57180.1476
0.052.04892.04890.2743
0.081.97251.97250.3954
0.101.67341.67340.5079
0.121.34281.34280.6090
0.151.07511.07510.6957
0.170.90130.90130.7653
0.200.81490.81490.8149
0.220.79030.79030.8422
0.250.79660.79660.8479
0.270.80610.80610.8330
0.300.79880.79880.7988
0.330.76500.76500.7472
0.350.70440.70440.6830
0.380.62480.62480.6119
0.400.53930.53930.5393
0.430.46320.46320.4720
0.450.41130.41130.4213
0.480.39480.39480.3997
0.500.41960.41960.4196
0.530.48510.48510.4881
0.550.58400.58400.5913
0.580.70290.70290.7097
0.600.82410.82410.8241
0.630.92830.92830.9172
0.650.99690.99690.9793
0.681.01641.01641.0031
0.700.98080.98080.9808
0.730.89470.89470.9094
0.750.77480.77480.8030
0.780.65000.65000.6801
0.800.55920.55920.5592
0.830.54520.54520.4552
0.850.64580.64580.3692
0.880.87890.87890.2985
0.901.22271.22270.2407
0.931.58741.58740.1930
0.951.78061.78060.1531
0.981.46171.46170.1183
1.000.08600.08600.0860

θ  π