Support Neumann conditions in p-FEM (#75)

* Support Neumann conditions in p-FEM

* Update test_odes.jl

* natural (Neumann) boundary conditions work!

* cos.(θ) .* F

* improve OP ldiv special cases

* Update test_odes.jl

* one-sided Neumann test

* LazyArrays v0.22

* Update test_odes.jl

* Update Project.toml

* improve cov

* Update test_legendre.jl

Author: dlfivefifty

Project.toml

@@ -1,7 +1,7 @@
 name = "ClassicalOrthogonalPolynomials"
 uuid = "b30e2e7b-c4ee-47da-9d5f-2c5c27239acd"
 authors = ["Sheehan Olver <[email protected]>"]
-version = "0.4.9"
+version = "0.4.10"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -22,6 +22,7 @@ LazyArrays = "5078a376-72f3-5289-bfd5-ec5146d43c02"
 LazyBandedMatrices = "d7e5e226-e90b-4449-9968-0f923699bf6f"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 QuasiArrays = "c4ea9172-b204-11e9-377d-29865faadc5c"
+SemiseparableMatrices = "f8ebbe35-cbfb-4060-bf7f-b10e4670cf57"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 
 [compat]
@@ -34,14 +35,14 @@ DomainSets = "0.5.6"
 FFTW = "1.1"
 FastGaussQuadrature = "0.4.3"
 FastTransforms = "0.12"
-FillArrays = "0.11, 0.12"
+FillArrays = "0.12"
 HypergeometricFunctions = "0.3.4"
 InfiniteArrays = "0.12"
 InfiniteLinearAlgebra = "0.6"
 IntervalSets = "0.5"
-LazyArrays = "0.21.20"
+LazyArrays = "0.22"
 LazyBandedMatrices = "0.7"
-QuasiArrays = "0.7, 0.8"
+QuasiArrays = "0.8"
 SpecialFunctions = "1.0"
 julia = "1.6"
 

src/ClassicalOrthogonalPolynomials.jl

@@ -28,7 +28,8 @@ import QuasiArrays: cardinality, checkindex, QuasiAdjoint, QuasiTranspose, Inclu
                     QuasiDiagonal, MulQuasiArray, MulQuasiMatrix, MulQuasiVector, QuasiMatMulMat,
                     ApplyQuasiArray, ApplyQuasiMatrix, LazyQuasiArrayApplyStyle, AbstractQuasiArrayApplyStyle,
                     LazyQuasiArray, LazyQuasiVector, LazyQuasiMatrix, LazyLayout, LazyQuasiArrayStyle,
-                    _getindex, layout_getindex, _factorize, AbstractQuasiArray, AbstractQuasiMatrix, AbstractQuasiVector
+                    _getindex, layout_getindex, _factorize, AbstractQuasiArray, AbstractQuasiMatrix, AbstractQuasiVector,
+                    AbstractQuasiFill, _dot
 
 import InfiniteArrays: OneToInf, InfAxes, Infinity, AbstractInfUnitRange, InfiniteCardinal, InfRanges
 import InfiniteLinearAlgebra: chop!, chop
@@ -351,6 +352,21 @@ function \(wA::WeightedOrthogonalPolynomial, wB::WeightedOrthogonalPolynomial)
     A \ ((A * (wA \ w_B)) .* B)
 end
 
+## special expansion routines for constants and x
+function _op_ldiv(P::AbstractQuasiMatrix{V}, f::AbstractQuasiFill{T,1}) where {T,V}
+    TV = promote_type(T,V)
+    Vcat(getindex_value(f)/_p0(P),Zeros{TV}(∞))
+end
+
+_op_ldiv(::AbstractQuasiMatrix{V}, ::QuasiZeros{T,1}) where {T,V} = Zeros{promote_type(T,V)}(∞)
+
+function _op_ldiv(P::AbstractQuasiMatrix{V}, f::Inclusion{T}) where {T,V}
+    axes(P,1) == f || throw(DimensionMismatch())
+    A,B,_ = recurrencecoefficients(P)
+    TV = promote_type(T,V)
+    c = inv(convert(TV,A[1]*_p0(P)))
+    Vcat(-B[1]c, c, Zeros{TV}(∞))
+end
 
 include("classical/hermite.jl")
 include("classical/jacobi.jl")

src/classical/fourier.jl

@@ -79,3 +79,8 @@ function broadcasted(::LazyQuasiArrayStyle{2}, ::typeof(*), s::BroadcastQuasiVec
                         Vcat([zeros(T,1,1)], Fill(Matrix(zero(T)*I,2,2),∞)),
                         Vcat([[one(T)/2 0]], Fill([0 -one(T)/2; one(T)/2 0],∞))))
 end
