Add orthogonal connection matrices

Project.toml

@@ -3,7 +3,6 @@ uuid = "b30e2e7b-c4ee-47da-9d5f-2c5c27239acd"
 authors = ["Sheehan Olver <[email protected]>"]
 version = "0.15.1"
 
-
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
 BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
@@ -22,6 +21,7 @@ IntervalSets = "8197267c-284f-5f27-9208-e0e47529a953"
 LazyArrays = "5078a376-72f3-5289-bfd5-ec5146d43c02"
 LazyBandedMatrices = "d7e5e226-e90b-4449-9968-0f923699bf6f"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+MatrixFactorizations = "a3b82374-2e81-5b9e-98ce-41277c0e4c87"
 QuasiArrays = "c4ea9172-b204-11e9-377d-29865faadc5c"
 RecurrenceRelationshipArrays = "b889d2dc-af3c-4820-88a8-238fa91d3518"
 RecurrenceRelationships = "807425ed-42ea-44d6-a357-6771516d7b2c"
@@ -51,6 +51,7 @@ InfiniteLinearAlgebra = "0.10"
 IntervalSets = "0.7"
 LazyArrays = "2.5.1"
 LazyBandedMatrices = "0.11"
+MatrixFactorizations = "3.1"
 MutableArithmetics = "1"
 QuasiArrays = "0.12"
 RecurrenceRelationshipArrays = "0.1.2"

src/ClassicalOrthogonalPolynomials.jl

@@ -48,6 +48,8 @@ import FastGaussQuadrature: jacobimoment
 import BlockArrays: blockedrange, _BlockedUnitRange, unblock, _BlockArray, block, blockindex, BlockSlice, blockvec
 import BandedMatrices: bandwidths
 
+using MatrixFactorizations: QRPackedQMatrix
+
 export OrthogonalPolynomial, Normalized, LanczosPolynomial,
             Hermite, Jacobi, Legendre, Chebyshev, ChebyshevT, ChebyshevU, ChebyshevInterval, Ultraspherical, Fourier, Laurent, Laguerre,
             HermiteWeight, JacobiWeight, ChebyshevWeight, ChebyshevGrid, ChebyshevTWeight, ChebyshevUWeight, UltrasphericalWeight, LegendreWeight, LaguerreWeight,

src/classical/jacobi.jl

