more conversions, with some bug fixes

Author: daanhb

Project.toml

@@ -4,6 +4,7 @@ authors = ["Daan Huybrechs <daan.huybrechs@kuleuven.be>"]
 version = "0.7.2"
 
 [deps]
+BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
 BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
 CompositeTypes = "b152e2b5-7a66-4b01-a709-34e65c35f657"
 DomainIntegrals = "cc6bae93-f070-4015-88fd-838f9505a86c"
@@ -33,10 +34,11 @@ PGFPlotsX = "8314cec4-20b6-5062-9cdb-752b83310925"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 
 [extensions]
-BasisFunctionsRecipesBaseExt = "RecipesBase"
 BasisFunctionsPGFPlotsXExt = "PGFPlotsX"
+BasisFunctionsRecipesBaseExt = "RecipesBase"
 
 [compat]
+BandedMatrices = "1.7.6"
 BlockArrays = "1"
 Calculus = "0.5"
 CompositeTypes = "0.1.3"

src/BasisFunctions.jl

@@ -4,7 +4,7 @@ module BasisFunctions
 
 ## Dependencies
 
-using StaticArrays, BlockArrays, SparseArrays, FillArrays
+using StaticArrays, BandedMatrices, BlockArrays, SparseArrays, FillArrays
 using ToeplitzMatrices, LinearAlgebra, GenericLinearAlgebra
 using FFTW, GenericFFT
 using DomainSets, DomainIntegrals

src/bases/poly/ops/connections.jl

@@ -1,10 +1,14 @@
 
-const JACOBI_OR_LEGENDRE = Union{Legendre,AbstractJacobi}
+isequaldict1(b1::Legendre, b2::AbstractJacobi) =
+    length(b1)==length(b2) && jacobi_α(b2) == 0 && jacobi_β(b2) == 0
+isequaldict2(b1::AbstractJacobi, b2::Legendre) = isequaldict1(b2, b1)
 
-isequaldict(b1::JACOBI_OR_LEGENDRE, b2::JACOBI_OR_LEGENDRE) =
-    length(b1)==length(b2) &&
-    jacobi_α(b1) == jacobi_α(b2) &&
-    jacobi_β(b1) == jacobi_β(b2)
+isequaldict1(b1::ChebyshevU, b2::Ultraspherical) = length(b1)==length(b2) && ultraspherical_λ(b2) == 1
+isequaldict2(b1::Ultraspherical, b2::ChebyshevU) = isequaldict1(b2, b1)
+
+isequaldict1(b1::Ultraspherical{T}, b2::Jacobi) where T =
+    length(b1)==length(b2) && ultraspherical_λ(b1) == one(T)/2 && (jacobi_α(b2) == jacobi_β(b2) == 0)
+isequaldict2(b1::Jacobi, b2::Ultraspherical) = isequaldict1(b2, b1)
 
 const RECURRENCEBASIS = Union{Monomials,OrthogonalPolynomials}
 
