Interpolations.jl implements (cardinal) B-splines, i.e. interpolating polynomial functions which are fit to the data and to each other.
We'll start with a simple 1D example - interpolating $f(x) = sin(x)$ over one period, using 10 grid points.
In [1]:
nx = 10
f(x) = sin(2pi/(nx-1) * (x-1))
xcoarse = 1:nx
ycoarse = f.(xcoarse);
xfine = 1:0.1:xcoarse[end]
yfine = f.(xfine);
An interpolation object is basically an array that supports evaluation with any real numbers.
To create an interpolation object, simply call the constructor interpolate, providing the array and some configuration for the interpolation:
In [2]:
using Interpolations
yitp = interpolate(ycoarse, BSpline(Linear()))
Out[2]:
There are a couple of noteworthy points here. First, note the extra argument to interpolate, BSpline(Linear()). This determines the behavior of the interpolating function inside the domain. Linear tells us that it will be a linear interpolation (duh...); in other words, for yitp(x) interpolation will be performed like this:
begin
ix = ifloor(x)
dx = x - ix
(1-dx)*ycoarse[ix] + dx*ycoarse[ix+1]
end
A B-spline interpolation is basically a piecewise polynomial function, where each grid cell is associated with its own polynomial. The coefficients are chosen such that the interpolating function passes through all the data points, while the transitions between polynomials are as smooth as possible given their degree. For example, each grid cell in a linear interpolation (such as yitp) is associated with a straight line, and the interpolation is continuous (with discontinuous derivative) over the entire domain.
Let's take a look at how the different interpolation degrees behave:
In [4]:
yitp_const = interpolate(ycoarse, BSpline(Constant()))
yconst = [yitp_const(x) for x in xfine]
yitp_linear = interpolate(ycoarse, BSpline(Linear()))
ylinear = [yitp_linear(x) for x in xfine]
yitp_quadratic = interpolate(ycoarse, BSpline(Quadratic(Line(OnCell()))))
yquadratic = [yitp_quadratic(x) for x in xfine];
You'll notice that when creating the quadratic interpolation, we had to give another input parameter to the interpolation type: the Line(OnCell()) argument which specifies the boundary condition. The first term (Line) specifies the behavior of the boundary condition, the second (OnCell) specifies where it applies. In this case yitp_quadratic will become linear for x < 0.5 or x > 10.5. The alternative, OnGrid, would apply the boundary condition for x < 1 or x > 10.
All interpolations of quadratic degree or higher require a prefiltering step, which entails solving a tridiagonal system of equations (details can be found for example in this paper), in order to make the interpolating function pass through the data points. Interpolations.jl takes care of solving this system for you, but in order to close the system a boundary condition is requred.
In [6]:
using Gadfly
plot(
layer(x=xcoarse,y=ycoarse,Geom.point),
layer(x=xfine,y=yconst,Geom.line,Theme(default_color=colorant"magenta")),
layer(x=xfine,y=ylinear,Geom.line,Theme(default_color=colorant"red")),
layer(x=xfine,y=yquadratic,Geom.line,Theme(default_color=colorant"green"))
)
Interpolations.jl supports interpolation in arbitrary dimensions - just give it an Array{T,N} where elements t::T support operations like a::Real * t + b::Real * t and it will figure it out. However, since it starts getting difficult to visualize higher-dimensional functions, we'll stick with an example in 2D, namely the function $g(x,y) = x^2 \cdot \sin(y)$.
In [7]:
xs = 1:5
ys = 1:8
g = Float64[x^2 * sin(y) for x in xs, y in ys]
gitp_quad2d = interpolate(g, BSpline(Quadratic(Line(OnCell()))))
display(plot(x=xs,y=ys,z=g,Geom.contour))
display(plot(x=1:.1:5, y=1:.1:8, z=[gitp_quad2d(x,y) for x in 1:.1:5, y in 1:.1:8], Geom.contour))
Simple, huh? =)
The behavior outside the domain is not carefully specified by interpolate, and most likely you should avoid using such objects outside the original domain. If you do want to control behavior outside the domain, you will likely be interested in extrapolate:
yextr = extrapolate(yitp, Flat)
The Flat parameter to extrapolate signifies that we want the value to remain constant (flat) outside of the domain. Other options are...(how the interpolation object will behave when an index outside of the domain $[1, 10]$ is given; with ExtrapError, a BoundsError is thrown, while with ExtrapNaN (the only other implemented extrapolation behavior at the moment), NaN is returned.)