Support 3-interval Hilbert (#74)

* Support 3-interval Hilbert

* test 3-interval ODE solve

Author: dlfivefifty

Project.toml

@@ -1,7 +1,7 @@
 name = "ClassicalOrthogonalPolynomials"
 uuid = "b30e2e7b-c4ee-47da-9d5f-2c5c27239acd"
 authors = ["Sheehan Olver <[email protected]>"]
-version = "0.4.8"
+version = "0.4.9"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -29,7 +29,7 @@ ArrayLayouts = "0.7.5"
 BandedMatrices = "0.16.11"
 BlockArrays = "0.16.6"
 BlockBandedMatrices = "0.11"
-ContinuumArrays = "0.9.1"
+ContinuumArrays = "0.9.4"
 DomainSets = "0.5.6"
 FFTW = "1.1"
 FastGaussQuadrature = "0.4.3"

src/stieltjes.jl

@@ -308,18 +308,17 @@ end
 ###
 
 
-@simplify function *(H::Hilbert, S::PiecewiseInterlace)
-    axes(H,2) == axes(S,1) || throw(DimensionMismatch())
-    @assert length(S.args) == 2
-    a,b = S.args
-    xa,xb = axes(a,1),axes(b,1)
-    Ha_a = inv.(xa .- xa') * a
-    Ha_b = inv.(xb .- xa') * a
-    Hb_a = inv.(xa .- xb') * b
-    Hb_b = inv.(xb .- xb') * b
-    c,d = Hb_a.args[1], Ha_b.args[1]
-    A,B,C,D = unitblocks(c \ Ha_a), unitblocks(c \ Hb_a), unitblocks(d \ Ha_b), unitblocks(d \ Hb_b)
-    PiecewiseInterlace(c,d)  * BlockBroadcastArray{promote_type(eltype(H),eltype(S))}(hvcat, 2, A, B, C, D)
+@simplify function *(H::Hilbert, W::PiecewiseInterlace)
+    axes(H,2) == axes(W,1) || throw(DimensionMismatch())
+    Hs = broadcast(function(a,b)
+                x,t = axes(a,1),axes(b,1)
+                H = inv.(x .- t') * b
+                H
+            end, [W.args...], permutedims([W.args...]))
+    N = length(W.args)
+    Ts = [broadcastbasis(+, broadcast(H -> H.args[1], Hs[k,:])...) for k=1:N]
+    Ms = broadcast((T,H) -> unitblocks(T\H), Ts, Hs)
+    PiecewiseInterlace(Ts...) * BlockBroadcastArray{promote_type(eltype(H),eltype(W))}(hvcat, N, permutedims(Ms)...)
 end
 
 

test/test_interlace.jl

@@ -58,6 +58,29 @@ import ClassicalOrthogonalPolynomials: PiecewiseInterlace, plotgrid
         # Δ \ (W'exp.(x))
     end
 
+    @testset "three-interval ODE" begin
+        d = (-1..0),(0..1),(1..2)
+        T = PiecewiseInterlace(chebyshevt.(d)...)
+        U = PiecewiseInterlace(chebyshevu.(d)...)
+
+        D = U \ (Derivative(axes(T,1))*T)
+        C = U \ T
+
+        A = BlockVcat(T[-1,:]',
+                        BlockBroadcastArray(vcat,unitblocks(T.args[1][end,:]),-unitblocks(T.args[2][begin,:]),Zeros((unitblocks(axes(T.args[3],2)),)))',
+                        BlockBroadcastArray(vcat,Zeros((unitblocks(axes(T.args[1],2)),)),unitblocks(T.args[2][end,:]),-unitblocks(T.args[3][begin,:]))',
+                        D-C)
+        N = 20
+        M = BlockArray(A[Block.(1:N+2), Block.(1:N)])
+
+        u = M \ [exp(-1); zeros(size(M,1)-1)]
+        x = axes(T,1)
+
+        F = factorize(T[:,Block.(Base.OneTo(N))])
+        @test F \ exp.(x) ≈ (T \ exp.(x))[Block.(1:N)] ≈ u
+    end
+
+
     @testset "plot" begin
         T1,T2 = chebyshevt(-1..0), chebyshevt(0..1)
         T = PiecewiseInterlace(T1, T2)

test/test_stieltjes.jl

@@ -266,7 +266,7 @@ end
 
                 @test Z * H[:,1] - H[:,2]/2 ≈ [sum(W[:,1]); zeros(∞)]
                 @test norm(-H[:,1]/2 + Z * H[:,2] - H[:,3]/2) ≤ 1E-12
-                
+
                 L = U \ ((x.^2 .- 1) .* Derivative(x) * T - x .* T)
                 c = T \ sqrt.(x.^2 .- 1)
                 @test [T[begin,:]'; L] \ [sqrt(2^2-1); zeros(∞)] ≈ c
@@ -351,6 +351,25 @@ end
         end
     end
 
+    @testset "three-interval" begin
+        d = (-2..(-1), 0..1, 2..3)
+        T = PiecewiseInterlace(chebyshevt.(d)...)
+        U = PiecewiseInterlace(chebyshevu.(d)...)
+        W = PiecewiseInterlace(Weighted.(U.args)...)
+        x = axes(W,1)
+        H = T \ inv.(x .- x') * W
+        c = W \ broadcast(x -> exp(x) *
+            if -2 ≤ x ≤ -1
+                sqrt(x+2)sqrt(-1-x)
+            elseif 0 ≤ x ≤ 1
+                sqrt(1-x)sqrt(x)
+            else
+                sqrt(x-2)sqrt(3-x)
+            end, x)
+        f = W * c
+        @test T[0.5,1:200]'*(H*c)[1:200] ≈ -3.0366466972156143
+    end
+
     #################################################
     # ∫f(x)g(x)(t-x)^a dx evaluation where f and g in Legendre
     #################################################