Quadratic Interpolation API

quadratic_interp(x, y, xi; bc=Left(QuadraticFit()), extrap=:none, deriv=0, search=Binary())

C1 piecewise quadratic spline interpolation at a single point.

Arguments

• x::AbstractVector: x-coordinates (sorted, length ≥ 2)
• y::AbstractVector: y-values (same length as x)
• xi::Real: Query point
• bc: Boundary condition (one of):
  • Left(QuadraticFit()): 3-point parabola fit at left (default, exact for polynomials)
  • Right(QuadraticFit()): 3-point parabola fit at right
  • Left(Deriv1(v)): First derivative = v at left endpoint
  • Left(Deriv2(v)): Second derivative = v at left endpoint
  • Right(Deriv1(v)): First derivative = v at right endpoint
  • Right(Deriv2(v)): Second derivative = v at right endpoint
  • MinCurvFit(): Minimize total curvature (globally smooth)
• extrap::Symbol: Extrapolation mode
  • :none (default): throws DomainError if outside domain
  • :constant: clamp to boundary values
  • :extension: extend the boundary polynomial
• deriv::Int: Derivative order (0, 1, or 2)
• search::AbstractSearchPolicy: Search algorithm for interval finding
  • Binary() (default): O(log n) binary search, stateless
  • HintedBinary(): O(1) if hint valid, O(log n) fallback
  • LinearBinary(linear_window=8): Linear search within window, then binary fallback

Returns

• Interpolated value (Float type)

Example

x = [0.0, 1.0, 2.0, 3.0]
y = x.^2  # [0, 1, 4, 9]

# Default: QuadraticFit (exact for quadratic polynomials)
quadratic_interp(x, y, 1.5)  # ≈ 2.25 (exact)

# With specific BC
quadratic_interp(x, y, 1.5; bc=Left(Deriv1(0.0)))  # zero slope at left
quadratic_interp(x, y, 1.5; bc=MinCurvFit())        # minimize curvature

# Optimized for sorted queries
sorted_queries = sort(rand(1000))
vals = quadratic_interp(x, y, sorted_queries; search=LinearBinary(linear_window=8))

quadratic_interp(x, y, x_targets; bc=Left(QuadraticFit()), extrap=:none, deriv=0, search=Binary())

Quadratic spline interpolation for multiple query points (allocating version).

Example

x = [0.0, 1.0, 2.0, 3.0]
y = x.^2
result = quadratic_interp(x, y, [0.5, 1.5, 2.5])
# result ≈ [0.25, 2.25, 6.25]

# Optimized for sorted queries
sorted_queries = sort(rand(1000))
vals = quadratic_interp(x, y, sorted_queries; search=LinearBinary(linear_window=8))

quadratic_interp(x, Y::AbstractMatrix; bc=Left(QuadraticFit()), extrap=:none)

Create a multi-Y quadratic series interpolant from a matrix where each column is a y-series.

Arguments

• x::AbstractVector: x-coordinates (length n)
• Y::AbstractMatrix: n×m matrix, each column is a y-series
• bc, extrap: Same as vector form

Example

x = collect(range(0.0, 1.0, 101))
Y = hcat(sin.(2π .* x), cos.(2π .* x))  # 101×2 matrix

sitp = quadratic_interp(x, Y)

quadratic_interp(grids::NTuple{N,AbstractVector}, data::AbstractArray{<:Any,N}; kwargs...)

Create an N-dimensional quadratic interpolant from grid vectors and data array.

Arguments

• grids::NTuple{N,AbstractVector}: Tuple of grid vectors for each dimension
• data::AbstractArray{<:Any,N}: Function values at grid points

Keywords

• bc=Left(QuadraticFit()): Boundary condition(s). Can be:
  • Single BC: Applied to all axes
  • NTuple{N}: Per-axis BCs
• extrap=:none: Extrapolation mode(s)
• search=Binary(): Search policy(s)

Returns

• QuadraticInterpolantND{Tg, Tv, N, ...}: Callable interpolant object

Examples

x = range(0.0, 2.0, 20)
y = range(0.0, 1.0, 15)
data = [xi^2 + yi^2 for xi in x, yi in y]
itp = quadratic_interp((x, y), data)
itp((1.0, 0.5))  # Evaluate
itp((1.0, 0.5); deriv=(1, 0))  # ∂f/∂x

quadratic_interp(grids, data, query; deriv=0, kwargs...)

One-shot ND quadratic interpolation at a single point.

quadratic_interp(grids, data, queries; deriv=0, kwargs...)

One-shot ND quadratic interpolation at multiple points (batch).

quadratic_interp!(output, x, y, x_targets; bc=Left(QuadraticFit()), extrap=:none, deriv=0, search=Binary())

In-place quadratic spline interpolation for multiple query points.

Arguments

• output: Pre-allocated output vector
• Other arguments same as quadratic_interp

Example

