Usability improvements for orthogonal polynomials

Author: daanhb

Project.toml

@@ -1,7 +1,7 @@
 name = "BasisFunctions"
 uuid = "4343a256-5453-507d-8aad-01a9d7189916"
 authors = ["Daan Huybrechs <daan.huybrechs@kuleuven.be>"]
-version = "0.6.5"
+version = "0.6.6"
 
 [deps]
 BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
@@ -38,25 +38,30 @@ BasisFunctionsPGFPlotsXExt = "PGFPlotsX"
 
 [compat]
 BlockArrays = "0.16"
+Calculus = "0.5"
 CompositeTypes = "0.1.3"
 DomainIntegrals = "0.4.5"
 DomainSets = "0.7"
+DoubleFloats = "1"
 FFTW = "1.6"
 FastGaussQuadrature = "0.4,0.5"
 FillArrays = "0.13,1"
 GaussQuadrature = "0.5.7"
 GenericFFT = "0.1.3"
 GenericLinearAlgebra = "0.3"
-GridArrays = "0.2.2"
+GridArrays = "0.2.3"
 IterativeSolvers = "0.9"
+LinearAlgebra = "1"
 MacroTools = "0.5"
 OrderedCollections = "1.2"
 PGFPlotsX = "1.4"
 QuadGK = "2.4"
+Random = "1"
 RecipesBase = "1.0"
 Reexport = "1.0"
 SpecialFunctions = "1.0, 2.0"
 StaticArrays = "1.4"
+Test = "1"
 ToeplitzMatrices = "0.7"
 julia = "1.9"
 

src/BasisFunctions.jl

@@ -300,7 +300,8 @@ export Legendre, Jacobi, Laguerre, Hermite,
     recurrence_eval, recurrence_eval_derivative, monic_recurrence_eval,
     monic_recurrence_coefficients,
     symmetric_jacobi_matrix, roots, gauss_rule, sorted_gauss_rule, first_moment,
-    leading_order_coefficient
+    leading_order_coefficient,
+    ops_roots
 
 # from specialOPS.jl
 export HalfRangeChebyshevIkind, HalfRangeChebyshevIIkind, WaveOPS,

src/bases/dictionary/basisfunction.jl

@@ -23,6 +23,25 @@ promote_rule(::Type{<:TypedFunction{S1,T1}}, ::Type{<:TypedFunction{S2,T2}}) whe
 "The supertype of functions that can be associated with a dictionary or a family of basis functions."
 abstract type AbstractBasisFunction{S,T} <: TypedFunction{S,T} end
 
+function expansion(φ::AbstractBasisFunction)
+    coef = zeros(dictionary(φ))
+    coef[index(φ)] = 1
+    expansion(dictionary(φ), coef)
+end
+
+roots(φ::AbstractBasisFunction) = roots(expansion(φ))
+
+function Base.:^(φ::AbstractBasisFunction, i::Int)
+    @assert i >= 0
+    if i == 1
+        φ
+    elseif i == 2
+        φ * φ
+    else
+        φ^(i-1) * φ
+    end
+end
+
 "A `BasisFunction` is one element of a dictionary."
 struct BasisFunction{S,T,D<:Dictionary{S,T},I} <: AbstractBasisFunction{S,T}
     dictionary  ::  D

src/bases/dictionary/expansions.jl

@@ -119,6 +119,17 @@ ei(dim,i, coefficients) = tuple((coefficients*Matrix{Int}(I, dim, dim)[:,i])...)
 differentiate(f::Expansion, var, order) = differentiate(f, ei(dimension(f), var, order))
 antidifferentiate(f::Expansion, var, order) = antidifferentiate(f, ei(dimension(f), var, order))
 
+function Base.:^(f::Expansion, i::Int)
+    @assert i >= 0
+    if i == 1
+        f
+    elseif i == 2
+        f * f
+    else
+        f^(i-1) * f
+    end
+end
+
 # To be implemented: Laplacian (needs multiplying functions)
 ## Δ(f::Expansion)
 
@@ -127,7 +138,9 @@ adjoint(f::Expansion) = differentiate(f)
 
 ∫(f::Expansion) = antidifferentiate(f)
 
-roots(f::Expansion) = roots(dictionary(f), coefficients(f))
+roots(f::Expansion) = expansion_roots(dictionary(f), coefficients(f))
+
+Base.@deprecate roots(dict::Dictionary, coef::AbstractVector) expansion_roots(dict, coef)
 
 show(io::IO, mime::MIME"text/plain", fun::Expansion) = composite_show(io, mime, fun)
 Display.displaystencil(fun::Expansion) = ["Expansion(", dictionary(fun), ", ", coefficients(fun), ")"]
