Ultraspherical conversion in Jacobi (#190)

* Ultraspherical conversion in Jacobi

* Specialize Chebyshev cases

* change equality to approx

* Add Legendre conversion tests

Author: jishnub

Project.toml

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

src/Spaces/Jacobi/JacobiOperators.jl

@@ -142,14 +142,22 @@ function Conversion(L::Jacobi,M::Jacobi)
     domain(L) == domain(M) || domain(L) == reverseorientation(domain(M)) ||
         throw(ArgumentError("Domains must be the same"))
 
-    if isapproxinteger(L.a-M.a) && isapproxinteger(L.b-M.b)
-        dm=domain(M)
-        if isapprox(M.a,L.a) && isapprox(M.b,L.b)
-            ConversionWrapper(Operator(I,L))
-        elseif (isapprox(M.b,L.b+1) && isapprox(M.a,L.a)) ||
-            (isapprox(M.b,L.b) && isapprox(M.a,L.a+1))
-            ConcreteConversion(L,M)
-        elseif M.b >= L.b && M.a >= L.a
+    dm=domain(M)
+    dl=domain(L)
+
+    if isapproxinteger(L.a-M.a) && isapproxinteger(L.b-M.b) && M.b >= L.b && M.a >= L.a
+        if isapprox(L.a,M.a) && isapprox(L.b,M.b)
+            return ConversionWrapper(Operator(I,L))
+        elseif (isapprox(L.b+1,M.b) && isapprox(L.a,M.a)) ||
+                (isapprox(L.b,M.b) && isapprox(L.a+1,M.a))
+            return ConcreteConversion(L,M)
+        elseif L.a ≈ L.b ≈ -0.5 && M.a ≈ M.b
+            return Conversion(L,Chebyshev(dl),Ultraspherical(M),M)
+        elseif L.a ≈ L.b && M.a ≈ M.b ≈ -0.5
+            return Conversion(L,Ultraspherical(L),Chebyshev(dm),M)
+        elseif L.a ≈ L.b && M.a ≈ M.b
+            return Conversion(L,Ultraspherical(L),Ultraspherical(M),M)
+        else
             # We split this into steps where a and b are changed by 1:
             # Define the intermediate space J = Jacobi(M.b, L.a, dm)
             # Conversion(L, M) == Conversion(J, M) * Conversion(L, J)
@@ -158,17 +166,19 @@ function Conversion(L::Jacobi,M::Jacobi)
             CLJ = [ConcreteConversion(Jacobi(b-1,L.a,dm), Jacobi(b, L.a, dm)) for b in M.b:-1:L.b+1]
             CJM = [ConcreteConversion(Jacobi(M.b,a-1,dm), Jacobi(M.b, a, dm)) for a in M.a:-1:L.a+1]
             C = [CJM; CLJ]
-            ConversionWrapper(TimesOperator(C))
-        else
-            error("Implement for $L → $M")
+            return ConversionWrapper(TimesOperator(C))
+        end
+    elseif isapproxinteger_addhalf(L.a - M.a) && isapproxinteger_addhalf(L.b - M.b)
+        if L.a ≈ L.b && M.a ≈ M.b ≈ -0.5
+            return Conversion(L,Ultraspherical(L),Chebyshev(dm),M)
+        elseif L.a ≈ L.b ≈ -0.5 && M.a ≈ M.b && M.a >= L.a
+            return Conversion(L,Chebyshev(dl),Ultraspherical(M),M)
+        elseif L.a ≈ L.b && M.a ≈ M.b && M.a >= L.a
+            return Conversion(L,Ultraspherical(L),Ultraspherical(M),M)
         end
-    elseif L.a ≈ L.b ≈ 0 && M.a ≈ M.b ≈ 0.5
-        Conversion(L,Ultraspherical(L),Ultraspherical(M),M)
-    elseif L.a ≈ L.b ≈ 0 && M.a ≈ M.b ≈ -0.5
-        Conversion(L,Ultraspherical(L),Chebyshev(M),M)
-    else # L.a - M.a ≈ L.b - M.b
-        error("Implement for $L → $M")
     end
+
+    error("Implement for $L → $M")
 end
 
 bandwidths(::ConcreteConversion{<:Jacobi,<:Jacobi}) = (0,1)

test/JacobiTest.jl

@@ -48,18 +48,64 @@ using Static
         testtransforms(Jacobi(-0.5,-0.5))
         @test norm(Fun(Fun(exp),Jacobi(-.5,-.5))-Fun(exp,Jacobi(-.5,-.5))) < 300eps()
 
-        @testset for d in [-1..1, 0..1]
-            f = Fun(x->x^2, Chebyshev(d))
-            C = space(f)
-            for J1 in (Jacobi(-0.5, -0.5, d), Legendre(d),
-                            Jacobi(0.5, 0.5, d), Jacobi(2.5, 1.5, d))
-                for J in (J1, NormalizedPolynomialSpace(J1))
-                    g = Fun(f, J)
-                    if !any(isnan, coefficients(g))
-                        @test Conversion(C, J) * f ≈ g
+        @testset for d in (-1..1, 0..1, ChebyshevInterval())
+            @testset "Chebyshev" begin
+                f = Fun(x->x^2, Chebyshev(d))
+                C = space(f)
+                for J1 in (Jacobi(-0.5, -0.5, d), Legendre(d),
+                                Jacobi(0.5, 0.5, d), Jacobi(2.5, 1.5, d))
+                    @testset for J in (J1, NormalizedPolynomialSpace(J1))
+                        g = Fun(f, J)
+                        if !any(isnan, coefficients(g))
+                            @test Conversion(C, J) * f ≈ g
+                        end
                     end
                 end
             end
+            @testset "legendre" begin
+                f = Fun(x->x^2, Legendre(d))
+                C = space(f)
+                for J1 in (Jacobi(-0.5, -0.5, d), Legendre(d),
+                                Jacobi(0.5, 0.5, d), Jacobi(1, 1, d))
+                    @testset for J in (J1, NormalizedPolynomialSpace(J1))
+                        g = Fun(f, J)
+                        if !any(isnan, coefficients(g))
+                            @test Conversion(C, J) * f ≈ g
+                        end
+                    end
+                end
+            end
+            @testset "Ultraspherical(0.5)" begin
+                f = Fun(x->x^2, Ultraspherical(0.5,d))
+                U = space(f)
+                for J1 in (Jacobi(-0.5, -0.5, d), Legendre(d),
+                                Jacobi(0.5, 0.5, d))
+                    @testset for J in (J1, NormalizedPolynomialSpace(J1))
+                        g = Fun(f, J)
+                        if !any(isnan, coefficients(g))
+                            @test Conversion(U, J) * f ≈ g
+                        end
+                    end
+                end
+            end
+            @testset "Ultraspherical(1)" begin
+                f = Fun(x->x^2, Ultraspherical(1,d))
+                U = space(f)
+                for J1 in (Legendre(d), Jacobi(0.5, 0.5, d), Jacobi(1.5, 0.5, d),
+                            Jacobi(0.5, 1.5, d))
+                    @testset for J in (J1,)
+                        g = Fun(f, J)
+                        if !any(isnan, coefficients(g))
+                            @test Conversion(U, J) * f ≈ g
+                        end
+                    end
+                end
+            end
+        end
+
+        @testset for (S1, S2) in ((Legendre(), Jacobi(-0.5, -0.5)), (Jacobi(-0.5, -0.5), Legendre()),
+                (Legendre(), Jacobi(1.5, 1.5)))
+            @test Conversion(S1, S2) * Fun(x->x^4, S1) ≈ Fun(x->x^4, S2)
         end
 
         @testset "inference tests" begin