Interpolations.jl icon indicating copy to clipboard operation
Interpolations.jl copied to clipboard

Higher order nonuniform interpolation

Open xtalax opened this issue 2 years ago • 6 comments

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.

xtalax avatar Sep 15 '22 12:09 xtalax

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.

mkitti avatar Sep 15 '22 14:09 mkitti

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

xtalax avatar Sep 15 '22 15:09 xtalax

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 avatar Sep 15 '22 17:09 mkitti

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

mkitti avatar Sep 15 '22 20:09 mkitti

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.

hyrodium avatar Jun 26 '23 14:06 hyrodium

Thanks, I'll take a look.

mkitti avatar Jun 26 '23 22:06 mkitti