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 and a bit of manual work. The basic steps are:
-
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.
-
Internally, BSplineKit performs (regular) spline 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 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.
-
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.
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):
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.
