Normalized jacobi conversion for alpha=beta=-0.5 (#227)

* Normalized jacobi conversion for alpha=beta=-0.5

* bump version to v0.6.22

Author: jishnub

Project.toml

@@ -1,6 +1,6 @@
 name = "ApproxFunOrthogonalPolynomials"
 uuid = "b70543e2-c0d9-56b8-a290-0d4d6d4de211"
-version = "0.6.21"
+version = "0.6.22"
 
 [deps]
 ApproxFunBase = "fbd15aa5-315a-5a7d-a8a4-24992e37be05"

src/Spaces/Jacobi/Jacobi.jl

@@ -22,7 +22,19 @@ const NormalizedJacobi{D<:Domain,R,T} = NormalizedPolynomialSpace{Jacobi{D,R,T},
 NormalizedJacobi(s...) = NormalizedPolynomialSpace(Jacobi(s...))
 NormalizedLegendre(d...) = NormalizedPolynomialSpace(Legendre(d...))
 
-normalization(::Type{T}, sp::Jacobi, k::Int) where T = FastTransforms.Anαβ(k, sp.a, sp.b)
+function normalization(::Type{T}, sp::Jacobi, k::Int) where T
+    if sp.a == sp.b == -0.5 && k == 0
+        # In this case, the expression for Anαβ has a division by zero, so we evaluate this using Mathematica
+        # In principle this may be generalized to arbitrary α + β = -1
+        # The exact expression from Mathematica in terms of the hypergeometric 2F1 function
+        # \fbox{$\frac{\, _2F_1(1,-\alpha ;\beta +2;-1)}{\beta +1}+\frac{\, _2F_1(1,-\beta ;\alpha +2;-1)}{\alpha +1}\text{ if }\alpha >-1\land \beta >-1$}
+        # or, by eliminating β and expressed in terms of the polygamma function,
+        # \fbox{$\frac{1}{2} \left(\psi ^{(0)}\left(\frac{1}{2}-\frac{\alpha }{2}\right)-\psi ^{(0)}\left(-\frac{\alpha }{2}\right)\right)+\frac{1}{2} \left(\psi ^{(0)}\left(\frac{\alpha }{2}+1\right)-\psi ^{(0)}\left(\frac{\alpha }{2}+\frac{1}{2}\right)\right)\text{ if }-1<\alpha <0$}
+        T(pi)
+    else
+        FastTransforms.Anαβ(T(k), T(sp.a), T(sp.b))
+    end
+end
 
 function Ultraspherical(J::Jacobi)
     if J.a == J.b

test/JacobiTest.jl

@@ -55,7 +55,7 @@ using Static
 
     @testset "Conversion" begin
         testtransforms(Jacobi(-0.5,-0.5))
-        @test norm(Fun(Fun(exp),Jacobi(-.5,-.5))-Fun(exp,Jacobi(-.5,-.5))) < 300eps()
+        @test norm(Fun(Fun(exp),Jacobi(-0.5,-0.5))-Fun(exp,Jacobi(-0.5,-0.5))) < 300eps()
 
         @testset for d in (-1..1, 0..1, ChebyshevInterval())
             @testset "Chebyshev" begin
@@ -70,6 +70,11 @@ using Static
                         end
                     end
                 end
+                # Jacobi(-0.5, -0.5) to NormalizedJacobi(-0.5, -0.5) can have a NaN in the (1,1) index
+                # if the normalization is not carefully implemented,
+                # so this test checks that this is not the case
+                g = Conversion(Jacobi(-0.5,-0.5), NormalizedJacobi(-0.5,-0.5)) * Fun(x->x^3, Jacobi(-0.5, -0.5))
+                @test coefficients(g) ≈ coefficients(Fun(x->x^3, NormalizedChebyshev()))
             end
             @testset "legendre" begin
                 f = Fun(x->x^2, Legendre(d))