ChebyshevApprox

ChebyshevApprox.jl is a Julia package for approximating continuous functions using Chebyshev polynomials. The package's focus is on multivariate functions that depend on an arbitrary number of variables. Both tensor-product polynomials and complete polynomials are implemented. Working with complete polynomials often leads to a considerable decrease in computation time with little loss of accuracy. The package allows the nodes to be either the roots of the Chebyshev polynomial (points of the first kind), the extrema of the Chebyshev polynomial (points of the second kind), or the Chebyshev extended points (Chebyshev roots scaled so that the boundry nodes equal -1.0 and 1.0). In addition to approximating functions, the package also uses the approximating polynomial to compute derivatives, gradients, and hessians. The package allows you to work with high precision floating point numbers, such as BigFloats, and is compatable with packages such at DoubleFloats.jl.

Installation

ChebyshevApprox.jl is a registered package. To install it simply type in the REPL:

using Pkg
Pkg.add("ChebyshevApprox")

Nodes

The package contains functions for computing several different types of approximating points. Which points you use may depend on your application.

To compute the Chebyshev roots within the [1.0, -1.0] interval use:

n = 11
points = nodes(n,:chebyshev_nodes)

where n, an integer, is the number of nodes and :chebyshev_nodes is a symbol indicating the type of nodes to be produced. Alternatives to :chebyshev_nodes are: chebyshev_extrema, and :chebyshev_extended. If you want these roots to be computed to high precision, then:

n = 11
points = nodes(BigFloat,n,:chebyshev_nodes)

or

using DoubleFloats
n = 11
points = nodes(Double64,n,:chebyshev_nodes)

analogously for the Chebyshev extrema and the extended Chebyshev points.

To compute nodes over bounded domains other than the [1.0,-1.0] interval, the nodes function accepts an optional third argument containing the domain in the form of a 1D array (a vector) containing two elements, where the first element is the upper bound on the interval and the second is the lower bound. For example,

n = 11
domain = [3.5,0.5]
points = nodes(n,:chebyshev_extrema,domain)
points = nodes(Double64,n,:chebyshev_extrema,domain)

would compute n extrema of the Chebyshev polynomial and scale those points to the [3.5,0.5] interval, using Double64 floating points number is the second case.

Grids

Once the approximating points have been computed, these points can be collected together in a Grid structure:

p1 = nodes(11,:chebyshev_nodes)
p2 = nodes(15,:chebyshev_nodes)

g = Grid((p1,p2))

Approximation Plans

The approximate a function the package makes use of an approximation plan, which is a stucture containing the information needed to produce an approximating function.

A_plan = CApproxPlan(g,order,domain)

where g is a Grid, order is an integer (complete polynomial approximation) or a tuple (tensor-product approximation), and domain is an array containing the upper and lower limits for each variable. For example,

dom_1 = [25.0,5.0]
dom_2 = [3.0,0.1]
dom = [dom_1 dom_2]

p1 = nodes(11,:chebyshev_nodes,dom_1)
p2 = nodes(15,:chebyshev_nodes,dom_2)

g = Grid((p1,p2))
order = (6,6)
A_plan = CApproxPlan(g,order,dom)

Function approximation

To approximate a function, you use the chebyshev_interp function, which itself returns a function. If the data for the function you wish to approximate, sampled on the Grid, g, are contained in the Array, y, then the approximating function is generated using:

y = [k^0.35*l^0.65 for k in p1.points, l in p2.points]
f_approx = chebyshev_interp(y,A_plan)

which can then be evaluated at a point, x, in the domain according to:

x = [9.0,0.7]
y_hat = f_approx(x)

Functions to approximate gradients and hessians are produced similarly:

f_grad = chebyshev_gradient(y,A_plan)
f_hess = chebyshev_hessian(y,A_plan)
grad_hat = f_grad(x)
hess_hat = f_hess(x)

There are multi-threaded versions of these functions: chebyshev_interp_threaded(), chebyshev_gradient_threaded(), and chebyshev_hessian_threaded().

Polynomials

Chebyshev polynomials are constructed using the chebyshev_polynomial() function, which takes two arguments. The first argument is an integer representing the order of the polynomial. The second argument is the point in the [1.0,-1.0] interval at which the polynominal is evaluated. This second argument can be a scalar or a 1D array. For example,

order = 5
x = 0.5
p = chebyshev_polynomial(order,x)

will return a 2D array (1 $\times$ 6) containing the Chebyshev polynomials of orders 0---5 evaluated at the point x. If x is a 1D array of points, as in:

