How to interpolate a scalar field?

General Usage
,
1

I want to get a polynomial scalar field that interpolate z as a function of x and y.

# consider an integer k,
# then the goal is that a called function f(x[k],y[k])
# would approximate z[k].

# abscissa example with dummy values
x = [1,2,3,4,1,2,3,4,1,2,3,4]
#ordinate example with dummy values
y = [1,1,1,1.1,2,2,2.1,2,3,3,3,3.2] 
# example with dummy values
# to be interpolated as a function of x and y
z = [10,11,12,13,10,11,12,13,10,11,12,13]

How to achieve this ?

For instance, getting the coefficients of the x^n, y^m, x^n y^m would be fine, where n and m are integers.

The goal is then that for any input x-value and y-value (not contained in x and y in the code above), an approximated z-value is obtained with the interpolation.

I have read several post about interpolation in Julia, but they seem to focus purely on linear interpolation of function (and not field).

2

If you want to fit a multivariate polynomial, it’s just linear equations, which are trivial to set up and solve in Julia. For example, if you want to fit to z = a_1 + a_2 x + a_3 y + a_4 x y + a_5 x^2 + a_6 y^2 you could do a least-square fit with:

a = [x.^0 x y x.*y x.^2 y.^2] \ z

(The matrix [ ... ] here is a generalization of the Vandermonde matrix that appears in univariate polynomial fitting: each column corresponds to one basis function.) If the number of basis functions (columns) equals length(z), then this becomes interpolation, i.e. it will exactly go through the given points.

You can also do an equivalent fit more easily and generally with FastChebInterp.jl:

using FastChebInterp, StaticArrays
c = chebregression(SVector.(x, y), z, (2,2))

which represents the polynomial (here of degree 2 in both x and y) in the basis of Chebyshev polynomials, and you can then just evaluate it at any given point with c([x,y]) or c(SVector(x,y)).

However, you have to be careful with polynomial fitting if you go to high degree compared to the number of data points (unless the function is smooth and the points are carefully selected, e.g. using Chebyshev points via chebpoints and chebinterp), since high-degree polynomials can go crazy in between the fitted points if you are not careful.

Other options include radial basis functions (e.g. Surrogates.jl or ScatteredInterpolation.jl, both of which implement several algorithms), splines (e.g. Dierckx.jl can do splines from irregular 2d points), dividing your points into a triangular mesh (e.g. Delaunay triangulation) and doing bi-linear interpolation within each triangle, and more.

See e.g. 2D Interpolation on an irregular grid

3 Likes