|
| 1 | +# # Half-range Chebyshev polynomials |
| 2 | +# In [this paper](https://doi.org/10.1137/090752456), [Daan Huybrechs](https://github.com/daanhb) introduced the so-called half-range Chebyshev polynomials |
| 3 | +# as the non-classical orthogonal polynomials with respect to the inner product: |
| 4 | +# ```math |
| 5 | +# \langle f, g \rangle = \int_0^1 f(x) g(x)\frac{{\rm d} x}{\sqrt{1-x^2}}. |
| 6 | +# ``` |
| 7 | +# By the variable transformation $y = 2x-1$, the resulting polynomials can be related to |
| 8 | +# orthogonal polynomials on $(-1,1)$ with the Jacobi weight $(1-y)^{-\frac{1}{2}}$ modified by the weight $(3+y)^{-\frac{1}{2}}$. |
| 9 | +# |
| 10 | +# We shall use the fact that: |
| 11 | +# ```math |
| 12 | +# \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, |
| 13 | +# ``` |
| 14 | +# and results from [this paper](https://arxiv.org/abs/2302.08448) to consider the half-range Chebyshev polynomials as |
| 15 | +# modifications of the Jacobi polynomials $P_n^{(-\frac{1}{2},0)}(y)$. |
| 16 | + |
| 17 | +using FastTransforms, LinearAlgebra, Plots, LaTeXStrings |
| 18 | +const GENFIGS = joinpath(pkgdir(FastTransforms), "docs/src/generated") |
| 19 | +!isdir(GENFIGS) && mkdir(GENFIGS) |
| 20 | +plotlyjs() |
| 21 | + |
| 22 | +# We truncate the generating function to ensure an absolute error of `eps()`: |
| 23 | +z = -1/(3+sqrt(8)) |
| 24 | +K = sqrt(-2z) |
| 25 | +N = log(abs(z), eps()*(1-abs(z))/K) - 1 |
| 26 | +d = K .* z .^(0:N) |
| 27 | + |
| 28 | +# Then, we convert this representation to the expansion in Jacobi polynomials $P_n^{(-\frac{1}{2}, 0)}(y)$: |
| 29 | +u = jac2jac(d, 0.0, 0.0, -0.5, 0.0; norm1 = false, norm2 = true) |
| 30 | + |
| 31 | +# Our working polynomial degree will be: |
| 32 | +n = 100 |
| 33 | + |
| 34 | +# We compute the connection coefficients between the modified orthogonal polynomials and the Jacobi polynomials: |
| 35 | +P = plan_modifiedjac2jac(Float64, n+1, -0.5, 0.0, u) |
| 36 | + |
| 37 | +# We store the connection to first kind Chebyshev polynomials: |
| 38 | +P1 = plan_jac2cheb(Float64, n+1, -0.5, 0.0; normjac = true) |
| 39 | + |
| 40 | +# We compute the Chebyshev series for the degree-$k\le n$ modified polynomial and its values at the Chebyshev points: |
| 41 | +q = k -> lmul!(P1, lmul!(P, [zeros(k); 1.0; zeros(n-k)])) |
| 42 | +qvals = k-> ichebyshevtransform(q(k)) |
| 43 | + |
| 44 | +# 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): |
| 45 | +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]) |
| 46 | +XQ = FastTransforms.modified_jacobi_matrix(P, XP) |
| 47 | + |
| 48 | +# And we plot: |
| 49 | +x = (chebyshevpoints(Float64, n+1, Val(1)) .+ 1 ) ./ 2 |
| 50 | +p = plot(x, qvals(0); linewidth=2.0, legend = false, xlim=(0,1), xlabel=L"x", |
| 51 | + ylabel=L"T^h_n(x)", title="Half-Range Chebyshev Polynomials and Their Roots", |
| 52 | + extra_plot_kwargs = KW(:include_mathjax => "cdn")) |
| 53 | +for k in 1:10 |
| 54 | + λ = (eigvals(SymTridiagonal(XQ.dv[1:k], XQ.ev[1:k-1])) .+ 1) ./ 2 |
| 55 | + plot!(x, qvals(k); linewidth=2.0, color=palette(:default)[k+1]) |
| 56 | + scatter!(λ, zero(λ); markersize=2.5, color=palette(:default)[k+1]) |
| 57 | +end |
| 58 | +p |
| 59 | +#savefig(joinpath(GENFIGS, "halfrange.html")) |
| 60 | +###```@raw html |
| 61 | +###<object type="text/html" data="../halfrange.html" style="width:100%;height:400px;"></object> |
| 62 | +###``` |
0 commit comments