JuliaMath / Interpolations.jl

Fast, continuous interpolation of discrete datasets in Julia
http://juliamath.github.io/Interpolations.jl/
Other
523 stars 110 forks source link

Higher order nonuniform interpolation #518

Open xtalax opened 1 year ago

xtalax commented 1 year ago

Hi,

I'm using Interpolations.jl over in MethodOfLines.jl, and am finding it very convenient for providing an interpolating interface to the solutions of PDES there. The solutions that it produces are most often nonuniform in at least the time dimension, and I would really like to be able to provide higher order interpolations for these solutions.

What would it take to implement higher than linear order nonuniform interpolations, for example a B-Spline interpolation? I'd be happy to implement this feature, if you could point me in the right direction.

mkitti commented 1 year ago

Hi. I we would want to extend the BSpline code to Gridded. The main challenge here is efficiently scaling the grids such that the derivatives match at the knots.

xtalax commented 1 year ago

can you expand on what you mean by "take the derivative match at the knots"?

mkitti commented 1 year ago

The main challenge here is efficiently scaling the grids such that the derivatives match at the knots.

I edited it to "The main challenge here is efficiently scaling the grids such that the derivatives match at the knots."

What we want is the interpolation to be "smooth". Effectively that means that the interpolation and its derivatives should be continuous at the known points, the knots.

This package usually calculate things on a unit grid and then uniformly scales everything. On a non-uniform grid the scaling is discontinuous at the knots. Thus we have rescale the derivatives as we move across the knots.

mkitti commented 1 year ago

Let's start with this example:

julia> using Interpolations, Plots

julia> knots = 1:0.5:5
1.0:0.5:5.0

julia> x = 1:0.1:5
1.0:0.1:5.0

julia> plot(knots, itp.(knots), line=:scatter; label = "Knots")

julia> plot!(x, itp.(x); label = "Interpolated")

image

One way we could introduce non-uniform spacing is just to rescale something to uniform spacing.

julia> rescale(x) = x > 2 ? (x-2)*0.1 + 2 : x
rescale (generic function with 1 method)

julia> plot(x, itp.(rescale.(x)))

image

Here we can see there is a sharp corner. The first derivative is discontinuous.

image

hyrodium commented 1 year ago

This package usually calculate things on a unit grid and then uniformly scales everything. On a non-uniform grid the scaling is discontinuous at the knots. Thus we have rescale the derivatives as we move across the knots.

One way we could introduce non-uniform spacing is just to rescale something to uniform spacing.

I'm not much familiar with the inside implementation of Interpolations.jl, but I believe rescaling is not the right strategy to interpolate on non-uniform knots. My package BasicBSpline.jl has documentation to interpolate values on non-uniform data points. I hope this doc would be helpful.

mkitti commented 1 year ago

Thanks, I'll take a look.