x = [0.0, 1.0, 2.0, 3.0]
y = x.^2
out = zeros(3)
quadratic_interp!(out, x, y, [0.5, 1.5, 2.5])
# out ≈ [0.25, 2.25, 6.25]

# Optimized for sorted queries
sorted_queries = sort(rand(1000))
output = zeros(1000)
quadratic_interp!(output, x, y, sorted_queries; search=LinearBinary(linear_window=8))

QuadraticInterpolant{Tg,Tv,X,Y,P}

Lightweight callable interpolant for quadratic spline interpolation. Returned by quadratic_interp(x, y) (2-argument form).

Type Parameters

• Tg<:AbstractFloat: Grid type (Float32, Float64) for x-coordinates
• Tv: Value type for y-values (can be Tg, Complex{Tg}, or other Number)
• X<:AbstractVector{Tg}: Type of x-coordinates
• Y<:AbstractVector{Tv}: Type of y-values
• P<:AbstractSearchPolicy: Search policy type

Fields

• x::X: x-coordinates (sorted)
• y::Y: y-values
• h::Vector{Tg}: Grid spacing (precomputed, geometry)
• a::Vector{Tv}: Quadratic coefficients (value-derived)
• d::Vector{Tv}: Slope coefficients (value-derived)
• extrap::ExtrapVal: Extrapolation mode
• search_policy::P: Default search policy for interval lookup

Usage

itp = quadratic_interp(x, y; bc=Right(Deriv1(6.0)))
val = itp(0.5)               # scalar evaluation
vals = itp.([0.5, 1.5])      # broadcast
vals = itp([0.5, 1.5])       # vector call

# Derivatives
d1 = itp(0.5; deriv=1)       # first derivative
d2 = itp(0.5; deriv=2)       # second derivative

# Complex values
x = [0.0, 1.0, 2.0, 3.0]
y = [1.0+2.0im, 3.0+4.0im, 5.0+6.0im, 7.0+8.0im]
itp = quadratic_interp(x, y)
val = itp(0.5)  # returns ComplexF64

# Custom search policy
itp = quadratic_interp(x, y; search=LinearBinary())
val = itp(0.5)               # uses LinearBinary() by default
val = itp(0.5; search=Binary())  # override with Binary()

QuadraticInterpolantND{Tg, Tv, N, NP1, G, S, B, E, P}

Generic N-dimensional quadratic interpolant with precomputed partial derivatives.

Stores function values AND all partial derivatives at grid nodes, enabling ultra-fast O(1) evaluation via tensor-product quadratic polynomials.

Type Parameters

• Tg: Grid/coordinate type (Float32 or Float64)
• Tv: Value type (can be Tg, Complex{Tg}, or other Number)
• N: Number of dimensions
• NP1: N + 1 (partials array dimensionality)
• G: Tuple type for grids
• S: Tuple type for spacings
• B: Tuple type for boundary conditions
• E: Tuple type for extrapolation modes
• P: Tuple type for search policies

Fields

• grids: N-tuple of grid vectors for each dimension
• spacings: N-tuple of grid spacing info (for O(1) h lookup)
• nodal_derivs: NodalDerivativesND containing partial derivatives at grid nodes
• bcs: N-tuple of boundary conditions used for construction
• extraps: N-tuple of extrapolation modes
• searches: N-tuple of search policies

Performance

• Construction: O(2^N x n^N) - computes all partial derivatives via recurrence
• Query: O(1) - tensor-product quadratic polynomial evaluation
• Memory: 2^N x n^N values

Thread-Safety

Immutable after construction; safe for concurrent read access.

Example

x = range(0.0, 2.0, 20)
y = range(0.0, 1.0, 15)
data = [xi^2 + yi^2 for xi in x, yi in y]

itp = quadratic_interp((x, y), data)
itp((1.0, 0.5))                  # Evaluate at (1.0, 0.5)
itp((1.0, 0.5); deriv=(1, 0))   # dfdx

MinCurvFit <: AbstractBC

Minimum curvature boundary condition for quadratic splines. Minimizes total curvature (∫S''² dx) by optimizing the initial slope d[1].

This is a singleton type (structure-only). The optimization is performed during interpolant construction based on the actual data values.

Mathematical Background

For quadratic splines, the slope d[i] depends on d[1] via recurrence:

• d[i] = α[i] * d[1] + β[i] where α[i] = (-1)^(i+1)

The optimal d[1] minimizes:

• ∫(S'')² dx = Σ 4a[i]²h[i] = Σ (s[i] - d[i])²/h[i]

Closed-form solution:

• d[1] = [Σ α[i]*(s[i] - β[i])/h[i]] / [Σ 1/h[i]]

Example

x = [0.0, 0.3, 0.8, 1.5, 2.5, 3.0, 4.0]
y = [0.0, 0.8, 1.2, 0.9, 0.3, 0.6, 1.0]

# Default BC uses QuadraticFit
itp_default = quadratic_interp(x, y)

# MinCurvFit gives globally smooth result via curvature minimization
itp_smooth = quadratic_interp(x, y; bc=MinCurvFit())