Skip to content

Commit 78c98fb

Browse files
include banded diff
1 parent 04d8e87 commit 78c98fb

File tree

1 file changed

+16
-3
lines changed

1 file changed

+16
-3
lines changed

examples/halfrange.jl

Lines changed: 16 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
# # Half-range Chebyshev polynomials
22
# 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:
3+
# as the semi-classical orthogonal polynomials with respect to the inner product:
44
# ```math
55
# \langle f, g \rangle = \int_0^1 f(x) g(x)\frac{{\rm d} x}{\sqrt{1-x^2}}.
66
# ```
@@ -19,10 +19,10 @@ const GENFIGS = joinpath(pkgdir(FastTransforms), "docs/src/generated")
1919
!isdir(GENFIGS) && mkdir(GENFIGS)
2020
plotlyjs()
2121

22-
# We truncate the generating function to ensure an absolute error of `eps()`:
22+
# We truncate the generating function to ensure a relative error less than `eps()` in the uniform norm on $(-1,1)$:
2323
z = -1/(3+sqrt(8))
2424
K = sqrt(-2z)
25-
N = log(abs(z), eps()*(1-abs(z))/K) - 1
25+
N = ceil(Int, log(abs(z), eps()/2*(1-abs(z))/K) - 1)
2626
d = K .* z .^(0:N)
2727

2828
# Then, we convert this representation to the expansion in Jacobi polynomials $P_n^{(-\frac{1}{2}, 0)}(y)$:
@@ -60,3 +60,16 @@ savefig(joinpath(GENFIGS, "halfrange.html"))
6060
###```@raw html
6161
###<object type="text/html" data="../halfrange.html" style="width:100%;height:400px;"></object>
6262
###```
63+
64+
# 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:
65+
= 3*u+XP[1:N+1,1:N]*u
66+
v = jac2jac(v̂, -0.5, 0.0, 0.5, 1.0; norm1 = true, norm2 = true)
67+
function threshold!(A::AbstractArray, ϵ)
68+
for i in eachindex(A)
69+
if abs(A[i]) < ϵ A[i] = 0 end
70+
end
71+
A
72+
end
73+
P′ = plan_modifiedjac2jac(Float64, n+1, 0.5, 1.0, v)
74+
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.
75+
DQ = UpperTriangular(threshold!(P′\(DP*(P*I)), 100eps())) # The semi-classical differentiation matrix representing 𝒟 Q(y) = Q̂(y) D_Q.

0 commit comments

Comments
 (0)