Skip to content

Commit 9a5fce8

Browse files
committed
expand orthogonal polynomials
1 parent c455ece commit 9a5fce8

File tree

13 files changed

+237
-153
lines changed

13 files changed

+237
-153
lines changed

Project.toml

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -49,7 +49,7 @@ FillArrays = "0.13,1"
4949
GaussQuadrature = "0.5.7"
5050
GenericFFT = "0.1.3"
5151
GenericLinearAlgebra = "0.3"
52-
GridArrays = "0.2.3"
52+
GridArrays = "0.2.4"
5353
IterativeSolvers = "0.9"
5454
LinearAlgebra = "1"
5555
MacroTools = "0.5"

src/BasisFunctions.jl

Lines changed: 6 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -92,7 +92,10 @@ import DomainIntegrals:
9292
isnormalized, isuniform,
9393
ismappedmeasure,
9494
jacobi_α, jacobi_β, laguerre_α
95-
using DomainIntegrals: ProductWeight, AbstractJacobiWeight
95+
using DomainIntegrals:
96+
ProductWeight,
97+
AbstractJacobiWeight,
98+
UltrasphericalWeight
9699
export integral
97100

98101
@reexport using GridArrays
@@ -295,7 +298,7 @@ export ChebyshevT, ChebyshevU,
295298
FastChebyshevTransformFFTW, InverseFastChebyshevTransformFFTW
296299

297300
# from bases/poly/orthopoly.jl and friends
298-
export Legendre, Jacobi, Laguerre, Hermite,
301+
export Legendre, Jacobi, Laguerre, Hermite, Ultraspherical,
299302
Monomials, GenericOPS,
300303
recurrence_eval, recurrence_eval_derivative, monic_recurrence_eval,
301304
monic_recurrence_coefficients,
@@ -406,6 +409,7 @@ include("bases/fourier/trig.jl")
406409
include("bases/poly/polynomials.jl")
407410
include("bases/poly/monomials.jl")
408411
include("bases/poly/quadrature.jl")
412+
include("bases/poly/ops/dlmf.jl")
409413
include("bases/poly/ops/orthopoly.jl")
410414
include("bases/poly/chebyshev/chebyshev.jl")
411415
include("bases/poly/ops/legendre.jl")

