Question on Derivative and Conversion on Chebyshev^2

[michakraus]

I am somewhat puzzled about what exactly is happening when computing derivatives and conversions on 2d Chebyshev product spaces. Some enlightenment would be greatly appreciated. This is also related to #725.

My question circles around the fact that Chebyshev(0..1)^2 is defined on Padua points, whereas derivatives and projections thereof, living in Ultraspherical(1, 0..1) ⊗ Chebyshev(0..1), Chebyshev(0..1) ⊗ Ultraspherical(1, 0..1) or Ultraspherical(1, 0..1)^2, are defined on tensor product meshes of Chebyshev points. I am wondering how exactly this works and how this is dealt with internally. I found that there are special methods for points(), transform() and itransform(), dispatching on TensorSpace{<:Tuple{<:Chebyshev{<:ChebyshevInterval},<:Chebyshev{<:ChebyshevInterval}}}, but I could not find similar methods for Derivative() or Conversion(), taking care of the special nature of the Padua points.

For context: I have a problem (namely computing Poincaré integral invariants), where

• I construct a Chebyshev product space,

C2 = Chebyshev(0..1)^2

• extract the Padua points and evaluate some scalar function on those points,

p = points(C2)
v = f.(p)

• transform those to coefficients and create a fun to work with, e.g., to compute derivatives,

plan = plan_transform(C2, v)
c = plan * v
f = Fun(C2, c)
# ...

In the end I need to extract all values, derivatives, etc., on a common grid, which is why everything that does not already live there is projected onto Ultraspherical(1, 0..1)^2.

I implemented this using both, ApproxFun and OrthogonalPolynomialsQuasi, where the former uses Padua points on the Chebyshev()^2 space and the latter uses a Chebyshev tensor grid. Not surprisingly, the ApproxFun implementation gives better results for a comparable number of points. However, even when using the full grid corresponding to the Padua grid, I cannot get comparable results in the OrthogonalPolynomialsQuasi implementation (accuracy is much worse). Trying to understand the reason for this and debugging my code lead to the above questions.

This also raises the question if using Padua points is possible together with OrthogonalPolynomialsQuasi (although this might not be the right place to ask this question)?

[dlfivefifty]

I’m actually surprised OrthogonalPolynomialsQuasi.jl works at all: I don’t remember implementing kronecker products with it.

Here you are touching some tricky points that haven’t really been thought through and the code is really just a thrown together prototype. I had a 4th year student 4 years ago look at Chebyshev U (sine series) and Fourier variants of Padua points and he worked out the details so it’s certainly possible to add something smarter, but I don’t have time to do this.

Happy to send the thesis and sketch out a software development plan if this is something you have time to contribute to. I think designing with 3D in mind would be a good idea.

[michakraus]

I’m actually surprised OrthogonalPolynomialsQuasi.jl works at all: I don’t remember implementing kronecker products with it.

Yes, in fact I patched this together in a very hacky way myself...

Happy to send the thesis and sketch out a software development plan if this is something you have time to contribute to. I think designing with 3D in mind would be a good idea.

Time is a rare resource for me as well, but I can see what I can do. I'd be happy to have a look at the thesis and try to get an idea how realistic it is for me to implement this.

[dlfivefifty]

There’s Julia code lying around for the transforms themselves, so you could probably roll something directly

[michakraus]

Sounds good. Is this available somewhere or can you send it to me (together with the aforementioned thesis)? Thanks!