matlabjuliapolynomialsfunction-approximationpolynomial-approximations

Create vector of complete Chebyshev polynomials in Julia (or MATLAB)?

Suppose we have a two dimensional function f(x,y) we want to approximate with a set of Chebyshev polynomials up to degree 2. Let the Chebyshev polynomial of degree j be Tj(x) or Tj(y). We usually approximate f(x,y) by constructing a function g(x,y) that is the tensor product of the one dimensional polynomials,

        

What I want to do is to generate a complete Chebyshev polynomial of level N. This is just the tensor product like above, but where the sum of the indices k+l must be less than or equal to N. So if N was 3, then we would have all the terms in the sum above except for T2(x)*T2(y) since 2+2=4 > 3. Many more terms get dropped as the dimensionality of the function increases.

Ultimately I'd like to do this with a flexible choice of the level and without having to write out a bunch of nested loops if I am working with more than 2 or 3 dimensions. It seems like @nloops would be the way to go but I can't quite figure it out.

For example, suppose I want to evaluate a 2 dimensional Chebyshev polynomial at (.5,.5). I can write an inline function that returns the one-dimensional Chebyshev polynomial of level N at point x.

cheb(d,x) = cos((d)*acos(x))
a = [.5, .5]
polys1d = [cheb(d,a[i]) for d = 0:2, i = 1:length(a)]

Creating the full tensor product polynomial in 2 dimensions (or even more) is easy. For example:

polys2d = kron(polys1d[:,1],polys1d[:,2])

But creating the complete polynomial in a general way is a bit trickier. I would like it do it in a way where I don't start off by constructing the full tensor product and then removing polynomials where the sum of the degrees is greater than the level. This will take up large amounts of memory if both the number of dimensions and the level are big.

Solution

  • If we define the following helper function:

    using Iterators   # install with Pkg.add("Iterators")

nlimitedkparts(n,k) = (diff([0;v]).-1 for v in subsets(1:(n+k),k))

    We can generate the following:

    julia> ["T_$(r[1])(x)*T_$(r[2])(y)" for r in nlimitedkparts(3,2)]
10-element Array{String,1}:
 "T_0(x)*T_0(y)"
 "T_0(x)*T_1(y)"
 "T_0(x)*T_2(y)"
 "T_0(x)*T_3(y)"
 "T_1(x)*T_0(y)"
 "T_1(x)*T_1(y)"
 "T_1(x)*T_2(y)"
 "T_2(x)*T_0(y)"
 "T_2(x)*T_1(y)"
 "T_3(x)*T_0(y)"

    Of course, you might want to do something else than generate these strings, but using the same nlimitedkparts function.