JuliaApproximation / AlgebraicCurveOrthogonalPolynomials.jl

MIT License
1 stars 3 forks source link

2D via Lanczos to Tensor product basis? #4

Open dlfivefifty opened 3 years ago

dlfivefifty commented 3 years ago

@MikaelSlevinsky's transform work by first expanding in a tensor product basis, which contains non-polynomial terms, before transforming to a polynomial basis. For example, to compute OPs on the triangle we expand in the basis

P_k(x) P_j(y/(1-x))

(I'm using OPs on 0..1 here) which includes non-polynomial terms y/(1-x)^j, before converting to the polynomial basis

P_{n-k}^(2k+1,0)(x)(1-x)^kP_k(y/(1-x))

The point is polynomials in x and y are a subspace of the non-polynomial basis.

Does this technique translate to OPs on the disk (a la Dunkl & Xu, not Zernike)

C_{n-k}^(k)(x)*ρ(x)^k*T_k(y/ρ(x))

where ρ(x) = sqrt(1-x^2)? A tensor basis like

T_k(x) T_j(y/ρ(x))

does not seem to work here... seems like we'd need a sum-space frame that also includes ρ(x)*T_k(x)*T_j(y/ρ(x))....

What about ρ(x) = (1-x^4)^(1/4)? Now we wouldn't know the OP basis but perhaps we can represent it by the conversion to a bigger basis a la LanczosPolynomial.

dlfivefifty commented 3 years ago

Hmm, for the disk what about something like

P_k(ρ(x)) * T_j(y/ρ(x)), x*P_k(ρ(x)) * T_j(y/ρ(x))

whee P_k are orthogonal on 0..1. I think this contains all polynomials a subspace, and we can use even/odd in x for fast transforms?

For ρ(x) = (1-x^4)^(1/4) we would need two more:

x^j * P_k(ρ(x)) * T_j(y/ρ(x)) for j = 0:3

so the expansion part is not so clear...

MikaelSlevinsky commented 3 years ago

Yeah I think it works on the disk since it's quite a lot like spherical harmonics without the z = cos(theta) variable transformation.

Don't know about the quartic weight, but perhaps instead of Givens rotations the isométries can be factored in terms of short Householder reflectors.

dlfivefifty commented 3 years ago

Works on the disk from what basis? I’m not talking about the radial coordinates basis

MikaelSlevinsky commented 3 years ago

I think the Dunkl-Xu will work. Just needs implementing

MikaelSlevinsky commented 3 years ago

It works like this. I assume we have no problem in the even k case. You write:

C_{n-k}^(k)(x)*ρ(x)^k*T_k(y/ρ(x))

as polynomials in x times T_k(y/ρ(x)) by unrolling C_{n-k}^(k)(x)*ρ(x)^k via Givens rotations.

For odd k, you leave one power of ρ(x) behind (and use the fact that T_k is odd). Discrete sine and cosine transforms follow.