L_SPLINE
=L_SPLINE(known_xs, known_ys, [x], [cond_start], [cond_end])
| Argument | Description | Example |
|---|---|---|
| known_xs | The x-values of the control points in increasing order | {2;4;8;10} |
| known_ys | The y-values of the control points | {3;1;5;5} |
| x | [optional] Values of x to interpolate | {2;2.1;...;9.9;10} |
| cond_start | [default="not"] Either "not", "free", or a specified slope at the start point | "not" |
| cond_end | [default="not"] Either "not", "free", or a specified slope at the end point | "not" |
Description
The L_SPLINE function creates an interpolating cubic spline that has a continuous first derivative and a continuous second derivative. It gives you the flexibility to choose either "not" (for the "not-a-knot" condition), "free" (for the "natural" condition), or to specify the value of the slope at each end point. The algorithm creates and solves a system of linear constraint equations and returns either the piecewise polynomial array structure, or the interpolated values corresponding to the optional x vector.
The default L_SPLINE(known_xs,known_ys) uses the "not-a-knot" end conditions.
L_SPLINE can be used to create a Natural Cubic Spline by specifying "free" for both cond_start and cond_end, but the L_NSPLINE function uses a more efficient algorithm.
Spline Constraints
I recommend watching the youtube series Splines in 5 Minutes by professor Steve Seitz. It does a very good job of describing the process of forming a system of constraint equations.
The 3 main constraints that define a C2 interpolating cubic spline are:
- The curve must pass through the control points
- The first derivative is continuous at the knots (a knot is the control point joining two polynomials)
- The second derivative is continuous at the knots
For a cubic spline created from n control points, these constraints produce 4(n-1)-2 equations, but we need two more constraints to have a fully constrained system (because we have 4(n-1) variables which are the coefficients of the n-1 polynomials). This is the purpose of the cond_start and cond_end parameters. These two parameters allow you to choose a constraint for the start and the end points. The constraints you choose will allow you to create the following common types of cubic splines:
- Natural Cubic Spline: Choose "free" for the start and end points. This sets the second derivative to zero at the ends. See L_NSPLINE for more information.
- Not-a-Knot Cubic Spline: Choose "not" for the start and end conditions. This forces the third derivative to be continuous at both the first knot and the last knot. It is the end condition used for the default spline function in Matlab. The first and second derivatives at the end points are not constrained.
- Clamped Cubic Spline: Enter the slope value for the start and end points (e.g. 0). Although "clamped" is a common name for this constraint, I do not like the term, because it implies that there is a direct relationship between an actual clamped beam and this end condition, which I don't think is accurate. In fact, the only constraints on the end are the slopes or first derivatives of the curve and the vertical position. There is no constraint that would restrict slipping (and therefore the overall length of the curve).
The Detailed Math
Each piece of the spline is a cubic polynomial \(p_i(x) = a_iu^3 + b_iu^2 + c_iu + d_i\) where \(u = (x-x_i)\) and is an index from 1 to n-1. This is a combination of nomenclature from the references below, but the main difference is the order of the coefficients. The system of constraint equations is simplified by scaling u to [0,1] by dividing by \({\Delta}_i=x_{i+1}-x_i\) to get \(t=(x-x_i)/{\Delta}_i\). This is shown in detail in the description of the L_CSPLINE function.
Scaling the x values in this way leads to the following definition for the polynomials and their derivatives:
$$ p_i(x) = q_i(t) = q_{i,1}t^3 + q_{i,2}t^2 + q_{i,3}t + q_{i,4} $$ $$ p'_i(x) = \frac{dq_i}{dt}\frac{dt}{dx} = \frac{1}{{\Delta}_i}\left(3q_{i,1}t^2 + 2q_{i,2}t + q_{i,3}\right) $$ $$ p''_i(x) = \frac{dq'_i}{dt}\frac{dt}{dx} = \frac{1}{{\Delta}_i^2}\left(6q_{i,1}t + 2q_{i,2}\right) $$ $$ p'''_i(x) = \frac{dq''_i}{dt}\frac{dt}{dx} = \frac{1}{{\Delta}_i^3}\left(6q_{i,1}\right) $$We will construct a system of linear equations to solve for the \(q_{i,j}\) coefficients, then we can calculate \(a_i=q_{i,1}/{\Delta}_i^3\) and \(b_i=q_{i,2}/{\Delta}_i^2\) and \(c_i=q_{i,1}/{\Delta}_i\).
The 3 main constraints above result in the following equations:
- (1a) Curve passes through \((x_i,y_i)\) : \( p_i(x_{i}) = q_i(0) = q_{i,4} = d_i = y_i \)
(1b) Curve passes through \((x_{i+1},y_{i+1})\) : \( p_i(x_{i+1}) = q_i(1) = q_{i,1}+q_{i,2}+q_{i,3}+q_{i,4} = y_{i+1} \) - (2) Continuous first derivative: \( p'_i(x_{i+1}) = p'_{i+1}(x_{i+1}) \) : \( q'_i(1) = q'_{i+1}(0) \) : \( 3q_{i,1}+2q_{i,2}+q_{i,3}-{\sigma}_iq_{i+1,3} = 0 \)
- (3) Continuous second derivative: \( p''_i(x_{i+1}) = p''_{i+1}(x_{i+1}) \) : \( q''_i(1) = q''_{i+1}(0) \) : \( 3q_{i,1}+q_{i,2}-{\sigma}_i^2q_{i+1,2} = 0 \)
where \({\sigma}_i = {\Delta}_i/{\Delta}_{i+1} = (x_{i+1}-x_i)/(x_{i+2}-x_{i+1}) \). The system of equations for 4 control points is shown below, with each row labeled with the reference to the corresponding constraint equation.
$$\begin{array}{c c} \begin{array} & (1a)_{i=1}: \\ (1b)_{i=1}: \\ (2)_{i=1}: \\ (3)_{i=1}: \\ (1a)_{i=2}: \\ (1b)_{i=2}: \\ (2)_{i=2}: \\ (3)_{i=2}: \\ (1a)_{i=3}: \\ (1b)_{i=3}: \\ ? \\ ? \end{array} \left[\begin{array}{cccc|cccc|cccc} 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 1 & 1 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 3 & 2 & 1 & 0 & 0 & 0 & -{\sigma}_1 & 0 & 0 & 0 & 0 & 0 \\ 3 & 1 & 0 & 0 & 0 & -({\sigma}_1)^2 & 0 & 0 & 0 & 0 & 0 & 0 \\ \hline 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 1 & 1 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 3 & 2 & 1 & 0 & 0 & 0 & -{\sigma}_2 & 0 \\ 0 & 0 & 0 & 0 & 3 & 1 & 0 & 0 & 0 & -({\sigma}_2)^2 & 0 & 0 \\ \hline 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 1 & 1 & 1 \\ ? & ? & ? & ? & ? & ? & ? & ? & ? & ? & ? & ? \\ ? & ? & ? & ? & ? & ? & ? & ? & ? & ? & ? & ? \end{array}\right] \end{array} \begin{bmatrix}q_{1,1} \\ q_{1,2} \\ q_{1,3} \\ q_{1,4} \\ \hline q_{2,1} \\ q_{2,2} \\ q_{2,3} \\ q_{2,4} \\ \hline q_{3,1} \\ q_{3,2} \\ q_{3,3} \\ q_{3,4} \end{bmatrix} = \begin{bmatrix}y_1 \\ y_2 \\ 0 \\ 0 \\ \hline y_2 \\ y_3 \\ 0 \\ 0 \\ \hline y_3 \\ y_4 \\ ? \\ ? \end{bmatrix}$$
Before we continue, we can first simplify the system by noting that constraint (1a) is trivial because \(q_{i,4}=y_i\). We can get rid of those rows and every 4th column by simple substitution. We end up with this:
$$\begin{array}{c c} \begin{array} & (1b)_{i=1}: \\ (2)_{i=1}: \\ (3)_{i=1}: \\ (1b)_{i=2}: \\ (2)_{i=2}: \\ (3)_{i=2}: \\ (1b)_{i=3}: \\ \text{cond_start} \\ \text{cond_end} \end{array} \left[\begin{array}{ccc|ccc|ccc} 1 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 3 & 2 & 1 & 0 & 0 & -{\sigma}_1 & 0 & 0 & 0 \\ 3 & 1 & 0 & 0 & -({\sigma}_1)^2 & 0 & 0 & 0 & 0 \\ \hline 0 & 0 & 0 & 1 & 1 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 3 & 2 & 1 & 0 & 0 & -{\sigma}_2 \\ 0 & 0 & 0 & 3 & 1 & 0 & 0 & -({\sigma}_2)^2 & 0 \\ \hline 0 & 0 & 0 & 0 & 0 & 0 & 1 & 1 & 1 \\ ? & ? & ? & ? & ? & ? & ? & ? & ? \\ ? & ? & ? & ? & ? & ? & ? & ? & ? \end{array}\right] \end{array} \begin{bmatrix}q_{1,1} \\ q_{1,2} \\ q_{1,3} \\ \hline q_{2,1} \\ q_{2,2} \\ q_{2,3} \\ \hline q_{3,1} \\ q_{3,2} \\ q_{3,3} \end{bmatrix} = \begin{bmatrix} y_2-y_1 \\ 0 \\ 0 \\ \hline y_3-y_2 \\ 0 \\ 0 \\ \hline y_4-y_3 \\ ? \\ ? \end{bmatrix}$$
The final two constraint equations come from our choice of end conditions. Below are the constraint equations for the three different options. The cond_start and cond_end parameters are independent, meaning that you could choose to define the slope for cond_start and set cond_end to "free" (or any other combination).
cond_start or cond_end = "free" (meaning "free to rotate", the second derivative is 0)
$$ \text{cond_start='free'} : p''_1(x_1)=0 \implies q''_1(0)=0 \implies q_{1,2} = 0 $$ $$ \text{cond_end='free'} : p''_{n-1}(x_n)=0 \implies q''_{n-1}(1)=0 \implies 3q_{n-1,1} + q_{n-1,2} = 0 $$cond_start or cond_end = "not" (referring to "not-a-knot", the third derivative is continuous at the first and/or last knots)
$$ \text{cond_start='not'} : p'''_1(x_2)=p'''_2(x_2) \implies q'''_1(1)=q'''_2(0) \implies q_{1,1} - {\sigma}_1^3q_{2,1} = 0 $$ $$ \text{cond_end='not'} : p'''_{n-2}(x_{n-1})=p'''_{n-1}(x_{n-1}) \implies q'''_{n-2}(1)=q'''_{n-1}(0) \implies q_{n-2,1} - {\sigma}_{n-2}^3q_{n-1,1} = 0 $$End Condition = slope specified (what some called "clamped")
$$ \text{cond_start}=m_1 : p'_1(x_1)=m_1 \implies q'_1(0)=m_1 \implies q_{1,3} = {\Delta}_1m_1 $$ $$ \text{cond_start}=m_n : p'_{n-1}(x_n)=m_n \implies q'_{n-1}(1)=m_2 \implies 3q_{n-1,1} + 2q_{n-1,2} + q_{n-1,3} = {\Delta}_{n-1}m_n $$After we add the desired combination of constraints, we solve the system Aq=B using q=MMULT(MINVERSE(A),B) and then convert the q coefficients to ai, bi, and ci. Using the MINVERSE function may not be an elegant solution given that the system of equations is already nearly tridiagonal, but it has worked so far for all the cases I have tested (which admittedly is not many). See reference [2] for another way to solve it.
Note: Please contact me if you see any typos or problems with the above equations.
Lambda Formula
This code for using L_SPLINE in Excel is provided under the License as part of the LAMBDA Library, but to use just this function, you may copy the following code directly into your spreadsheet.
Code for AFE Workbook Module (Excel Labs Add-in)
/**
* Creates a C2 Interpolating Spline with Specified End Conditions
* Must use at least 3 control points
*/
L_SPLINE = LAMBDA(known_xs,known_ys,[x],[cond_start],[cond_end],[debug],
IF(MAX(ROWS(known_ys),COLUMNS(known_ys))<3,NA(),
LET(doc,"https://www.vertex42.com/lambda/spline.html",
cond_start,IF(ISOMITTED(cond_start),"not",cond_start),
cond_end,IF(ISOMITTED(cond_end),"not",cond_end),
xs,IF(AND(ROWS(known_xs)=1,COLUMNS(known_xs)>1),TRANSPOSE(known_xs),known_xs),
ys,IF(AND(ROWS(known_ys)=1,COLUMNS(known_ys)>1),TRANSPOSE(known_ys),known_ys),
h,DROP(xs,1,0)-DROP(xs,-1,0),
si,DROP(h,-1,0)/DROP(h,1,0),
n,ROWS(xs),
dmat,{1,1,1;3,2,1;3,1,0},
constraints,REDUCE(SEQUENCE(3,3,0,0),SEQUENCE(n-2),LAMBDA(acc,i,LET(
smat,VSTACK({0,0,0},HSTACK(0,0,-INDEX(si,i,1)),HSTACK(0,-(INDEX(si,i,1)^2),0)),
IF(i=1,IF(n>3,HSTACK(dmat,smat,SEQUENCE(3,3*(n-3),0,0)),HSTACK(dmat,smat)),
IF(i=n-2,VSTACK(acc,HSTACK(SEQUENCE(3,3*(n-3),0,0),dmat,smat)),
VSTACK(acc,HSTACK(SEQUENCE(3,3*(i-1),0,0),dmat,smat,SEQUENCE(3,3*(n-2-i),0,0)))
))
))),
amatrix,VSTACK(constraints,
HSTACK(SEQUENCE(1,3*(n-2),0,0),{1,1,1}),
SWITCH(cond_start,
"free",HSTACK(0,1,0,SEQUENCE(1,3*(n-2),0,0)),
"not",HSTACK(1,0,0,-(INDEX(si,1,1)^3),SEQUENCE(1,3*(n-2)-1,0,0)),
HSTACK(0,0,1,SEQUENCE(1,3*(n-2),0,0))
),
SWITCH(cond_end,
"free",HSTACK(SEQUENCE(1,3*(n-2),0,0),3,1,0),
"not",IF(n>3,
HSTACK(SEQUENCE(1,3*(n-3),0,0),1,0,0,-(INDEX(si,n-2,1)^3),0,0),
HSTACK(0,0,0,3,1,0)
),
HSTACK(SEQUENCE(1,3*(n-2),0,0),3,2,1)
)
),
bmatrix,VSTACK(
REDUCE({0;0;0},SEQUENCE(n-2),LAMBDA(acc,i,IF(i=1,
VSTACK(INDEX(ys,2,1)-INDEX(ys,1,1),0,0),
VSTACK(acc,INDEX(ys,i+1,1)-INDEX(ys,i,1),0,0)
))),
INDEX(ys,n,1)-INDEX(ys,n-1,1),
SWITCH(cond_start,"free",0,"not",0,cond_start*INDEX(h,1,1)),
SWITCH(cond_end,"free",0,"not",0,cond_end*INDEX(h,n-1,1))
),
coeffs,HSTACK(WRAPROWS(MMULT(MINVERSE(amatrix),bmatrix),3)/(h^{3,2,1}),DROP(ys,-1,0)),
pp,MAP(HSTACK(xs,coeffs),LAMBDA(cell,IFERROR(cell,""))),
res,IF(debug=TRUE,HSTACK(amatrix,bmatrix),
IF(ISOMITTED(x),pp,
L_PPVAL(pp,x)
)
),
IF(AND(n=3,cond_start="not",cond_end="not"),"Error: redundant knot constraint",res)
)));
Named Function for Google Sheets
Name: L_SPLINE
Description: C2 Interpolating Spline with Specified End Conditions
Arguments:
Function:
[in the works]
L_SPLINE Examples
See Also
- [1] Wikipedia - Spline
- [2] Natural Spline Interpolation (PDF), Eberly, D., 2008.
- [3] Wikipedia - Spline Interpolation
- Matlab - spline function
