The interpolant is a sum of polynomial pieces, so its derivatives are polynomial pieces too, and both packages evaluate them analytically: the derivative keyword gives an interpolant of a derivative instead of the data, built from the differentiated kernel. There is no finite difference and no step size to tune, and the result is exact to rounding for the interpolant that is being differentiated.
x = range(0, 2π, length = 40) # 40 uniformly spaced samples
y = sin.(x) # the values at the samples
d1 = convolution_interpolation(x, y; derivative = 1) # the first derivative of the interpolant
d2 = convolution_interpolation(x, y; derivative = 2) # the second derivative
d1(1.0), d2(1.0) # cos(1) and −sin(1), approximately(0.5403022965038063, -0.8414703517312542)import numpy as np
from convinterp import convolution_interpolation
x = np.linspace(0.0, 2 * np.pi, 40) # 40 uniformly spaced samples
y = np.sin(x) # the values at the samples
d1 = convolution_interpolation(x, y, derivative=1) # the first derivative of the interpolant
d2 = convolution_interpolation(x, y, derivative=2) # the second derivative
d1(1.0), d2(1.0) # cos(1) and −sin(1), approximately| f′(1) | f″(1) | |
|---|---|---|
| interpolant | 0.540302296504 | -0.841470351731 |
| exact | 0.540302305868 | -0.841470984808 |
Each kernel supports derivatives up to a fixed order (the table on the interpolation page): one for :a3 to :a7, three for :b5, five for :b7, six for :b9 and seven for :b11 and :b13. Asking for more is an error. The limit is set by accuracy, not by the polynomials: a kernel of degree p has non-trivial derivatives up to order p, but each order multiplies the rounding error by 1/h, and beyond the listed order that amplification exceeds what the derivative is worth in double precision.
Differentiating the interpolant costs about one order of convergence per derivative order: a 7th-order kernel gives a 6th-order first derivative and a 5th-order second derivative. Runge’s function again, with :b7, and the largest error of the interpolant and its first three derivatives against the spacing:
runge(x) = 1 / (1 + 25x^2) # Runge's function
drunge = (x -> -50x / (1 + 25x^2)^2, # its first derivative
x -> (3750x^2 - 50) / (1 + 25x^2)^3, # its second derivative
x -> (15000x - 375000x^3) / (1 + 25x^2)^4) # its third derivative
xr = range(-1, 1, length = 2001) # fine points for measuring the error
ns = [41, 81, 161, 321, 641, 1281] # the numbers of samples
hs = 2 ./ (ns .- 1) # the corresponding spacings
orders = 0:3 # the interpolant and three derivatives
exact = (runge, drunge...) # the function and its derivatives
errors = [begin
xn = range(-1, 1, length = n)
itp_m = convolution_interpolation(xn, runge.(xn); kernel = :b7, derivative = m)
maximum(abs.(itp_m.(xr) .- exact[m+1].(xr)))
end for n in ns, m in orders]
plot(xscale = :log10, yscale = :log10, xlabel = "spacing h", ylabel = "largest error",
legend = :bottomleft, xflip = true)
for m in orders
plot!(hs, errors[:, m+1], marker = :circle, label = "derivative $m")
end| derivative 0 | derivative 1 | derivative 2 | derivative 3 | |
|---|---|---|---|---|
| observed convergence order | 7.62 | 6.04 | 4.99 | 3.94 |
In more dimensions derivative takes one order per axis. A single number applies to every axis, which gives the mixed partial derivative; a tuple gives each axis its own order, with zero for the axes that are only interpolated:
xg = range(0, 2π, length = 60) # knots along x
yg = range(0, 1, length = 30) # knots along y
data = [sin(xi) * exp(yi) for xi in xg, yi in yg] # sin(x)·exp(y) on the grid
dx = convolution_interpolation((xg, yg), data; derivative = (1, 0)) # ∂/∂x
dxy = convolution_interpolation((xg, yg), data; derivative = (1, 1)) # ∂²/∂x∂y
dxx = convolution_interpolation((xg, yg), data; derivative = (2, 0)) # ∂²/∂x²
dx(1.0, 0.5), dxy(1.0, 0.5), dxx(1.0, 0.5)(0.8908079058783364, 0.8908079058756695, -1.387351006247964)xg = np.linspace(0.0, 2 * np.pi, 60) # knots along x
yg = np.linspace(0.0, 1.0, 30) # knots along y
data = np.sin(xg)[:, None] * np.exp(yg)[None, :] # sin(x)·exp(y) on the grid
dx = convolution_interpolation((xg, yg), data, derivative=(1, 0)) # ∂/∂x
dxy = convolution_interpolation((xg, yg), data, derivative=(1, 1)) # ∂²/∂x∂y
dxx = convolution_interpolation((xg, yg), data, derivative=(2, 0)) # ∂²/∂x²
dx(1.0, 0.5), dxy(1.0, 0.5), dxx(1.0, 0.5)| at (1, ½) | ∂f/∂x | ∂²f/∂x∂y | ∂²f/∂x² |
|---|---|---|---|
| interpolant | 0.890807905878 | 0.890807905876 | -1.387351006248 |
| exact | 0.890807904293 | 0.890807904293 | -1.387351111330 |
The kernel can also be chosen per axis, kernel = (:b9, :b5), when one axis needs higher derivatives than another and the wider stencil is only worth paying for there.
Derivatives work on nonuniform grids too, in the Julia package, with the nonuniform kernels’ own differentiated weights; the nonuniform page has the details.