Fast multidimensional Chebyshev interpolation on a hypercube (Cartesian-product) domain, using a separable (tensor-product) grid of Chebyshev interpolation points.

For domain upper and lower bounds lb and ub, and a given order tuple, you would create an interpolator for a function f via:

using FastChebInterp
x = chebpoints(order, lb, ub) # an array of StaticVector
c = chebinterp(f.(x), lb, ub)

and then evaluate the interpolant for a point y (a vector) via c(y).

We also provide a function chebgradient(c,y) that returns a tuple (c(y), ∇c(y)) of the interpolant and its gradient at a point y.

The FastChebInterp package also supports complex and vector-valued functions f. In this case, c(y) returns a vector of interpolants, and one can use chebjacobian(c,y) to compute the tuple (c(y), J(y)) of the interpolant and its Jacobian matrix at y.

ChebPoly

A multidimensional Chebyshev-polynomial interpolation object. Given a c::ChebPoly, you can evaluate it at a point x with c(x), where x is a vector (or a scalar if c is 1d).

Jevaluate(x, c::Array{T,N}, ::Val{dim}, i1, len)

Similar to evaluate above, but returns a tuple (v,J) of the evaluated value v and the Jacobian J with respect to x[1:dim].

chebgradient(c::ChebPoly, x)

Return a tuple (v, ∇v) where v is the value c(x) of the Chebyshev polynomial c at x, and ∇v is the gradient of this value with respect to x.

(Requires c to be a scalar-valued polynomial; if c is vector-valued, you should use chebjacobian instead.)

If x is a scalar, returns a scalar ∇v, i.e. `(v, ∂v/∂x).

chebinterp(vals, lb, ub; tol=eps)

Given a multidimensional array vals of function values (at points corresponding to the coordinates returned by chebpoints), and arrays lb and ub of the lower and upper coordinate bounds of the domain in each direction, returns a Chebyshev interpolation object (a ChebPoly).

This object c = chebinterp(vals, lb, ub) can be used to evaluate the interpolating polynomial at a point x via c(x).

The elements of vals can be vectors as well as numbers, in order to interpolate vector-valued functions (i.e. to interpolate several functions at once).

The tol argument specifies a relative tolerance below which Chebyshev coefficients are dropped; it defaults to machine precision for the precision of float(vals). Passing tol=0 will keep all coefficients up to the order passed to chebpoints.

chebinterp_v1(vals, lb, ub; tol=eps)

Like chebinterp(vals, lb, ub), but slices off the first dimension of vals and treats it as a vector of values to interpolate.

For example, if vals is a 2×31×32 array of numbers, then it is treated equivalently to calling chebinterp with a 31×32 array of 2-component vectors.

(This function is mainly useful for calling from Python, where arrays of vectors are painful to construct.)

chebjacobian(c::ChebPoly, x)

Return a tuple (v, J) where v is the value c(x) of the Chebyshev polynomial c at x, and J is the Jacobian of this value with respect to x.

That is, if v is a vector, then J is a matrix of length(v) × length(x) giving the derivatives of each component. If v is a scalar, then J is a 1-row matrix; in this case you may wish to call chebgradient instead.

chebpoints(order, lb, ub)

Return an array of Chebyshev points (as SVector values) for the given order (an array or tuple of polynomial degrees), and hypercube lower-and upper-bound vectors lb and ub. If ub and lb are numbers, returns an array of numbers.

These are the points where you should evaluate a function in order to create a Chebyshev interpolant with chebinterp.

(Note that the number of points along each dimension is 1 + the order in that direction.)

chebregression(x, y, [lb, ub,] order)

Return a Chebyshev polynomial (ChebPoly) constructed by performing a least-square fit of Chebyshev polnomials of the given order, where x are the coordinates of the data points y. lb and ub are the lower and upper bounds, respectively, of the Chebyshev domain; these should normally enclose all of the points in x, and default to the minimum and maximum coordinates in x if they are omitted.

In the 1d case, x is an array of scalars, lb < ub are scalars, and order is an integers. In the N-dimensional case, order is an N-tuple of integers (the order in each dimension), lb and ub are N-component vectors, and x is an array of N-component vectors (or a matrix with N columns, interpreted as the vector components).

y can be a vector or numbers or a vector of vectors (for vector- valued Chebyshev fits). The latter case can also be input as a matrix whose columns are the vector componnents. size(x,1) and size(y,1) must match, and must exceed prod(order .+ 1).

evaluate(x, c::Array{T,N}, ::Val{dim}, i1, len)

Low-level polynomial-evaluation function, which performs a multidimensional Clenshaw recurrence by recursing on the coefficient (c) array dimension dim (from N to 1). The current dimension (via column-major order) is accessed by c[i1 + (i-1)*Δi], i.e. i1 is the starting index and Δi = len ÷ size(c,dim) is the stride. len is the product of size(c)[1:dim]. The interpolation point x should lie within [-1,+1] in each coordinate.