+
+
+# support cos.(θ) .* F
+Base.broadcasted(::LazyQuasiArrayStyle{2}, ::typeof(*), c::Broadcasted{<:Any,<:Any,typeof(cos),<:Tuple{<:Inclusion{<:Any,RealNumbers}}}, F::Fourier) = materialize(c) .* F
+Base.broadcasted(::LazyQuasiArrayStyle{2}, ::typeof(*), c::Broadcasted{<:Any,<:Any,typeof(sin),<:Tuple{<:Inclusion{<:Any,RealNumbers}}}, F::Fourier) = materialize(c) .* F

src/classical/jacobi.jl

@@ -33,49 +33,16 @@ summary(io::IO, w::JacobiWeight) = print(io, "(1-x)^$(w.a) * (1+x)^$(w.b) on -1.
 
 sum(P::JacobiWeight) = jacobimoment(P.a, P.b)
 
-struct LegendreWeight{T} <: AbstractJacobiWeight{T} end
-LegendreWeight() = LegendreWeight{Float64}()
-legendreweight(d::AbstractInterval{T}) where T = LegendreWeight{float(T)}()[affine(d,ChebyshevInterval{T}())]
-
-function getindex(w::LegendreWeight{T}, x::Number) where T
-    x ∈ axes(w,1) || throw(BoundsError())
-    one(T)
-end
-
-getproperty(w::LegendreWeight{T}, ::Symbol) where T = zero(T)
-
-sum(::LegendreWeight{T}) where T = 2one(T)
-
-_weighted(::LegendreWeight, P) = P
-_weighted(::SubQuasiArray{<:Any,1,<:LegendreWeight}, P) = P
-
-broadcasted(::LazyQuasiArrayStyle{1}, ::typeof(*), ::LegendreWeight{T}, ::LegendreWeight{V}) where {T,V} =
-    LegendreWeight{promote_type(T,V)}()
-
-broadcasted(::LazyQuasiArrayStyle{1}, ::typeof(sqrt), w::LegendreWeight{T}) where T = w
-
-broadcasted(::LazyQuasiArrayStyle{1}, ::typeof(Base.literal_pow), ::Base.RefValue{typeof(^)}, w::LegendreWeight, ::Base.RefValue{Val{k}}) where k = w
 
 # support auto-basis determination
 
 singularities(a::AbstractAffineQuasiVector) = singularities(a.x)
-singularitiesbroadcast(_, L::LegendreWeight) = L # Assume we stay smooth
-singularitiesbroadcast(::typeof(exp), L::LegendreWeight) = L
-singularitiesbroadcast(::typeof(Base.literal_pow), ::typeof(^), L::LegendreWeight, ::Val) = L
 
 
 for op in (:+, :*)
     @eval singularitiesbroadcast(::typeof($op), A, B, C, D...) = singularitiesbroadcast(*, singularitiesbroadcast(*, A, B), C, D...)
 end
-for op in (:+, :-, :*)
-    @eval begin
-        singularitiesbroadcast(::typeof($op), A::LegendreWeight, B::LegendreWeight) = LegendreWeight{promote_type(eltype(A), eltype(B))}()
-        singularitiesbroadcast(::typeof($op), L::LegendreWeight, ::NoSingularities) = L
-        singularitiesbroadcast(::typeof($op), ::NoSingularities, L::LegendreWeight) = L
-    end
-end
-singularitiesbroadcast(::typeof(^), L::LegendreWeight, ::NoSingularities) = L
-singularitiesbroadcast(::typeof(/), ::NoSingularities, L::LegendreWeight) = L # can't find roots
+
 
 _parent(::NoSingularities) = NoSingularities()
 _parent(a) = parent(a)
@@ -87,38 +54,13 @@ singularitiesbroadcast(F, V::Union{NoSingularities,SubQuasiArray}...) = singular
 singularitiesbroadcast(::typeof(*), V::Union{NoSingularities,SubQuasiArray}...) = singularitiesbroadcast(*, map(_parent,V)...)[_parentindices(V...)...]
 
 
-singularitiesbroadcast(::typeof(*), ::LegendreWeight, b::AbstractJacobiWeight) = b
-singularitiesbroadcast(::typeof(*), a::AbstractJacobiWeight, ::LegendreWeight) = a
-
 abstract type AbstractJacobi{T} <: OrthogonalPolynomial{T} end
 
-singularities(::AbstractJacobi{T}) where T = LegendreWeight{T}()
-singularities(::Inclusion{T,<:AbstractInterval}) where T = LegendreWeight{T}()
-singularities(d::Inclusion{T,<:Interval}) where T = LegendreWeight{T}()[affine(d,ChebyshevInterval{T}())]
-
-struct Legendre{T} <: AbstractJacobi{T} end
-Legendre() = Legendre{Float64}()
-
-legendre() = Legendre()
-legendre(d::AbstractInterval{T}) where T = Legendre{float(T)}()[affine(d,ChebyshevInterval{T}()), :]
-
-"""
-     legendrep(n, z)
-
-computes the `n`-th Legendre polynomial at `z`.
-"""
-legendrep(n::Integer, z::Number) = Base.unsafe_getindex(Legendre{typeof(z)}(), z, n+1)
+include("legendre.jl")
 
+singularitiesbroadcast(::typeof(*), ::LegendreWeight, b::AbstractJacobiWeight) = b
+singularitiesbroadcast(::typeof(*), a::AbstractJacobiWeight, ::LegendreWeight) = a
 
-==(::Legendre, ::Legendre) = true
-
-OrthogonalPolynomial(w::LegendreWeight{T}) where {T} = Legendre{T}()
-orthogonalityweight(::Legendre{T}) where T = LegendreWeight{T}()
-
-function qr(P::Legendre)
-    Q = Normalized(P)
-    QuasiQR(Q, Diagonal(Q.scaling))
-end
 
 struct Jacobi{T} <: AbstractJacobi{T}
     a::T
@@ -221,12 +163,8 @@ function plotgrid(Pn::SubQuasiArray{T,2,<:AbstractJacobi,<:Tuple{Inclusion,Any}}
     ChebyshevGrid{2,T}(40maximum(jr))
 end
 
-
-function ldiv(::Legendre{V}, f::AbstractQuasiVector) where V
-    T = ChebyshevT{V}()
-    [cheb2leg(paddeddata(T \ f)); zeros(V,∞)]
-end
-
+ldiv(P::Jacobi{V}, f::Inclusion{T}) where {T,V} = _op_ldiv(P, f)
+ldiv(P::Jacobi{V}, f::AbstractQuasiFill{T,1}) where {T,V} = _op_ldiv(P, f)
 function ldiv(P::Jacobi{V}, f::AbstractQuasiVector) where V
     T = ChebyshevT{V}()
     [cheb2jac(paddeddata(T \ f), P.a, P.b); zeros(V,∞)]
@@ -237,23 +175,7 @@ end
 # Mass Matrix
 #########
 
-# massmatrix(P) = Weighted(P)'P
-massmatrix(P::Legendre{T}) where T = Diagonal(convert(T,2) ./ (2(0:∞) .+ 1))
-
-
-
-"""
-    legendre_massmatrix
-
-computes the massmatrix by first re-expanding in Legendre
-"""
-function legendre_massmatrix(Ac, B)
-    A = parent(Ac)
-    P = Legendre{eltype(B)}()
-    (P\A)'*massmatrix(P)*(P\B)
-end
 
-@simplify *(Ac::QuasiAdjoint{<:Any,<:Legendre}, B::Legendre) = massmatrix(Legendre{promote_type(eltype(Ac), eltype(B))}())
 @simplify *(Ac::QuasiAdjoint{<:Any,<:AbstractJacobi}, B::AbstractJacobi) = legendre_massmatrix(Ac,B)
 @simplify *(Ac::QuasiAdjoint{<:Any,<:AbstractJacobi}, B::Weighted{<:Any,<:AbstractJacobi}) = legendre_massmatrix(Ac,B)
 
@@ -287,14 +209,6 @@ end
 # Jacobi Matrix
 ########
 
-jacobimatrix(::Legendre{T}) where T =  Tridiagonal((one(T):∞)./(1:2:∞), Zeros{T}(∞), (one(T):∞)./(3:2:∞))
-
-# These return vectors A[k], B[k], C[k] are from DLMF. Cause of MikaelSlevinsky we need an extra entry in C ... for now.
-function recurrencecoefficients(::Legendre{T}) where T
-    n = zero(T):∞
-    ((2n .+ 1) ./ (n .+ 1), Zeros{T}(∞), n ./ (n .+ 1))
-end
-
 function jacobimatrix(J::Jacobi)
     b,a = J.b,J.a
     n = 0:∞
@@ -316,19 +230,6 @@ function recurrencecoefficients(P::Jacobi)
     (A,B,C)
 end
 
-# explicit special case for normalized Legendre
-# todo: do we want these explicit constructors for normalized Legendre?
-# function jacobimatrix(::Normalized{<:Any,<:Legendre{T}}) where T
-#     b = (one(T):∞) ./sqrt.(4 .*(one(T):∞).^2 .-1)
-#     Symmetric(_BandedMatrix(Vcat(zeros(∞)', (b)'), ∞, 1, 0), :L)
-# end
-# function recurrencecoefficients(::Normalized{<:Any,<:Legendre{T}}) where T
-#     n = zero(T):∞
-#     nn = one(T):∞
-#     ((2n .+ 1) ./ (n .+ 1) ./ sqrt.(1 .-2 ./(3 .+2n)), Zeros{T}(∞), Vcat(zero(T),nn ./ (nn .+ 1) ./ sqrt.(1 .-4 ./(3 .+2nn))))
-# end
-
-@simplify *(X::Identity, P::Legendre) = ApplyQuasiMatrix(*, P, P\(X*P))
 
 
 
@@ -519,25 +420,11 @@ end
 end
 
 
-###
-# Splines
-###
-
-function \(A::Legendre, B::HeavisideSpline)
-    @assert B.points == -1:2:1
-    Vcat(1, Zeros(∞,1))
-end
 
 ###
 # sum
 ###
 
-function _sum(P::Legendre{T}, dims) where T
-    @assert dims == 1
-    Hcat(convert(T, 2), Zeros{T}(1,∞))
-end
-
-_sum(p::SubQuasiArray{T,1,Legendre{T},<:Tuple{Inclusion,Int}}, ::Colon) where T = parentindices(p)[2] == 1 ? convert(T, 2) : zero(T)
 _sum(P::AbstractJacobi{T}, dims) where T = 2 * (Legendre{T}() \ P)[1:1,:]