src/bases/dictionary/gram.jl

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -21,14 +21,15 @@ dict_innerproduct(Φ1::Dictionary, i::Int, Φ2::Dictionary, j, measure; options.
2121
dict_innerproduct(Φ1::Dictionary, i, Φ2::Dictionary, j, measure; options...) =
2222
dict_innerproduct_native(Φ1, i, Φ2, j, measure; options...)
2323

24-
# - innerproduct_native: if not specialized, called innerproduct1
24+
# - innerproduct_native: if not specialized, calls innerproduct1
2525
dict_innerproduct_native(Φ1::Dictionary, i, Φ2::Dictionary, j, measure; options...) =
2626
dict_innerproduct1(Φ1, i, Φ2, j, measure; options...)
2727
# - innerproduct1: possibility to dispatch on the first dictionary without ambiguity.
2828
# If not specialized, we call innerproduct2
2929
dict_innerproduct1(Φ1::Dictionary, i, Φ2, j, measure; options...) =
3030
dict_innerproduct2(Φ1, i, Φ2, j, measure; options...)
3131
# - innerproduct2: possibility to dispatch on the second dictionary without ambiguity
32+
# the warning might indicate that an expensive numerical integration takes place
3233
function dict_innerproduct2(Φ1, i, Φ2::Dictionary, j, measure; verbose=false, options...)
3334
verbose && println("Invoking default inner product")
3435
default_dict_innerproduct(Φ1, i, Φ2, j, measure; options...)

src/bases/poly/chebyshev/ChebyshevT.jl

Lines changed: 12 additions & 22 deletions
Original file line numberDiff line numberDiff line change
@@ -62,16 +62,13 @@ function transform_dict(s::ChebyshevT{T}; chebyshevpoints=:nodes, options...) wh
6262
end
6363

6464
# Parameters alpha and beta of the corresponding Jacobi polynomial
65-
jacobi_α(b::ChebyshevT{T}) where {T} = -one(T)/2
66-
jacobi_β(b::ChebyshevT{T}) where {T} = -one(T)/2
65+
jacobi_α(b::ChebyshevT{T}) where T = -one(T)/2
66+
jacobi_β(b::ChebyshevT{T}) where T = -one(T)/2
6767

68-
69-
70-
# See DLMF, Table 18.9.1
71-
# http://dlmf.nist.gov/18.9#i
72-
rec_An(b::ChebyshevT{T}, n::Int) where {T} = n==0 ? one(T) : 2one(T)
73-
rec_Bn(b::ChebyshevT{T}, n::Int) where {T} = zero(T)
74-
rec_Cn(b::ChebyshevT{T}, n::Int) where {T} = one(T)
68+
# Recurrence relation
69+
rec_An(b::ChebyshevT{T}, n::Int) where T = chebyshev_1st_rec_An(T, n)
70+
rec_Bn(b::ChebyshevT{T}, n::Int) where T = chebyshev_1st_rec_Bn(T, n)
71+
rec_Cn(b::ChebyshevT{T}, n::Int) where T = chebyshev_1st_rec_Cn(T, n)
7572

7673
# We can define this O(1) evaluation method, but only for points that are
7774
# real and lie in [-1,1]
@@ -114,22 +111,15 @@ hasmeasure(dict::ChebyshevT) = true
114111
measure(dict::ChebyshevT{T}) where T = ChebyshevTWeight{T}()
115112
issymmetric(::ChebyshevT) = true
116113

117-
dict_innerproduct_native(b1::ChebyshevT, i::PolynomialDegree, b2::ChebyshevT, j::PolynomialDegree, m::ChebyshevTWeight;
118-
T = coefficienttype(b1), options...) =
114+
function dict_innerproduct_native(b1::ChebyshevT, i::PolynomialDegree,
115+
b2::ChebyshevT, j::PolynomialDegree, m::ChebyshevTWeight; options...)
116+
T = promote_type(domaintype(b1),domaintype(b2))
119117
innerproduct_chebyshev_full(i, j, T)
120-
121-
function innerproduct_chebyshev_full(i, j, T)
122-
if i == j
123-
if i == PolynomialDegree(0)
124-
convert(T, pi)
125-
else
126-
convert(T, pi)/2
127-
end
128-
else
129-
zero(T)
130-
end
131118
end
132119

120+
innerproduct_chebyshev_full(i, j, T) =
121+
i == j ? chebyshev_1st_hn(T, value(i)) : zero(T)
122+
133123
function dict_innerproduct_native(b1::ChebyshevT, i::PolynomialDegree,
134124
b2::ChebyshevT, j::PolynomialDegree, measure::Union{LegendreWeight,Lebesgue}; options...)
135125
n1 = degree(i)

src/bases/poly/chebyshev/ChebyshevU.jl

Lines changed: 10 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -30,28 +30,22 @@ issymmetric(::ChebyshevU) = true
3030
measure(dict::ChebyshevU{T}) where {T} = ChebyshevUWeight{T}()
3131
hasmeasure(::ChebyshevU) = true
3232

33-
function dict_innerproduct_native(b1::ChebyshevU, i::PolynomialDegree, b2::ChebyshevU, j::PolynomialDegree, m::ChebyshevUWeight;
34-
T = coefficienttype(b1), options...)
35-
if i == j
36-
convert(T, pi)/2
37-
else
38-
zero(T)
39-
end
33+
function dict_innerproduct_native(b1::ChebyshevU, i::PolynomialDegree,
34+
b2::ChebyshevU, j::PolynomialDegree, m::ChebyshevUWeight; options...)
35+
T = promote_type(domaintype(b1), domaintype(b2))
36+
i == j ? chebyshev_2nd_hn(T, value(i)) : zero(T)
4037
end
4138

4239

4340
# Parameters alpha and beta of the corresponding Jacobi polynomial
44-
jacobi_α(b::ChebyshevU{T}) where {T} = one(T)/2
45-
jacobi_β(b::ChebyshevU{T}) where {T} = one(T)/2
41+
jacobi_α(b::ChebyshevU{T}) where T = one(T)/2
42+
jacobi_β(b::ChebyshevU{T}) where T = one(T)/2
4643

4744

48-
# See DLMF, Table 18.9.1
49-
# http://dlmf.nist.gov/18.9#i
50-
rec_An(b::ChebyshevU{T}, n::Int) where {T} = convert(T, 2)
51-
52-
rec_Bn(b::ChebyshevU{T}, n::Int) where {T} = zero(T)
53-
54-
rec_Cn(b::ChebyshevU{T}, n::Int) where {T} = one(T)
45+
# Recurrence relation
46+
rec_An(b::ChebyshevU{T}, n::Int) where T = chebyshev_2nd_rec_An(T, n)
47+
rec_Bn(b::ChebyshevU{T}, n::Int) where T = chebyshev_2nd_rec_Bn(T, n)
48+
rec_Cn(b::ChebyshevU{T}, n::Int) where T = chebyshev_2nd_rec_Cn(T, n)
5549

5650

5751
"A Chebyshev polynomial of the second kind"

src/bases/poly/ops/connections.jl

Lines changed: 1 addition & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -4,7 +4,6 @@ const JACOBI_OR_LEGENDRE = Union{Legendre,AbstractJacobi}
44
isequaldict(b1::JACOBI_OR_LEGENDRE, b2::JACOBI_OR_LEGENDRE) =
55
length(b1)==length(b2) &&
66
jacobi_α(b1) == jacobi_α(b2) &&
7-
jacobi_β(b2) == jacobi_β(b2)
7+
jacobi_β(b1) == jacobi_β(b2)
88

99
# TODO: implement conversions
10-

src/bases/poly/ops/dlmf.jl

Lines changed: 104 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,104 @@
1+
# Formulas from the Digital Library of Mathematical Functions.
2+
3+
"""
4+
Pochhammer's Symbol.
5+
6+
See [DLMF](https://dlmf.nist.gov/5.2#iii) formula (5.2.4).
7+
"""
8+
pochhammer(a, n) = n == 0 ? one(a) : prod(a+i for i in 0:n-1)
9+
# or in terms of the gamma function (DLMF, 5.2.5): gamma(n+a) / gamma(a)
10+
11+
##############################
12+
# DLMF Table 18.3.1
13+
# https://dlmf.nist.gov/18.3
14+
##############################
15+
# Formulas for squared norm hn and leading order coefficient kn
16+
17+
# The formula for An and A0
18+
dlmf_18d3d1_A(α, β, n) = 2^+β+1)*gamma(n+α+1)*gamma(n+β+1) / (factorial(n)*(2n+α+β+1)*gamma(n+α+β+1))
19+
dlmf_18d3d1_A0(α, β) = 2^+β+1)*gamma+1)*gamma+1) / gamma+β+2)
20+
21+
jacobi_hn(α, β, n) = n == 0 ? dlmf_18d3d1_A0(α, β) : dlmf_18d3d1_A(α, β, n)
22+
jacobi_kn(α, β, n) = pochhammer(n+α+β+1, n) / (2^n * factorial(n))
23+
24+
ultraspherical_hn::T, n) where T = 2^(1-2λ)*T(π)*gamma(n+2λ) / ((n+λ)*gamma(λ)^2*factorial(n))
25+
ultraspherical_kn(λ, n) = 2^n*pochhammer(λ, n) / factorial(n)
26+
27+
chebyshev_1st_hn(::Type{T}, n) where T = n == 0 ? T(π) : T(π)/2
28+
chebyshev_1st_kn(::Type{T}, n) where T = n == 0 ? one(T) : T(2)^(n-1)
29+
30+
chebyshev_2nd_hn(::Type{T}, n) where T = T(π)/2
31+
chebyshev_2nd_kn(::Type{T}, n) where T = T(2)^n
32+
33+
chebyshev_3rd_hn(::Type{T}, n) where T = T(π)
34+
chebyshev_3rd_kn(::Type{T}, n) where T = T(2)^n
35+
36+
chebyshev_4rd_hn(::Type{T}, n) where T = T(π)
37+
chebyshev_4rd_kn(::Type{T}, n) where T = T(2)^n
38+
39+
legendre_hn(::Type{T}, n) where T = 2 / T(2n+1)
40+
legendre_kn(::Type{T}, n) where T = 2^n*pochhammer(one(T)/2, n) / factorial(n)
41+
42+
laguerre_hn(α, n) = gamma(n+α+1) / factorial(n)
43+
laguerre_kn::T, n) where T = (-one(T))^n / factorial(n)
44+
45+
hermite_hn(::Type{T}, n) where T = sqrt(T(π)) * 2^n * factorial(n)
46+
hermite_kn(::Type{T}, n) where T = T(2)^n
47+
48+
49+
##############################
50+
# DLMF Table 18.9.1
51+
# http://dlmf.nist.gov/18.9#i
52+
##############################
53+
# Recurrence coefficients in the form
54+
# p_{n+1}(x) = (A_n x + B_n) p_n(x) - C_n p_{n-1}(x)
55+
# with initial values
56+
# p_0(x) = 1 and p_1(x) = A_0x + B_0
57+
58+
59+
# Jacobi, (18.9.2)
60+
function jacobi_rec_An(α, β, n)
61+
if (n == 0) &&+ β + 1 == 0)
62+
+β)/2+1
63+
else
64+
(2*n + α + β + 1) * (2n + α + β + 2) / (2 * (n+1) * (n + α + β + 1))
65+
end
66+
end
67+
function jacobi_rec_Bn(α, β, n)
68+
if (n == 0) && ((α + β + 1 == 0) ||+β == 0))
69+
-β)/2
70+
else
71+
^2 - β^2) * (2*n + α + β + 1) / (2 * (n+1) * (n + α + β + 1) * (2*n + α + β))
72+
end
73+
end
74+
function jacobi_rec_Cn(α, β, n)
75+
(n + α) * (n + β) * (2*n + α + β + 2) / ((n+1) * (n + α + β + 1) * (2*n + α + β))
76+
end
77+
78+
# Legendre
79+
legendre_rec_An(::Type{T}, n) where T = (2*n+1)/T(n+1)
80+
legendre_rec_Bn(::Type{T}, n) where T = zero(T)
81+
legendre_rec_Cn(::Type{T}, n) where T = n/T(n+1)
82+
83+
# Ultraspherical
84+
ultraspherical_rec_An(λ, n) = 2(n+λ) / (n+1)
85+
ultraspherical_rec_Bn(λ, n) = zero(λ)
86+
ultraspherical_rec_Cn(λ, n) = (n+2λ-1) / (n+1)
87+
88+
# Chebyshev
89+
chebyshev_1st_rec_An(::Type{T}, n) where T = n == 0 ? one(T) : T(2)
90+
chebyshev_1st_rec_Bn(::Type{T}, n) where T = zero(T)
91+
chebyshev_1st_rec_Cn(::Type{T}, n) where T = one(T)
92+
chebyshev_2nd_rec_An(::Type{T}, n) where T = T(2)
93+
chebyshev_2nd_rec_Bn(::Type{T}, n) where T = zero(T)
94+
chebyshev_2nd_rec_Cn(::Type{T}, n) where T = one(T)
95+
96+
# Laguerre
97+
laguerre_rec_An::T, n) where T = -one(T)/(n+1)
98+
laguerre_rec_Bn(α, n) = (2n+α+1) / (n+1)
99+
laguerre_rec_Cn(α, n) = (n+α) / (n+1)
100+
101+
# Hermite
102+
hermite_rec_An(::Type{T}, n) where T = T(2)
103+
hermite_rec_Bn(::Type{T}, n) where T = zero(T)
104+
hermite_rec_Cn(::Type{T}, n) where T = T(2n)

