Weighted(Jacobi(1,1)) -> Ultraspherical(-1/2) for simpler operators (#60)

* Weighted(Jacobi(1,1)) -> Ultraspherical(-1/2) for simpler operators

* tests pass

* tests pass

* fix dirichlet

* Update test_dirichlet.jl

Author: dlfivefifty

Project.toml

@@ -1,7 +1,7 @@
 name = "PiecewiseOrthogonalPolynomials"
 uuid = "4461d12d-4663-4550-8580-cb764c85e20f"
 authors = ["Sheehan Olver <solver@mac.com>"]
-version = "0.3"
+version = "0.4"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"

src/PiecewiseOrthogonalPolynomials.jl

@@ -9,8 +9,10 @@ import ClassicalOrthogonalPolynomials: grid, plotgrid, ldiv, pad, adaptivetransf
 import ContinuumArrays: @simplify, factorize, TransformFactorization, AbstractBasisLayout, MemoryLayout, layout_broadcasted, ExpansionLayout, basis, plan_transform, plan_grid_transform, grammatrix, weaklaplacian, MulPlan, InvPlan
 import LazyArrays: paddeddata
 import LazyBandedMatrices: BlockBroadcastMatrix, BlockVec, BandedLazyLayouts, AbstractLazyBandedBlockBandedLayout, UpperOrLowerTriangular
-import Base: size, axes, getindex, +, -, *, /, ==, \, OneTo, oneto, replace_in_print_matrix, copy, diff, getproperty, adjoint, transpose, tail, _sum, inv
+import Base: size, axes, getindex, +, -, *, /, ==, \, OneTo, oneto, replace_in_print_matrix, copy, diff, getproperty, adjoint, transpose, tail, _sum, inv, show, summary
 import LinearAlgebra: BlasInt
+import InfiniteArrays: OneToInf
+import FillArrays: AbstractFill
 import MatrixFactorizations: reversecholcopy
 
 export PiecewisePolynomial, ContinuousPolynomial, DirichletPolynomial, Derivative, Block, weaklaplacian, grammatrix

src/arrowhead.jl

@@ -112,6 +112,8 @@ arrowheadlayout(::BandedLazyLayouts) = LazyArrowheadLayout()
 symmetriclayout(lay::ArrowheadLayouts) = SymmetricLayout{typeof(lay)}()
 
 MemoryLayout(::Type{<:ArrowheadMatrix{<:Any,<:Any,<:Any,<:Any,<:AbstractVector{D}}}) where D = arrowheadlayout(MemoryLayout(D))
+# Hack since ∞ diagonal fill is only DiagonalLayout{FillLayout}
+MemoryLayout(::Type{<:ArrowheadMatrix{<:Any,<:Any,<:Any,<:Any,<:AbstractVector{<:Diagonal{<:Any,<:AbstractFill{<:Any,1,<:Tuple{OneToInf}}}}}}) = LazyArrowheadLayout()
 
 sublayout(::ArrowheadLayouts,
           ::Type{<:NTuple{2,BlockSlice{<:BlockRange{1, Tuple{OneTo{Int}}}}}}) = ArrowheadLayout()

src/continuouspolynomial.jl

@@ -14,6 +14,9 @@ axes(B::ContinuousPolynomial{0}) = axes(PiecewisePolynomial(B))
 axes(B::ContinuousPolynomial{1}) =
     (Inclusion(first(B.points) .. last(B.points)), blockedrange(Vcat(length(B.points), Fill(length(B.points) - 1, ∞))))
 
+show(io::IO, Q::ContinuousPolynomial{λ}) where λ = summary(io, Q)
+summary(io::IO, Q::ContinuousPolynomial{λ}) where λ = print(io, "ContinuousPolynomial{$λ}($(Q.points))")
+
 ==(P::PiecewisePolynomial, C::ContinuousPolynomial{0}) = P == PiecewisePolynomial(C)
 ==(C::ContinuousPolynomial{0}, P::PiecewisePolynomial) = PiecewisePolynomial(C) == P
 ==(::PiecewisePolynomial, ::ContinuousPolynomial{1}) = false
@@ -25,6 +28,8 @@ axes(B::ContinuousPolynomial{1}) =
 
 getindex(P::ContinuousPolynomial{0,T}, x::Number, Kk::BlockIndex{1}) where {T} = PiecewisePolynomial(P)[x, Kk]
 
+bubblebasis(::Type{T}) where T = Ultraspherical{T}(-one(real(T))/2)
+
 function getindex(P::ContinuousPolynomial{1,T}, x::Number, Kk::BlockIndex{1}) where {T}
     K, k = block(Kk), blockindex(Kk)
     if K == Block(1)
@@ -33,7 +38,7 @@ function getindex(P::ContinuousPolynomial{1,T}, x::Number, Kk::BlockIndex{1}) wh
         b = searchsortedlast(P.points, x)
         if b == k
             α, β = convert(T, P.points[b]), convert(T, P.points[b+1])
-            Weighted(Jacobi{T}(1, 1))[affine(α.. β, ChebyshevInterval{real(T)}())[x], Int(K)-1]
+            bubblebasis(T)[affine(α.. β, ChebyshevInterval{real(T)}())[x], Int(K)+1]
         else
             zero(T)
         end
@@ -66,11 +71,11 @@ end
 
 function plan_transform(C::ContinuousPolynomial{1,T}, Ns::NTuple{N,Block{1}}, dims=ntuple(identity,Val(N))) where {N,T}
     P = Legendre{T}()
-    W = Weighted(Jacobi{T}(1,1))
+    W = bubblebasis(T)
     # TODO: this is unnecessarily not triangular which makes it much slower than necessary and prevents in-place.
     # However, the speed of the Legendre transform will far exceed this so its not high priority.
     # Probably using Ultraspherical(-1/2) would be better.
-    R̃ = [[T[1 1; -1 1]/2; Zeros{T}(∞,2)] (P \ W)]
+    R̃ = [[T[1 1; -1 1]/2; Zeros{T}(∞,2)] (P \ W)[:,3:end]]
     ContinuousPolynomialTransform(plan_transform(ContinuousPolynomial{0}(C), Ns .+ 1, dims).F, 
                                   InvPlan(map(N -> R̃[1:Int(N),1:Int(N)], Ns .+ 1), _doubledims(dims...)),
                                   dims)
@@ -167,11 +172,11 @@ function adaptivetransform_ldiv(Q::ContinuousPolynomial{1,V}, f::AbstractQuasiVe
     c̃ = paddeddata(c)
     N = max(2,div(length(c̃), M, RoundUp)) # degree
     P = Legendre{T}()
-    W = Weighted(Jacobi{T}(1,1))
+    W = bubblebasis(T)
 
     # Restrict hat function to each element, add in bubble functions and compute connection
     # matrix to Legendre. [1 1; -1 1]/2 are the Legendre coefficients of the hat functions.
-    R̃ = [[T[1 1; -1 1]/2; Zeros{T}(∞,2)] (P \ W)]
+    R̃ = [[T[1 1; -1 1]/2; Zeros{T}(∞,2)] (P \ W)[:,3:end]]
 
     # convert from Legendre to piecewise restricted hat + Bubble
     dat = R̃[1:N,1:N] \ reshape(pad(c̃, M*N), M, N)'
@@ -200,7 +205,7 @@ grid(P::ContinuousPolynomial{1}, M::Block{1}) = grid(ContinuousPolynomial{0}(P),
 function \(P::ContinuousPolynomial{0}, C::ContinuousPolynomial{1})
     T = promote_type(eltype(P), eltype(C))
     @assert P.points == C.points
-    v = (convert(T, 2):2:∞) ./ (3:2:∞)
+    v = one(T) ./ (3:2:∞)
     N = length(P.points)
     ArrowheadMatrix(_BandedMatrix(Ones{T}(2, N)/2, oneto(N-1), 0, 1),
         (_BandedMatrix(Fill(v[1], 1, N-1), oneto(N-1), 0, 0),),
@@ -234,15 +239,15 @@ function grammatrix(C::ContinuousPolynomial{1, T, <:AbstractRange}) where T
 
     N = length(r) - 1
     h = step(r) # 2/N
-    a = ((convert(T,4):4:∞) .* (convert(T,-2):2:∞)) ./ ((1:2:∞) .* (3:2:∞) .* (-1:2:∞))
-    b = (((convert(T,2):2:∞) ./ (3:2:∞)).^2 .* (convert(T,2) ./ (1:2:∞) .+ convert(T,2) ./ (5:2:∞)))
+    a = [-2 /((2k+1)*(2k+3)*(2k+5)) for k = -1:∞]
+    b = [4 /((2k+1)*(2k+3)*(2k+5)) for k = 0:∞]
 
     a11 = LazyBandedMatrices.Bidiagonal(Vcat(h/3, Fill(2h/3, N-1), h/3), Fill(h/6, N), :U)
-    a21 = _BandedMatrix(Fill(h/3, 2, N), N+1, 1, 0)
-    a31 = _BandedMatrix(Vcat(Fill(-2h/15, 1, N), Fill(2h/15, 1, N)), N+1, 1, 0)
+    a21 = _BandedMatrix(Fill(h/6, 2, N), N+1, 1, 0)
+    a31 = _BandedMatrix(Vcat(Fill(-h/30, 1, N), Fill(h/30, 1, N)), N+1, 1, 0)
 
     Symmetric(ArrowheadMatrix(a11, (a21, a31), (),
-                Fill(_BandedMatrix(Vcat((-h*a/2)',
+                Fill(_BandedMatrix(Vcat((h*a/2)',
                 Zeros{T}(1,∞),
                 (h*b/2)'), ∞, 0, 2), N)))
 end
@@ -278,12 +283,17 @@ function diff(C::ContinuousPolynomial{1,T}; dims=1) where T
     r = C.points
     N = length(r)
     s = one(T) ./ (r[2:end]-r[1:end-1])
-    v = mortar(Fill(T(2) * s, ∞)) .* mortar(Fill.((-convert(T, 2):-2:-∞), N - 1))
-    z = Zeros{T}(axes(v))
-    H = BlockBroadcastArray(hcat, z, v)
-    M = BlockVcat(Hcat(Ones{T}(N) .* [zero(T); s] , -Ones{T}(N) .* [s; zero(T)] ), H)
-    P = ContinuousPolynomial{0}(C)
-    ApplyQuasiMatrix(*, P, _BandedBlockBandedMatrix(M', axes(P, 2), (0, 0), (0, 1)))
+    ContinuousPolynomial{0}(r) * ArrowheadMatrix(_BandedMatrix(Vcat([0; s]', [-s; 0]'), length(s), 0, 1), (), (),
+                                                 (-2s) .* Ref(Eye{T}(∞)))
+end
+
+function diff(C::ContinuousPolynomial{1,T,<:AbstractRange}; dims=1) where T
+    # Legendre() \ (D*Weighted(Jacobi(1,1)))
+    r = C.points
+    N = length(r)
+    s = 1 ./ step(r)
+    ContinuousPolynomial{0}(r) * ArrowheadMatrix(_BandedMatrix(Vcat(Fill(s, 1, N), Fill(-s, 1, N)), N-1, 0, 1), (), (),
+                                                 Fill(-2s*Eye{T}(∞), N-1))
 end
 
 function weaklaplacian(C::ContinuousPolynomial{1,T,<:AbstractRange}) where T
@@ -294,7 +304,7 @@ function weaklaplacian(C::ContinuousPolynomial{1,T,<:AbstractRange}) where T
     t1 = Vcat(-si, Fill(-2si, N-2), -si)
     t2 = Fill(si, N-1)
     Symmetric(ArrowheadMatrix(LazyBandedMatrices.Bidiagonal(t1, t2, :U), (), (),
-        Fill(Diagonal(convert(T, -16) .* (1:∞) .^ 2 ./ (s .* ((2:2:∞) .+ 1))), N-1)))
+        Fill(Diagonal(convert(T,-4) ./ (s*(convert(T,3):2:∞))), N-1)))
 end
 
 

src/dirichlet.jl

@@ -11,6 +11,10 @@ PiecewisePolynomial(P::DirichletPolynomial{T}) where {T} = PiecewisePolynomial(C
 
 axes(B::DirichletPolynomial) = (Inclusion(first(B.points) .. last(B.points)), blockedrange(Vcat(length(B.points)-2, Fill(length(B.points) - 1, ∞))))
 
+show(io::IO, Q::DirichletPolynomial) = summary(io, Q)
+summary(io::IO, Q::DirichletPolynomial) = print(io, "DirichletPolynomial($(Q.points))")
+
+
 ==(P::DirichletPolynomial, C::DirichletPolynomial) = P.points == C.points
 ==(P::PiecewisePolynomial, C::DirichletPolynomial) = false
 ==(C::DirichletPolynomial, P::PiecewisePolynomial) = false
@@ -170,23 +174,23 @@ function weaklaplacian(C::DirichletPolynomial{T,<:AbstractRange}) where T
     t1 = Fill(-2si, N-2)
     t2 = Fill(si, N-3)
     Symmetric(ArrowheadMatrix(LazyBandedMatrices.Bidiagonal(t1, t2, :U), (), (),
-        Fill(Diagonal(convert(T, -16) .* (1:∞) .^ 2 ./ (s .* ((2:2:∞) .+ 1))), N-1)))
+        Fill(Diagonal(convert(T, -4) ./ (s*(convert(T,3):2:∞))), N-1)))
 end
 
 function grammatrix(Q::DirichletPolynomial{T, <:AbstractRange}) where T
     r = Q.points
 
     N = length(r) - 1
     h = step(r) # 2/N
-    a = ((convert(T,4):4:∞) .* (convert(T,-2):2:∞)) ./ ((1:2:∞) .* (3:2:∞) .* (-1:2:∞))
-    b = (((convert(T,2):2:∞) ./ (3:2:∞)).^2 .* (convert(T,2) ./ (1:2:∞) .+ convert(T,2) ./ (5:2:∞)))
+    a = [-2 /((2k+1)*(2k+3)*(2k+5)) for k = -1:∞]
+    b = [4 /((2k+1)*(2k+3)*(2k+5)) for k = 0:∞]
 
     a11 = LazyBandedMatrices.Bidiagonal(Fill(2h/3, N-1), Fill(h/6, N-2), :U)
-    a21 = _BandedMatrix(Fill(h/3, 2, N), N-1, 0, 1)
-    a31 = _BandedMatrix(Vcat(Fill(-2h/15, 1, N), Fill(2h/15, 1, N)), N-1, 0, 1)
+    a21 = _BandedMatrix(Fill(h/6, 2, N), N-1, 0, 1)
+    a31 = _BandedMatrix(Vcat(Fill(-h/30, 1, N), Fill(h/30, 1, N)), N-1, 0, 1)
 
     Symmetric(ArrowheadMatrix(a11, (a21, a31), (),
-                Fill(_BandedMatrix(Vcat((-h*a/2)',
+                Fill(_BandedMatrix(Vcat((h*a/2)',
                 Zeros{T}(1,∞),
                 (h*b/2)'), ∞, 0, 2), N)))
 end

test/test_arrowhead.jl

@@ -8,7 +8,7 @@ import Base: oneto, OneTo
 @testset "ArrowheadMatrix" begin
     @testset "Algebra" begin
         n = 4; p = 5;
-        A = ArrowheadMatrix(BandedMatrix(0 => randn(n) .+ 10, 1 => randn(n-1), -1 => randn(n-1)), 
+        A = ArrowheadMatrix(BandedMatrix(0 => randn(n) .+ 10, 1 => randn(n-1), -1 => randn(n-1)),
                                 ntuple(_ -> BandedMatrix((0 => randn(n-1), -1 => randn(n-1)), (n,n-1)), 2),
                                 ntuple(_ -> BandedMatrix((0 => randn(n), 1 => randn(n-1)), (n-1,n)), 3),
                             fill(BandedMatrix((0 => randn(p) .+ 10, 2 => randn(p-2), -1=> randn(p-1)), (p, p)), n-1))
@@ -30,7 +30,7 @@ import Base: oneto, OneTo
 
     @testset "mul" begin
         n = 4; p = 5;
-        A = ArrowheadMatrix(BandedMatrix(0 => randn(n) .+ 10, 1 => randn(n-1), -1 => randn(n-1)), 
+        A = ArrowheadMatrix(BandedMatrix(0 => randn(n) .+ 10, 1 => randn(n-1), -1 => randn(n-1)),
                                 ntuple(_ -> BandedMatrix((0 => randn(n-1), -1 => randn(n-1)), (n,n-1)), 2),
                                 ntuple(_ -> BandedMatrix((0 => randn(n), 1 => randn(n-1)), (n-1,n)), 3),
                             fill(BandedMatrix((0 => randn(p) .+ 10, 2 => randn(p-2), -1=> randn(p-1)), (p, p)), n-1))
@@ -47,7 +47,7 @@ import Base: oneto, OneTo
 
     @testset "UpperTriangular" begin
         n = 4; p = 5;
-        A = ArrowheadMatrix(BandedMatrix(0 => randn(n) .+ 10, 1 => randn(n-1), -1 => randn(n-1)), 
+        A = ArrowheadMatrix(BandedMatrix(0 => randn(n) .+ 10, 1 => randn(n-1), -1 => randn(n-1)),
                                 [BandedMatrix((0 => randn(n-1), -1 => randn(n-1)), (n,n-1)) for _=1:2],
                                 [BandedMatrix((0 => randn(n), 1 => randn(n-1)), (n-1,n)) for _=1:2],
                             fill(BandedMatrix((0 => randn(p) .+ 10, 2 => randn(p-2), -1=> randn(p-1)), (p, p)), n-1))
@@ -74,15 +74,15 @@ import Base: oneto, OneTo
 
     @testset "reversecholesky" begin
         n = 3; p = 5;
-        A = ArrowheadMatrix(BandedMatrix(0 => randn(n) .+ 10, 1 => randn(n-1), -1 => randn(n-1)), 
+        A = ArrowheadMatrix(BandedMatrix(0 => randn(n) .+ 10, 1 => randn(n-1), -1 => randn(n-1)),
                                 [BandedMatrix((0 => randn(n-1), -1 => randn(n-1)), (n,n-1)) for _=1:2],
                                 BandedMatrix{Float64, Matrix{Float64}, OneTo{Int}}[],
                             fill(BandedMatrix((0 => randn(p) .+ 10, 2 => randn(p-2)), (p, p)), n-1))
 
         @test reversecholesky(Symmetric(Matrix(A))).U ≈ reversecholesky!(Symmetric(copy(A))).U
 
 
-        A = ArrowheadMatrix(BandedMatrix(0 => randn(n) .+ 10, 1 => randn(n-1), -1 => randn(n-1)), 
+        A = ArrowheadMatrix(BandedMatrix(0 => randn(n) .+ 10, 1 => randn(n-1), -1 => randn(n-1)),
                                 [BandedMatrix((0 => randn(n-1), -1 => randn(n-1)), (n,n-1)) for _=1:2],
                                 BandedMatrix{Float64, Matrix{Float64}, OneTo{Int}}[],
                             fill(BandedMatrix((0 => randn(p) .+ 10, 1 => randn(p-1), 2 => randn(p-2)), (p, p)), n-1))
@@ -103,7 +103,7 @@ import Base: oneto, OneTo
         @test blockbandwidths(L) == (1,1)
         @test subblockbandwidths(L) == (1,1)
 
-        @test stringmime("text/plain", L[Block.(OneTo(3)), Block.(OneTo(3))]) == "3×3-blocked 9×10 ArrowheadMatrix{Float64}:\n  0.5   0.5    ⋅    ⋅   │   0.666667    ⋅          ⋅        │   ⋅    ⋅    ⋅ \n   ⋅    0.5   0.5   ⋅   │    ⋅         0.666667    ⋅        │   ⋅    ⋅    ⋅ \n   ⋅     ⋅    0.5  0.5  │    ⋅          ⋅         0.666667  │   ⋅    ⋅    ⋅ \n ───────────────────────┼───────────────────────────────────┼───────────────\n -0.5   0.5    ⋅    ⋅   │   0.0         ⋅          ⋅        │  0.8   ⋅    ⋅ \n   ⋅   -0.5   0.5   ⋅   │    ⋅         0.0         ⋅        │   ⋅   0.8   ⋅ \n   ⋅     ⋅   -0.5  0.5  │    ⋅          ⋅         0.0       │   ⋅    ⋅   0.8\n ───────────────────────┼───────────────────────────────────┼───────────────\n   ⋅     ⋅     ⋅    ⋅   │  -0.666667    ⋅          ⋅        │  0.0   ⋅    ⋅ \n   ⋅     ⋅     ⋅    ⋅   │    ⋅        -0.666667    ⋅        │   ⋅   0.0   ⋅ \n   ⋅     ⋅     ⋅    ⋅   │    ⋅          ⋅        -0.666667  │   ⋅    ⋅   0.0"
+        @test stringmime("text/plain", L[Block.(OneTo(3)), Block.(OneTo(3))]) == "3×3-blocked 9×10 ArrowheadMatrix{Float64}:\n  0.5   0.5    ⋅    ⋅   │   0.333333    ⋅          ⋅        │   ⋅    ⋅    ⋅ \n   ⋅    0.5   0.5   ⋅   │    ⋅         0.333333    ⋅        │   ⋅    ⋅    ⋅ \n   ⋅     ⋅    0.5  0.5  │    ⋅          ⋅         0.333333  │   ⋅    ⋅    ⋅ \n ───────────────────────┼───────────────────────────────────┼───────────────\n -0.5   0.5    ⋅    ⋅   │   0.0         ⋅          ⋅        │  0.2   ⋅    ⋅ \n   ⋅   -0.5   0.5   ⋅   │    ⋅         0.0         ⋅        │   ⋅   0.2   ⋅ \n   ⋅     ⋅   -0.5  0.5  │    ⋅          ⋅         0.0       │   ⋅    ⋅   0.2\n ───────────────────────┼───────────────────────────────────┼───────────────\n   ⋅     ⋅     ⋅    ⋅   │  -0.333333    ⋅          ⋅        │  0.0   ⋅    ⋅ \n   ⋅     ⋅     ⋅    ⋅   │    ⋅        -0.333333    ⋅        │   ⋅   0.0   ⋅ \n   ⋅     ⋅     ⋅    ⋅   │    ⋅          ⋅        -0.333333  │   ⋅    ⋅   0.0"
     end
 
     @testset "Helmholtz solve" begin
@@ -128,7 +128,7 @@ import Base: oneto, OneTo
 
     @testset "Dirichlet" begin
         n = 5; p = 5;
-        A = ArrowheadMatrix(BandedMatrix(0 => randn(n-2) .+ 10, 1 => randn(n-3), -1 => randn(n-3)), 
+        A = ArrowheadMatrix(BandedMatrix(0 => randn(n-2) .+ 10, 1 => randn(n-3), -1 => randn(n-3)),
                                 [BandedMatrix((0 => randn(n-2), 1 => randn(n-2)), (n-2,n-1)) for _=1:2],
                                 BandedMatrix{Float64, Matrix{Float64}, OneTo{Int}}[],
                             fill(BandedMatrix((0 => randn(p) .+ 10, 2 => randn(p-2)), (p, p)), n-1))

test/test_continuouspolynomial.jl

@@ -61,6 +61,7 @@ using ContinuumArrays: plan_grid_transform
             JR = Block.(1:10)
             KR = Block.(1:11)
             R = P\C
+            @test P[0.1,1:10]'*R[1:10,1:10] ≈ C[0.1,1:10]'
             @test R[KR,JR]'*((P'P)[KR,KR]*R[KR,JR]) ≈ (C'C)[JR,JR]
             @test (P'C)[JR,JR] ≈ (C'P)[JR,JR]'
         end
@@ -87,6 +88,8 @@ using ContinuumArrays: plan_grid_transform
             D = Derivative(x)
             A = P\D*C
 
+            @test (D*expand(C, exp))[0.1] ≈ exp(0.1)
+
             xx = rand(5)
             c = (x, p) -> ContinuousPolynomial{1, eltype(x)}(r)[x, p]
 
@@ -147,7 +150,7 @@ using ContinuumArrays: plan_grid_transform
         x = SVector(0.1, 0.2)
         @test view(C, x, 1) .* view(C, x, 1) ≈ C[x,1] .* C[x,1]
 
-        @test C \ (C[:,1] .* C[:,1]) ≈ [1; zeros(3); -1/4; zeros(∞)]
+        @test C \ (C[:,1] .* C[:,1]) ≈ [1; zeros(3); -1/2; zeros(∞)]
         @test C \ (C[:,1] .* C[:,3]) ≈ zeros(∞)
     end
 
@@ -158,8 +161,8 @@ using ContinuumArrays: plan_grid_transform
         x = axes(P,1)
         D = Derivative(x)
 
-        a = expand(P, x -> abs(x) ≤ 1/3 ? 2 : 3)    
-        L = (D*C)'* (a .* (D*C)) 
+        a = expand(P, x -> abs(x) ≤ 1/3 ? 2 : 3)
+        L = (D*C)'* (a .* (D*C))
     end
 
     @testset "expand" begin

test/test_dirichlet.jl

@@ -69,4 +69,14 @@ using PiecewiseOrthogonalPolynomials: ArrowheadMatrix, plan_grid_transform
         @test Q[0.1, Block.(1:20)]' * V * Q[0.2,Block.(1:21)] ≈ f(0.1,0.2)
         @test F \ V ≈ vals
     end
+
+    @testset "show" begin
+        Q = DirichletPolynomial(-1:1)
+        C = ContinuousPolynomial{1}(-1:1)
+        P = ContinuousPolynomial{0}(-1:1)
+
+        @test sprint(show, Q) == "DirichletPolynomial(-1:1)"
+        @test sprint(show, C) == "ContinuousPolynomial{1}(-1:1)"
+        @test sprint(show, P) == "ContinuousPolynomial{0}(-1:1)"
+    end
 end