@@ -18,7 +22,11 @@ function conversion1(::Type{T}, src::RECURRENCEBASIS, dest::RECURRENCEBASIS; opt
         E = extension_operator(T, dest1, dest)
         E * conversion(T, src, dest1; options...)
     else
-        conversion_using_recurrences(T, src, dest; options...)
+        if isequaldict(src, dest)
+            IdentityOperator(src)
+        else
+            conversion_using_recurrences(T, src, dest; options...)
+        end
     end
 end
 
@@ -82,9 +90,121 @@ function conversion_using_recurrences(::Type{T}, src::Jacobi, dest::Ultraspheric
         diag = T[jacobi_to_ultraspherical(k-1, jacobi_α(src)) for k in 1:length(src)]
         ArrayOperator(Diagonal(diag), src, dest)
     else
-        generic_conversion_using_recurrences(T, src, dest)
+        generic_conversion_using_recurrences(T, src, dest; options...)
+    end
+end
+
+function conversion_using_recurrences(::Type{T}, src::Ultraspherical, dest::Jacobi; options...) where T
+    if (jacobi_α(dest) == jacobi_β(dest)) && (jacobi_α(src) ≈ jacobi_α(dest))
+        inv(conversion_using_recurrences(T, dest, src); options...)
+    else
+        generic_conversion_using_recurrences(T, src, dest; options...)
+    end
+end
+
+
+# Conversion from one ultraspherical basis to the ultraspherical basis with lambda+1
+# Implementation of formula (3.3) of Olver and Townsend, A fast and well-conditioned spectral method,
+# SIAM Review, 2013.
+function banded_ultraspherical_to_ultraspherical(::Type{T}, n, λ) where T
+    A = BandedMatrix{T}(undef, (n,n), (0,2))
+    if n > 0
+        A[1,1] = 1
+    end
+    if n > 1
+        A[2,2] = λ / (λ+1)
+        A[1,2] = 0
+    end
+    for k in 2:n-1
+        K = k+1
+        A[K,K] = λ / (λ+k)
+        A[K-1,K] = 0
+        A[K-2,K] = -λ / (λ+k)
+    end
+    A
+end
+
+# Convert from ChebyshevT expansion to ultraspherical expansion with lambda=1
+function banded_chebyshevt_to_ultraspherical(::Type{T}, n) where T
+    A = BandedMatrix{T}(undef, (n,n), (0,2))
+    if n > 0
+        A[1,1] = 1
+    end
+    onehalf = one(T)/2
+    if n > 1
+        A[1,2] = 0
+        A[2,2] = onehalf
+    end
+    for k in 2:n-1
+        A[k+1,k+1] = onehalf
+        A[k-1,k+1] = -onehalf
+        A[k,k+1] = 0
+    end
+    A
+end
+
+function conversion_using_recurrences(::Type{T}, src::Ultraspherical, dest::Ultraspherical; options...) where T
+    @assert length(src) == length(dest)
+    if ultraspherical_λ(src) ≈ ultraspherical_λ(dest)-1
+        ArrayOperator(banded_ultraspherical_to_ultraspherical(T, length(src), ultraspherical_λ(src)), src, dest)
+    else
+        generic_conversion_using_recurrences(T, src, dest; options...)
+    end
+end
+
+function conversion_using_recurrences(::Type{T}, src::ChebyshevT, dest::Ultraspherical; options...) where T
+    @assert length(src) == length(dest)
+    if ultraspherical_λ(dest) == 1
+        ArrayOperator(banded_chebyshevt_to_ultraspherical(T, length(src)), src, dest)
+    else
+        generic_conversion_using_recurrences(T, src, dest; options...)
+    end
+end
+
+function diagonal_chebyshevt_to_jacobi(::Type{T}, n::Int) where T
+    A = T[chebyshevt_to_jacobi(k, T) for k in 0:n-1]
+    Diagonal(A)
+end
+
+function conversion_using_recurrences(::Type{T}, src::ChebyshevT, dest::Jacobi; options...) where T
+    @assert length(src) == length(dest)
+    onehalf = one(T)/2
+    if jacobi_α(dest) ≈ -onehalf && jacobi_β(dest) ≈ -onehalf
+        ArrayOperator(diagonal_chebyshevt_to_jacobi(T, length(src)), src, dest)
+    else
+        generic_conversion_using_recurrences(T, src, dest; options...)
     end
 end
 
-conversion_using_recurrences(::Type{T}, src::Ultraspherical, dest::Jacobi; options...) where T =
-    inv(conversion_using_recurrences(T, dest, src); options...)
+function conversion_using_recurrences(::Type{T}, src::Jacobi, dest::ChebyshevT; options...) where T
+    onehalf = one(T)/2
+    if jacobi_α(src) ≈ -onehalf && jacobi_β(src) ≈ -onehalf
+        inv(conversion_using_recurrences(T, dest, src; options...))
+    else
+        generic_conversion_using_recurrences(T, src, dest; options...)
+    end
+end
+
+function diagonal_chebyshevu_to_jacobi(::Type{T}, n::Int) where T
+    A = T[chebyshevu_to_jacobi(k, T) for k in 0:n-1]
+    Diagonal(A)
+end
+
+function conversion_using_recurrences(::Type{T}, src::ChebyshevU, dest::Jacobi; options...) where T
+    @assert length(src) == length(dest)
+    onehalf = one(T)/2
+    if jacobi_α(dest) ≈ onehalf && jacobi_β(dest) ≈ onehalf
+        ArrayOperator(diagonal_chebyshevu_to_jacobi(T, length(src)), src, dest)
+    else
+        generic_conversion_using_recurrences(T, src, dest; options...)
+    end
+end
+
+function conversion_using_recurrences(::Type{T}, src::Jacobi, dest::ChebyshevU; options...) where T
+    onehalf = one(T)/2
+    if jacobi_α(src) ≈ onehalf && jacobi_β(src) ≈ onehalf
+        inv(conversion_using_recurrences(T, dest, src; options...))
+    else
+        generic_conversion_using_recurrences(T, src, dest; options...)
+    end
+end

src/bases/poly/ops/dlmf.jl

@@ -87,7 +87,7 @@ function jacobi_eval(n::Int, x, α, β)
     end
     z = z1
     for i in 2:n
-        z = (jacobi_rec_An(i, α, β)*x + jacobi_rec_Bn(i, α, β))*z1 - jacobi_rec_Cn(i, α, β)*z0
+        z = (jacobi_rec_An(i-1, α, β)*x + jacobi_rec_Bn(i-1, α, β))*z1 - jacobi_rec_Cn(i-1, α, β)*z0
         z0 = z1
         z1 = z
     end

src/bases/poly/ops/jacobi.jl

@@ -2,11 +2,6 @@
 "Abstract supertype of Jacobi polynomials."
 abstract type AbstractJacobi{T} <: IntervalOPS{T} end
 
-iscompatible(b1::AbstractJacobi, b2::AbstractJacobi) =
-	jacobi_α(b1) == jacobi_α(b2) && jacobi_β(b1) == jacobi_β(b2)
-isequaldict(b1::AbstractJacobi, b2::AbstractJacobi) =
-	length(b1)==length(b2) && iscompatible(b1,b2)
-
 isorthogonal(b::AbstractJacobi, μ::DomainIntegrals.AbstractJacobiWeight) =
 	jacobi_α(b) == jacobi_α(μ) && jacobi_β(b) == jacobi_β(μ)
 
@@ -62,6 +57,9 @@ jacobi_β(b::Jacobi) = b.β
 
 measure(b::Jacobi) = JacobiWeight(b.α, b.β)
 
+iscompatible(b1::Jacobi, b2::Jacobi) = jacobi_α(b1) == jacobi_α(b2) && jacobi_β(b1) == jacobi_β(b2)
+
+
 rec_An(b::Jacobi, n::Int) = jacobi_rec_An(n, jacobi_α(b), jacobi_β(b))
 rec_Bn(b::Jacobi, n::Int) = jacobi_rec_Bn(n, jacobi_α(b), jacobi_β(b))
 rec_Cn(b::Jacobi, n::Int) = jacobi_rec_Cn(n, jacobi_α(b), jacobi_β(b))
@@ -103,6 +101,10 @@ const Gegenbauer = Ultraspherical
 jacobi_α(b::Ultraspherical{T}) where T = b.λ - one(T)/2
 jacobi_β(b::Ultraspherical{T}) where T = b.λ - one(T)/2
 
+ultraspherical_λ(b::Ultraspherical) = b.λ
+ultraspherical_λ(b::Jacobi{T}) where T =
+	jacobi_α(b) == jacobi_β(b) ? jacobi_α(b)+one(T)/2 : throw(ArgumentError("Jacobi polynomial is not ultraspherical."))
+
 similar(b::Ultraspherical, ::Type{T}, n::Int) where T = Ultraspherical{T}(n, b.λ)
 
 measure(b::Ultraspherical{T}) where T = UltrasphericalWeight(b.λ)
@@ -113,6 +115,9 @@ rec_Cn(b::Ultraspherical, n::Int) = ultraspherical_rec_Cn(n, b.λ)
 
 same_ops_family(b1::Ultraspherical, b2::Ultraspherical) = b1.λ == b2.λ
 
+iscompatible(b1::Ultraspherical, b2::Ultraspherical) = ultraspherical_λ(b1) == ultraspherical_λ(b2)
+
+
 ## Printing
 function show(io::IO, b::Jacobi{Float64})
 	if jacobi_α(b) == 0 && jacobi_β(b) == 0
@@ -132,4 +137,4 @@ end
 
 show(io::IO, b::Ultraspherical{Float64}) = print(io, "Ultraspherical($(length(b)), $(b.λ))")
 show(io::IO, b::Ultraspherical{T}) where T =
-	print(io, "Ultraspherical{$(T)}($(length(b)), $(b.λ))")
+	print(io, "Ultraspherical{$(T)}($(length(b)), $(ultraspherical_λ(b)))")

src/computations/conversion.jl

@@ -143,7 +143,7 @@ projection3(T, src, dest, μ::Measure; options...) =
 default_projection(src::Dictionary, dest::Dictionary, μ = measure(dest); options...) =
     default_projection(operatoreltype(src,dest), src, dest, μ; options...)
 function default_projection(T, src, dest, μ = measure(dest); options...)
-    if src == dest
+    if isequaldict(src, dest)
         IdentityOperator{T}(src, dest)
     else
         G = mixedgram(T, dest, src, μ; options...)
@@ -168,7 +168,7 @@ conversion(::Type{T}, src::D, dest::D; options...) where {T,D} =
     conversion_same_type(T, src, dest; options...)
 
 function conversion_same_type(::Type{T}, src, dest; options...) where {T}
-    if src == dest
+    if isequaldict(src, dest)
         IdentityOperator{T}(src)
     elseif extensible_dictionaries(src, dest)
         extension(T, src, dest; options...)

src/test/Test.jl

@@ -15,7 +15,7 @@ export random_index
 export test_tolerance
 include("test_util_functions.jl")
 
-export test_generic_dict_interface
+export test_generic_dict_interface, test_generic_conversion
 export supports_approximation, supports_interpolation
 export suitable_interpolation_grid, suitable_function
 include("test_dictionary.jl")

src/test/test_dictionary.jl

@@ -456,7 +456,22 @@ function test_generic_dict_derivative(basis)
     # end
 end
 
-
+# Test the conversion from b1 to b2
+function test_generic_conversion(b1, b2)
+    f = random_expansion(b1)
+    x = fixed_point_in_domain(b1)
+    C = conversion(b1, b2)
+    tolerance = sqrt(eps(prectype(b1)))
+    @test src(C) == b1
+    @test dest(C) == b2
+    g = C*f
+    @test g(x) ≈ f(x)
+    if length(b1) == length(b2)
+        Cinv = conversion(b2, b1)
+        f2 = Cinv * (C*f)
+        @test norm(coefficients(f2-f2)) < tolerance
+    end
+end
 
 
 function test_generic_dict_interface(@nospecialize basis)

test/dictionaries/test_dictionaries_functionality.jl

@@ -1,20 +1,4 @@
 
-function test_generic_conversion(b1, b2)
-    f = random_expansion(b1)
-    x = fixed_point_in_domain(b1)
-    C = conversion(b1, b2)
-    tolerance = sqrt(eps(prectype(b1)))
-    @test src(C) == b1
-    @test dest(C) == b2
-    g = C*f
-    @test g(x) ≈ f(x)
-    if length(b1) == length(b2)
-        Cinv = conversion(b2, b1)
-        f2 = Cinv * (C*f)
-        @test norm(coefficients(f2-f2)) < tolerance
-    end
-end
-
 
 function test_dictionary_conversions()
     # OPS are tested separately in test_ops

test/runtests.jl

@@ -9,6 +9,7 @@ using DoubleFloats,
     LinearAlgebra,
     Random,
     StaticArrays,
+    BandedMatrices,
     DomainSets,
     DomainIntegrals