# Building a cubic spline interpolator using points and specifying derivatives

**URL:** <https://discourse.julialang.org/t/building-a-cubic-spline-interpolator-using-points-and-specifying-derivatives/64238>\
**Category:** General Usage\
**Tags:** interpolations, splines\
**Created:** [July 7, 2021, 10:45pm UTC](https://discourse.julialang.org/t/building-a-cubic-spline-interpolator-using-points-and-specifying-derivatives/64238 "2021-07-07T22:45:39Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![donkeysaddle](https://avatars.discourse-cdn.com/v4/letter/d/3e96dc/32.png) [@donkeysaddle](https://discourse.julialang.org/u/donkeysaddle)\
**Post date:** [July 7, 2021, 10:45pm UTC](https://discourse.julialang.org/t/building-a-cubic-spline-interpolator-using-points-and-specifying-derivatives/64238/1 "2021-07-07T22:45:39Z")

</div>

Is there a way to build a cubic spline in either Interpolations.jl or Dierckx.jl (or somewhere else) where you can set what the derivatives should be at specified points. Ideally, these derivatives would not have to be at the same locations as the values you are using (in fact, having the derivatives and the values correspond to the same points might not allow for continuity of the second derivative?).

pseudo code:

```julia
#f: a function defined elsewhere
val_times = [1,7,10]
vals = [f(i) for i in val_times]

#f_dev: a function defined elsewhere that gives the derivative of f
dev_times = [3, 5, 10, 11]
devs = [f_dev(i) for i in dev_times]

#build interpolator using the two sets of points 
int = CubicSpline((val_times,vals), (dev_times,devs))

```

---

<div class="post-metadata">

**Author:** ![zdenek\_hurak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/zdenek_hurak/32/53118_2.png) [@zdenek\_hurak](https://discourse.julialang.org/u/zdenek_hurak)\
**Post date:** [July 7, 2021, 11:12pm UTC](https://discourse.julialang.org/t/building-a-cubic-spline-interpolator-using-points-and-specifying-derivatives/64238/2 "2021-07-07T23:12:29Z")

</div>

Although I will not recommend a Julia package, I will at least mention that the functionality you are after could be better search for using the keyword(s) “collocation methods”.

---

<div class="post-metadata">

**Author:** ![donkeysaddle](https://avatars.discourse-cdn.com/v4/letter/d/3e96dc/32.png) [@donkeysaddle](https://discourse.julialang.org/u/donkeysaddle)\
**Post date:** [July 7, 2021, 11:16pm UTC](https://discourse.julialang.org/t/building-a-cubic-spline-interpolator-using-points-and-specifying-derivatives/64238/3 "2021-07-07T23:16:50Z")

</div>

Hopefully your comment will attract the right person…

Thanks!

---

<div class="post-metadata">

**Author:** ![jipolanco](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jipolanco/32/12129_2.png) [@jipolanco](https://discourse.julialang.org/u/jipolanco)\
**Post date:** [July 8, 2021, 8:31am UTC](https://discourse.julialang.org/t/building-a-cubic-spline-interpolator-using-points-and-specifying-derivatives/64238/4 "2021-07-08T08:31:25Z")

</div>

> Although I will not recommend a Julia package, I will at least mention that the functionality you are after could be better search for using the keyword(s) “collocation methods”.

~~I’m a bit lost, I don’t see the issue with mentioning a different package.~~  
Apologies, I misunderstood your comment!

> Is there a way to build a cubic spline in either Interpolations.jl or Dierckx.jl (or somewhere else) where you can set what the derivatives should be at specified points…

You can indeed do this with [BSplineKit.jl](https://github.com/jipolanco/BSplineKit.jl) and a bit of manual work. The basic steps are:

1. Construct a B-spline basis \{b\_j\}\_{j = 1}^N using B-splines of order k = 4 (corresponding to cubic splines). In BSplineKit, the basis can be directly constructed using the locations of all data points. In this basis, a spline can then be written as S(x) = \sum\_{j = 1}^N \alpha\_j b\_j(x), where the \alpha\_j are the B-spline coefficients.

2. Internally, BSplineKit performs [(regular) spline interpolation](https://jipolanco.github.io/BSplineKit.jl/dev/interpolation/) by solving a linear system C \boldsymbol{\alpha} = \boldsymbol{y} to determine the coefficients \alpha\_j from values y\_i at points x\_i. The matrix C, called a [collocation matrix](https://jipolanco.github.io/BSplineKit.jl/dev/collocation/#Matrices) in BSplineKit, contains the evaluations of all B-splines at every data point x\_i, C\_{ij} \equiv b\_j(x\_i). Due to the local support of each B-spline, this matrix is usually banded for a proper choice of the x\_i's.

3. In your case, you want to solve a similar linear system of the form \left( \begin{array}{c} C \\ C' \end{array} \right) \boldsymbol{\alpha} = \left( \begin{array}{c} \boldsymbol{y} \\ \boldsymbol{y}' \end{array} \right), where \boldsymbol{y}' contains the known derivatives at points x'\_i. Similarly to the collocation matrix C, the submatrix C' contains the B-spline derivatives at the points x'\_i, C\_{ij}' = b\_j'(x'\_i).

Below is a small example showing how to do this with BSplineKit.

```julia
using BSplineKit
using SparseArrays
using CairoMakie

f(x) = cos(2x)
f_dev(x) = -2sin(2x)

val_times = [1, 7, 10]
vals = f.(val_times)

dev_times = [3, 5, 10, 11]
devs = f_dev.(dev_times)

# Create B-spline basis, determining knots from interpolation points
# The function `make_knots` is internally used to determine knots from data
# points when doing interpolations.
spline_order = BSplineOrder(4) # corresponds to cubic splines
times = sort(vcat(val_times, dev_times))
ts = SplineInterpolations.make_knots(times, order(spline_order))
B = BSplineBasis(spline_order, ts; augment = Val(false))

# Create collocation matrices for values and derivatives
# Note that we create sparse matrices instead of the default banded matrices
Cval = collocation_matrix(B, val_times, Derivative(0), SparseMatrixCSC{Float64})
Cdev = collocation_matrix(B, dev_times, Derivative(1), SparseMatrixCSC{Float64})

# Determine B-spline coefficients from linear system
A = vcat(Cval, Cdev)
y = vcat(vals, devs)
coefs = A \ y

S = Spline(B, coefs)
S′ = diff(S)
S″ = diff(S, Derivative(2))

fig = Figure()
ax = Axis(fig[1, 1]; xlabel = "Time")
plot!(ax, 1..11, x -> S(x); color = :blue, label = "f(t)")
plot!(ax, 1..11, x -> S′(x); color = :orange, label = "f′(t)")
plot!(ax, 1..11, x -> S″(x); color = :gray, linestyle = :dash, label = "f″(t)")
scatter!(ax, val_times, vals; color = :blue)
scatter!(ax, dev_times, devs; color = :orange)
scatter!(ax, knots(B), zero.(knots(B)); marker = :x, color = :black, markersize = 14)
axislegend(ax; position = :ct, orientation = :horizontal)
save("interp.png", fig)

```

As you can see in the figure, the spline and its derivative properly pass through the input values (markers):

 ![interp](https://global.discourse-cdn.com/julialang/original/3X/1/7/17b1fc2c5f404737ed782e482f9b6ef251fb4f6c.png)

> in fact, having the derivatives and the values correspond to the same points might not allow for continuity of the second derivative?

Not if the B-spline knots (represented by crosses) are not repeated in the interior of the domain. Note that, in the interior, the knot locations never match the data points. As you can see in the figure, the second derivative (dashed line) is continuous even though both the value and the derivative were specified at t = 10.

* * *

EDIT 1: use internal `make_knots` function (used for interpolations) to determine B-spline knots from data points. This is to make sure that we’re solving a square linear system.

EDIT 2: plot knot locations and second derivative.

---

<div class="post-metadata">

**Author:** ![zdenek\_hurak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/zdenek_hurak/32/53118_2.png) [@zdenek\_hurak](https://discourse.julialang.org/u/zdenek_hurak)\
**Post date:** [July 8, 2021, 9:56am UTC](https://discourse.julialang.org/t/building-a-cubic-spline-interpolator-using-points-and-specifying-derivatives/64238/5 "2021-07-08T09:56:17Z")

</div>

> [@jipolanco](#):
>
> > Although I will not recommend a Julia package, I will at least mention that the functionality you are after could be better search for using the keyword(s) “collocation methods”.
> 
> I’m a bit lost, I don’t see the issue with mentioning a different package.

It is just that I was **not able** to give a tip on a suitable package because I am still not familiar with them.

---

<div class="post-metadata">

**Author:** ![rafael.guerra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rafael.guerra/32/216610_2.png) [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)\
**Post date:** [July 8, 2021, 1:34pm UTC](https://discourse.julialang.org/t/building-a-cubic-spline-interpolator-using-points-and-specifying-derivatives/64238/6 "2021-07-08T13:34:36Z")

</div>

@jipolanco thank you for this.  
It is a real gem.

---

<div class="post-metadata">

**Author:** ![donkeysaddle](https://avatars.discourse-cdn.com/v4/letter/d/3e96dc/32.png) [@donkeysaddle](https://discourse.julialang.org/u/donkeysaddle)\
**Post date:** [July 8, 2021, 7:49pm UTC](https://discourse.julialang.org/t/building-a-cubic-spline-interpolator-using-points-and-specifying-derivatives/64238/7 "2021-07-08T19:49:01Z")

</div>

@jipolanco This is so awesome! What a great explanation. Will likely be following up with questions…

Thanks so much, really appreciate it.

- DS