Specifying the coefficients for a 2D polynomial in ApproxFun.jl

Specific Domains Numerics
1

I am trying to understand how ApproxFun.jl constructs 2D functions from their coefficients. In 1D, it all makes sense, for example:

julia> using ApproxFun

julia> f = Fun(Interval(-1, 1), [1,2,3])
Fun(Chebyshev(【-1.0,1.0】),[1.0,2.0,3.0])

julia> f(0.1)
-1.74

julia> 1 + 2*cos(acos(0.1)) + 3*cos(2acos(0.1))
-1.74

So, ApproxFun just uses the coefficients f_k = [1,2,3] in the formula f(x) \approx \sum_k f_k \cos(k \acos(x)) (sorry about the bad formatting). You can see a similar discussion in this part of the docs, but what they don’t seem to discuss is what happens in 2D:

julia> f = Fun(Interval(-1, 1)^2, [1,2,3])
Fun(Chebyshev(【-1.0,1.0】)⊗Chebyshev(【-1.0,1.0】),[1.0,2.0,3.0])

julia> f(0.1, 0.1)
1.50

So how does ApproxFun calculate f(x, y) from x, y and the coefficients [1,2,3]?

2

It goes in order of polynomial degree lexigraphically, that is:

T_0(x)T_0(y) 
T_0(x)T_1(y) 
T_1(x)T_0(y)
T_0(x)T_2(y)
T_1(x)T_1(y)
T_2(x)T_0(y)
…

For example:

julia> f=Fun(Chebyshev()^2,[1.,2.,3.])
Fun(Chebyshev(【-1.0,1.0】)⊗Chebyshev(【-1.0,1.0】),[1.0, 2.0, 3.0])

julia> f(0.1,0.2)
1.7000000000000002

julia> 1 + 2*0.2 + 3*0.1
1.7
1 Like
3

Thanks @dlfivefifty !

Just let me use one more example to make sure that I understand…

julia> using ApproxFun

julia> f = Fun((x, y) -> 1 + 1.2cos(2acos(x)) + 1.3cos(acos(y)), Interval(-1,1)^2)
Fun(Chebyshev(【-1.0,1.0】)⊗Chebyshev(【-1.0,1.0】),[1.0,1.3,3.31642e-17,4.11708e-17,0.0,1.2])

julia> f.coefficients
6-element Array{Float64,1}:
 1.0        
 1.3        
 3.31642e-17
 4.11708e-17
 0.0        
 1.2

The first coefficient is for T_0(x) T_0(y) == cos(0.0acos(x))*cos(0.0acos(y)), the second coefficient is for T_0(x) T_1(y) == cos(0.0acos(x))*cos(1.0acos(y)) and the last coefficient is for T_2(x) T_0(y) == cos(2.0acos(x))*cos(0.0acos(y)).

So a function that I could use to compute the polynomial from the coefficients is:

function compute_poly(c, x, y)
    ret = 0.0
    indeces = [[0,0],[0,1],[1,0],[0,2],[1,1],[2,0]]
    for (i, coeff) in enumerate(c)
        ret += coeff*cos(indeces[i][1]*acos(x))*cos(indeces[i][2]*acos(y))
    end
    return ret
end

I know that ApproxFun does things in a much more clever way, but I am just trying to understand the mathematics.

One thing that I noticed on my adventure trying to implement compute_poly() above is that the ordering of indeces is different from what you get with Combinatorics.jl:

julia> using Combinatorics

julia> a = multiset_combinations([0,0,1,1,2,2], 2) |> collect
  6-element Array{Array{Int64,1},1}:
 [0,0]
 [0,1]
 [0,2]
 [1,1]
 [1,2]
 [2,2]

So apparently people in combinatorics have different conventions from people in numerics?

One last question. Although your explanation makes sense now, it was not excactly intuitive to me… Do you think that this would be worth adding to the documentation?

4

So apparently people in combinatorics have different conventions from people in numerics?

I wouldn’t say this is a standard convention, it’s probably more common to just a matrix. The reason for this convention is that it means every coefficient is hit in a unique location (allowing most of the 1D code in ApproxFun to work also for 2D). The ordering from Combinatorics.jl wouldn’t work since the coefficients would be in different places if you went up to 3.

Do you think that this would be worth adding to the documentation?

Yes of course :wink: Feel free to make a pull request, or even just an Issue.

1 Like
5

Great, thanks for all of your help!

I will try and make a PR asap.