-
Notifications
You must be signed in to change notification settings - Fork 27
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Merge branch 'master' into pochhammervec
- Loading branch information
Showing
19 changed files
with
1,565 additions
and
513 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,77 @@ | ||
| # # Half-range Chebyshev polynomials | ||
| # In [this paper](https://doi.org/10.1137/090752456), [Daan Huybrechs](https://github.com/daanhb) introduced the so-called half-range Chebyshev polynomials | ||
| # as the semi-classical orthogonal polynomials with respect to the inner product: | ||
| # ```math | ||
| # \langle f, g \rangle = \int_0^1 f(x) g(x)\frac{{\rm d} x}{\sqrt{1-x^2}}. | ||
| # ``` | ||
| # By the variable transformation $y = 2x-1$, the resulting polynomials can be related to | ||
| # orthogonal polynomials on $(-1,1)$ with the Jacobi weight $(1-y)^{-\frac{1}{2}}$ modified by the weight $(3+y)^{-\frac{1}{2}}$. | ||
| # | ||
| # We shall use the fact that: | ||
| # ```math | ||
| # \frac{1}{\sqrt{3+y}} = \sqrt{\frac{2}{3+\sqrt{8}}}\sum_{n=0}^\infty P_n(y) \left(\frac{-1}{3+\sqrt{8}}\right)^n, | ||
| # ``` | ||
| # and results from [this paper](https://arxiv.org/abs/2302.08448) to consider the half-range Chebyshev polynomials as | ||
| # modifications of the Jacobi polynomials $P_n^{(-\frac{1}{2},0)}(y)$. | ||
|
|
||
| using FastTransforms, LinearAlgebra, Plots, LaTeXStrings | ||
| const GENFIGS = joinpath(pkgdir(FastTransforms), "docs/src/generated") | ||
| !isdir(GENFIGS) && mkdir(GENFIGS) | ||
| plotlyjs() | ||
|
|
||
| # We truncate the generating function to ensure a relative error less than `eps()` in the uniform norm on $(-1,1)$: | ||
| z = -1/(3+sqrt(8)) | ||
| K = sqrt(-2z) | ||
| N = ceil(Int, log(abs(z), eps()/2*(1-abs(z))/K) - 1) | ||
| d = K .* z .^(0:N) | ||
|
|
||
| # Then, we convert this representation to the expansion in Jacobi polynomials $P_n^{(-\frac{1}{2}, 0)}(y)$: | ||
| u = jac2jac(d, 0.0, 0.0, -0.5, 0.0; norm1 = false, norm2 = true) | ||
|
|
||
| # Our working polynomial degree will be: | ||
| n = 100 | ||
|
|
||
| # We compute the connection coefficients between the modified orthogonal polynomials and the Jacobi polynomials: | ||
| P = plan_modifiedjac2jac(Float64, n+1, -0.5, 0.0, u) | ||
|
|
||
| # We store the connection to first kind Chebyshev polynomials: | ||
| P1 = plan_jac2cheb(Float64, n+1, -0.5, 0.0; normjac = true) | ||
|
|
||
| # We compute the Chebyshev series for the degree-$k\le n$ modified polynomial and its values at the Chebyshev points: | ||
| 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^{(-\frac{1}{2}, 0)}(y)$ and the modified plan, we may compute the modified Jacobi matrix and the corresponding roots (as eigenvalues): | ||
| XP = SymTridiagonal([-inv((4n-1)*(4n-5)) for n in 1:n+1], [4n*(2n-1)/(4n-1)/sqrt((4n-3)*(4n+1)) for n in 1:n]) | ||
| XQ = FastTransforms.modified_jacobi_matrix(P, XP) | ||
| SymTridiagonal(XQ.dv[1:10], XQ.ev[1:9]) | ||
|
|
||
| # And we plot: | ||
| x = (chebyshevpoints(Float64, n+1, Val(1)) .+ 1 ) ./ 2 | ||
| p = plot(x, qvals(0); linewidth=2.0, legend = false, xlim=(0,1), xlabel=L"x", | ||
| ylabel=L"T^h_n(x)", title="Half-Range Chebyshev Polynomials and Their Roots", | ||
| extra_plot_kwargs = KW(:include_mathjax => "cdn")) | ||
| for k in 1:10 | ||
| λ = (eigvals(SymTridiagonal(XQ.dv[1:k], XQ.ev[1:k-1])) .+ 1) ./ 2 | ||
| plot!(x, qvals(k); linewidth=2.0, color=palette(:default)[k+1]) | ||
| scatter!(λ, zero(λ); markersize=2.5, color=palette(:default)[k+1]) | ||
| end | ||
| p | ||
| savefig(joinpath(GENFIGS, "halfrange.html")) | ||
| ###```@raw html | ||
| ###<object type="text/html" data="../halfrange.html" style="width:100%;height:400px;"></object> | ||
| ###``` | ||
|
|
||
| # By [Theorem 2.20](https://arxiv.org/abs/2302.08448) it turns out that the *derivatives* of the half-range Chebyshev polynomials are a linear combination of at most two polynomials orthogonal with respect to $\sqrt{(3+y)(1-y)}(1+y)$ on $(-1,1)$. This fact enables us to compute the banded differentiation matrix: | ||
| v̂ = 3*[u; 0]+XP[1:N+2, 1:N+1]*u | ||
| v = jac2jac(v̂, -0.5, 0.0, 0.5, 1.0; norm1 = true, norm2 = true) | ||
| function threshold!(A::AbstractArray, ϵ) | ||
| for i in eachindex(A) | ||
| if abs(A[i]) < ϵ A[i] = 0 end | ||
| end | ||
| A | ||
| end | ||
| P′ = plan_modifiedjac2jac(Float64, n+1, 0.5, 1.0, v) | ||
| DP = UpperTriangular(diagm(1=>[sqrt(n*(n+1/2)) for n in 1:n])) # The classical differentiation matrix representing 𝒟 P^{(-1/2,0)}(y) = P^{(1/2,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]) |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,135 @@ | ||
| struct SymmetricToeplitzPlusHankel{T} <: AbstractMatrix{T} | ||
| v::Vector{T} | ||
| n::Int | ||
| end | ||
|
|
||
| function SymmetricToeplitzPlusHankel(v::Vector{T}) where T | ||
| n = (length(v)+1)÷2 | ||
| SymmetricToeplitzPlusHankel{T}(v, n) | ||
| end | ||
|
|
||
| size(A::SymmetricToeplitzPlusHankel{T}) where T = (A.n, A.n) | ||
| getindex(A::SymmetricToeplitzPlusHankel{T}, i::Integer, j::Integer) where T = A.v[abs(i-j)+1] + A.v[i+j-1] | ||
|
|
||
| struct SymmetricBandedToeplitzPlusHankel{T} <: BandedMatrices.AbstractBandedMatrix{T} | ||
| v::Vector{T} | ||
| n::Int | ||
| b::Int | ||
| end | ||
|
|
||
| function SymmetricBandedToeplitzPlusHankel(v::Vector{T}, n::Integer) where T | ||
| SymmetricBandedToeplitzPlusHankel{T}(v, n, length(v)-1) | ||
| end | ||
|
|
||
| size(A::SymmetricBandedToeplitzPlusHankel{T}) where T = (A.n, A.n) | ||
| function getindex(A::SymmetricBandedToeplitzPlusHankel{T}, i::Integer, j::Integer) where T | ||
| v = A.v | ||
| if abs(i-j) < length(v) | ||
| if i+j-1 ≤ length(v) | ||
| v[abs(i-j)+1] + v[i+j-1] | ||
| else | ||
| v[abs(i-j)+1] | ||
| end | ||
| else | ||
| zero(T) | ||
| end | ||
| end | ||
| bandwidths(A::SymmetricBandedToeplitzPlusHankel{T}) where T = (A.b, A.b) | ||
|
|
||
| # | ||
| # Jac*W-W*Jac' = G*J*G' | ||
| # This returns G and J, where J = [0 I; -I 0], respecting the skew-symmetry of the right-hand side. | ||
| # | ||
| function compute_skew_generators(A::SymmetricToeplitzPlusHankel{T}) where T | ||
| v = A.v | ||
| n = size(A, 1) | ||
| J = [zero(T) one(T); -one(T) zero(T)] | ||
| G = zeros(T, n, 2) | ||
| G[n, 1] = one(T) | ||
| u2 = reverse(v[2:n+1]) | ||
| u2[1:n-1] .+= v[n+1:2n-1] | ||
| G[:, 2] .= -u2 | ||
| G, J | ||
| end | ||
|
|
||
| function cholesky(A::SymmetricToeplitzPlusHankel{T}) where T | ||
| n = size(A, 1) | ||
| G, J = compute_skew_generators(A) | ||
| L = zeros(T, n, n) | ||
| r = A[:, 1] | ||
| r2 = zeros(T, n) | ||
| l = zeros(T, n) | ||
| v = zeros(T, n) | ||
| col1 = zeros(T, n) | ||
| STPHcholesky!(L, G, r, r2, l, v, col1, n) | ||
| return Cholesky(L, 'L', 0) | ||
| end | ||
|
|
||
| function STPHcholesky!(L::Matrix{T}, G, r, r2, l, v, col1, n) where T | ||
| @inbounds @simd for k in 1:n-1 | ||
| x = sqrt(r[1]) | ||
| for j in 1:n-k+1 | ||
| L[j+k-1, k] = l[j] = r[j]/x | ||
| end | ||
| for j in 1:n-k+1 | ||
| v[j] = G[j, 1]*G[1,2]-G[j,2]*G[1,1] | ||
| end | ||
| F12 = k == 1 ? T(2) : T(1) | ||
| r2[1] = (r[2] - v[1])/F12 | ||
| for j in 2:n-k | ||
| r2[j] = (r[j+1]+r[j-1] + r[1]*col1[j] - col1[1]*r[j] - v[j])/F12 | ||
| end | ||
| r2[n-k+1] = (r[n-k] + r[1]*col1[n-k+1] - col1[1]*r[n-k+1] - v[n-k+1])/F12 | ||
| cst = r[2]/x | ||
| for j in 1:n-k | ||
| r[j] = r2[j+1] - cst*l[j+1] | ||
| end | ||
| for j in 1:n-k | ||
| col1[j] = -F12/x*l[j+1] | ||
| end | ||
| c1 = G[1, 1] | ||
| c2 = G[1, 2] | ||
| for j in 1:n-k | ||
| G[j, 1] = G[j+1, 1] - l[j+1]*c1/x | ||
| G[j, 2] = G[j+1, 2] - l[j+1]*c2/x | ||
| end | ||
| end | ||
| L[n, n] = sqrt(r[1]) | ||
| end | ||
|
|
||
| function cholesky(A::SymmetricBandedToeplitzPlusHankel{T}) where T | ||
| n = size(A, 1) | ||
| b = A.b | ||
| R = BandedMatrix{T}(undef, (n, n), (0, bandwidth(A, 2))) | ||
| r = A[1:b+2, 1] | ||
| r2 = zeros(T, b+3) | ||
| l = zeros(T, b+3) | ||
| col1 = zeros(T, b+2) | ||
| SBTPHcholesky!(R, r, r2, l, col1, n, b) | ||
| return Cholesky(R, 'U', 0) | ||
| end | ||
|
|
||
| function SBTPHcholesky!(R::BandedMatrix{T}, r, r2, l, col1, n, b) where T | ||
| @inbounds @simd for k in 1:n | ||
| x = sqrt(r[1]) | ||
| for j in 1:b+1 | ||
| l[j] = r[j]/x | ||
| end | ||
| for j in 1:min(n-k+1, b+1) | ||
| R[k, j+k-1] = l[j] | ||
| end | ||
| F12 = k == 1 ? T(2) : T(1) | ||
| r2[1] = r[2]/F12 | ||
| for j in 2:b+1 | ||
| r2[j] = (r[j+1]+r[j-1] + r[1]*col1[j] - col1[1]*r[j])/F12 | ||
| end | ||
| r2[b+2] = (r[b+1] + r[1]*col1[b+2] - col1[1]*r[b+2])/F12 | ||
| cst = r[2]/x | ||
| for j in 1:b+2 | ||
| r[j] = r2[j+1] - cst*l[j+1] | ||
| end | ||
| for j in 1:b+2 | ||
| col1[j] = -F12/x*l[j+1] | ||
| end | ||
| end | ||
| end |
Oops, something went wrong.