@@ -146,7 +159,10 @@ iscompatible(d1::Dictionary, d2::Dictionary) = false
 
 iscompatible(e1::Expansion, e2::Expansion) = iscompatible(dictionary(e1),dictionary(e2))
 
-function compatible_dictionary(d1::Dictionary, d2::Dictionary)
+compatible_dictionary(d1::Dictionary, d2::Dictionary) =
+    default_compatible_dictionary(d1, d2)
+
+function default_compatible_dictionary(d1::Dictionary, d2::Dictionary)
     if d1 == d2
         d1
     else
@@ -184,7 +200,6 @@ function to_compatible_expansions(e1::Expansion, e2::Expansion)
 end
 
 function (+)(f::Expansion, g::Expansion)
-    @assert iscompatible(f, g)
     dict, coef = expansion_sum(dictionary(f),dictionary(g),coefficients(f),coefficients(g))
     Expansion(dict, coef)
 end
@@ -208,7 +223,7 @@ function expansion_sum2(dict1, dict2, coef1, coef2)
     if iscompatible(dict1, dict2)
         default_expansion_sum(dict1, dict2, coef1, coef2)
     else
-        error("Multiplication of expansions in these dictionaries is not implemented.")
+        error("Summation of expansions in these dictionaries is not implemented.")
     end
 end
 
@@ -265,3 +280,23 @@ iterate(e::Expansion, state) = iterate(coefficients(e), state)
 Base.collect(e::Expansion) = coefficients(e)
 
 Base.BroadcastStyle(e::Expansion) = Base.Broadcast.DefaultArrayStyle{dimension(e)}()
+
+# Interoperate with basis functions.
+# TODO: implement cleaner using promotion mechanism
+
+Base.:*(φ1::AbstractBasisFunction, φ2::AbstractBasisFunction) = expansion(φ1) * expansion(φ2)
+Base.:*(f::Expansion, φ2::AbstractBasisFunction) = f * expansion(φ2)
+Base.:*(φ1::AbstractBasisFunction, f::Expansion) = expansion(φ1) * f
+
+Base.:*(A::Number, φ::AbstractBasisFunction) = A * expansion(φ)
+Base.:*(φ::AbstractBasisFunction, A::Number) = A * expansion(φ)
+
+Base.:+(φ1::AbstractBasisFunction, φ2::AbstractBasisFunction) = expansion(φ1) + expansion(φ2)
+Base.:+(f::Expansion, φ2::AbstractBasisFunction) = f + expansion(φ2)
+Base.:+(φ1::AbstractBasisFunction, f::Expansion) = expansion(φ1) + f
+
+Base.:-(φ1::AbstractBasisFunction, φ2::AbstractBasisFunction) = expansion(φ1) - expansion(φ2)
+Base.:-(f::Expansion, φ2::AbstractBasisFunction) = f - expansion(φ2)
+Base.:-(φ1::AbstractBasisFunction, f::Expansion) = expansion(φ1) - f
+
+

src/bases/modifiers/mapped_dict.jl

