I'm building a package for B-spline interpolation in Julia, and I've come across a boundary condition that I want to implement but can't wrap my head around how to do it (mathematically).
Basically, it boils down to a more general problem, which I suppose comes up in discretizing PDE's and such:
Given a system of equations you usually close by introducing "ghost cells" and boundary conditions, how do you close it instead by considering the outermost data points your boundary conditions (i.e. only solve the problem on the interior domain)?
A more formal definition of the question follows after…
Some background and examples
Introduction to B-Splines (the way I've treated them)
A b-spline interpolant on a uniform grid $i\in [1,N] \cap \mathbb N$ allows us to evaluate the interpolant at any $x\in[1,N]$, by defining a piecewise function $f_i(x) = \sum_{k}c_{i+k}p_k(x-i)$. As an example, for a cubic b-spline it turns out that $p_k(x – i) \equiv 0\,\,\, \forall\,\, x : x-i\notin[-1,2]$, so each piece is actually just a sum over four terms:
$$
f_i(x) = c_{i-1} p_{-1}(x-i) + c_i p_0(x-i) + c_{i+1}p_1(x-i) + c_{i+2}p_2(x-i)
$$
and each piece $f_i$ is valid on $x\in[i,i+1]$. To find $c_i$ and $p_k$ one requires continuity in function value and a number of derivatives (first and second order for cubic b-splines) as well as that the basis functions are symmetric (i.e. that $p_{-1}(x-i) = p_2(1-(x-i))$ and $p_0(x-i) = p_1(1-(x-i))$). For my question, the actual functions $p_k$ are irrelevant; I'm currently interested in finding $c_i$, which is done by applying these continuity requirements on each "hand-over" from $f_i$ to $f_{i+1}$ (e.g. $f_2(3) = f_3(3)$, $f'_2(3) = f'_3(3)$; $f''_2(3) = f''_3(3)$ for cubic b-splines).
It turns out that for quadratic or higher b-splines, it is necessary to solve an equation system on this form [1] to find them
$$
M_{ij} c_j = v_i + b_i
$$
$c_j$ are the sought coefficients, $v_i$ are the data points that we wish our interpolant to pass through, and $M_{ij}$ and $b_i$ are determined by the type of interpolation and the boundary conditions we are imposing.
Example: Cubic b-spline with flat boundary conditions
For a cubic b-spline with flat boundary conditions (i.e. $f'_1(1) = 0$ ${}^\dagger$), this system looks like
$$
\left[\begin{array}{}
-1 & 1 & & \\
\frac{1}{6} & \frac{2}{3} & \frac{1}{6} & \\
& \frac{1}{6} & \frac{2}{3} & \frac{1}{6} \\
& & \ddots & \ddots & \ddots
\end{array}\right]
\left[\begin{array}{}c_0\\c_1\\c_2\\\vdots\end{array}\right]
=
\left[\begin{array}{}0\\v_1\\v_2\\\vdots\end{array}\right]
$$
where the top row of $M_{ij}$, as well as $b_i = \mathbf 0$, are determined by the flat boundary condition.
Question by example: Cubic b-spline for interior domain only
I also want my library to support not having to specify a boundary conditions, but instead use the outer end points to close the system (and thus have the interpolant valid for only $x\in[2,N-1]$). Intuitively it feels like this should be posible, since the boundary condition in effect lets me define a ghost point $c_0$ and solve for it, but whichever way I try I just come up with the same underspecified equation system,
$$
\left[\begin{array}{}
\frac{1}{6} & \frac{2}{3} & \frac{1}{6} & \\
& \frac{1}{6} & \frac{2}{3} & \frac{1}{6} \\
& & \ddots & \ddots & \ddots
\end{array}\right]
\left[\begin{array}{}c_2\\c_3\\\vdots\end{array}\right]
=
\left[\begin{array}{}v_2\\c_3\\\vdots\end{array}\right].
$$
What do I do about $c_1$ and $v_1$ to express the fact that $v_1$ is considered "outside" the domain, and its utility is just to close ?
$\dagger$: In this post I've only cared to mention one end of the domain; for the moment we can consider the symmetric problem where we apply exactly the same conditions in the other end. It was just too verbose to write it out all the time 🙂 Also, the interval starts at $1$ rather than $0$ because Julia is 1-indexed so it makes sense for my implementation; it is really only a matter of notation.
[1]: Thévenaz, Philippe, Thierry Blu, and Michael Unser. "Interpolation revisited [medical images application]." Medical Imaging, IEEE Transactions on 19.7 (2000): 739-758. (Available online here at the time of writing.)
Update 18/2 2016
Ok, I've spent some time with this and tried to figure out how to implement the "not-a-knot" condition suggested by bubba below.
Continuing with the cubic example, the interpolant between $i$ and $i+1$ is given by a linear combination of basis functions and coefficients, that touches the coefficients ranging from $i-1$ to $i+2$;
$$
f_i(x) = \sum_{k=i-1}^{i+2} c_k p_{k-i}(x-i)
$$
If I had $N$ data points and a boundary condition, I would invent an additional data point (as in your suggestion 4.) by applying the boundary condition. Instead, I want to make do without that ghost cell (for implementational reasons) which makes it impossible to evaluate the b-spline on $x\in[1,2)$ (for now), since $f_1$ is dependent on $c_0$ which doesn't then exist. Thus, I take the C3 continuity to mean that $f'''_2(3) = f'''_3(3)$, yielding
$$
-c_1 + 3 c_2 – 3 c_3 + c_4 = – c_2 + 3 c_3 – 3 c_4 + c_5 \\
\\
-c_1 + 4 c_2 – 6 c_3 + 4 c_4 – c_5 = 0
$$
and then I let the $f_2$ piece go all the way to $x=1$, but applying this doesn't yield sensible results. Where am I going wrong?
Best Answer
I don't understand all the details of your question, but there are a few different ways to avoid the need for boundary conditions when constructing splines.
The fundamental problem is that a cubic spline with $n$ segments typically has $n+1$ data points but $n+3$ control points. You solve a linear system of equations to get the control points, but you need two additional pieces of data (the boundary conditions) in order to get $n+3$ equations. If you don't want the user to specify these boundary conditions, you need to balance things some other way. Here are some approaches:
The "not-a-knot" condition. Here you assume that the spline is C3 continuous at the second and last-but-one knots. In other words, the first two segments are actually a single cubic (and similarly the last two). So, the spline has only $n-2$ segments, rather than $n$.
The natural spline condition. You assume that second derivatives are zero at the start and end of the spline.
Estimate first derivatives. Construct a quadratic polynomial interpolating the first three points, and use its first derivative as the first derivative at the start of the spline. Do the same sort of thing at the end-point.
Fabricate two extra data points. As in #3 above, construct a quadratic polynomial from the first three points, and use this to calculate an additional data point to the left of the first one. Do the same sort of thing at the end-point.