Spline

Soit X_i et Y_i les abscisses et ordonnées de n nœuds, les abscisses étant classées par ordre croissant.

1. Calculer la matrice colonne h des largeurs des n–1 intervalles entre nœuds :

${\displaystyle \mathbf {h} _{i}=X_{i+1}-X_{i}}$, pour i de 1 à n–1.

2. Calculer la matrice colonne F des n–2 dérivées secondes aux nœuds :

2.1. Pour les nœuds 2 à n–1 :

${\displaystyle \mathbf {F} _{i}={\frac {Y_{i+1}-Y_{i}}{\mathbf {h} _{i}}}-{\frac {Y_{i}-Y_{i-1}}{\mathbf {h} _{i-1}}}}$, pour i de 2 à n–1.

2.2. Il faut de plus valoriser les dérivées secondes aux extrémités, F₁ et F_n.

Si l'on pose F₁ = 0 et F_n = 0, on obtient une spline naturelle ;

si l'on possède des informations sur les extrémités, on affecte à F₁ et/ou à F_n la valeur ad hoc.

3. Construire la matrice carrée R_i,j, de dimensions n×n, telle que, i désigne l'indice des lignes et j celui des colonnes :

3.1 Mettre tous les éléments de la matrice R à 0.

3.2 Sur la première ligne, seul le premier élément n'est pas nul : R_1,1 = 1.

3.3 Sur la dernière ligne, seul le dernier élément n'est pas nul : R_n,n = 1.

3.4 La sous-matrice de R de dimension (n–2) × n en excluant la première et la dernière ligne est tridiagonale. Tous ses éléments sont nuls, sauf :

• La diagonale centrale de la sous-matrice vaut :

${\displaystyle \mathbf {R} _{i,i}={\frac {h_{i-1}+h_{i}}{3}}}$ pour i = 2, ..., n–1.

• La diagonale supérieure vaut :

${\displaystyle \mathbf {R} _{i,{i+1}}={\frac {h_{i}}{6}}}$ pour i = 2, ..., n–1.

• La diagonale inférieure vaut :

${\displaystyle \mathbf {R} _{i,{i-1}}={\frac {h_{i-1}}{6}}}$ pour i = 2, ..., n–1.

4. Calculer la matrice colonne M de dimension n, telle que :

${\displaystyle \mathbf {M} =\mathbf {R} ^{-1}\times \mathbf {F} }$.

5. Calculer les matrices colonnes C et C', de dimensions n–1, contenant respectivement les constantes d'intégration d'ordre 1 et 2 des n–1 polynômes, de la façon suivante :

Pour i de 1 à n–1 :

${\displaystyle \mathbf {C} _{i}={\frac {Y_{i+1}-Y_{i}}{\mathbf {h} _{i}}}-{\frac {\mathbf {h} _{i}}{6}}(\mathbf {M} _{i+1}-\mathbf {M} _{i})}$

