From 0af4a9e3e48295b24b0beaa21557f34bb261b482 Mon Sep 17 00:00:00 2001 From: MikaelSlevinsky Date: Tue, 26 Nov 2024 12:12:59 -0600 Subject: [PATCH] add two new constructors for the GramMatrix based on modified OP moments --- docs/Project.toml | 1 + examples/semiclassical.jl | 40 ++++++++++++++------------ src/GramMatrix.jl | 59 ++++++++++++++++++++++++++++++++++++++- 3 files changed, 81 insertions(+), 19 deletions(-) diff --git a/docs/Project.toml b/docs/Project.toml index 969efec..ed41228 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -3,6 +3,7 @@ ApproxFun = "28f2ccd6-bb30-5033-b560-165f7b14dc2f" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" FastTransforms = "057dd010-8810-581a-b7be-e3fc3b93f78c" LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f" +LazyArrays = "5078a376-72f3-5289-bfd5-ec5146d43c02" Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306" PlotlyJS = "f0f68f2c-4968-5e81-91da-67840de0976a" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" diff --git a/examples/semiclassical.jl b/examples/semiclassical.jl index 0b39f7b..c8c26cc 100644 --- a/examples/semiclassical.jl +++ b/examples/semiclassical.jl @@ -4,9 +4,9 @@ # \langle f, g \rangle = \int_{-1}^1 f(x) g(x) w(x){\rm d} x, # ``` # where $w(x) = w^{(\alpha,\beta,\gamma,\delta,\epsilon)}(x) = (1-x)^\alpha(1+x)^\beta(2+x)^\gamma(3+x)^\delta(5-x)^\epsilon$ is a modification of the Jacobi weight. -# We shall use results from [this paper](https://arxiv.org/abs/2302.08448) to consider these semi-classical orthogonal polynomials as modifications of the Jacobi polynomials $P_n^{(\alpha,\beta)}(x)$. +# We shall use results from [this paper](https://arxiv.org/abs/2302.08448) to consider these semi-classical orthogonal polynomials as modifications of the orthonormalized Jacobi polynomials $\tilde{P}_n^{(\alpha,\beta)}(x)$. -using ApproxFun, FastTransforms, LinearAlgebra, Plots, LaTeXStrings +using ApproxFun, FastTransforms, LazyArrays, LinearAlgebra, Plots, LaTeXStrings const GENFIGS = joinpath(pkgdir(FastTransforms), "docs/src/generated") !isdir(GENFIGS) && mkdir(GENFIGS) plotlyjs() @@ -30,7 +30,7 @@ P1 = plan_jac2cheb(Float64, n+1, α, β; normjac = true) q = k -> lmul!(P1, lmul!(P, [zeros(k); 1.0; zeros(n-k)])) qvals = k -> ichebyshevtransform(q(k)) -# With the symmetric Jacobi matrix for $P_n^{(\alpha, \beta)}(x)$ and the modified plan, we may compute the modified Jacobi matrix and the corresponding roots (as eigenvalues): +# With the symmetric Jacobi matrix for $\tilde{P}_n^{(\alpha, \beta)}(x)$ and the modified plan, we may compute the modified Jacobi matrix and the corresponding roots (as eigenvalues): x = Fun(x->x, NormalizedJacobi(β, α)) XP = SymTridiagonal(Symmetric(Multiplication(x, space(x))[1:n+1, 1:n+1])) XQ = FastTransforms.modified_jacobi_matrix(P, XP) @@ -61,18 +61,21 @@ function threshold!(A::AbstractArray, ϵ) A end P′ = plan_modifiedjac2jac(Float64, n+1, α+1, β+1, v.coefficients) -DP = UpperTriangular(diagm(1=>[sqrt(n*(n+α+β+1)) for n in 1:n])) # The classical differentiation matrix representing 𝒟 P^{(α,β)}(y) = P^{(α+1,β+1)}(y) D_P. -DQ = UpperTriangular(threshold!(P′\(DP*(P*I)), 100eps())) # The semi-classical differentiation matrix representing 𝒟 Q(y) = Q̂(y) D_Q. -UpperTriangular(DQ[1:10,1:10]) +DP = UpperTriangular(diagm(1=>[sqrt(n*(n+α+β+1)) for n in 1:n])) # The classical differentiation matrix representing 𝒟 P^{(α,β)}(x) = P^{(α+1,β+1)}(x) D_P. +DQ = UpperTriangular(threshold!(P′\(DP*(P*I)), 100eps())) # The semi-classical differentiation matrix representing 𝒟 Q(x) = Q̂(x) D_Q. +UpperTriangular(DQ[1:10, 1:10]) -# A faster method now exists via the `GramMatrix` architecture and its associated displacement equation. We compute `U`: -U = Symmetric(Multiplication(u, space(u))[1:n+1, 1:n+1]) - -# Then we form a `GramMatrix` together with the Jacobi matrix: -G = GramMatrix(U, XP) +# A faster method now exists via the `GramMatrix` architecture and its associated displacement equation. Given the modified orthogonal polynomial moments implied by the normalized Jacobi series for $u(x)$, we pad this vector to the necessary size and construct the `GramMatrix` with these moments, the multiplication operator, and the constant $\tilde{P}_0^{(\alpha,\beta)}(x)$: +μ = PaddedVector(u.coefficients, 2n+1) +x = Fun(x->x, NormalizedJacobi(β, α)) +XP2 = SymTridiagonal(Symmetric(Multiplication(x, space(x))[1:2n+1, 1:2n+1])) +p0 = Fun(NormalizedJacobi(β, α), [1])(0) +G = GramMatrix(μ, XP2, p0) +G[1:10, 1:10] # And compute its cholesky factorization. The upper-triangular Cholesky factor represents the connection between original Jacobi and semi-classical Jacobi as ${\bf P}^{(\alpha,\beta)}(x) = {\bf Q}(x) R$. R = cholesky(G).U +R[1:10, 1:10] # Every else works almost as before, including evaluation on a Chebyshev grid: q = k -> lmul!(P1, ldiv!(R, [zeros(k); 1.0; zeros(n-k)])) @@ -99,11 +102,12 @@ savefig(joinpath(GENFIGS, "semiclassical1.html")) ###``` # And banded differentiation: -V = Symmetric(Multiplication(v, space(v))[1:n+1, 1:n+1]) -x = Fun(x->x, NormalizedJacobi(β+1, α+1)) -XP′ = SymTridiagonal(Symmetric(Multiplication(x, space(x))[1:n+1, 1:n+1])) -G′ = GramMatrix(V, XP′) +μ′ = PaddedVector(v.coefficients, 2n+1) +x′ = Fun(x->x, NormalizedJacobi(β+1, α+1)) +XP′ = SymTridiagonal(Symmetric(Multiplication(x′, space(x′))[1:2n+1, 1:2n+1])) +p0′ = Fun(NormalizedJacobi(β+1, α+1), [1])(0) +G′ = GramMatrix(μ′, XP′, p0′) R′ = cholesky(G′).U -DP = UpperTriangular(diagm(1=>[sqrt(n*(n+α+β+1)) for n in 1:n])) # The classical differentiation matrix representing 𝒟 P^{(α,β)}(y) = P^{(α+1,β+1)}(y) D_P. -DQ = UpperTriangular(threshold!(R′*(DP*(R\I)), 100eps())) # The semi-classical differentiation matrix representing 𝒟 Q(y) = Q̂(y) D_Q. -UpperTriangular(DQ[1:10,1:10]) +DP = UpperTriangular(diagm(1=>[sqrt(n*(n+α+β+1)) for n in 1:n])) # The classical differentiation matrix representing 𝒟 P^{(α,β)}(x) = P^{(α+1,β+1)}(x) D_P. +DQ = UpperTriangular(threshold!(R′*(DP*(R\I)), 100eps())) # The semi-classical differentiation matrix representing 𝒟 Q(x) = Q̂(x) D_Q. +UpperTriangular(DQ[1:10, 1:10]) diff --git a/src/GramMatrix.jl b/src/GramMatrix.jl index 5321005..e7e8a70 100644 --- a/src/GramMatrix.jl +++ b/src/GramMatrix.jl @@ -50,6 +50,64 @@ GramMatrix(W::WT, X::XT) where {T, WT <: AbstractMatrix{T}, XT <: AbstractMatrix @inline bandwidths(G::GramMatrix) = bandwidths(G.W) @inline MemoryLayout(G::GramMatrix) = MemoryLayout(G.W) +""" + GramMatrix(μ::AbstractVector, X::AbstractMatrix) + +Construct a GramMatrix from modified orthogonal polynomial moments and the multiplication operator. +In the standard (classical) normalization, ``p_0(x) = 1``, so that the moments +``\\mu_n = \\langle p_{n-1}, 1\\rangle`` are in fact the first column of the Gram matrix. +The recurrence is built from ``X^\\top W = WX``. +""" +GramMatrix(μ::AbstractVector{T}, X::XT) where {T, XT <: AbstractMatrix{T}} = GramMatrix(μ, X, one(T)) +function GramMatrix(μ::AbstractVector{T}, X::XT, p0::T) where {T, XT <: AbstractMatrix{T}} + N = length(μ) + n = (N+1)÷2 + @assert N == size(X, 1) == size(X, 2) + @assert bandwidths(X) == (1, 1) + W = Matrix{T}(undef, N, N) + if n > 0 + @inbounds for m in 1:N + W[m, 1] = p0*μ[m] + end + end + if n > 1 + @inbounds for m in 2:N-1 + W[m, 2] = (X[m-1, m]*W[m-1, 1] + (X[m, m]-X[1, 1])*W[m, 1] + X[m+1, m]*W[m+1, 1])/X[2, 1] + end + end + @inbounds @simd for n in 3:n + for m in n:N-n+1 + W[m, n] = (X[m-1, m]*W[m-1, n-1] + (X[m, m]-X[n-1, n-1])*W[m, n-1] + X[m+1, m]*W[m+1, n-1] - X[n-2, n-1]*W[m, n-2])/X[n, n-1] + end + end + return GramMatrix(Symmetric(W[1:n, 1:n], :L), eval(XT.name.name)(view(X, 1:n, 1:n))) +end + +function GramMatrix(μ::PaddedVector{T}, X::XT, p0::T) where {T, XT <: AbstractMatrix{T}} + N = length(μ) + b = length(μ.args[2])-1 + n = (N+1)÷2 + @assert N == size(X, 1) == size(X, 2) + @assert bandwidths(X) == (1, 1) + W = BandedMatrix{T}(undef, (N, N), (b, 0)) + if n > 0 + @inbounds for m in 1:min(N, b+1) + W[m, 1] = p0*μ[m] + end + end + if n > 1 + @inbounds for m in 2:min(N-1, b+2) + W[m, 2] = (X[m-1, m]*W[m-1, 1] + (X[m, m]-X[1, 1])*W[m, 1] + X[m+1, m]*W[m+1, 1])/X[2, 1] + end + end + @inbounds @simd for n in 3:n + for m in n:min(N-n+1, b+n) + W[m, n] = (X[m-1, m]*W[m-1, n-1] + (X[m, m]-X[n-1, n-1])*W[m, n-1] + X[m+1, m]*W[m+1, n-1] - X[n-2, n-1]*W[m, n-2])/X[n, n-1] + end + end + return GramMatrix(Symmetric(W[1:n, 1:n], :L), eval(XT.name.name)(view(X, 1:n, 1:n))) +end + # # X'W-W*X = G*J*G' # This returns G, where J = [0 1; -1 0], respecting the skew-symmetry of the right-hand side. @@ -201,7 +259,6 @@ function compute_skew_generators(W::ChebyshevGramMatrix{T}) where T @inbounds @simd for j in 1:n-1 G[j, 2] = -(μ[n+2-j] + μ[n+j])/2 end - G[n, 2] = -μ[2]/2 G end