# FastChebInterp

Fast multidimensional Chebyshev interpolation on a hypercube (Cartesian-product) domain, using a separable (tensor-product) grid of Chebyshev interpolation points, as well as Chebyshev regression (least-square fits) from an arbitrary set of points. In both cases we support arbitrary dimensionality, complex and vector-valued functions, and fast derivative and Jacobian computation.

## Usage

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 `SVector` from [StaticArrays.jl](https://github.com/JuliaArrays/StaticArrays.jl), or scalars in 1d
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`

. (You can also use automatic differentiation, e.g. via the ForwardDiff.jl package,
but `chebgradient`

is slightly faster and also supports derivatives of complex-valued functions, unlike ForwardDiff. ChainRules are also implemented in this package to speed up derivative computations using AD tools like Zygote.jl.)

The FastChebInterp package also supports complex and vector-valued functions `f`

. In
the latter 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`

.

### Regression from arbitrary points

We also have a function `chebregression(x, y, [lb, ub,], order)`

that
can perform multidimensional Chebyshev least-square fitting. It
returns a Chebyshev polynomial of a given `order`

(tuple) fit
to a set of points `x[i]`

and values `y[i]`

, optionally in a box
with bounds `lb, ub`

(which default to bounding box for `x`

).

### 1d Example

Here is an example interpolating the (highly oscillatory) 1d function `f(x) = sin(2x + 3cos(4x))`

for `0 ≤ x ≤ 10`

, with a degree-200 Chebyshev polynomial

f(x) = sin(2x + 3cos(4x))
x = chebpoints(200, 0, 10)
c = chebinterp(f.(x), 0, 10)

We can then compare the "exact" function and its derivative at a set of points:

julia> xx = 0:0.1:10;
julia> maximum(@. abs(c(xx) - f(xx)))
2.6336643132285342e-5

and we see that the interpolant `c`

matches `f`

to about five decimal digits.

The function `chebgradient`

returns both the interpolant and its derivative, e.g. at `x = 1.234`

, and we can compare it to the exact values

julia> chebgradient(c, x)
(0.008334535719968672, -13.700695443638699)
julia> f(x) # exact function value
0.008336024670192028
julia> cos(2x + 3cos(4x)) * (2 - 12sin(4x)) # exact derivative
-13.700760631142602

Interpolation is most efficient and accurate if we evaluate our function at the points given by `chebpoints`

. However, we can also perform least-square polynomial fitting (in the Chebyshev basis, which is well behaved even at high degree) from an *arbitrary* set of points — this is useful if the points were specified externally, or if we want to "smooth" the data by fitting to a polynomial of lower degree than for interpolation. For example, we can fit the same function above, again to a degree-200 Chebyshev polynomial, using 10000 *random* points in the domain:

xr = rand(10000) * 10 # 10000 uniform random points in [0, 10]
c = chebregression(xr, f.(xr), 0, 10, 200) # fit to a degree-200 polynomial
```jl
julia> maximum(@. abs(c(xx) - f(xx)))
1.4655330320523241e-5

### 2d Example

Here is a 2d example, interpolating the function `g(x) = sin(x₁ + cos(x₂))`

for `1 ≤ x₁ ≤ 2`

and `3 ≤ x₂ ≤ 4`

, using order 10 in the `x₁`

direction and order 20 in the `x₂`

direction:

g(x) = sin(x[1] + cos(x[2]))
lb, ub = [1,3], [2, 4] # lower and upper bounds of the domain, respectively
x = chebpoints((10,20), lb, ub)
c = chebinterp(g.(x), lb, ub)

Let's evaluate the interpolant at an arbitrary point `(1.2, 3.4)`

and compare it to the exact value:

julia> c([1.2, 3.4]) # polynomial interpolant
0.23109384193446084
julia> g([1.2, 3.4]) # exact value
0.23109384193445792
julia> g([1.2, 3.4]) - c([1.2, 3.4])
-2.914335439641036e-15

In this case, because the function is smooth and not very wiggly in the domain, our low-degree polynomial suffices to interpolate to nearly machine (`Float64`

) precision.

Note that FastChebInterp uses StaticArrays.jl internally to work with vectors, and e.g. a `Vector`

like `[1.2, 3.4]`

gets converted to an `SVector`

internally. If you are working with many such coordinate vectors, it is often advisable to use StaticArrays yourself, in which case you can pass `SVector`

coordinates directly to FastChebInterp functions.

If we inspect the `c`

object, we see that it is actually of a lower degree than we requested:

julia> c
Chebyshev order (10, 16) interpolator on [1,2] × [3,4]

What happened is that `chebinterp`

computed the order-`(10,20)`

polynomial as requested, but found that the order `> 16`

terms in the `x₂`

direction were all smaller than machine precision, so it discarded them (since lower-degree polynomials are cheaper to work with). You can control this behavior by passing the `tol`

keyword argument to `chebinterp`

: passing `tol=0`

prevents it from discarding any terms, and a larger `tol`

can be used to obtain an even lower-degree polynomial:

julia> chebinterp(g.(x), lb, ub, tol=0) # prevent terms from being dropped
Chebyshev order (10, 20) interpolator on [1,2] × [3,4]
julia> chebinterp(g.(x), lb, ub, tol=0.01) # request only about 1% accuracy
Chebyshev order (2, 2) interpolator on [1,2] × [3,4]

The `chebgradient`

function now returns the function and its *gradient* (partial derivatives with respect to `x₁`

and `x₂`

) vector, which we can compare to the analytical gradient:

julia> chebgradient(c, [1.2, 3.4])
(0.23109384193446084, [0.9729314653248615, 0.248623978845799])
g(x) = sin(x[1] + cos(x[2]))
julia> [cos(1.2 + cos(3.4)), cos(1.2 + cos(3.4)) * -sin(3.4)] # exact gradient
2-element Vector{Float64}:
0.9729314653252673
0.24862397884579854

and we see that the derivative matches to high accuracy.

## Related packages

This package was inspired by functionality in ChebyshevApprox.jl, but was rewritten in order to get more performance and other features. The ApproxFun.jl package also performs Chebyshev interpolation and many other tasks. BasicInterpolators.jl also provides Chebyshev interpolation in 1d and 2d, and Surrogates.jl provides some other interpolation schemes.