${\displaystyle \mathbf {C'} _{i}=Y_{i}-\mathrm {M} _{i}{\frac {\mathbf {h} _{i}^{2}}{6}}}$

6. Calcul d'une valeur d’interpolation pour une abscisse x située dans l'intervalle k, entre les nœuds d'abscisses X_k et X_k+1 (avec k de 1 à n–1) :

${\displaystyle {\begin{array}{lcl}y&=&\mathbf {M} _{k}{\frac {(X_{k+1}-x)^{3}}{6\mathbf {h} _{k}}}\\&+&\mathbf {M} _{k+1}{\frac {(x-X_{k})^{3}}{6\mathbf {h} _{k}}}\\&+&\mathbf {C} _{k}(x-X_{k})\\&+&\mathbf {C'} _{k}\end{array}}}$

Exemple

Cet exemple présente l'approximation de la courbe y = sin x par une spline cubique naturelle pour laquelle on donne les cinq valeurs suivantes (5 nœuds) de la fonction :

${\displaystyle n=5}$

${\displaystyle {\begin{array}{rrr}i&X_{i}&Y_{i}\\1&0{,}0&0{,}0000\\2&0{,}5&0{,}4794\\3&1{,}0&0{,}8415\\4&1{,}5&0{,}9975\\5&2{,}0&0{,}9093\\\end{array}}}$

${\displaystyle \mathbf {h} ={\begin{bmatrix}0{,}5\\0{,}5\\0{,}5\\0{,}5\\\end{bmatrix}}{\text{ ; }}}$${\displaystyle \mathbf {F} =\left[{\begin{array}{r}0{,}0000\\-0{,}2348\\-0{,}4120\\-0{,}4884\\0{,}0000\\\end{array}}\right]{\text{ ; }}}$${\displaystyle \mathbf {R} =\left[{\begin{array}{rrrrr}1{,}0000&0&0&0&0\\0{,}0833&0{,}3333&0{,}0833&0&0\\0&0{,}0833&0{,}3333&0{,}0833&0\\0&0&0{,}0833&0{,}3333&0{,}0833\\0&0&0&0&1{,}0000\\\end{array}}\right]{\text{ ; }}}$

${\displaystyle \mathbf {M} =\left[{\begin{array}{r}0{,}0000\\-0{,}5061\\-0{,}7928\\-1{,}2671\\0{,}0000\\\end{array}}\right]{\text{ ; }}}$${\displaystyle \mathbf {C} =\left[{\begin{array}{r}1{,}0010\\0{,}7480\\0{,}3516\\-0{,}2820\\\end{array}}\right]{\text{ ; }}}$${\displaystyle \mathbf {C'} =\left[{\begin{array}{r}0{,}0000\\0{,}5005\\0{,}8745\\1{,}0503\\\end{array}}\right]}$

Commentaire sur l'exemple : La fonction ${\displaystyle sin(i)}$ est bien représentée par son approximation par une spline naturelle à gauche du graphique parce que la dérivée seconde de ${\displaystyle sin(i)}$ est bien nulle pour ${\displaystyle i=0}$, comme le suppose le calcul de la spline naturelle. En revanche, à droite du graphique, l'approximation par une spline naturelle est moins exacte, en effet, la dérivée seconde de ${\displaystyle sin(i)}$ n'est pas nulle pour ${\displaystyle X_{5}=2{,}0}$ (elle vaut ${\displaystyle -0,035}$), valeur différente de l'hypothèse de nullité posée pour la détermination de la spline naturelle.

Algorithme de calcul d'une spline cubique de lissage^[2].

Soit X_i et Y_i les abscisses et ordonnées de n nœuds, les abscisses étant classées par ordre croissant (n > 2).

1. Calculer la matrice colonne h des largeurs des n–1 intervalles entre nœuds :

${\displaystyle \mathbf {h} _{i}=X_{i+1}-X_{i}}$, pour i de 1 à n–1.

2. Construire la matrice tridiagonale Q (n × (n–2)) telle que toutes les valeurs de la matrice sont nulles, sauf :

• La diagonale centrale :

${\displaystyle \mathbf {Q} _{i,{i-1}}={\frac {1}{h_{i-1}}}+{\frac {1}{h_{i}}}}$, pour i = 2 à n–1 ;

• La diagonale supérieure :

${\displaystyle \mathbf {Q} _{i,i}=-{\frac {1}{h_{i}}}}$, pour i = 1 à n–2 ;

• La diagonale inférieure :

${\displaystyle \mathbf {Q} _{i,{i-2}}=-{\frac {1}{h_{i-1}}}}$, pour i = 3 à n ;

3. Construire la matrice carrée tridiagonale R ((n–2)×(n–2)) telle que toutes les valeurs de la matrice sont nulles, sauf :

• La diagonale principale :

${\displaystyle \mathbf {R} _{i,i}={\frac {h_{i}+h_{i+1}}{3}}}$, pour i = 1 à n–2 ;

• La diagonale supérieure :

${\displaystyle \mathbf {R} {i,{i+1}}=-{\frac {h_{i}}{6}}}$ pour i = 1 à n–3 ;

• La diagonale inférieure :

${\displaystyle \mathbf {R} {i,{i-1}}=-{\frac {h_{i-1}}{6}}}$ pour i = 2 à n–2 ;

4. Calculer la matrice carrée K (n × n) :

${\displaystyle \mathbf {K} =\mathbf {Q} \mathbf {R} ^{-1}\mathbf {Q} ^{t}}$;

5. Soit ${\displaystyle \mathbf {I} }$ la matrice unitaire n × n et ${\displaystyle \lambda }$ le paramètre de lissage.

On obtient les ordonnées s de la spline de lissage en chacun des n nœud d'abscisse X et d'ordonnée Y pour la valeur ${\displaystyle \lambda }$ du paramètre de lissage par :

${\displaystyle \mathbf {s} =(\mathbf {I} +\lambda \mathbf {K} )^{-1}\mathbf {Y} }$.

6. On peut interpoler la valeur de la spline de lissage pour toute abscisse comprise entre les nœuds extrêmes X₁ et X_n en appliquant l'algorithme d'interpolation des splines aux points d'abscisse X et d'ordonnée s.

Exemple :

Cet exemple présente le lissage des quatre nœuds de coordonnées (X_i, Y_i), figurés en rouge dans le graphique ci-dessus, par une spline cubique naturelle avec le paramètre de lissage λ = 1,5.

${\displaystyle {\begin{array}{rrr}i&X_{i}&Y_{i}\\1&0{,}0&0{,}0\\2&1{,}0&6{,}0\\3&3{,}0&5{,}0\\4&5{,}0&0{,}0\\\end{array}}}$ ;

${\displaystyle \mathbf {h} ={\begin{bmatrix}1{,}0\\2{,}0\\2{,}0\\\end{bmatrix}}{\text{ ; }}}$${\displaystyle \mathbf {Q} ={\begin{bmatrix}{\begin{array}{rr}-1{,}0000&0{,}0000\\1{,}5000&-0{,}5000\\-0{,}5000&1{,}0000\\0{,}0000&-0{,}5000\\\end{array}}\end{bmatrix}}{\text{ ; }}}$${\displaystyle \mathbf {R} ={\begin{bmatrix}{\begin{array}{rr}1{,}0000&-0{,}1667\\-0{,}1667&1{,}3333\\\end{array}}\end{bmatrix}}}$ ;

${\displaystyle \mathbf {K} =\mathbf {Q} \mathbf {R} ^{-1}\mathbf {Q} ^{t}={\begin{bmatrix}{\begin{array}{rrrr}1{,}0213&-1{,}4681&0{,}3830&0{,}0638\\-1{,}4681&2{,}2979&-0{,}9255&0{,}0957\\0{,}3830&-0{,}9255&0{,}8936&-0{,}3511\\0{,}0638&0{,}0957&-0{,}3511&0{,}1915\\\end{array}}\end{bmatrix}}}$ ;

I étant la matrice unité n × n et avec le paramètre λ, on calcule :

${\displaystyle (\mathbf {I} +\lambda \mathbf {K} )={\begin{bmatrix}{\begin{array}{rrrr}2.53191&-2.20213&0.57447&0.09574\\-2.20213&4.44681&-1.38830&0.14362\\0.57447&-1.38830&2.34043&-0.52660\\0.09574&0.14362&-0.52660&1.28723\\\end{array}}\end{bmatrix}}{\text{, et : }}}$${\displaystyle (\mathbf {I} +\lambda \mathbf {K} )^{-1}={\begin{bmatrix}{\begin{array}{rrrr}0.7051&0.3583&0.0206&-0.0840\\0.3583&0.4599&0.1843&-0.0026\\0.0206&0.1843&0.5800&0.2152\\-0.0840&-0.0026&0.2152&0.8714\\\end{array}}\end{bmatrix}}}$.

d'où l'on calcule ${\displaystyle \mathbf {s} =(\mathbf {I} +\lambda \mathbf {K} )^{-1}\mathbf {Y} }$