Question

Matlab 三次样条插值(有条件)

x=[0 1 2 3 4 5 6 7 8 9 10];
y=[0.0 0.79 1.53 2.19 2.71 3.03 3.27 2.89 3.06 3.19 3.29];
在x=0处y的二阶导数为0.8,在x=10处y的二阶导数为0.2,请问应该怎样进行三次样条插值？？

Answer 1

这个代码改一下：

function [yi,yxi,yxxi] = csinterp( x, y, xi )

% CSINTERP 1-D piecewise cubic spline interpolation of tabulated
% data and calculation of first and second derivatives.
% [YI,YXI,YXXI] = CSINTERP(X,Y,XI) interpolates function values,
% first and second derivatives using the rows/columns (X,Y) at
% the nodes specified in XI. X and Y should be row or column
% vectors, each containing at least three values, otherwise the
% function will return an empty output.
% The vector X specifies the points at which the data Y is given.
% Output data YI, YXI and YXXI are same size as XI.
% X and XI can be non-uniformly spaced, but repeated values in X
% will be removed and the resulting vector will be sorted.
% This function also extrapolates the input data (no NaNs are
% returned); use extrapolated values at your own risk.
%
% Requires bindex.m (available at Matlab FEX).
%
% Examples:
% x = [0.1:0.1:0.3 1:3]; y = sin( x );
% xi = linspace(0,3,200);
% [yi,yxi,yxxi] = csinterp(x,y,xi);
% figure(1)
% plot(x,y,'o',xi,yi,'r-'), box on, grid on
% xlabel('x')
% ylabel('闹信y(x)')
% title('Continuity test')
% legend('exact','csinterp'拍和)
%
% n = 101; m = n;
% x = linspace(0,2*pi,n);
% y = sin( x ); yx = cos( x ); yxx = -y;
% xi = 2*pi*rand(1,m);
% [yi,yxi,yxxi] = csinterp(x,y,xi);
% figure(2)
% subplot(311)
% plot(x,y,xi,yi,'o'), box on, grid on
% title('y(x)')
% legend('exact','csinterp')
% subplot(312)
% plot(x,yx,xi,yxi,'o'), box on, grid on
% title('dy/dx')
% legend('exact','csinterp'袭弯盯)
% subplot(313)
% plot(x,yxx,xi,yxxi,'o'), box on, grid on
% title('d^2y/dx^2')
% legend('exact','csinterp')
%
% n = 2001; m = 5001;
% x = linspace(0,2*pi,n);
% y = sin( x ); yx = cos( x ); yxx = -y;
% xi = sort( 2*pi*rand(1,m) );
% [yi,yxi,yxxi] = csinterp(x,y,xi);
% figure(3)
% subplot(311)
% plot(x,y,'o',xi,yi,'r-'), box on, grid on
% ylabel('y(x)')
% title('Large dataset test')
% legend('exact','csinterp')
% subplot(312)
% plot(x,yx,'o',xi,yxi,'r-'), box on, grid on
% ylabel('dy/dx')
% legend('exact','csinterp')
% subplot(313)
% plot(x,yxx,'o',xi,yxxi,'r-'), box on, grid on
% ylabel('d^2y/dx^2')
% legend('exact','csinterp')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% First version: 20/11/2007
%
% Contact: orodrig@ualg.pt
%
% Any suggestions to improve the performance of this
% code will be greatly appreciated.
%
% Reference:
% C. Moler, Numerical Computing with Matlab (NCM)
% http://www.mathworks.com/moler
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Let's go:

% Declare an empty output just in case calculations won't be
% performed:

yi = [];
yxi = [];
yxxi = [];

% Error checking:

sx = size(x);
sy = size(y);
sxi = size(xi);

nx = length(x);
ny = length(y);

if ( nargin ~= 3 )
error('This function requires three input arguments...');
end

if nx ~= ny
error('The length of x must match the length of y...')
end

if min(sx) ~= 1
error('Input x should be a row or column vector...')
end

if min(sy) ~= 1
error('Input y should be a row or column vector...')
end

% Convert x and y to column vectors to avoid problems
% with mixed row/column input data:

x = x(:);
y = y(:);

% Remove repeated elements:

[x,ui] = unique(x);

y = y(ui);

nx = length(x);

% Proceed with calculations only when the # of function points
% is greater than 2:

if nx > 2

yi = zeros( sxi );
yxi = zeros( sxi );
yxxi = zeros( sxi );

dy = diff( y );

h = diff( x ); hsq = h.*h;

h1 = h(1);
h2 = h(2);
h3 = h(nx-2);
h4 = h(nx-1);
hb = h(1:nx-2);
hf = h(2:nx-1);

% Bracketting:

ix = bindex(xi,x,1);

% Calculate the local variable:

if sxi(1) == 1
s = xi(:) - x(ix);
else
s = xi - x( ix );
end

% Calculate the slopes:

slopes = dy./h;

s1 = slopes(1);
s2 = slopes(2);
s3 = slopes(nx-2);
s4 = slopes(nx-1);
sb = slopes(1:nx-2);
sf = slopes(2:nx-1);

% Calculate the vectors filling the matrix of the system:

lodiag(1:nx-2) = hf;
lodiag( nx-1 ) = h3 + h4; l4 = lodiag( nx-1 );

cendiag( 1 ) = h2;
cendiag(2:nx-1) = 2*( hb + hf );
cendiag( nx ) = h3;

updiag( 1 ) = h1 + h2; u1 = updiag(1);
updiag(2:nx-1) = hb;

% Right side of the system (not-a-knot end conditions):

r1 = ( ( h1 + 2*u1 )*h2*s1 + h1^2*s2 )/u1;
rn = ( ( h4 + 2*l4 )*h3*s4 + h4^2*s3 )/l4;
r = 3*( hf.*sb + hb.*sf );
r = [r1;r;rn];

% Fill the matrix of the system (it is tridiagonal, so let us use a sparse matrix):

M = sparse( diag( cendiag ) + diag( lodiag, -1 ) + diag( updiag, 1 ) );

% Solve the linear system:

d = M\r;

db = d(1:nx-1);
df = d(2: nx );

% Chapter 3, page 14:

% Calculate the remaining spline coefficients

c = ( 3*slopes - 2*db - df )./h ;
b = ( db - 2*slopes + df )./hsq;

% Interpolate the function and its derivatives

s2 = s.*s;
s3 = s2.*s;

yi = y( ix ) + s.*d( ix ) + s2.*c( ix ) + s3.*b( ix );
yxi = d( ix ) + 2*s.*c( ix ) + 3*s2.*b( ix );
yxxi = 2*c( ix ) + 6*s.*b( ix );

if sxi(1) == 1
yi = yi(:)';
yxi = yxi(:)';
yxxi = yxxi(:)';
end

end