src/bases/poly/ops/hermite.jl

Lines changed: 6 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -29,22 +29,15 @@ issymmetric(::Hermite) = true
2929

3030
gauss_rule(dict::Hermite{T}) where T = GaussHermite{T}(length(dict))
3131

32-
# See DLMF, Table 18.9.1
33-
# http://dlmf.nist.gov/18.9#i
34-
rec_An(b::Hermite, n::Int) = 2
35-
rec_Bn(b::Hermite, n::Int) = 0
36-
rec_Cn(b::Hermite, n::Int) = 2*n
32+
# recurrence relation
33+
rec_An(b::Hermite{T}, n::Int) where T = hermite_rec_An(T, n)
34+
rec_Bn(b::Hermite{T}, n::Int) where T = hermite_rec_Bn(T, n)
35+
rec_Cn(b::Hermite{T}, n::Int) where T = hermite_rec_Cn(T, n)
3736

3837
coefficients_of_x(b::Hermite{T}) where {T} = (c=zeros(b); c[2]=one(T)/2; c)
3938

40-
function dict_innerproduct_native(d1::Hermite, i::PolynomialDegree, d2::Hermite, j::PolynomialDegree, measure::HermiteWeight; options...)
41-
T = coefficienttype(d1)
42-
if i == j
43-
sqrt(convert(T, pi)) * (1<<value(i)) * factorial(value(i))
44-
else
45-
zero(T)
46-
end
47-
end
39+
dict_innerproduct_native(d1::Hermite{T}, i::PolynomialDegree, d2::Hermite, j::PolynomialDegree, measure::HermiteWeight; options...) where T =
40+
i == j ? hermite_hn(T, value(i)) : zero(T)
4841

4942
show(io::IO, b::Hermite{Float64}) = print(io, "Hermite($(length(b)))")
5043
show(io::IO, b::Hermite{T}) where T = print(io, "Hermite{$(T)}($(length(b)))")

0 commit comments

Comments
 (0)