x=linspace(0,1,11);
y=sin(x);
xx=linspace(0,1,101);
yy = spline(x,y,xx);
plot(x,y,'*',xx,yy)
pp = spline(x,y);
pp
pp =
form: 'pp'
breaks: [1x11 double]
coefs: [10x4 double]
pieces: 10
order: 4
dim: 1
pp.coefs
ans =
-0.1652 -0.0003 1.0000 0
-0.1652 -0.0499 0.9950 0.0998
-0.1615 -0.0994 0.9801 0.1987
-0.1567 -0.1479 0.9553 0.2955
-0.1501 -0.1949 0.9211 0.3894
-0.1422 -0.2399 0.8776 0.4794
-0.1325 -0.2826 0.8253 0.5646
-0.1227 -0.3223 0.7648 0.6442
-0.1073 -0.3592 0.6967 0.7174
-0.1073 -0.3913 0.6216 0.7833
pp.breaks
ans =
Columns 1 through 5
0 0.1000 0.2000 0.3000 0.4000
Columns 6 through 10
0.5000 0.6000 0.7000 0.8000 0.9000
Column 11
1.0000
ppval(pp,xx)
ans =
Columns 1 through 5
0 0.0100 0.0200 0.0300 0.0400
Columns 6 through 10
0.0500 0.0600 0.0699 0.0799 0.0899
Columns 11 through 15
0.0998 0.1098 0.1197 0.1296 0.1395
Columns 16 through 20
0.1494 0.1593 0.1692 0.1790 0.1889
Columns 21 through 25
0.1987 0.2085 0.2182 0.2280 0.2377
Columns 26 through 30
0.2474 0.2571 0.2667 0.2764 0.2860
Columns 31 through 35
0.2955 0.3051 0.3146 0.3240 0.3335
Columns 36 through 40
0.3429 0.3523 0.3616 0.3709 0.3802
Columns 41 through 45
0.3894 0.3986 0.4078 0.4169 0.4259
Columns 46 through 50
0.4350 0.4439 0.4529 0.4618 0.4706
Columns 51 through 55
0.4794 0.4882 0.4969 0.5055 0.5141
Columns 56 through 60
0.5227 0.5312 0.5396 0.5480 0.5564
Columns 61 through 65
0.5646 0.5729 0.5810 0.5891 0.5972
Columns 66 through 70
0.6052 0.6131 0.6210 0.6288 0.6365
Columns 71 through 75
0.6442 0.6518 0.6594 0.6669 0.6743
Columns 76 through 80
0.6816 0.6889 0.6961 0.7033 0.7104
Columns 81 through 85
0.7174 0.7243 0.7311 0.7379 0.7446
Columns 86 through 90
0.7513 0.7578 0.7643 0.7707 0.7771
Columns 91 through 95
0.7833 0.7895 0.7956 0.8016 0.8076
Columns 96 through 100
0.8134 0.8192 0.8249 0.8305 0.8360
Column 101
0.8415
unmkpp(pp)
ans =
Columns 1 through 5
0 0.1000 0.2000 0.3000 0.4000
Columns 6 through 10
0.5000 0.6000 0.7000 0.8000 0.9000
Column 11
1.0000
[x,P]=unmkpp(pp)
x =
Columns 1 through 5
0 0.1000 0.2000 0.3000 0.4000
Columns 6 through 10
0.5000 0.6000 0.7000 0.8000 0.9000
Column 11
1.0000
P =
-0.1652 -0.0003 1.0000 0
-0.1652 -0.0499 0.9950 0.0998
-0.1615 -0.0994 0.9801 0.1987
-0.1567 -0.1479 0.9553 0.2955
-0.1501 -0.1949 0.9211 0.3894
-0.1422 -0.2399 0.8776 0.4794
-0.1325 -0.2826 0.8253 0.5646
-0.1227 -0.3223 0.7648 0.6442
-0.1073 -0.3592 0.6967 0.7174
-0.1073 -0.3913 0.6216 0.7833
help csape
CSAPE Cubic spline interpolation with various end conditions.
PP = CSAPE(X,Y) returns the cubic spline interpolant (in ppform) to the
given data (X,Y) using Lagrange end conditions (the default in table below).
The interpolant matches, at the data site X(j), the given data value
Y(:,j), j=1:length(X). The data values may be scalars, vectors, matrices,
or even ND-arrays. Data points with the same site are averaged.
For interpolation to gridded data, see below.
PP = CSAPE(X,Y,CONDS) uses the end conditions specified by CONDS, with
corresponding end condition values endcondvals .
If there are two more data values than data sites, then the first (last)
data value is taken as the value for the left (right) end condition, i.e.,
endcondvals = Y(:,[1 end]).
Otherwise, default values are used.
CONDS may be a *string* whose first character matches one of the
following: 'complete' or 'clamped', 'not-a-knot', 'periodic',
'second', 'variational', with the following meanings:
'complete' : match endslopes (as given, with
default as under *default*)
'not-a-knot' : make spline C^3 across first and last interior
break (ignoring given end condition values if any)
'periodic' : match first and second derivatives at first data
point with those at last data point
(ignoring given end condition values if any)
'second' : match end second derivatives (as given,
with default [0 0], i.e., as in variational)
'variational' : set end second derivatives equal to zero
(ignoring given end condition values if any)
The *default* : match endslopes to the slope of the cubic that
matches the first four data at the respective end.
By giving CONDS as a 1-by-2 matrix instead, it is possible to
specify *different* conditions at the two endpoints, namely
CONDS(i) with value endcondvals(:,i), with i=1 (i=2) referring to the
left (right) endpoint.
CONDS(i)=j means that the j-th derivative is being specified to
be endcondvals(:,i) , j=1,2. CONDS(1)=0=CONDS(2) means periodic end
conditions.
If CONDS(i) is not specified or is different from 0, 1 or 2, then
the default value for CONDS(i) is 1 and the default value of
endcondvals(:,i) is taken. If no end condition values are specified,
then the default value for endcondvals(:,i) is taken to be
deriv. of cubic interpolant to nearest four points, if CONDS(i)=1;
0 if CONDS(i)=2.
For backward compatibility, it is also possible to specify endcondvals as
an optional fourth input argument, but only for univariate data.
It is also possible to handle gridded data, by having X be a cell array
containing m univariate meshes and, correspondingly, having Y be an
m-dimensional array (or an (m+length(d))-dimensional array if the function
is to have values of size d ).
Correspondingly, CONDS is a cell array with m entries, and,
consistent with the tensor product procedure used, also the tensor product
of the univariate end conditions is enforced and must be supplied with
values; see the example below.
For example,
fnplt(csape( [0:4], [1 0 -1 0 1;0 1 0 -1 0], 'periodic')), axis equal
plots a circle, while
x = linspace(0,2*pi,21); pp = csape( x, [1 sin(x) 0], [1 2] );
gives a good approximation to the sine function on the interval [0 .. 2*pi]
(matching its slope 1 at the left endpoint, x(1) = 0, and its second
derivative 0 at the right endpoint, x(21) = 2*pi, in addition to its value
at every x(i), i=1:21).
As an illustration of the specification of end conditions in a multivariate
setting, here is complete bicubic interpolation, with the data explicitly
derived from the bicubic polynomial g(x,y) = x^3y^3, to make it easy for
you to see exactly where the slopes and slopes of slopes (i.e., cross
derivatives) must be placed in the data values supplied. Since our g is a
bicubic polynomial, its interpolant, f , must be g itself. We test this.
sites = {[0 1],[0 2]}; coefs = zeros(4,4); coefs(1,1) = 1;
g = ppmak(sites,coefs);
Dxg = fnval(fnder(g,[1 0]),sites);
Dyg = fnval(fnder(g,[0 1]),sites);
Dxyg = fnval(fnder(g,[1 1]),sites);
f = csape(sites,[Dxyg(1,1), Dxg(1,:), Dxyg(1,2); ...
Dyg(:,1), fnval(g,sites), Dyg(:,2) ; ...
Dxyg(2,1), Dxg(2,:), Dxyg(2,2)], ...
{'complete','complete'});
if any(squeeze(fnbrk(f,'c'))-coefs), 'something went wrong', end
As a multivariate vector-valued example, here is a sphere, done as a
parametric bicubic spline, using prescribed slopes in one direction and
periodic end conditions in the other:
x = 0:4; y=-2:2; s2 = 1/sqrt(2);
clear v
v(3,:,:) = [0 1 s2 0 -s2 -1 0].'*[1 1 1 1 1];
v(2,:,:) = [1 0 s2 1 s2 0 -1].'*[0 1 0 -1 0];
v(1,:,:) = [1 0 s2 1 s2 0 -1].'*[1 0 -1 0 1];
sph = csape({x,y},v,{'clamped','periodic'});
values = fnval(sph,{0:.1:4,-2:.1:2});
surf(squeeze(values(1,:,:)),squeeze(values(2,:,:)),squeeze(values(3,:,:)))
% the previous two lines could have been replaced by: fnplt(sph)
axis equal, axis off
See also csapi, spapi, spline.
Reference page in Help browser
doc csape
x=[1,2,3,4];
y=[2,4,5,7];
xx=linspace(1,4,100);
yy=spline(x,y,xx);
plot(xx,yy,x,y,'*')
t=[1,2,3,4];
x=[2,3,4,6,4,3];
y=[1,0,-1,0,2,4];
plot(x,y,'*')
tt=linspace(1,4,100);
xx=spline(t,x,tt);
yy=spline(t,y,tt);
plot(xx,yy)
plot(xx,yy,x,y,'*')
x=[2,3,4,6,4,3,2];
y=[1,0,-1,0,2,4,1];
xx=spline(t,x,tt);
??? Error using ==> chckxy at 89
The number of sites, 4, is incompatible with the number
of values, 7.
Error in ==> spline at 55
[x,y,sizey,endslopes] = chckxy(x,y);
t=[1,2,3,4,5,6];
x=[2,3,4,6,4,3,2];
t=[1,2,3,4,5,6,7];
x=[2,3,4,6,4,3,2];
y=[1,0,-1,0,2,4,1];
xx=spline(t,x,tt);
yy=spline(t,y,tt);
plot(xx,yy,x,y,'*')
tt=linspace(1,7,100);
xx=spline(t,x,tt);
yy=spline(t,y,tt);
plot(xx,yy,x,y,'*')
x = [1,2,3,4];
y = sin(x);
pp = spline(x,y);
f = @(x) ppval(pp,x);
quad(f,1,4)
ans =
1.2135
diary off