Best way to take derivatives of unevenly spaced data (with interpolations? discrete derivatives?)

CosmoProf · January 28, 2021, 4:57am

Hi, I have a function f(t), where f is numerically sampled at discrete points in time t that are NOT evenly sampled in time. (This is unavoidable.)

I need to take the first and second derivatives of f(t), compute a complicated function (call it D(t)) of f' and f'' over that time interval, and integrate the function D(t) over the interval.

Does anyone know how best to do this? There are ways to estimate derivatives of discrete data, but it gets messy for unevenly spaced intervals. There are probably ways to calculate and manipulate derivatives of interpolations, but I’m not sure how to do that.

Thanks for any suggestions…!

rafael.guerra · January 28, 2021, 6:00am

Please check Interpolations.jl, in the General usage section.

CosmoProf · January 28, 2021, 6:26am

I was hoping for more detailed suggestions than this.

macd · January 28, 2021, 6:42am

I think your best bet would be to use ApproxFun. Once you have your data as a “Fun”, it is straightforward to take derivatives and integrate. Unfortunately, you need your data to be on a Chebshev grid for ApproxFun. If you are on an even grid, or and uneven one, a rational approximation can help. I wrote a small package (https://github.com/macd/BaryRational.jl) which was recently registered. The readme has an example on how to use the AAA algorithm to get such data into ApproxFun as a Fun.

dlfivefifty · January 28, 2021, 7:45am

No you can use ApproxFun with an uneven grid:

https://juliaapproximation.github.io/ApproxFun.jl/latest/faq/

You may want to make the matrix rectangular to get a least squares fit to avoid over fitting

jipolanco · January 28, 2021, 8:48am

In my opinion splines are great for this. The idea is to construct a spline interpolating your unevenly spaced data, and then to evaluate the spline derivatives at any point in space.

This can be done with my own BSplineKit.jl package, but there are other very good alternatives that can do the same.

This is an example using BSplineKit (I’ll add it soon to the docs):

using BSplineKit
using LaTeXStrings
using Plots
using Random

# Actual function and its derivative
f(x) = sin(2π * x)
f′(x) = 2π * cos(2π * x)

# Generate uneven grid of ``N`` points in ``[0, 1]``
N = 20
rng = MersenneTwister(42)
x = sort!(rand(rng, N))
x[1] = 0; x[end] = 1;  # make sure endpoints are included
y = f.(x)              # actual function at grid points

# Interpolate data with B-splines of order 6 (polynomial degree 5)
yint = interpolate(x, y, BSplineOrder(6))

# Extract resulting spline and generate its first derivative
S = spline(yint)             # this is the spline interpolation of f
S′ = diff(S, Derivative(1))  # this is a spline approximating f′

# Evaluate splines and plot
xplot = 0:0.01:1
plt = plot(f, 0, 1; label = L"f(x)", xlabel = L"x", title = "Function interpolation")
plot!(plt, xplot, S.(xplot); label = "Interpolation")
scatter!(plt, x, y; label = "Input data")

plt = plot(f′, 0, 1; label = L"f'(x)", xlabel = L"x", legend = :top, title = "Derivative approximation")
plot!(plt, xplot, S′.(xplot); label = "Interpolation")
scatter!(plt, x, S′.(x); label = "Interpolation at grid points")

This gives the following approximations for f and its derivative:

interpolation
derivative

rafael.guerra · January 28, 2021, 10:25am

@CosmoProf, fair enough. Maybe Dierckx.jl would be worth your attention. The points can be unevenly spaced.
PS: code edited to use ApproxFun for numerical integration and used consistent spline orders as per @jipolanco

using Dierckx, Random, Plots
f(t) = sin(t)/t + (t^2)/10
f′(t) = cos(t)/t - sin(t)/t^2 + 2*t/10
f′′(t) = -2*cos(t)/(t^2) + 2*sin(t)/(t^3) - sin(t)/t + 2/10  # zero @ t0 = 1.20608
RCurv(t)::Float64 = (1.0 + f′(t)^2)^(3/2) / abs(f′′(t))  # division by 0 @ t0 = 1.20608

# spline degree=k, order=k (Dierckx), order=k+1 (BSplineKit)
k = 4;  str=" (degree=" * string(k) * ")" 
t0 = 1.21   # singularity at t0 = 1.20608 
tr = range(t0,20,length=100)
tir = [tr[1]; sort(tr[shuffle(2:length(tr)-1)[1:55]]); tr[end]]
spl = Spline1D(tir, f.(tir), k=k)  # irregular input
dfir = derivative(spl, tir)  # 1st derivative in irregular grid
spl2 = Spline1D(tir, dfir, k=k) 
d2fir = derivative(spl2, tir)  # 2nd derivative

spl = Spline1D(tir, RCurv.(tir), k=k)  # irregular input
# iRCir = [Dierckx.integrate(spl, t0, t) for t in tir]  # builtin integration not accurate near t0
IntegrateStep = (tir[end]-tir[1])/(10000-1)
IntegrateRange = tir[1]:IntegrateStep:tir[end]
RCurvD = cumsum(spl.(IntegrateRange))*IntegrateStep

using ApproxFun
RCr = cumsum(Fun(RCurv, t0..tr[end])) # numerical integration using ApproxFun

using BSplineKit
S = BSplineKit.spline(BSplineKit.interpolate(tir, f.(tir), BSplineOrder(k+1)))
Sd1 = BSplineKit.diff(S, BSplineKit.Derivative(1))
Sd2 = diff(Sd1, BSplineKit.Derivative(1))

RCurvS(t) = (1.0 + Sd1(t)^2)^(3/2)/abs(Sd2(t))  # function of BSplines
RCurvDel = cumsum(RCurvS.(IntegrateRange))*IntegrateStep

h = Array{Plots.Plot{Plots.GRBackend}}(undef,4)
plot(tr,f, label="f(t)", legend=:topleft)
h[1] = scatter!(tir, f.(tir), mc=:green, label="Uneven samples")

plot(tr,f′, label="f′(t)", legend=:topleft)
scatter!(tir, Sd1.(tir), mc=:red, label="BSplineKit"*str)
h[2] = scatter!(tir, dfir, mc=:green, label="Dierckx"*str)

plot(tr,f′′, label="f′′(t)", legend=:topright)
scatter!(tir, Sd2.(tir), mc=:red, label="BSplineKit"*str)
h[3] = scatter!(tir, d2fir, mc=:green,label="Dierckx"*str)

plot(tr,RCr.(tr),lc=:blue, label="ApproxFun integral{ RCurv(t)=(1 + f′(t)^2)^(3/2)/|f′′(t)| }", legend=:topleft)
plot!(IntegrateRange, RCurvDel, lc=:red, ls=:dash, label="BSplineKit"*str*" + user integration")
h[4] = plot!(IntegrateRange, RCurvD, lc=:green, ls=:dash, label="Dierckx"*str*" + user integration")

plot(h..., ms=4, ma=0.5, msw =0, legendfontsize=10, layout=(2,2), size=(1600,1600))
savefig("Dierckx_BSplineKit_ApproxFun_interpolations.png")

rafael.guerra · January 28, 2021, 12:24pm

@jipolanco, your package seems to have most of Dierckx.jl’s features and to go far beyond. Would you be able to comment if it can completely replace Dierckx? Thank you.

jipolanco · January 28, 2021, 12:40pm

Some of the features of Dierckx.jl, such as parametric splines and roots of splines, are not yet available in BSplineKit.jl, even though it shouldn’t be too difficult to add this.

Conversely, there are features, such as splines of arbitrary order, that are not available in Dierckx. This also includes a large number of tools for applying B-splines for the solution of boundary value problems, which was my main motivation for this package.

macd · January 28, 2021, 4:23pm

That is precisely my point. You use the AAA algorithm to generate a model from an uneven grid, evaluate that model on the Chebyshev points. Then use ApproxFun.transform to transform those values into the Chebyshev coefficients and then construct a Fun. I believe this is better way to get this type of data into ApproxFun than using a Vandermonde approach.
But ApproxFun is a heavy dependency and splines are good at this too. I’ve used Dierckx.jl and it’s quite nice. Didn’t know about BSplineKit.jl, will have to take a look!

dlfivefifty · January 28, 2021, 4:34pm

But once you have a rational approximation you can just work with that directly…

we once had a basic version called RatFun.jl implemented but never did AAA and it’s gone a bit stale

macd · January 28, 2021, 7:11pm

You’re right, of course, that you could work with the rational approximation directly, if all of that functionality was implemented there, but it is not. So getting it into ApproxFun is the next best (well maybe not best, but certainly easy) thing if you want integrals of complicated functions of the first and second derivatives of the numerically defined function.

At the risk of beating a dead horse and borrowing @jipolanco 's example and then extending to defining a ‘complicated’ function of the function and its first and second derivatives, and then integrating that, we see

using ApproxFun
using BaryRational
using Plots
using Random

# Actual function and its derivatives
f(x) = sin(2π * x)
f′(x) = 2π * cos(2π * x)
f′′(x) = -4π^2 * sin(2π * x)

# Generate uneven grid of ``N`` points in ``[0, 1]``
N = 20
rng = MersenneTwister(42)
x = sort!(rand(rng, N))
x[1] = 0; x[end] = 1;  # make sure endpoints are included
y = f.(x)              # actual function at grid points

modl = aaa(x, y)

M = 129
pts = chebyshevpoints(M)
S = Chebyshev()

# The Fun and its first two derivatives
nfun = Fun(S, ApproxFun.transform(S, modl.(pts)))
D′  = nfun'
D′′ = nfun''

xplot = 0:0.01:1
plt = plot(f, 0, 1; label = "f(x) exact", xlabel = "x", title = "Function interpolation")
plot!(plt, xplot, nfun.(xplot); label = "f(x) Interpolation")
scatter!(plt, x, y; label = "Input data")
savefig("function")

plt = plot(f′, 0, 1; label = "f'(x) exact",  xlabel = "x", legend = :top, title = "1st derivative approximation")
plot!(plt, D′, 0, 1; label = "f'(x) approx")
scatter!(plt, x, D′.(x); label = "Approximation at grid points")
savefig("1_derivative")

plt = plot(f′′, 0, 1; label = "f''(x) exact", xlabel = "x", legend = :top, title = "2nd derivative approximation")
plot!(plt, xplot, D′′.(xplot); label = "f''(x) approx")
scatter!(plt, x, D′′.(x); label = "Approximation at grid points")
savefig("2_derivative")

G = (nfun + 4*D′)^2 + sin(D′′)

Gint = cumsum(G)

plt = plot(xplot, Gint.(xplot); label = "Integral approximation", title="Integral Approximation of (f(x) + 4f'(x))^2 + sin(f''(x))")
scatter!(plt, x, Gint.(x); label = "Approximation at grid points")
savefig("integral")

function
1_derivative
2_derivative
integral

dlfivefifty · January 28, 2021, 7:15pm

Awesome example! Mind adding it to ApproxFunExamples?

macd · January 28, 2021, 7:19pm

Thanks and sure thing.

CosmoProf · January 29, 2021, 2:09am

Hi folks, thanks for all the suggestions; I have more comments coming soon, but for now I cannot add BaryRational, I get the error:

julia> Pkg.add("BaryRational")
ERROR: The following package names could not be resolved:
 * BaryRational (not found in project, manifest or registry)

CosmoProf · January 29, 2021, 2:49am

Hi again everyone, thanks again for the all the suggestions so far.
Currently I am basing a trial calc on the example from jipolanco.
I am substituting very fine sampling of the splines for an “integral”.

My concerns are that the first and last couple of points of the spline interpolation seem to be bad fits; and more importantly, each time I run the program (with different randomized time points), the final integrated number is different, even though the underlying function is the same.

I’ve tried pumping up the spline interpolation order, and the number of sampled time points (which I couldn’t do with my real data). That makes the variation less, but still very large.
If anyone has any ideas on how to get a more consistent result, please let me know!

Here is my code:

using Plots
# using Interpolations
using BSplineKit


# Actual function and its derivatives:
f(x) = sin(x)/x + ((x^2)/10.0)
fP(x) = (cos(x)/x) - (sin(x)/x^2) + (2.0*x/10.0)
fPP(x) = -(2*cos(x)/(x^2)) + (2*sin(x)/(x^3)) - (sin(x)/x) + (2.0/10.0)

t = sort([((19.0*rand(Float64))+1.0) for i in 1:57])

S = spline(BSplineKit.interpolate(t, f.(t), BSplineOrder(12)))
Sd1 = diff(S, Derivative(1))
Sd2 = diff(Sd1, Derivative(1))

RCurv(x) = (((1.0+(Sd1(x)^2))^(3/2))/abs(Sd2(x)))
IntegrateStep = (t[end-2]-t[3])/(10000-1)
IntegrateRange = t[3]:IntegrateStep:t[end-2]
RCurvDel = RCurv.(IntegrateRange)
RCurvIntegratedAverage = sum(RCurvDel)/length(RCurvDel)
println("\n","RCurvIntegratedAverage = ",RCurvIntegratedAverage,"\n")

And here are some plotting commands to visualize it:

plot(IntegrateRange,RCurvDel)

PlotRng = 1.0:0.001:20
RealVsSplineFunctionPlot = plot(PlotRng, f.(PlotRng), color="aqua", linewidth=7.0);
            scatter!(t, f.(t), markersize=9, color = "blue", legend = false);
            plot!(PlotRng, S.(PlotRng), color="red", linewidth=2.0);
            scatter!(t, S.(t), color="red", markersize=5); display(RealVsSplineFunctionPlot)
#
RealVsSplineFIRSTDerivPlot = plot(PlotRng, fP.(PlotRng), color="aqua", linewidth=7.0);
            scatter!(t, fP.(t), markersize=9, color = "blue", legend = false);
            plot!(PlotRng, Sd1.(PlotRng), color="red", linewidth=2.0);
            scatter!(t, Sd1.(t), color="red", markersize=5); display(RealVsSplineFIRSTDerivPlot)
#
RealVsSplineSECONDDerivPlot = plot(PlotRng, fPP.(PlotRng), color="aqua", linewidth=7.0);
            scatter!(t, fPP.(t), markersize=9, color = "blue", legend = false);
            plot!(PlotRng, Sd2.(PlotRng), color="red", linewidth=2.0);
            scatter!(t, Sd2.(t), color="red", markersize=5); display(RealVsSplineSECONDDerivPlot)

Thanks!

macd · January 29, 2021, 7:27am

I did double check and it is in the general registry. Perhaps updating the registry would help? From the pkg repl (use ] to enter) type registry update and see if that makes a difference.

jipolanco · January 29, 2021, 10:14am

I ran your script, and indeed you obtain bad fits near the boundaries. This is what you get for the second derivative:

RealVsSplineSECONDDerivPlot

(For next time, it would be easier to help you if you included your results in the post.)

What you observe is expected, since near the boundaries you’re extrapolating, that is, you’re evaluating the splines outside of their domain of definition. This is because the first and last points of your t vector, which set the domain of the spline, are random, instead of being exactly located at 1 and 20. That’s precisely why, in my example, I took care of including the endpoints in the generated grid.

rafael.guerra · January 29, 2021, 5:09pm

It would be useful to compare the different methods used here. Is there a simple standard numeric test that is used to compare interpolation code accuracy?

dlfivefifty · January 29, 2021, 6:43pm

If you know the true answer just take the norm of difference on a dense grid

If you can increase the density of samples then you can investigate convergence rate (here AAA + ApproxFun should have spectral convergence so will beat algebraic convergence of splines )

If all you have is the data of a single sample it’s likely an ill-posed question : there’s no notion of “correct”