@@ -28,13 +28,13 @@ broadcasted(::LazyQuasiArrayStyle{1}, ::typeof(Base.literal_pow), ::Base.RefValu
 The quasi-vector representing the Jacobi weight function ``(1-x)^a (1+x)^b`` on ``[-1,1]``. See also [`jacobiweight`](@ref) and [`Jacobi`](@ref).
 # Examples
 ```jldoctest
-julia> J=JacobiWeight(1.0,1.0)
+julia> w = JacobiWeight(1.0,1.0)
 (1-x)^1.0 * (1+x)^1.0 on -1..1
 
-julia> J[0.5]
+julia> w[0.5]
 0.75
 
-julia> axes(J)
+julia> axes(w)
 (Inclusion(-1.0 .. 1.0 (Chebyshev)),)
 ```
 """
@@ -56,13 +56,13 @@ The [`JacobiWeight`](@ref) affine-mapped to interval `d`.
 
 # Examples
 ```jldoctest
-julia> J = jacobiweight(1, 1, 0..1)
+julia> w = jacobiweight(1, 1, 0..1)
 (1-x)^1 * (1+x)^1 on -1..1 affine mapped to 0 .. 1
 
-julia> axes(J)
+julia> axes(w)
 (Inclusion(0 .. 1),)
 
-julia> J[0.5]
+julia> w[0.5]
 1.0
 ```
 """

src/classical/ultraspherical.jl

@@ -1,10 +1,23 @@
 
 
 
-##
-# Ultraspherical
-##
-
+"""
+    UltrasphericalWeight{T}(a,b)
+    UltrasphericalWeight(a,b)
+
+The quasi-vector representing the Ultraspherical weight function ``(1-x^2)^{λ-1/2}`` on ``[-1,1]``. See also [`Ultraspherical`](@ref).
+# Examples
+```jldoctest
+julia> w = UltrasphericalWeight(1.0)
+(1-x^2)^0.5 on -1..1
+
+julia> w[0.5]
+0.8660254037844386
+
+julia> axes(w)
+(Inclusion(-1.0 .. 1.0 (Chebyshev)),)
+```
+"""
 struct UltrasphericalWeight{T,Λ} <: AbstractJacobiWeight{T}
     λ::Λ
 end
@@ -17,7 +30,7 @@ AbstractQuasiArray{T}(w::UltrasphericalWeight) where T = UltrasphericalWeight{T}
 AbstractQuasiVector{T}(w::UltrasphericalWeight) where T = UltrasphericalWeight{T}(w.λ)
 
 show(io::IO, w::UltrasphericalWeight) = summary(io, w)
-summary(io::IO, w::UltrasphericalWeight) = print(io, "UltrasphericalWeight($(w.λ))")
+summary(io::IO, w::UltrasphericalWeight) = print(io, "(1-x^2)^$(w.λ-1/2) on -1..1")
 
 ==(a::UltrasphericalWeight, b::UltrasphericalWeight) = a.λ == b.λ
 
@@ -30,6 +43,8 @@ sum(w::UltrasphericalWeight{T}) where T = sqrt(convert(T,π))*exp(loggamma(one(T
 
 hasboundedendpoints(w::UltrasphericalWeight) = 2w.λ ≥ 1
 
+broadcasted(::LazyQuasiArrayStyle{1}, ::typeof(sqrt), w::UltrasphericalWeight) = UltrasphericalWeight(w.λ/2 + one(w.λ)/4)
+
 
 struct Ultraspherical{T,Λ} <: AbstractJacobi{T}
     λ::Λ
@@ -308,3 +323,10 @@ end
 \(A::Ultraspherical, w_B::WeightedUltraspherical) = (UltrasphericalWeight(one(A.λ)/2) .* A) \ w_B
 \(A::Legendre, wB::WeightedUltraspherical) = Ultraspherical(A) \ wB
 
+function \(P::OrthonormalWeighted{<:Any,<:Ultraspherical}, Q::OrthonormalWeighted{<:Any,<:Ultraspherical})
+    @assert isone(2P.P.P.λ) # only Legendre is implemented
+    @assert 2Q.P.P.λ == 5
+    c = sqrt.(2*(5:2:∞) ./ ((3:∞) .* (4:∞) ))
+    s = sqrt.(((1:∞) .* (2:∞)) ./ ((3:∞) .* (4:∞) ))
+    -QRPackedQMatrix(BandedMatrix(-2 => -(s ./ (c .+ 1))), c .+ 1)
+end

test/test_ultraspherical.jl

@@ -1,7 +1,7 @@
 using ClassicalOrthogonalPolynomials, ContinuumArrays, BandedMatrices, LazyArrays, Test
 using ForwardDiff
 using LazyArrays: rowsupport, colsupport
-using ClassicalOrthogonalPolynomials: grammatrix
+using ClassicalOrthogonalPolynomials: grammatrix, OrthonormalWeighted
 
 @testset "Ultraspherical" begin
     @testset "Conversion" begin
@@ -204,4 +204,36 @@ using ClassicalOrthogonalPolynomials: grammatrix
         @test (C \ diff(U,1))[1:10,1:10] == (C \ diff(U))[1:10,1:10]
         @test (C³ \ diff(U,2))[1:10,1:10] == (C³ \ diff(diff(U)))[1:10,1:10]
     end
+
+    @testset "orthonormal" begin
+        P = OrthonormalWeighted(Ultraspherical(1/2))
+        W = OrthonormalWeighted(Ultraspherical(5/2))
+        @test P'P == W'W == I
+        @test P[0.1,1:10] ≈ Normalized(Legendre())[0.1,1:10]
+        @test W[0.1,1:10] ≈ (1 - 0.1.^2) * Normalized(Jacobi(2,2))[0.1,1:10]
+        Q = P\W
+
+        @test (Normalized(Legendre()) \ (JacobiWeight(1,1) .* Normalized(Jacobi(2,2))))[1:10,1:10] ≈ Q[1:10,1:10]
+        @test Q[1:12,1:10]'Q[1:12,1:10] ≈ I
+
+        X = jacobimatrix(Normalized(Legendre()))
+        @test Q[1:10,1:10] ≈ -qr(X^2 - I).Q[1:10,1:10]
+
+        @testset "derivation" begin
+            c = sqrt.(2*(5:2:∞) ./ ((3:∞) .* (4:∞) ))
+            s = sqrt.(((1:∞) .* (2:∞)) ./ ((3:∞) .* (4:∞) ))
+
+            n = 10
+            @test (c.^2 + s.^2)[1:n] ≈ ones(n)
+            G = [Matrix(1.0I,n,n) for k=1:n-2]
+            for k = 1:n-2
+                G[k][[k,k+2],[k,k+2]] = [c[k] s[k]; -s[k] c[k]] 
+            end
+
+            @test Q[1:n,1:n-2] ≈  *(G...)[:,1:n-2]
+
+            @test qr(I - X^2).τ[1:10] ≈ c[1:10] .+ 1
+            @test qr(I - X^2).factors[band(-2)][1:10] ≈ -(s ./ (c .+ 1))[1:10]
+        end
+    end
 end