[ANN] ConvolutionInterpolations.jl: Smooth multi-dimensional high order of accuracy interpolation

Package Announcements
,
81

3D Vector Line Integral: From Discrete Samples to Machine Precision

Note (edited): The initial results I posted was made with a corrupted version of the package, giving only 3rd order convergence. With the correct released version, 7th order convergence is recovered.

This is a quick demonstration of application of ConvolutionInterpolations.jl (v0.15.0) to compute a 3D vector line integral entirely from discrete samples, using interpolation, differentiation, and integration from the same package.

The problem

The task is to compute the line integral:

\int_C \mathbf{F} \cdot d\mathbf{r}

where \mathbf{F}(x,y,z) = (yz,\, xz,\, xy) is a conservative vector field with potential \varphi = xyz, and the curve is a helix

\mathbf{r}(t) = \bigl(\cos t,\, \sin t,\, t/(2\pi)\bigr), \quad t \in [0, 2\pi].

Since \mathbf{F} = \nabla\varphi, the exact answer is

\int_C \mathbf{F} \cdot d\mathbf{r} = \varphi\bigl(\mathbf{r}(2\pi)\bigr) - \varphi\bigl(\mathbf{r}(0)\bigr) = (1)(0)(1) - (1)(0)(0) = 0.

This makes the absolute error trivial to measure: any nonzero result is pure numerical error.

The pipeline

The calculation proceeds in five steps, all handled by ConvolutionInterpolations.jl:

  1. 3D interpolation: Sample \mathbf{F} on a uniform 3D grid with M points per dimension, build three 3D interpolants (one per field component) using :b7 kernels.
  2. 1D interpolation: Sample the curve \mathbf{r}(t) at N discrete points, build 1D interpolants for x(t), y(t), z(t).
  3. 1D differentiation: Build derivative interpolants for the tangent vector \mathbf{r}'(t) = (x'(t),\, y'(t),\, z'(t)) using derivative=1.
  4. Evaluation: At each point on a fine t grid, evaluate \mathbf{F}\bigl(\mathbf{r}(t)\bigr) \cdot \mathbf{r}'(t) by composing the 3D and 1D interpolants.
  5. 1D integration: Build an antiderivative interpolant of the integrand using derivative=-1, evaluate at the endpoints.

The code

using ConvolutionInterpolations
using Printf

# Sample the vector field F(x,y,z) = (yz, xz, xy) on a 3D grid
function sample_field(M)
    xs = range(-1.5, 1.5, length=M)
    ys = range(-1.5, 1.5, length=M)
    zs = range(-0.5, 1.5, length=M)
    F1 = [y*z for x in xs, y in ys, z in zs]
    F2 = [x*z for x in xs, y in ys, z in zs]
    F3 = [x*y for x in xs, y in ys, z in zs]
    return xs, ys, zs, F1, F2, F3
end

# Compute the line integral for given N (curve samples) and M (field grid)
function line_integral(N, M; kernel_3d=:b7, kernel_1d=:b7)
    # Step 1: 3D interpolants for each field component
    xs, ys, zs, F1, F2, F3 = sample_field(M)
    itp_F1 = convolution_interpolation((xs, ys, zs), F1; kernel=kernel_3d)
    itp_F2 = convolution_interpolation((xs, ys, zs), F2; kernel=kernel_3d)
    itp_F3 = convolution_interpolation((xs, ys, zs), F3; kernel=kernel_3d)

    # Step 2: 1D interpolants for the curve r(t) = (cos(t), sin(t), t/(2pi))
    t_grid = range(0.0, 2pi, length=N)
    x_data = cos.(t_grid)
    y_data = sin.(t_grid)
    z_data = collect(t_grid) ./ (2pi)

    itp_x = convolution_interpolation(t_grid, x_data; kernel=kernel_1d)
    itp_y = convolution_interpolation(t_grid, y_data; kernel=kernel_1d)
    itp_z = convolution_interpolation(t_grid, z_data; kernel=kernel_1d)

    # Step 3: 1D derivative interpolants for the tangent vector r'(t)
    itp_dx = convolution_interpolation(t_grid, x_data; kernel=kernel_1d, derivative=1)
    itp_dy = convolution_interpolation(t_grid, y_data; kernel=kernel_1d, derivative=1)
    itp_dz = convolution_interpolation(t_grid, z_data; kernel=kernel_1d, derivative=1)

    # Step 4: Evaluate integrand F(r(t)) . r'(t) on a fine grid
    N_fine = 4 * N
    t_fine = range(0.0, 2pi, length=N_fine)
    integrand = zeros(N_fine)
    for (i, t) in enumerate(t_fine)
        cx, cy, cz = itp_x(t), itp_y(t), itp_z(t)
        f1, f2, f3 = itp_F1(cx, cy, cz), itp_F2(cx, cy, cz), itp_F3(cx, cy, cz)
        dx, dy, dz = itp_dx(t), itp_dy(t), itp_dz(t)
        integrand[i] = f1*dx + f2*dy + f3*dz
    end

    # Step 5: Integrate via antiderivative
    itp_int = convolution_interpolation(t_fine, integrand; kernel=kernel_1d, derivative=-1)
    return itp_int(2pi) - itp_int(0.0)
