JuliaApproximation / SemiclassicalOrthogonalPolynomials.jl

A Julia repository for semiclassical orthogonal polynomials
MIT License
7 stars 3 forks source link

Associated 2-band OPs Golub Welsch? #29

Closed dlfivefifty closed 3 months ago

dlfivefifty commented 3 years ago

The two-band analogue of Chebyshev T is now implemented as T = TwoBandJacobi(ρ, -1/2, -1/2, 1/2), orthogonal w.r.t. |x|/(sqrt(1-x^2)sqrt(x^2-ρ^2)). It has the nice property that the Hilbert transform of its weight is 0, so the Hilbert transform is a simple map from T to Q = Associated(T), the associated orthogonal polynomials.

The surprising thing is Q has an orthogonality measure with a delta function at 0, I believe its something like sqrt(1-x^2)sqrt(x^2-ρ^2)/|x| + ρ δ(x) (where the first term is only valid for ρ < |x| <1). This has the annoying feature that Golub--Welsch gives extra points associated with the δ:

julia> ρ = 0.5; T = TwoBandJacobi(ρ, -1/2, -1/2, 1/2);

julia> Q = Associated(T);

julia> golubwelsch(Q[:,Base.OneTo(7)])[1] # grid point at 0 which I don't want
7-element Vector{Float64}:
 -0.9434855817366542
 -0.7905694150420907
 -0.5998624484455107
  7.771561172376096e-16
  0.5998624484455127
  0.790569415042095
  0.9434855817366555

Anyone seen anything like this before and know if there's a modification to Gauss quadrature, i.e. an anti-Gauss–Lobatto that prevents certain points?

I just need this for the expansion in Q. I can work around it by using a Vandermonde matrix. But it's an interesting problem.

@MikaelSlevinsky @TSGut @ajt60gaibb @marcusdavidwebb

TSGut commented 3 years ago

There are "anti-Gaussian" quadrature methods which as far as I know interlace the points of Gaussian quadrature. Not sure if that is why you said "anti-Gauss–Lobatto" or if you just meant "the opposite of what Gauss–Lobatto does".

dlfivefifty commented 3 years ago

I meant "opposite of what Gauss-Lobatto" (or Radau): instead of insisting on a specific point, I want to forbid it.

Note its a bit delicate as for even points we get two points tending to (but not exactly) zero

MikaelSlevinsky commented 3 years ago

Is there a characterization of the non-negative weights that live in the kernel of the Hilbert transform?

In the even case, does the pair of Gauss points breach |x| = ρ?

dlfivefifty commented 3 years ago

Is there a characterization of the non-negative weights that live in the kernel of the Hilbert transform?

I think the kernel is only |x|/(sqrt(1-x^2)sqrt(ρ^2-x^2)).

In the even case, does the pair of Gauss points breach |x| = ρ?

Not even close:

julia> golubwelsch(Q[:,Base.OneTo(2)])[1]
2-element Vector{Float64}:
 -0.3354101966249685
  0.3354101966249685

julia> golubwelsch(Q[:,Base.OneTo(4)])[1]
4-element Vector{Float64}:
 -0.8260064989841912
 -0.10635236062957643
  0.10635236062957698
  0.8260064989841945

julia> golubwelsch(Q[:,Base.OneTo(6)])[1]
6-element Vector{Float64}:
 -0.9146586174230956
 -0.6894624123565193
 -0.03501609669916683
  0.0350160966991655
  0.6894624123565204
  0.9146586174230971
MikaelSlevinsky commented 3 years ago

That's interesting: I can see why the odd symmetry forces a root at the non-negligible Dirac point mass, but had no idea that it would be "strong" enough to pull other points into zero-weight territory (in the even case).

dlfivefifty commented 3 years ago

It sort of makes sense; since one typically expects finite section to converge

TSGut commented 3 months ago

Since two bands were moved out of the package I guess it's fine to close this for now.