Julia code to approximate continuous functions using Chebyshev polynomials.
In many applications it is important to be able to approximate accurately continuous functions of several variables. This code can approximate functions that depend upon an arbitrary number of variables.
Chebyshev nodes in the [1.0, -1.0]
interval are computed using the command
nodes = chebyshev_nodes(n)
where n
, an integer, is the number of nodes. Similarly, Chebyshev nodes scaled to a bounded interval are computed using the command
nodes = chebyshev_nodes(n,domain)
where domain
is a 1D array with two elements, the first representing the upper bound on the interval and the second the lower bound.
To approximate a 1D function using a Chebyshev polynomial one first computes the Chebyshev weights. There are several ways of doing this. The central command is
w = chebyshev_weights(y,nodes,order,domain)
where y
is a 1D array containing function evaluations at nodes and order is an integer (or 1D integer array containing one element) representing the order of the polynomial. nodes
can be either a 1D array or a tuple containing a 1D array. An alternative is the command
w = chebyshev_weights(y,poly,order,domain)
where poly
is an array containing a 2D array containing the Chebyshev polynomials evaluated at each node. Regardless of which command is used, the Chebyshev weights are computed using Chebyshev regression.
Given the Chebyshev weights, the function can be approximated at a given point, x
, within domain using either
y_approx = chebyshev_evaluate(w,[x],order,domain)
or
y_approx = clenshaw_evaluate(w,[x],order,domain)
The latter command evaluates the polynominal using Clenshaw's recursion. For functions where the domain of x
coincides with the [1.0, -1.0]
interval the domain argument can be omitted, both when computing the weights and when evaluating the polynomial.
The codes can approximate functions of an arbitrary number of variables. The only real difference is that in the higher dimension cases the codes allow the polynomials to be either tensor products of 1D polynomials or complete polynomials. When the function depends on six or fewer variables, and in the tensor-product case, to approximate a function (of three variables, say) the relevant commands are
w = chebyshev_weights(y,nodes_1,nodes_2,nodes_3,order,domain) --- (1)
where order
is a 1D integer array whose elements represents the polynomial orders in each successive dimension, or
w = chebyshev_weights(y,nodes,order,domain) --- (2)
where nodes
is a tuple containing 1D arrays, or
w = chebyshev_weights(y,poly,order,domain) --- (3)
where poly
is a tuple containing 2D arrays of Chebyshev polynomials evaluated at each node. To go beyond functions of six variables, commands (2) or (3) can be used, but not (1).
Then to evaluate the polynomial at a point x
, the relevant commands are either
y_approx = chebyshev_evaluate(w,x,order,domain)
or
y_approx = clenshaw_evaluate(w,x,order,domain)
In each of these commands, the variable domain
is a 2D array (2*3 matrix for the 3-variable case) whose columns summarize the upper- and lower-bound on each variable.
For the case where a complete polynomial rather than a tensor-product polynomial is desired, the commands are the same as above, but the order
variable is now simply an integer rather than a 1D array of integers.
As with the 1D case, if the domain of each variable is [1.0, -1.0]
, then the domain variable can be omitted.
Working with complete polynomials rather than tensor-product polynomials often leads to a considerable decrease in computation time with little loss of accuracy.
Have fun.
10/15/2014
about 2 months ago
110 commits