From dea45045222c0b92976bb094fb150c7252fd9027 Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Tue, 25 Feb 2025 16:52:34 +0000 Subject: [PATCH 1/3] Add orthogonal connection matrices --- test/test_ultraspherical.jl | 20 +++++++++++++++++++- 1 file changed, 19 insertions(+), 1 deletion(-) diff --git a/test/test_ultraspherical.jl b/test/test_ultraspherical.jl index b58a7bec..d368c395 100644 --- a/test/test_ultraspherical.jl +++ b/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,22 @@ 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(3/2)) + @test P'P == W'W == I + P\W + + Q = Normalized(Legendre()) \ (JacobiWeight(1,1) .* Normalized(Jacobi(2,2))) + @test Q[1:12,1:10]'Q[1:12,1:10] ≈ I + + c = sqrt.(2*(5:2:∞) ./ ((3:∞) .* (4:∞) )) + s = sqrt.(((1:∞) .* (2:∞)) ./ ((3:∞) .* (4:∞) )) + + @test (c.^2 + s.^2)[1:10] ≈ ones(10) + G₁ + + + end end \ No newline at end of file From 9b7a77b611efa760fbc2eb3202409baf554e2146 Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Tue, 25 Feb 2025 19:51:49 +0000 Subject: [PATCH 2/3] Update test_ultraspherical.jl --- test/test_ultraspherical.jl | 18 +++++++++++++++--- 1 file changed, 15 insertions(+), 3 deletions(-) diff --git a/test/test_ultraspherical.jl b/test/test_ultraspherical.jl index d368c395..1a7b6f86 100644 --- a/test/test_ultraspherical.jl +++ b/test/test_ultraspherical.jl @@ -214,12 +214,24 @@ using ClassicalOrthogonalPolynomials: grammatrix, OrthonormalWeighted Q = Normalized(Legendre()) \ (JacobiWeight(1,1) .* Normalized(Jacobi(2,2))) @test Q[1:12,1:10]'Q[1:12,1:10] ≈ I + X = jacobimatrix(Normalized(Legendre())) + @test Q[1:10,1:10] ≈ -qr(I - X^2).Q[1:10,1:10] + c = sqrt.(2*(5:2:∞) ./ ((3:∞) .* (4:∞) )) s = sqrt.(((1:∞) .* (2:∞)) ./ ((3:∞) .* (4:∞) )) - @test (c.^2 + s.^2)[1:10] ≈ ones(10) - G₁ + 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] - + MatrixFactorizations.QRPackedQ(BandedMatrix(-2 => -(s ./ (c .+ 1))), c .+ 1)[1:10,1:10] end end \ No newline at end of file From 6ca0f88a0e26050ec3afcd4ea330cac17844b48f Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Mon, 17 Mar 2025 16:13:46 +0000 Subject: [PATCH 3/3] Orthogonal connection matrix for Ultraspherical --- Project.toml | 3 ++- src/ClassicalOrthogonalPolynomials.jl | 2 ++ src/classical/jacobi.jl | 12 ++++----- src/classical/ultraspherical.jl | 32 ++++++++++++++++++---- test/test_ultraspherical.jl | 38 ++++++++++++++------------- 5 files changed, 57 insertions(+), 30 deletions(-) diff --git a/Project.toml b/Project.toml index aa709430..83afb2a6 100644 --- a/Project.toml +++ b/Project.toml @@ -3,7 +3,6 @@ uuid = "b30e2e7b-c4ee-47da-9d5f-2c5c27239acd" authors = ["Sheehan Olver "] 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" diff --git a/src/ClassicalOrthogonalPolynomials.jl b/src/ClassicalOrthogonalPolynomials.jl index c2c5d590..2ae883a3 100644 --- a/src/ClassicalOrthogonalPolynomials.jl +++ b/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, diff --git a/src/classical/jacobi.jl b/src/classical/jacobi.jl index 70673124..9406d291 100644 --- a/src/classical/jacobi.jl +++ b/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 ``` """ diff --git a/src/classical/ultraspherical.jl b/src/classical/ultraspherical.jl index 7b75f0ee..cecf7edc 100644 --- a/src/classical/ultraspherical.jl +++ b/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 \ No newline at end of file diff --git a/test/test_ultraspherical.jl b/test/test_ultraspherical.jl index 1a7b6f86..75d57e68 100644 --- a/test/test_ultraspherical.jl +++ b/test/test_ultraspherical.jl @@ -207,31 +207,33 @@ using ClassicalOrthogonalPolynomials: grammatrix, OrthonormalWeighted @testset "orthonormal" begin P = OrthonormalWeighted(Ultraspherical(1/2)) - W = OrthonormalWeighted(Ultraspherical(3/2)) + W = OrthonormalWeighted(Ultraspherical(5/2)) @test P'P == W'W == I - P\W + @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 - Q = Normalized(Legendre()) \ (JacobiWeight(1,1) .* Normalized(Jacobi(2,2))) + @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(I - X^2).Q[1:10,1:10] + @test Q[1:10,1:10] ≈ -qr(X^2 - I).Q[1:10,1:10] - c = sqrt.(2*(5:2:∞) ./ ((3:∞) .* (4:∞) )) - s = sqrt.(((1:∞) .* (2:∞)) ./ ((3:∞) .* (4:∞) )) + @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] + 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 - MatrixFactorizations.QRPackedQ(BandedMatrix(-2 => -(s ./ (c .+ 1))), c .+ 1)[1:10,1:10] + @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 \ No newline at end of file