@@ -362,8 +362,8 @@ function expansion_multiply(dict1::MappedDict, dict2::MappedDict, coef_src1, coe
     mapped_dict(mset, forward_map(dict1)), mcoef
 end
 
-function roots(dict::MappedDict, coef)
-    r = roots(superdict(dict), coef)
+function expansion_roots(dict::MappedDict, coef)
+    r = expansion_roots(superdict(dict), coef)
     forward_map(dict).(r)
 end
 

src/bases/modifiers/operated_dict.jl

@@ -308,7 +308,7 @@ function (*)(op::DictionaryOperator, dict::OperatedDict)
     OperatedDict(operator(dict) ∘ op)
 end
 
-roots(dict::OperatedDict, coef) = roots(superdict(dict), operator(dict)*coef)
+expansion_roots(dict::OperatedDict, coef) = expansion_roots(superdict(dict), operator(dict)*coef)
 
 expansion_multiply1(dict1::OperatedDict, dict2, coef1, coef2) =
 	expansion_multiply(superdict(dict1), dict2, operator(dict1)*coef1, coef2)

src/bases/poly/chebyshev/ChebyshevT.jl

@@ -2,7 +2,7 @@
 """
 A basis of Chebyshev polynomials of the first kind on the interval `[-1,1]`.
 """
-struct ChebyshevT{T} <: OPS{T}
+struct ChebyshevT{T} <: IntervalOPS{T}
     n			::	Int
 end
 
@@ -11,8 +11,17 @@ show(io::IO, d::ChebyshevT{T}) where T = print(io, "ChebyshevT{$(T)}($(length(d)
 
 ChebyshevT(n::Int) = ChebyshevT{Float64}(n)
 
+ChebyshevT(d::PolynomialDegree) = ChebyshevT(value(d)+1)
+ChebyshevT{T}(d::PolynomialDegree) where T = ChebyshevT{T}(value(d)+1)
+
 similar(b::ChebyshevT, ::Type{T}, n::Int) where {T} = ChebyshevT{T}(n)
 
+to_chebyshev_dict(dict::IntervalOPS{T}) where T = ChebyshevT{T}(length(dict))
+to_chebyshev(f) = to_chebyshev(expansion(f))
+to_chebyshev(f::Expansion) = to_chebyshev(dictionary(f), coefficients(f))
+to_chebyshev(dict::Dictionary, coef) =
+	conversion(dict, to_chebyshev_dict(dict)) * expansion(dict, coef)
+
 convert(::Type{ChebyshevT{T}}, d::ChebyshevT) where {T} = ChebyshevT{T}(d.n)
 
 "Mapped Chebyshev polynomials."
@@ -63,8 +72,6 @@ rec_An(b::ChebyshevT{T}, n::Int) where {T} = n==0 ? one(T) : 2one(T)
 rec_Bn(b::ChebyshevT{T}, n::Int) where {T} = zero(T)
 rec_Cn(b::ChebyshevT{T}, n::Int) where {T} = one(T)
 
-support(b::ChebyshevT{T}) where {T} = ChebyshevInterval{T}()
-
 # We can define this O(1) evaluation method, but only for points that are
 # real and lie in [-1,1]
 # Note that if x is not Real, recurrence_eval will be called by the OPS supertype
@@ -178,7 +185,8 @@ end
 transform_from_grid(T, src::GridBasis, dest::ChebyshevT, grid; options...) =
 	inv(transform_to_grid(T, dest, src, grid; options...))
 
-function roots(dict::ChebyshevT, coefficients::AbstractVector)
+function expansion_roots(dict::ChebyshevT, coefficients::AbstractVector)
+	@assert length(dict) == length(coefficients)
 	T = eltype(coefficients)
 	n = length(dict)-1
 	# construct the colleague matrix (according to ATAP)

src/bases/poly/chebyshev/ChebyshevU.jl

@@ -1,11 +1,14 @@
 
 "A basis of Chebyshev polynomials of the second kind on the interval `[-1,1]`."
-struct ChebyshevU{T} <: OPS{T}
+struct ChebyshevU{T} <: IntervalOPS{T}
     n			::	Int
 end
 
 ChebyshevU(n::Int) = ChebyshevU{Float64}(n)
 
+ChebyshevU(d::PolynomialDegree) = ChebyshevU(value(d)+1)
+ChebyshevU{T}(d::PolynomialDegree) where T = ChebyshevU{T}(value(d)+1)
+
 similar(b::ChebyshevU, ::Type{T}, n::Int) where {T} = ChebyshevU{T}(n)
 
 show(io::IO, d::ChebyshevU{Float64}) = print(io, "ChebyshevU($(length(d)))")
@@ -49,8 +52,6 @@ rec_Bn(b::ChebyshevU{T}, n::Int) where {T} = zero(T)
 
 rec_Cn(b::ChebyshevU{T}, n::Int) where {T} = one(T)
 
-support(b::ChebyshevU{T}) where {T} = ChebyshevInterval{T}()
-
 
 "A Chebyshev polynomial of the second kind"
 struct ChebyshevUPolynomial{T} <: OrthogonalPolynomial{T}

src/bases/poly/monomials.jl

@@ -12,6 +12,11 @@ Monomials(n) = Monomials{Float64}(n)
 show(io::IO, b::Monomials{Float64}) = print(io, "Monomials($(length(b)))")
 show(io::IO, b::Monomials{T}) where T = print(io, "Monomials{$(T)}($(length(b)))")
 
+to_monomials_dict(dict::PolynomialBasis{T}) where T = Monomials{T}(length(dict))
+to_monomials(f) = to_monomials(expansion(f))
+to_monomials(f::Expansion) = to_monomials(dictionary(f), coefficients(f))
+to_monomials(dict::Dictionary, coef) = conversion(dict, to_monomials_dict(dict)) * expansion(dict, coef)
+
 support(dict::Monomials{T}) where {T} = DomainSets.FullSpace{T}()
 
 size(dict::Monomials) = (dict.n,)

src/bases/poly/ops/hermite.jl

@@ -10,6 +10,9 @@ end
 
 Hermite(n::Int) = Hermite{Float64}(n)
 
+Hermite(d::PolynomialDegree) = Hermite(value(d)+1)
+Hermite{T}(d::PolynomialDegree) where T = Hermite{T}(value(d)+1)
+
 similar(b::Hermite, ::Type{T}, n::Int) where {T} = Hermite{T}(n)
 
 support(b::Hermite{T}) where {T} = DomainSets.FullSpace(T)