In this notebook, we show how to use Julia to interpolate a function using cubic splines. The code and this notebook are available online on github
In [1]:
# we need to import the library first
using splines
In [2]:
function f(x::Float64,y::Float64)
if x==0 && y==0
return 1.0
else
t = sqrt( x^2+y^2 )
return sin(pi*t)/(pi*t)
end
end
f(x::Array{Float64,1}, y::Array{Float64,1}) = [f(x[i],y[i]) for i=1:length(x)]
Out[2]:
The grid over which we approximate the function is defined by its boundaries $[a_i, b_i]$ and the number of points $n_i$ in each dimension $i \in [1,d]$.
In each dimension $i\in[1,d]$ the $n_i$ points $(a_i=s^1_i, s^2_i, ..., s^{n_j}_i=b_i)$ are regularly spaced, i.e. $s^j_i = a_i + j\frac{b_i-a_i}{n_i}$. These points are called nodes.
For the two dimensional function defined above, we define $[a_1,b_1]=[a_2,b_2]=[-5,5]$ and we use 50 points in each dimension.
In [3]:
# define the approximation space
a = [-5.0, -5.0] # a_1, a_2
b = [5.0, 5.0] # b_1, b_2
orders = [15,15] # n_2, n_2
Out[3]:
The cartesian product of these nodes is used to approximate the compact state $\mathcal{S}=[a_1, b_1], ... [a_d, b_d]$, by a finite set of points $S=(s_{i_1,...i_d})_{i_1\in[1,n_1], ..., i_d\in[1,n_d]}$. By definition, the point indexed by i_1,...id is simply the tuple $(s{i1},...s{i_d})$.
In general, this multivariate index is sufficient and it is not needed to list these points one after all. In case it is needed, the convention in Julia is that the first index varies faster meaning that points are listed like in $s_{11}, s_{21}, s_{31}, s_{12}, s_{22}, s_{32}$
In [4]:
# optional plot the grid
using PyPlot
for x=linspace(a[1],b[1],orders[1]), y=linspace(a[2],b[2],orders[2])
PyPlot.plot(x,y,"o", color="red")
end
# vectorized version
#all_points = [ [x y] for x=linspace(a[1],b[1],orders[1]), y=linspace(a[2],b[2],orders[2])]
#all_points = vcat(all_points...)
#PyPlot.plot(all_points[:,1], all_points[:,2], "x", color="blue")
The interpolating data is supposed to fit the exact values $v_{i_1, ..., i_d}=f(s_{i_1, ..., i_d})$ of the function on the grid. We compute them below.
In [5]:
vals = Float64[f(x,y) for x=linspace(a[1],b[1],orders[1]),
y=linspace(a[2],b[2],orders[2])]
display(size(vals)) # the result is a 50x50 array
Before we can interpolate the function, we need to precompute some parameters of the interpolant (here the b-spline coefficients, more of it below). This can be done once for several interpolating calls.
The function interpolant_cspline does the required work provided with information describing the interpolation space and the values on the grid.
In [8]:
interp = interpolant_cspline(a, b, orders, vals)
Out[8]:
This function can readily be readily evaluated on multiple arguments.
In [9]:
interp(0.0, 0.1)
Out[9]:
When evaluating the interpolant on N multiple points, it is more efficient to prepare an array of $N$ point with size $N \times 2$.
In [10]:
points = [[0.0 z] for z = linspace(-6.0, 6.0, 100)]
points = vcat(points...)
interpolated_values = interp(points)
display("Array of points")
display((points))
display("Array of interpolated values")
display(interpolated_values)
In [11]:
plot(points[:,2], interpolated_values,"-",label="interpolated values")
plot(points[:,2], f(points[:,1], points[:,2]),":", color="black",label="true values")
nodes = linspace(a[2], b[2], orders[2])
plot(nodes, f(0*nodes,nodes),"o", color="red", label="nodes" )
legend()
Out[11]:
TODO: explain how basis splines are defined on knots and how to directly deal with the spline coefficients.
In [12]:
# fit coefficients of the spline (note the size of the matrix ( 15+2 x 15+2 )
coeffs = filter_coeffs(a,b,orders,vals)
display(coeffs)
In [12]:
# cartesian product with 50 points along each dimension
fine_grid = vcat( [ [x y] for x = linspace(0,3,10), y = linspace(0,3,10) ]... );
In [12]:
# evaluate the spline
interpolated_values_2 = eval_UC_spline( a, b, orders, coeffs, fine_grid );
# evaluate with gradient
V, dV = eval_UC_spline_G( a, b, orders, coeffs, fine_grid );
In [13]:
xm = reshape(fine_grid[:,1], 10,10)
ym = reshape(fine_grid[:,2], 10,10)
V_mesh = reshape(V, 10,10)
contourf(xm,ym,V_mesh,100)
colorbar()
# add gradient field
# see api of quiver at http://matplotlib.org/api/pyplot_api.html
dV_m_1 = reshape(dV[:,1],10,10)
dV_m_2 = reshape(dV[:,1],10,10)
#streamplot( linspace(a[1], b[1], 50), linspace(a[1], b[2], 50), dV_m_1, dV_m_2)
quiver( linspace(0,3, 10), linspace(0,3, 10), dV_m_1, dV_m_2)
Out[13]: