98
3 Approximation de fonctions et de données
Programme 3.1. cubicspline : spline d’interpolation cubique
function s= c u b ic splin e(x ,y ,zi , type , der)
% C U B I CSP LINE calcule une spline cubique
% S = C U B I CSPLI NE(X ,Y , ZI ) calcule la valeur aux a b s c isses
% ZI de la spline d ’ i n t e r pola tion cubique n a t ur elle qui
% i n t erpo le les valeurs Y aux noeuds X.
% S = C U B I CSPLI NE(X ,Y , ZI , TYPE , DER) si TYPE =0 calcule la
% valeur aux a b s cisse s ZI de la spline cubique
% i n t e rpola nt les valeurs Y et dont la dérivée
% première aux e x t rém ités vaut DER (1) et DER (2).
% Si TYPE =1 alors DER (1) et DER (2) sont les valeurs de
% la dérivée seconde aux e x t r émités.
[n , m ]= size ( x );
if n == 1
x = x’;
y = y’;
n = m;
end
if nargin == 3
der0 = 0; dern = 0; type = 1;
else
der0 = der (1); dern = der (2);
end
h = x (2: end ) -x (1: end -1);
e = 2*[ h (1); h (1: end -1)+ h (2: end ); h ( end )];
A = spdiags ([[h ; 0] e [0; h ]] , -1:1 ,n , n );
d = ( y (2: end) -y (1: end -1))./ h ;
rhs = 3*( d (2: end ) -d (1: end -1));
if type == 0
A (1 ,1) = 2*h (1);
A (1 ,2) = h (1);
A (n ,n ) = 2*h ( end ); A ( end , end -1) = h ( end );
rhs = [3*(d (1) -der0 ); rhs; 3*( dern - d( end ))];
else
A (1 ,:) = 0; A (1 ,1) = 1;
A (n ,:) = 0; A(n , n) = 1;
rhs = [ der0 ; rhs ; dern ];
end
S = zeros (n ,4);
S (: ,3) = A\ rhs;
for m = 1:n -1
S (m ,4) = (S ( m +1 ,3) -S (m ,3))/3/ h( m );
S (m ,2) = d( m ) - h( m )/3*( S( m + 1 ,3)+2*S(m ,3));
S (m ,1) = y( m );
end
S = S (1:n -1 , 4: -1:1);
pp = mkpp (x ,S ); s = ppval (pp , zi );
return
La commande MATLAB spline (voir aussi la toolbox splines)
spline
force la dérivée troisième de s 3 à être continue en x 1 et x n−1 . On donne
à cette condition le nom curieux de condition not-a-knot. Les paramètres
d’entrée sont les vecteurs x, y et le vecteur zi (ayant la même signification que précédemment). Les commandes mkpp et ppval utilisées dans le
mkpp
ppval Programme 3.1 servent à construire et évaluer un polynôme composite.
Exemple 3.8 Considérons à nouveau les données de la Table 3.1 correspondant à la colonne K = 0.67 et calculons la spline cubique associée s3. Les
3 Approximation de fonctions et de données
Programme 3.1. cubicspline : spline d’interpolation cubique
function s= c u b ic splin e(x ,y ,zi , type , der)
% C U B I CSP LINE calcule une spline cubique
% S = C U B I CSPLI NE(X ,Y , ZI ) calcule la valeur aux a b s c isses
% ZI de la spline d ’ i n t e r pola tion cubique n a t ur elle qui
% i n t erpo le les valeurs Y aux noeuds X.
% S = C U B I CSPLI NE(X ,Y , ZI , TYPE , DER) si TYPE =0 calcule la
% valeur aux a b s cisse s ZI de la spline cubique
% i n t e rpola nt les valeurs Y et dont la dérivée
% première aux e x t rém ités vaut DER (1) et DER (2).
% Si TYPE =1 alors DER (1) et DER (2) sont les valeurs de
% la dérivée seconde aux e x t r émités.
[n , m ]= size ( x );
if n == 1
x = x’;
y = y’;
n = m;
end
if nargin == 3
der0 = 0; dern = 0; type = 1;
else
der0 = der (1); dern = der (2);
end
h = x (2: end ) -x (1: end -1);
e = 2*[ h (1); h (1: end -1)+ h (2: end ); h ( end )];
A = spdiags ([[h ; 0] e [0; h ]] , -1:1 ,n , n );
d = ( y (2: end) -y (1: end -1))./ h ;
rhs = 3*( d (2: end ) -d (1: end -1));
if type == 0
A (1 ,1) = 2*h (1);
A (1 ,2) = h (1);
A (n ,n ) = 2*h ( end ); A ( end , end -1) = h ( end );
rhs = [3*(d (1) -der0 ); rhs; 3*( dern - d( end ))];
else
A (1 ,:) = 0; A (1 ,1) = 1;
A (n ,:) = 0; A(n , n) = 1;
rhs = [ der0 ; rhs ; dern ];
end
S = zeros (n ,4);
S (: ,3) = A\ rhs;
for m = 1:n -1
S (m ,4) = (S ( m +1 ,3) -S (m ,3))/3/ h( m );
S (m ,2) = d( m ) - h( m )/3*( S( m + 1 ,3)+2*S(m ,3));
S (m ,1) = y( m );
end
S = S (1:n -1 , 4: -1:1);
pp = mkpp (x ,S ); s = ppval (pp , zi );
return
La commande MATLAB spline (voir aussi la toolbox splines)
spline
force la dérivée troisième de s 3 à être continue en x 1 et x n−1 . On donne
à cette condition le nom curieux de condition not-a-knot. Les paramètres
d’entrée sont les vecteurs x, y et le vecteur zi (ayant la même signification que précédemment). Les commandes mkpp et ppval utilisées dans le
mkpp
ppval Programme 3.1 servent à construire et évaluer un polynôme composite.
Exemple 3.8 Considérons à nouveau les données de la Table 3.1 correspondant à la colonne K = 0.67 et calculons la spline cubique associée s3. Les