Unweighted Laplacian (#99)

* Update test_disk.jl

* Fix JacobiTriangle evaluation

* Add Rz

* Update test_triangle.jl

* redefine WeightedTriangle

* Update triangle.jl

* some support for general Weighted

* add ==, etc., higher order L

* increase coverage

* fix tests

* Update test_triangle.jl

* add unweighted Laplacian

* v0.2

* More Zernike conversions

Author: dlfivefifty

Project.toml

@@ -1,6 +1,6 @@
 name = "MultivariateOrthogonalPolynomials"
 uuid = "4f6956fd-4f93-5457-9149-7bfc4b2ce06d"
-version = "0.1.4"
+version = "0.2"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -44,8 +44,9 @@ StaticArrays = "1"
 julia = "1.6"
 
 [extras]
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["RecipesBase", "Test"]
+test = ["ForwardDiff", "RecipesBase", "Test"]

src/ModalInterlace.jl

@@ -14,6 +14,7 @@ axes(Z::ModalInterlace) = blockedrange.(oneto.(Z.MN))
 
 blockbandwidths(R::ModalInterlace) = R.bandwidths
 subblockbandwidths(::ModalInterlace) = (0,0)
+copy(M::ModalInterlace) = M
 
 
 function Base.view(R::ModalInterlace{T}, KJ::Block{2}) where T

src/MultivariateOrthogonalPolynomials.jl

@@ -1,10 +1,12 @@
 module MultivariateOrthogonalPolynomials
+using StaticArrays: iszero
+using QuasiArrays: AbstractVector
 using ClassicalOrthogonalPolynomials, FastTransforms, BlockBandedMatrices, BlockArrays, DomainSets,
       QuasiArrays, StaticArrays, ContinuumArrays, InfiniteArrays, InfiniteLinearAlgebra,
       LazyArrays, SpecialFunctions, LinearAlgebra, BandedMatrices, LazyBandedMatrices, ArrayLayouts,
       HarmonicOrthogonalPolynomials
 
-import Base: axes, in, ==, *, ^, \, copy, OneTo, getindex, size, oneto
+import Base: axes, in, ==, *, ^, \, copy, OneTo, getindex, size, oneto, all
 import DomainSets: boundary
 
 import QuasiArrays: LazyQuasiMatrix, LazyQuasiArrayStyle

src/disk.jl

@@ -130,6 +130,8 @@ axes(P::Zernike{T}) where T = (Inclusion(UnitDisk{T}()),blockedrange(oneto(∞))
 
 copy(A::Zernike) = A
 
+Base.summary(io::IO, P::Zernike) = print(io, "Zernike($(P.a), $(P.b))")
+
 orthogonalityweight(Z::Zernike) = ZernikeWeight(Z.a, Z.b)
 
 zerniker(ℓ, m, a, b, r::T) where T = sqrt(convert(T,2)^(m+a+b+2-iszero(m))/π) * r^m * normalizedjacobip((ℓ-m) ÷ 2, b, m+a, 2r^2-1)
@@ -248,6 +250,16 @@ getindex(W::WeightedZernikeLaplacianDiag, k::Integer) = W[findblockindex(axes(W,
     WZ.P * Diagonal(WeightedZernikeLaplacianDiag{eltype(eltype(WZ))}())
 end
 
+@simplify function *(Δ::Laplacian, Z::Zernike)
+    a,b = Z.a,Z.b
+    @assert a == 0
+    T = promote_type(eltype(eltype(Δ)),eltype(Z)) # TODO: remove extra eltype
+    D = Derivative(Inclusion(ChebyshevInterval{T}())) 
+    Δs = BroadcastVector{AbstractMatrix{T}}((C,B,A) -> 4(HalfWeighted{:b}(C)\(D*HalfWeighted{:b}(B)))*(B\(D*A)), Normalized.(Jacobi.(b+2,a:∞)), Normalized.(Jacobi.(b+1,(a+1):∞)), Normalized.(Jacobi.(b,a:∞)))
+    Δ = ModalInterlace(Δs, (ℵ₀,ℵ₀), (-2,2))
+    Zernike(a,b+2) * Δ
+end
+
 ###
 # Fractional Laplacian
 ###
@@ -290,20 +302,16 @@ fractionalcfs2d(l::Integer, m::Integer, β) = fractionalcfs(l,m,β,2)
 
 function \(A::Zernike{T}, B::Zernike{V}) where {T,V}
     TV = promote_type(T,V)
-    A.a == B.a && A.b == B.b && return Eye{TV}(∞)
-    @assert A.a == 0 && A.b == 1
-    @assert B.a == 0 && B.b == 0
-    ModalInterlace{TV}((Normalized.(Jacobi{TV}.(1,0:∞)) .\ Normalized.(Jacobi{TV}.(0,0:∞))) ./ sqrt(convert(TV, 2)), (ℵ₀,ℵ₀), (0,2))
+    A.a == B.a && A.b == B.b && return Eye{TV}((axes(A,2),))
+    st = Int(A.a - B.a + A.b - B.b)
+    ModalInterlace{TV}((Normalized.(Jacobi{TV}.(A.b,A.a:∞)) .\ Normalized.(Jacobi{TV}.(B.b,B.a:∞))) .* convert(TV, 2)^(-st/2), (ℵ₀,ℵ₀), (0,2st))
 end
 
-function \(A::Zernike{T}, B::Weighted{V,Zernike{V}}) where {T,V}
+function \(A::Zernike{T}, wB::Weighted{V,Zernike{V}}) where {T,V}
     TV = promote_type(T,V)
-    A.a == B.P.a == A.b == B.P.b == 0 && return Eye{TV}(∞)
-    if A.a == A.b == 0
-        @assert B.P.a == 0 && B.P.b == 1
-        ModalInterlace{TV}((Normalized.(Jacobi{TV}.(0, 0:∞)) .\ HalfWeighted{:a}.(Normalized.(Jacobi{TV}.(1, 0:∞)))) ./ sqrt(convert(TV, 2)), (ℵ₀,ℵ₀), (2,0))
-    else
-        Z = Zernike{TV}() 
-        (A \ Z) * (Z \ B)
-    end
+    B = wB.P
+    A.a == B.a == A.b == B.b == 0 && return Eye{TV}((axes(A,2),))
+    c = Int(B.a - A.a + B.b - A.b)
+    @assert iszero(B.a)
+    ModalInterlace{TV}((Normalized.(Jacobi{TV}.(A.b, A.a:∞)) .\ HalfWeighted{:a}.(Normalized.(Jacobi{TV}.(B.b, B.a:∞)))) .* convert(TV, 2)^(-c/2), (ℵ₀,ℵ₀), (2Int(B.b), 2Int(A.a+A.b)))
 end

src/triangle.jl

@@ -39,14 +39,42 @@ end
 
 TriangleWeight(a::T, b::T, c::T) where T = TriangleWeight{float(T),T}(a, b, c)
 
-const WeightedTriangle{T} = WeightedBasis{T,<:TriangleWeight,<:JacobiTriangle}
+const WeightedTriangle{T,V} = Weighted{T,JacobiTriangle{T,V}}
 
-WeightedTriangle(a, b, c) = TriangleWeight(a,b,c) .* JacobiTriangle(a,b,c)
+WeightedTriangle(a, b, c) = Weighted(JacobiTriangle(a,b,c))
 
 axes(P::TriangleWeight{T}) where T = (Inclusion(UnitTriangle{T}()),)
+function getindex(P::TriangleWeight, xy::StaticVector{2})
+    x,y = xy
+    x^P.a * y^P.b * (1-x-y)^P.c
+end
+
+all(::typeof(isone), w::TriangleWeight) = iszero(w.a) && iszero(w.b) && iszero(w.c)
+==(w::TriangleWeight, v::TriangleWeight) = w.a == v.a && w.b == v.b && w.c == v.c
+
+==(wA::WeightedTriangle, B::JacobiTriangle) = wA.P == B == JacobiTriangle(0,0,0)
+==(B::JacobiTriangle, wA::WeightedTriangle) = wA == B
+
+==(w_A::WeightedBasis{<:Any,<:TriangleWeight,<:JacobiTriangle}, w_B::WeightedBasis{<:Any,<:TriangleWeight,<:JacobiTriangle}) = arguments(w_A) == arguments(w_B)
+
+function ==(w_A::WeightedBasis{<:Any,<:TriangleWeight,<:JacobiTriangle}, B::JacobiTriangle)
+    w,A = arguments(w_A)
+    all(isone,w) && A == B
+end
+
+==(B::JacobiTriangle, w_A::WeightedBasis{<:Any,<:TriangleWeight,<:JacobiTriangle}) = w_A == B
+
+function ==(w_A::WeightedBasis{<:Any,<:TriangleWeight,<:JacobiTriangle}, wB::WeightedTriangle)
+    w,A = arguments(w_A)
+    w.a == A.a && w.b == A.b && w.c == A.c && A == wB.P
+end
+
+==(wB::WeightedTriangle, w_A::WeightedBasis{<:Any,<:TriangleWeight,<:JacobiTriangle}) = w_A == wB
 
 Base.summary(io::IO, P::TriangleWeight) = print(io, "x^$(P.a)*y^$(P.b)*(1-x-y)^$(P.c) on the unit triangle")
 
+orthogonalityweight(P::JacobiTriangle) = TriangleWeight(P.a, P.b, P.c)
+
 @simplify function *(Dx::PartialDerivative{1}, P::JacobiTriangle)
     a,b,c = P.a,P.b,P.c
     n = mortar(Fill.(oneto(∞),oneto(∞)))
@@ -117,6 +145,18 @@ function Ry(a,b,c)
     _BandedBlockBandedMatrix(dat, axes(k,1), (0,1), (0,1))
 end
 
+function Rz(a,b,c)
+    n = mortar(Fill.(oneto(∞),oneto(∞)))
+    k = mortar(Base.OneTo.(oneto(∞)))
+    dat = PseudoBlockArray(Vcat(
+        (-(k .+ (b-1) ) .* (n .+ k .+ (b+c-1)) ./ ((2n .+ (a+b+c)) .* (2k .+ (b+c-1) )))',
+        ((k .- n .- a ) .* (k .+ (b+c)) ./ ((2n .+ (a+b+c)) .* (2k .+ (b+c-1) )))',
+        (-(k .+ (b-1) ) .* (k .- n .- 1) ./ ((2n .+ (a+b+c)) .* (2k .+ (b+c-1) )))',
+        ((n .+ k .+ (a+b+c) ) .* (k .+ (b+c)) ./ ((2n .+ (a+b+c)) .* (2k .+ (b+c-1) )))'
+        ), (blockedrange(Fill(2,2)), axes(n,1)))
+    _BandedBlockBandedMatrix(dat, axes(k,1), (0,1), (0,1))
+end
+
 function Lx(a,b,c)
     n = mortar(Fill.(oneto(∞),oneto(∞)))
     k = mortar(Base.OneTo.(oneto(∞)))
@@ -151,35 +191,94 @@ end
 
 
 function \(w_A::WeightedTriangle, w_B::WeightedTriangle)
-    wA,A = w_A.args
-    wB,B = w_B.args
-
-    @assert wA.a == A.a && wA.b == A.b && wA.c == A.c
-    @assert wB.a == B.a && wB.b == B.b && wB.c == B.c
+    A,B = w_A.P,w_B.P
 
-    if A.a + 1 == B.a && A.b == B.b && A.c == B.c
+    if A.a == B.a && A.b == B.b && A.c == B.c
+        Eye{promote_type(eltype(A),eltype(B))}((axes(B,2),))
+    elseif A.a + 1 == B.a && A.b == B.b && A.c == B.c
         Lx(B.a, B.b, B.c)
     elseif A.a == B.a && A.b + 1 == B.b && A.c == B.c
         Ly(B.a, B.b, B.c)
     elseif A.a == B.a && A.b == B.b && A.c + 1 == B.c
         Lz(B.a, B.b, B.c)
+    elseif A.a < B.a
+        w_C = WeightedTriangle(A.a+1,A.b,A.c)
+        (w_A \ w_C) * (w_C \ w_B)
+    elseif A.b < B.b
+        w_C = WeightedTriangle(A.a,A.b+1,A.c)
+        (w_A \ w_C) * (w_C \ w_B)
+    elseif A.c < B.c
+        w_C = WeightedTriangle(A.a,A.b,A.c+1)
+        (w_A \ w_C) * (w_C \ w_B)
     else
-        error("not implemented for $A and $wB")
+        error("not implemented for $w_A and $w_B")
     end
 end
 
-\(w_A::JacobiTriangle, w_B::WeightedTriangle) =
-    (TriangleWeight(0,0,0) .* w_A) \ w_B
+function \(w_A::WeightedBasis{<:Any,<:TriangleWeight,<:JacobiTriangle}, w_B::WeightedTriangle)
+    wA,A = w_A.args
+    w_A == Weighted(A) && Weighted(A) \ w_B
+    all(isone,wA) && return A \ w_B
+    error("Not implemented")
+end
+
+function \(w_A::WeightedTriangle, w_B::WeightedBasis{<:Any,<:TriangleWeight,<:JacobiTriangle})
+    wB,B = w_B.args
+    w_B == Weighted(B) && w_A \ Weighted(B)
+    all(isone,wB) && return w_A \ B
+    error("Not implemented")
+end
+
+
+function \(w_A::WeightedBasis{<:Any,<:TriangleWeight,<:JacobiTriangle}, w_B::WeightedBasis{<:Any,<:TriangleWeight,<:JacobiTriangle})
+    wA,A = w_A.args
+    w_A == Weighted(A) && Weighted(A) \ w_B
+    all(isone,wA) && return A \ w_B
+    error("Not implemented")
+end
+
+function \(w_A::WeightedBasis{<:Any,<:TriangleWeight,<:JacobiTriangle}, B::JacobiTriangle)
+    wA,A = w_A.args
+    w_A == Weighted(A) && return Weighted(A) \ B
+    all(isone,wA) && return A \ B
+    error("Not implemented")
+end
+
+function \(A::JacobiTriangle, w_B::WeightedBasis{<:Any,<:TriangleWeight,<:JacobiTriangle})
+    wB,B = w_B.args
+    w_B == Weighted(B) && return A \ Weighted(B)
+    all(isone,wB) && return A \ B
+    error("Not implemented")
+end
+
+function \(A::JacobiTriangle, w_B::WeightedTriangle)
+    w_B.P == JacobiTriangle() && return A \ w_B.P
+    A == JacobiTriangle() && return Weighted(A) \ w_B
+    (TriangleWeight(0,0,0) .* A) \ w_B
+end
+function \(w_A::WeightedTriangle, B::JacobiTriangle)
+    w_A.P == JacobiTriangle() && return w_A.P \ B
+    w_A \ (TriangleWeight(0,0,0) .* B)
+end
 
 function \(A::JacobiTriangle, B::JacobiTriangle)
-    if A.a == B.a && A.b == B.b && A.c == B.c
-        Eye((axes(B,2),))
+    if A == B
+        Eye{promote_type(eltype(A),eltype(B))}((axes(B,2),))
     elseif A.a == B.a + 1 && A.b == B.b && A.c == B.c
         Rx(B.a, B.b, B.c)
     elseif A.a == B.a && A.b == B.b + 1 && A.c == B.c
         Ry(B.a, B.b, B.c)
-    # elseif A.a == B.a && A.b == B.b && A.c == B.c + 1
-    #     Rz(B.a, B.b, B.c)
+    elseif A.a == B.a && A.b == B.b && A.c == B.c + 1
+        Rz(B.a, B.b, B.c)
+    elseif A.a > B.a
+        C = JacobiTriangle(B.a+1,B.b,B.c)
+        (A \ C) * (C \ B)
+    elseif A.b > B.b
+        C = JacobiTriangle(B.a,B.b+1,B.c)
+        (A \ C) * (C \ B)
+    elseif A.c > B.c
+        C = JacobiTriangle(B.a,B.b,B.c+1)
+        (A \ C) * (C \ B)
     else
         error("not implemented for $A and $B")
     end
@@ -431,14 +530,18 @@ function tri_forwardrecurrence(N::Int, X, Y, x, y)
     N < 1 && return ret
     ret[1] = 1
     N < 2 && return ret
-    ret[2] = x/X[1,1]-1
-    ret[3] = -Y[2,1]/(Y[3,1]*X[1,1]) * (x-X[1,1]) + (y-Y[1,1])/Y[3,1]
 
     X_N = X[Block.(1:N), Block.(1:N-1)]
     Y_N = Y[Block.(1:N), Block.(1:N-1)]
 
+    let n = 1
+        A = TriangleRecurrenceA(n, X_N, Y_N)
+        B = TriangleRecurrenceB(n, X_N, Y_N)
+        ret[Block(2)] .= A*[x; y] + vec(B)
+    end
+
     for n = 2:N-1
-        # P[n+1,xy] == (A[n]*[x*Eye(n); y*Eye(n)] + B[n])*P[n,xy] - C[n]*P[n-1,xy]
+        # P[n+1,xy] == (A*[x*Eye(n); y*Eye(n)] + B)*P[n,xy] - C*P[n-1,xy]
         A = TriangleRecurrenceA(n, X_N, Y_N)
         B = TriangleRecurrenceB(n, X_N, Y_N)
         C = TriangleRecurrenceC(n, X_N, Y_N)