Arc length for a spline

Given some spline (e.g. Dierckx.jl, Interpolations.jl), is there a good way to get its arc length?
Right now I’m using this:

using LinearAlgebra, Dierckx
t = range(0, stop = π, length = 100)
x = cos.(t)
y = sin.(t)
spl = ParametricSpline(t, hcat(x, y)')
tl = range(0, stop = π, length = 10000)
sum(norm(derivative(spl, i))*step(tl) for i in tl) # should be π

Just pinging this cause I have a feeling that posting things late CET gets swallowed by newer, more timely, posts…

Arc length of a cubic is an elliptic integral, so there is no nice closed-form expression for it. The best way is to use numerical integration. The function is only smooth in between the knots, so you will want to use numerical integration divided into intervals at each knot.

3 Likes

Wow, ok. I naively thought this was “super easy, barely an inconvenience™”. Ok… Thanks!!!

For information, a direct calculation of the distances between successive interpolated points provides here a much more accurate estimate of the total length of the arc than the summation of the norms of the interpolated derivatives.

using LinearAlgebra, Dierckx
t = LinRange(0, π, 100)
x, y = cos.(t), sin.(t)
spl = ParametricSpline(t, hcat(x, y)')

# Integral of norm of derivatives:
tl = LinRange(0, π, 10_000)
π1 = sum(norm.(eachcol(derivative(spl,tl)))*step(tl))
# Sum of the distances:
xy = spl(tl)
Δxy = xy[:,2:end] - xy[:,1:end-1]
π2 = sum(norm.(eachcol(Δxy)))

# π1 = 3.141906841...
# π2 = 3.141592637...
# π  = 3.1415926536...

I’ve since resorted to using numerical integration:

using QuadGK
π3, _ = quadgk(t -> norm(derivative(spl, t)), 0, π)
(3.1415926535893455, 1.318944953254686e-13)

fast and most accurate :slight_smile:

1 Like

@yakir12, yes that is much more accurate! Thanks :slight_smile:

However, in terms of execution time, it may not be that much better?
The simple method above seems to get the arc length π2 with ~1e-7 accuracy, much faster (x200?) than using QuadGK.
Could you confirm?