order = 5
x = chebyshev_nodes(11)
p = chebyshev_polynomial(order,x)

then p will be a 2D array (11 $\times$ 6) containing the Chebyshev polynomials of orders 0---5 evaluated at each element in x.

Weights

ChebyshevApprox.jl uses Chebyshev regression to compute the weights in the Chebyshev polynomial. The central function for computing Chebyshev weights is the following:

w = chebyshev_weights(y,nodes,order,domain)

where y is a nD array containing the function evaluations at nodes, nodes is a tuple of 1D arrays containing Chebyshev-roots, order is a tuple (tensor-product polynomial) or an integer (complete polynomial) specifying the order of the polynomial in each dimension, and domain is a 2D array containing the upper and lower bounds on the approximating interval in each dimension. So,

order_x1  = 5
domain_x1 = [25.0,5.0]
nodes_x1  = chebyshev_nodes(11,domain_x1)

order_x2  = 7
domain_x2 = [3.0,0.1]
nodes_x2  = chebyshev_nodes(15,domain_x2)

order = (order_x1,order_x2)
grid_points = (nodes_x1,nodes_x2)
domain = [domain_x1 domain_x2]

y = [k^0.35*l^0.65 for k in nodes_x1, l in nodes_x2]

w = chebyshev_weights(y,grid_points,order,domain)

would compute the weights, w, (a 2D array in this example) in a tensor-product polynomial. The domain-argument is optional, needed only if one or more variable does not have domain [1.0,-1.0]. The nodes-argument can be an array-of-arrays (instead of a tuple). Alternatively, the polynominals can be computed and entered directly into the chebyshev_weights() function:

poly_1 = chebyshev_polynomial(order_x1,nodes_x1)
poly_2 = chebyshev_polynomial(order_x2,nodes_x2)
poly   = (poly_1,poly_2)

w = chebyshev_weights(y,poly,order)

The poly-argument can be an array-of-arrays (instead of a tuple). The weights, w are returned in a (multi-dimensional) array.

If the solution nodes are instead the Chebyshev-extrema, then the analogue to the above is the use the chebyshev_weights_extrema() function. For example:

order_x1  = 5
domain_x1 = [25.0,5.0]
nodes_x1  = chebyshev_extrema(11,domain_x1)

order_x2  = 7
domain_x2 = [3.0,0.1]
nodes_x2  = chebyshev_extrema(15,domain_x2)

order = [order_x1,order_x2]
grid_points = (nodes_x1,nodes_x2)
domain = [domain_x1 domain_x2]

y = [k^0.35*l^0.65 for k in nodes_x1, l in nodes_x2]

w = chebyshev_weights_extrema(y,grid_points,order,domain)

The third possibility is to use chebyshev_weights_extended().

Function evaluation

ChebyshevApprox uses the chebyshev_evaluate() function, which accommodates several methods, to evaluate Chebyshev polynomials. If x is a 1D array representing the point at which the polynomial is to be evaluated, then:

yhat = chebyshev_evaluate(w,x,order,domain)

For the case where a complete polynomial rather than a tensor-product polynomial is to be evaluated, the order variable is now simply an integer rather than a tuple of integers.

Derivatives

The chebyshev_derivative() function can be used to approximate the partial derivative of a function with respect to a designated variable. For example, the partial derivative with respect to the 3'rd variable evaluated at point x can be computed by:

deriv = chebyshev_derivative(w,x,3,order,domain)

where deriv is a floating point number.

Gradients

Gradients are computed using the chebyshev_gradient() function.

grad = chebyshev_gradient(w,x,order,domain)

where grad is a 2D array with one row.

Multi-threading

Computing the weights in a multivariate Chebyshev polynomial can be time-consuming for functions whose dimensions are large, or where the number of nodes and/or the order of the polynomals is large. For this reason, multi-threaded functions for computing the weights are provided. For example, if the nodes are Chebyshev-roots:

w = chebyshev_weights_threaded(y,nodes,order,domain)

etc.

For completeness, there are also:

w = chebyshev_weights_extrema_threaded(y,nodes,order,domain)
w = chebyshev_weights_extended_threaded(y,nodes,order,domain)

Related packages

If you are looking to approximate a function of one variable, then there is:

• ApproxFun.jl

For multivariate functions, there is:

• SmolyakApprox.jl
• HyperbolicCrossApprox.jl
• PiecewiseLinearApprox.jl

Finally, if you are approximating complex-valued multivariate functions, then a package insipired by ChebyshevApprox.jl and offering similar functionality is:

• FastChebInterp.jl