end

Results

With M = 20 (field grid) held fixed and sweeping the curve resolution N:

N_curve  M_field  |  |error|       rate
  -------  -------  |  -----------  ------
       8      20  |   7.131e-05      -
      10      20  |   2.364e-05    4.95
      12      20  |   9.937e-06    4.75
      15      20  |   1.398e-05   -1.53
      20      20  |   8.758e-07    9.63
      25      20  |   4.564e-07    2.92
      30      20  |   3.160e-11   52.53
      40      20  |   2.231e-11    1.21
      50      20  |   6.192e-12    5.74
      60      20  |   1.917e-12    6.43
      80      20  |   2.732e-13    6.77
     100      20  |   5.793e-14    6.95
     160      20  |   1.659e-15    7.56

Clean 7th-order convergence from N = 60 through N = 160, reaching 1.7 \times 10^{-15} at N = 160, essentially machine precision. The jump around N = 25-30 is where the :detect boundary condition polynomial fit gains enough points to start working at full accuracy.

Note on M: the vector field \mathbf{F} = (yz, xz, xy) is polynomial (degree 2), so the :b7 kernel reproduces it exactly regardless of grid resolution. The 3D interpolation is exact to machine precision even at M = 10. All convergence above comes purely from the 1D curve interpolation, differentiation, and integration steps.

Discussion

The observed 7th-order convergence is consistent with the :b7 kernel’s interpolation order. Differentiation loses one order (to 6th), but the antiderivative step recovers it, giving 7th order overall for the integrated result.

The pipeline composes naturally: 3D interpolation for the field, 1D interpolation + differentiation for the curve geometry, and 1D integration for the final scalar. All from the same package, using the same kernel family and the same API (convolution_interpolation with kernel, derivative).

More broadly, vector calculus is a natural application of the per-dimension derivative support. Given a 3D field sampled on a grid, all the standard differential operators (\nabla\mathbf{F}, \nabla\cdot\mathbf{F}, \nabla\times\mathbf{F}, \nabla^2\mathbf{F}) follow directly:

# gradient of a scalar field f(x,y,z)
df_dx = convolution_interpolation((xs, ys, zs), f; kernel=:b7, derivative=(1,0,0))
df_dy = convolution_interpolation((xs, ys, zs), f; kernel=:b7, derivative=(0,1,0))
df_dz = convolution_interpolation((xs, ys, zs), f; kernel=:b7, derivative=(0,0,1))
# grad(f) = [df_dx, df_dy, df_dz]

# divergence of a vector field F = (F1, F2, F3)
dF1_dx = convolution_interpolation((xs, ys, zs), F1; kernel=:b7, derivative=(1,0,0))
dF2_dy = convolution_interpolation((xs, ys, zs), F2; kernel=:b7, derivative=(0,1,0))
dF3_dz = convolution_interpolation((xs, ys, zs), F3; kernel=:b7, derivative=(0,0,1))
# div(F) = dF1_dx(x,y,z) + dF2_dy(x,y,z) + dF3_dz(x,y,z)

# curl of a vector field F = (F1, F2, F3)
dF3_dy = convolution_interpolation((xs, ys, zs), F3; kernel=:b7, derivative=(0,1,0))
dF2_dz = convolution_interpolation((xs, ys, zs), F2; kernel=:b7, derivative=(0,0,1))
dF1_dz = convolution_interpolation((xs, ys, zs), F1; kernel=:b7, derivative=(0,0,1))
dF3_dx = convolution_interpolation((xs, ys, zs), F3; kernel=:b7, derivative=(1,0,0))
dF2_dx = convolution_interpolation((xs, ys, zs), F2; kernel=:b7, derivative=(1,0,0))
dF1_dy = convolution_interpolation((xs, ys, zs), F1; kernel=:b7, derivative=(0,1,0))
# curl(F) = (dF3_dy - dF2_dz, dF1_dz - dF3_dx, dF2_dx - dF1_dy)

# laplacian of a scalar field f(x,y,z)
d2f_dx2 = convolution_interpolation((xs, ys, zs), f; kernel=:b7, derivative=(2,0,0))
d2f_dy2 = convolution_interpolation((xs, ys, zs), f; kernel=:b7, derivative=(0,2,0))
d2f_dz2 = convolution_interpolation((xs, ys, zs), f; kernel=:b7, derivative=(0,0,2))
# laplacian(f) = d2f_dx2(x,y,z) + d2f_dy2(x,y,z) + d2f_dz2(x,y,z)

Each of these returns a callable functor that can be evaluated at any point in the domain.

1 Like