Dear all,
I am trying to interpolate a function f: R^2 → R with Interpolations.jl.
But as far as I understand the documentation this is only possible for functions from R^n to R^n.
Are you aware of a way to build an interpolation?
Thank you so much!
Hi! There are multiple packages than can do multivariate interpolations. You can have a look at the following thread: 3D interpolation in Julia. To name some you might want to look at: Dierckx.jl, GridInterpolations.jl, GMT.jl, DIVAnd.jl, RadialBasisFunctions.jl, ScatteredInterpolation.jl, KernelInterpolation.jl. Which of these is suited the best for your use case depends on how your data looks like, especially the location of your data points in \mathbb{R}^2 (e.g., grid data or scattered data). If you have any questions about one (or multiple) of these packages, feel free to ask.
2 Likes
If your data is scattered, NaturalNeighbours.jl may be of interest.
Yes, but I think the person asking the question wants to interpolate 2D data.
Sorry, I posted a similar question and I thought you were replying to that…
1 Like
Please check this simple example:
# R² -> R
using Interpolations
f(x, y) = 2x - 5y
x = y = -10:5.0:10
v = f.(x, y')
itp = interpolate((x,y), v, Gridded(Linear()))
itp(1.0, 2.0) ≈ f(1.0, 2.0) # true (= -8.0)
1 Like
Dear @JoshuaLampert, @rafael.guerra thank you for your answers! Maybe this minimal example helps:
# x coord; dim = n x 1
X = [0.0 1.0 3.0 0.0 1.0 3.0 0.0 1.0 3.0]
# y coord; dim = n x 1
Y = [0.0 0.0 0.0 1.0 1.0 1.0 3.0 3.0 3.0]
# PDE sol; dim = n x 1
sol = [0.5 1.0 1.5 1.0 0.5 1.0 1.5. 2.0 0.5]
(X[i], Y[i]) is a node in a grid and sol[i] is the solution at (X[i], Y[i]).
I want to build now an interpolation of the vector sol on the whole domain.
AFAIU for the answers and the other threads, I need a a n x n matrix for sol?
If your data is from a grid, you will get the most efficiency from providing sol as a matrix that conforms to the grid. Do you have data on the entire grid, or do you need a scattered interpolant?
If you do actually want to interpolate with data like from your example, you can do e.g.
julia> using NaturalNeighbours
julia> X = [0.0 1.0 3.0 0.0 1.0 3.0 0.0 1.0 3.0] |> vec;
julia> Y = [0.0 0.0 0.0 1.0 1.0 1.0 3.0 3.0 3.0] |> vec;
julia> sol = [0.5 1.0 1.5 1.0 0.5 1.0 1.5 2.0 0.5] |> vec;
julia> itp = interpolate(X, Y, sol); # also see ?interpolate for some keywords
julia> itp(0.2, 0.3)
0.69
julia> itp([0.5, 0.2, 0.3], [0.1, 0.2, 0.9]; method = Laplace()) # see ?interpolate
3-element Vector{Float64}:
0.75
0.66
0.8300000000000001
For your data example, you could reshape the input as follows:
using Interpolations
X = [0.0 1.0 3.0 0.0 1.0 3.0 0.0 1.0 3.0]
Y = [0.0 0.0 0.0 1.0 1.0 1.0 3.0 3.0 3.0]
sol = [0.5 1.0 1.5 1.0 0.5 1.0 1.5 2.0 0.5]
x, y = unique(X), unique(Y)
v = reshape(sol, length(x),length(y))
itp = interpolate((x,y), v, Gridded(Linear()))
itp(0.5, 2.0) # 1.25
NB:
The x axis is along the vertical columns of v and the y-axis along the rows.
2 Likes
If you are trying to interpolate a smooth function f(x,y), e.g. the solution of a PDE as a function of some parameters, and you have the freedom to choose the x,y points, then you will probably be much better off choosing the x,y points to lie on a Chebyshev grid and then use Chebyshev interpolation, which converges exponentially fast for smooth functions. See FastChebInterp.jl.
On the other hand if your x,y points come from how you solve the PDE, e.g. a finite-difference grid, then you should choose your interpolation based on what computational method you used. e.g. if you used 2nd-order centered differences then you should just use bilinear interpolation, and if you used finite elements you should use the FEM basis functions.
4 Likes