Splines is one of the topics in my numerical methods class this semester, so I was curious about the SimmSpline.
It is a cubic spline, which means that it is a piecewise cubic polynomial; a different polynomial in each interval between two knots.
A spline is smooth, which means that the function, and the 1st and 2nd derivatives are required to be continuous at the knots. SimmSpline is an interpolating spline (as opposed to a smoothing spline), so it is also required to pass exactly through the knots. These requirements are not sufficient to fully specify the polynomial coefficients. Two additional boundary conditions are needed.
The code suggested that SimmSpline estimates the 3rd derivative at the first and last knot, from finite differences on the knot data, and then requires that the first and last polynomial have these 3rd derivative values. I ran some tests and that is indeed the case. I had not seen this type of boundary condition before, but it makes some sense.
More common are the following types of boundary conditions:
- natural: second derivatives at the first and last knot are required to be zero
- clamped: first derivatives at the first and last knot are required to be zero (or some specified non-zero value)
- not-a-knot: third derivative is required to be continuous at the second knot and the second last knot
Equations for the natural spline can be found here:
https://en.wikipedia.org/wiki/Spline_interpolation.
The word "spline" comes from ship building where thin pieces of wood were used to draw smooth shapes. The natural spline is the mathematical representation of this. The endpoint second derivative is zero in a wooden spline because this minimizes the elastic energy due to bending. Here are some pictures:
https://www.boatdesign.net/threads/spli ... 244/page-2
I ran some tests on a simple example to compare all four splines. The extrapolation is linear in all cases.
Here is the code:
Code: Select all
% define knots
x = [-1 0 1 2 3 4 5 6];
y = [0 0 0 1 1 0 0 0];
% construct the SimmSpline
import org.opensim.modeling.*;
sspline = SimmSpline();
for i = 1:numel(x)
sspline.addPoint(x(i),y(i));
end
% evaluate the SimmSpline on 1000 points and extrapolate somewhat
xx = linspace(min(x)-1, max(x)+1, 1000)';
yy = zeros(numel(x),4); % we will store 4 spline curves
for i = 1:numel(xx)
yy(i,1) = sspline.calcValue(Vector(StdVectorDouble(xx(i))));
end
% use my own spline code, three different boundary conditions
p2 = makespline(x,y,'natural');
yy(:,2) = evalspline(p2,x,y,xx);
p3 = makespline(x,y,'clamped');
yy(:,3) = evalspline(p3,x,y,xx);
p4 = makespline(x,y,'not-a-knot');
yy(:,4) = evalspline(p4,x,y,xx);
% plot them all to compare
subplot(2,2,1)
plot(xx, yy, x,y,'o');
title('knots and splines')
legend('SimmSpline','natural','clamped','not-a-knot', 'data')
% use finite differences to get the derivatives of the splines)
% (we could do this analytically but this is OK)
yd = [NaN(1,4) ; diff(yy)./diff(xx)];
ydd = [diff(yd)./diff(xx) ; NaN(1,4)];
yddd = [diff(ydd)./diff(xx) ; NaN(1,4)];
% plot the derivatives
subplot(2,2,2)
plot(xx,yd);
title('1st derivative')
legend('SimmSpline','natural','clamped','not-a-knot')
subplot(2,2,3)
plot(xx,ydd);
title('2nd derivative')
subplot(2,2,4)
plot(xx,yddd);
title('3rd derivative')
ylim([-10 10])
% this code will confirm that SimmSpline uses 3rd derivative estimated from the data as a boundary condition
yd_data = [diff(y)./diff(x) NaN];
ydd_data = [NaN diff(yd_data)./diff(x)];
yddd_data = diff(ydd_data)./diff(x);
fprintf('3rd derivative estimated from first 4 data points: %6.2f\n', yddd_data(2))
fprintf('3rd derivative estimated from last 4 data points: %6.2f\n', yddd_data(end-1))
%=====================================================================
function p = makespline(x,y,bc)
% Makes an interpolating natural spline for data points x,y
%
% Inputs:
% x x values for the N knots
% y y values for the N knots
% bc boundary conditions: 'natural', 'clamped', 'not-a-knot'
%
% Outputs:
% p 3*(N-1) spline coefficients, stored as [p11 p21 p31 p12 p22 p32 ...]
% pij is the ith coefficient of the polynomial for the jth
% time interval
% The polynomial for time interval j (from x(j) to x(j+1) is:
% y = y(j) + p1j*(x-x(j)) + p2j*(x-x(j))^2 + p3j*(x-x(j))^3
N = numel(x); % determine how many knots there are
% make space to store the system of linear equations A*x = b
A = zeros(3*(N-1),3*(N-1)); % matrix which will be used to solve p
b = zeros(3*(N-1),1); % right hand side
% set up the 3*(N-1) linear equations
row = 0;
for i = 1:N-1
col = 3*(i-1)+1; % index of p1 for interval x(i) to x(i+1)
% equation 1: spline must go through x(i+1),y(i+1)
row = row+1;
A(row, col) = x(i+1)-x(i);
A(row, col+1) = (x(i+1)-x(i))^2;
A(row, col+2) = (x(i+1)-x(i))^3;
b(row) = y(i+1)-y(i);
% equations 2&3: 1st and 2nd derivative must be continous at x(i)
if (i > 1) % we can only do this at node 2 and higher
row = row+1;
% continuity of first derivative
A(row, col-3) = 1;
A(row, col-2) = 2*(x(i)-x(i-1));
A(row, col-1) = 3*(x(i)-x(i-1))^2;
A(row, col) = -1;
b(row) = 0;
% continuity of second derivative
row = row+1;
A(row, col-2) = 2;
A(row, col-1) = 6*(x(i)-x(i-1));
A(row, col+1) = -2;
b(row) = 0;
end
% for first or last interval, use boundary conditions specified in
% the input bc
if strcmp(bc, 'natural')
% at first interval, second derivative must be zero at x(i)
if (i == 1)
row = row+1;
A(row, col+1) = 2; % this means p2 must be zero
b(row) = 0;
end
% at last interval, second derivative must be zero at x(i+1)
if (i == (N-1))
row = row+1;
A(row, col+1) = 2; % this means p2+6*p3*(x(i+1)-x(i)) must be zero
A(row, col+2) = 6*(x(i+1)-x(i));
b(row) = 0;
end
elseif strcmp(bc, 'clamped')
% at first interval, first derivative must be zero at x(i)
if (i == 1)
row = row+1;
% this row says p1 must be zero for this interval
A(row, col) = 1;
b(row) = 0;
end
% at last interval, first derivative must be zero at x(i+1)
if (i == (N-1))
row = row+1;
% this row says p1 + 2*p2*((x(i+1)-x(i))) +
% 3*p3*(x(i+1)-x(i))^2 must be zero
A(row, col) = 1;
A(row, col+1) = 2*(x(i+1)-x(i));
A(row, col+2) = 3*(x(i+1)-x(i))^2;
b(row) = 0;
end
elseif strcmp(bc, 'not-a-knot')
% third derivative must be continuous at x(2)
if (i == 1)
row = row+1;
% this row says 6*p3 = 6*p3(from next interval)
A(row, col+2) = 1;
A(row, col+5) = -1;
b(row) = 0;
end
% third derivative must be continuous at x(N-1)
if (i == (N-1))
row = row+1;
% this row says 6*p3 = 6*p3(from next interval)
A(row, col-1) = 1;
A(row, col+2) = -1;
b(row) = 0;
end
end
end
% solve the polynomial coefficients
p = A\b;
end
%=======================================================================
function j = locate(x,xq)
% find the integer j, such that x(j) <= xq < x(j+1)
%
% Inputs:
% x............Array of knots (must be increasing)
% xq...........Query point
%
% Output:
% j............Index of knot interval where xq is found
% check that x values are increasing
if min(diff(x)) < 0
error('locate: knot array x is not increasing');
end
% check to see if xq is within the range defined by x
N = numel(x);
if xq < x(1)
j = 0;
elseif xq > x(N)
j = N;
else
% find the two surrounding knots by bisection
ilow = 1;
ihigh = N;
while (ihigh - ilow) > 1
imid = round((ihigh+ilow)/2);
if xq < x(imid)
ihigh = imid;
else
ilow = imid;
end
end
% we now know that xq is between x(ilow) and x(ilow+1)
j = ilow;
end
end
%=========================================================================
function yq = evalspline(p,x,y,xq)
% Evaluates a cubic spline at a series of query points xq
%
% Inputs:
% p 3*(N-1) spline coefficients, stored as [p11 p21 p31 p12 p22 p32 ...]
% pij is the ith coefficient of the polynomial for the jth
% time interval
% The polynomial for time interval j (from x(j) to x(j+1) is:
% y = y(j) + p1j*(x-x(j)) + p2j*(x-x(j))^2 + p3j*(x-x(j))^3
% x the N x-values at the knots (must be increasing)
% y the N y-values at the knots
% xq a list of query points where the spline must be evaluated
%
% Outputs:
% yq spline function values, evaluated at all elements of xq
N = numel(x); % the number of knots
Nq = numel(xq); % determine how many elements in xq
yq = zeros(size(xq)); % create space, same size as xq, to store yq
% go through all elements of xq, and compute yq for each of them
for iq = 1:Nq
% find the knot x(j) that is just below xq(iq)
j = locate(x,xq(iq));
% if j is outside of the intervals 1..N-1, do linear extrapolation
if j==0
slope = p(1);
yq(iq) = y(1) + slope*(xq(iq)-x(1));
elseif j==N
slope = p(3*(N-2)+1) + 2*p(3*(N-2)+2)*(x(N)-x(N-1)) + 3*p(3*(N-2)+3)*(x(N)-x(N-1))^2;
yq(iq) = y(N) + slope*(xq(iq)-x(N));
else
% use the cubic polynomial for interval j
% from p, extract the 3 polynomial coefficients for the
% interval x(i) to x(i+1)
p1 = p(3*(j-1) + 1);
p2 = p(3*(j-1) + 2);
p3 = p(3*(j-1) + 3);
% compute the polynomial value at xq(iq)
t = xq(iq) - x(j);
yq(iq) = y(j) + p1*t + p2*t^2 + p3*t^3;
end
end
end
And the results:
![untitled.png](./download/file.php?id=6685&sid=ab95af02c399dbbec73c40432754acb4)
- untitled.png (44.11 KiB) Viewed 753 times
The last plot shows that the SimmSpline has 3rd derivative values of +1 and -1 at the ends, which are the same as the result of a finite difference approximation on the knot data.
Conclusions:
- if you want to extrapolate, it's probably best to use a natural spline
- if you want to use the spline in the first or last knot interval, natural spline is safest also
- elsewhere, it does not matter much which spline is used
In OpenSim, SimmSpline is typically used to define a relationship between two variables. In OpenSim models, you often see extra points added before or after the region of interest. This is to reduce the influence of the boundary conditions, as the examples show.