Start Julia and run the following commands:
julia> using Pkg
julia> Pkg.add("NormalHermiteSplines")
To use this package, begin your code with
using NormalHermiteSplines
The first example is the function φ(x,y)=sin(4.0*sqrt(x^2+y^2)) defined in the region Ω=[-1,1]x[-1,1].
We'll construct an interpolating normal spline using this function and its gradient values sampled on set of 1000 Halton nodes ([1]) distributed in the [-1,1]x[-1,1] square.
using NormalHermiteSplines
nodes = get_2D_halton_nodes(1000) # generates Halton data set in [0, 1] x [0, 1]
n_1 = size(nodes, 2)
u = Vector{Float64}(undef, n_1) # function values
d_nodes = Matrix{Float64}(undef, 2, n_1) # directional derivative nodes
es = Matrix{Float64}(undef, 2, n_1) # derivative directions
du = Vector{Float64}(undef, n_1) # directional derivative values
grad = [0.0; 0.0]
for i = 1:n_1
nodes[1, i] = nodes[1, i] * 2.0 - 1.0 # transforming Halton nodes to [-1, 1] x [-1, 1]
nodes[2, i] = nodes[2, i] * 2.0 - 1.0
d_nodes[1,i] = nodes[1,i]
d_nodes[2,i] = nodes[2,i]
d = sqrt(nodes[1, i]^2 + nodes[2, i]^2) # here d > 0
u[i] = sin(4.0 * d)
grad[1] = 4.0 * nodes[1, i] * cos(4.0 * d) / d
grad[2] = 4.0 * nodes[2, i] * cos(4.0 * d) / d
es[1,i] = grad[1] # no need to normalize 'es' vectors
es[2,i] = grad[2]
du[i] = norm(grad)
end
# Hermite spline must be constructed with RK_H1 or RK_H2 kernel.
# Here value of the 'scaling parameter' ε is estimated in the interpolate procedure.
rk = RK_H1()
#
spline = interpolate(nodes, u, d_nodes, es, du, rk)
grid = get_2D_grid2(100) # regular grid covering [-1, 1] x [-1, 1]
σ = evaluate(spline, grid) # spline values on the grid
σ1 = evaluate_one(spline, [0.5; 0.5])
# ≈ 0.308071
#
# exact function value at [0.5; 0.5]:
# ≈ 0.308072
g1 = evaluate_gradient(spline, [0.5; 0.5])
# ≈ -2.69080
# ≈ -2.69090
#
# exact function gradient at [0.5; 0.5]:
# ≈ -2.69086
# ≈ -2.69086
The spline surface and filled 2-D contour plots:
Approximation error plots:
Spline was evaluated on a uniform Cartesian grid of size 101x101. Accuracy of the interpolation was measured by calculating the Root Mean Square Error (RMSE) and the Maximum Absolute Error (MAE). For this case
RMSE: 1.6E-03, MAE: 1.1E-01, estimated value of the scaling parameter ε is 8.8, estimation of the Gram matrix condition number is 1.0E+11.
The second example is the function Ψ(x,y,z)=cos(π*x)*cos(y-0.5)*sin(π*(z-0.5)) defined in the region Ω=[0,1]x[0,1]x[0,1].
We'll construct an interpolating normal spline using function Ψ values sampled on set of 1000 non-uniform random nodes distributed in the unit cube [0,1]x[0,1]x[0,1].
using NormalHermiteSplines
nodes = get_3D_random_grid(9) # generates 1000 non-uniform random nodes
n_1 = size(nodes, 2)
u = Vector{Float64}(undef, n_1) # function values
for i = 1:n_1
x = nodes[1,i]
y = nodes[2,i]
z = nodes[3,i]
u[i] = cos(π*x)*cos(y-0.5)*sin(π*(z-0.5))
end
# Here spline is being constructed with RK_H2 kernel,
# the 'scaling parameter' ε is defined explicitly.
rk = RK_H2(5.0)
#
spline = interpolate(nodes, u, rk)
grid = get_3D_grid(50) # creates the uniform Cartesian grid of size 51x51x51 in [0, 1]x[0, 1]x[0, 1]
σ = evaluate(spline, grid)
σ1 = evaluate_one(spline, [0.8; 0.6; 0.8])
# ≈ -0.65122
#
# exact function value at [0.8; 0.6; 0.8]:
# ≈ -0.65124
g1 = evaluate_gradient(spline, [0.8; 0.6; 0.8])
# ≈ -1.4862
# ≈ 0.0653
# ≈ -1.4863
#
# exact function gradient at [0.8; 0.6; 0.8]:
# ≈ -1.4865
# ≈ 0.0653
# ≈ -1.4865
The spline plots:
Approximation error plots:
Spline was evaluated on a uniform Cartesian grid of size 51x51x51. For this case RMSE: 2.9E-03, MAE: 5.1E-02, value of the scaling parameter ε is 5.0, estimation of the Gram matrix condition number is 1.0E+14.
Further examples are given in documentation.
The normal splines method for one-dimensional function interpolation and linear ordinary differential and integral equations was proposed in [2]. An idea of the multivariate splines in Sobolev space was initially formulated in [8], however it was not well-suited to solving real-world problems. Using that idea the multivariate generalization of the normal splines method was developed for two-dimensional problem of low-range computerized tomography in [3] and applied for solving a mathematical economics problem in [4]. At the same time an interpolation scheme with Matérn kernels was developed in [9], this scheme coincides with interpolating normal splines method. Further results related to applications of the normal splines method were reported at the seminars and conferences [5,6,7].
For more information and explanation see Documentation.
References
[1] Halton sequence
[2] V. Gorbunov, The method of normal spline collocation. USSR Computational Mathematics and Mathematical Physics, Vol. 29, No. 1, 1989
[3] I. Kohanovsky, Normal Splines in Computing Tomography (Нормальные сплайны в вычислительной томографии). Avtometriya, No.2, 1995
[4] V. Gorbunov, I. Kohanovsky, K. Makedonsky, Normal splines in reconstruction of multi-dimensional dependencies. Papers of WSEAS International Conference on Applied Mathematics, Numerical Analysis Symposium, Corfu, 2004
[5] I. Kohanovsky, Multidimensional Normal Splines and Problem of Physical Field Approximation, International Conference on Fourier Analysis and its Applications, Kuwait, 1998.
[6] I. Kohanovsky, Inequality-Constrained Multivariate Normal Splines with Some Applications in Finance. 27th GAMM-Seminar on Approximation of Multiparametric functions, Max-Planck-Institute for Mathematics in the Sciences, Leipzig, Germany, 2011.
[7] V. Gorbunov, I. Kohanovsky, Heterogeneous Parallel Method for the Construction of Multi-dimensional Smoothing Splines. ESCO 2014 4th European Seminar on Computing, 2014
[8] A. Imamov, M. Dzhurabaev, Splines in S.L. Sobolev spaces (Сплайны в пространствах С.Л.Соболева). Deposited manuscript. Dep. UzNIINTI, No 880, 1989.
[9] J. Dix, R. Ogden, An Interpolation Scheme with Radial Basis in Sobolev Spaces H^s(R^n). Rocky Mountain J. Math. Vol. 24, No.4, 1994.