Skip to content

Commit c455ece

Browse files
daanhbdaanhb
committed
add ultraspherical polynomials
1 parent e920168 commit c455ece

File tree

5 files changed

+116
-83
lines changed

5 files changed

+116
-83
lines changed

‎Project.toml‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -40,7 +40,7 @@ BasisFunctionsPGFPlotsXExt = "PGFPlotsX"
4040
BlockArrays = "0.16"
4141
Calculus = "0.5"
4242
CompositeTypes = "0.1.3"
43-
DomainIntegrals = "0.5"
43+
DomainIntegrals = "0.5.1"
4444
DomainSets = "0.7"
4545
DoubleFloats = "1"
4646
FFTW = "1.6"

‎src/BasisFunctions.jl‎

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -83,7 +83,7 @@ import DomainSets:
8383

8484
import DomainIntegrals:
8585
integral, integrate,
86-
support,
86+
support, moment,
8787
weightfun, weightfunction,
8888
unsafe_weightfun, unsafe_weightfunction,
8989
weight, unsafe_weight,
@@ -92,7 +92,7 @@ import DomainIntegrals:
9292
isnormalized, isuniform,
9393
ismappedmeasure,
9494
jacobi_α, jacobi_β, laguerre_α
95-
using DomainIntegrals: ProductWeight
95+
using DomainIntegrals: ProductWeight, AbstractJacobiWeight
9696
export integral
9797

9898
@reexport using GridArrays
@@ -414,6 +414,7 @@ include("bases/poly/ops/laguerre.jl")
414414
include("bases/poly/ops/hermite.jl")
415415
include("bases/poly/ops/generic_op.jl")
416416
include("bases/poly/ops/specialOPS.jl")
417+
include("bases/poly/ops/connections.jl")
417418

418419
include("util/plot.jl")
419420

‎src/bases/poly/ops/connections.jl‎

Lines changed: 10 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,10 @@
1+
2+
const JACOBI_OR_LEGENDRE = Union{Legendre,AbstractJacobi}
3+
4+
isequaldict(b1::JACOBI_OR_LEGENDRE, b2::JACOBI_OR_LEGENDRE) =
5+
length(b1)==length(b2) &&
6+
jacobi_α(b1) == jacobi_α(b2) &&
7+
jacobi_β(b2) == jacobi_β(b2)
8+
9+
# TODO: implement conversions
10+

‎src/bases/poly/ops/jacobi.jl‎

Lines changed: 99 additions & 77 deletions
Original file line numberDiff line numberDiff line change
@@ -1,17 +1,97 @@
11

2+
# Supporting functions
3+
4+
# See DLMF (18.9.2)
5+
# http://dlmf.nist.gov/18.9#i
6+
function jacobi_rec_An(α::T, β::T, n::Int) where T
7+
if (n == 0) && (α + β + 1 == 0)
8+
one(T)/2*(α+β)+1
9+
else
10+
T(2*n + α + β + 1) * (2n + α + β + 2) / T(2 * (n+1) * (n + α + β + 1))
11+
end
12+
end
13+
function jacobi_rec_Bn(α::T, β::T, n::Int) where T
14+
if (n == 0) && ((α + β + 1 == 0) || (α+β == 0))
15+
one(T)/2*(α-β)
16+
else
17+
T(α^2 - β^2) * (2*n + α + β + 1) / T(2 * (n+1) * (n + α + β + 1) * (2*n + α + β))
18+
end
19+
end
20+
function jacobi_rec_Cn(α::T, β::T, n::Int) where T
21+
T(n + α) * (n + β) * (2*n + α + β + 2) / T((n+1) * (n + α + β + 1) * (2*n + α + β))
22+
end
23+
24+
# TODO: move these elsewhere. But where?
25+
# The packages GridArrays (which defines the nodes) and DomainIntegrals (which define the
26+
# jacobi_x functions) do not depend on each other.
27+
jacobi_α(x::GridArrays.JacobiNodes) = x.α
28+
jacobi_β(x::GridArrays.JacobiNodes) = x.β
29+
jacobi_α(w::GridArrays.JacobiWeights) = w.α
30+
jacobi_β(w::GridArrays.JacobiWeights) = w.β
31+
32+
33+
"Abstract supertype of Jacobi polynomials."
34+
abstract type AbstractJacobi{T} <: IntervalOPS{T} end
35+
36+
iscompatible(b1::AbstractJacobi, b2::AbstractJacobi) =
37+
jacobi_α(b1) == jacobi_α(b2) && jacobi_β(b1) == jacobi_β(b2)
38+
isequaldict(b1::AbstractJacobi, b2::AbstractJacobi) =
39+
length(b1)==length(b2) && iscompatible(b1,b2)
40+
41+
isorthogonal(b::AbstractJacobi, μ::DomainIntegrals.AbstractJacobiWeight) =
42+
jacobi_α(b) == jacobi_α(μ) && jacobi_β(b) == jacobi_β(μ)
43+
44+
issymmetric(dict::AbstractJacobi) = jacobi_α(dict)≈jacobi_β(dict)
45+
46+
rec_An(b::AbstractJacobi, n::Int) = jacobi_rec_An(jacobi_α(b), jacobi_β(b), n)
47+
rec_Bn(b::AbstractJacobi, n::Int) = jacobi_rec_Bn(jacobi_α(b), jacobi_β(b), n)
48+
rec_Cn(b::AbstractJacobi, n::Int) = jacobi_rec_Cn(jacobi_α(b), jacobi_β(b), n)
49+
50+
interpolation_grid(b::AbstractJacobi) = JacobiNodes(length(b), jacobi_α(b), jacobi_β(b))
51+
iscompatiblegrid(b::AbstractJacobi, grid::JacobiNodes) =
52+
length(b) == length(grid) &&
53+
jacobi_α(b) ≈ jacobi_α(grid) && jacobi_β(b) ≈ jacobi_β(grid)
54+
isorthogonal(b::AbstractJacobi, μ::GaussJacobi) =
55+
jacobi_α(b) ≈ jacobi_α(μ) && jacobi_β(b) ≈ jacobi_β(μ) && opsorthogonal(b, μ)
56+
57+
gauss_rule(b::AbstractJacobi) = GaussJacobi(length(b), jacobi_α(b), jacobi_β(b))
58+
59+
first_moment(b::AbstractJacobi) = moment(measure(b))
60+
61+
function dict_innerproduct_native(d1::AbstractJacobi, i::PolynomialDegree, d2::AbstractJacobi, j::PolynomialDegree, measure::AbstractJacobiWeight; options...)
62+
T = coefficienttype(d1)
63+
if iscompatible(d1, d2) && isorthogonal(d1, measure)
64+
if i == j
65+
a = jacobi_α(d1)
66+
b = jacobi_β(d1)
67+
n = value(i)
68+
2^(a+b+1)/(2n+a+b+1) * gamma(n+a+1)*gamma(n+b+1)/factorial(n)/gamma(n+a+b+1)
69+
else
70+
zero(T)
71+
end
72+
else
73+
dict_innerproduct1(d1, i, d2, j, measure; options...)
74+
end
75+
end
76+
77+
78+
279
"""
380
A basis of the classical Jacobi polynomials on the interval `[-1,1]`.
481
These polynomials are orthogonal with respect to the weight function
582
```
683
w(x) = (1-x)^α (1+x)^β.
784
```
885
"""
9-
struct Jacobi{T} <: IntervalOPS{T}
86+
struct Jacobi{T} <: AbstractJacobi{T}
1087
n :: Int
1188
α :: T
1289
β :: T
1390

14-
Jacobi{T}(n, α = 0, β = 0) where {T} = new{T}(n, α, β)
91+
function Jacobi{T}(n, α = 0, β = 0) where T
92+
@assert α > -1 && β > -1
93+
new{T}(n, α, β)
94+
end
1595
end
1696

1797
Jacobi(n::Int; α = 0, β = 0) = Jacobi(n, α, β)
@@ -28,83 +108,33 @@ const JacobiExpansion{T,C} = Expansion{T,T,Jacobi{T},C}
28108

29109
similar(b::Jacobi, ::Type{T}, n::Int) where {T} = Jacobi{T}(n, b.α, b.β)
30110

31-
iscompatible(d1::Jacobi, d2::Jacobi) = d1.α == d2.α && d1.β == d2.β
32-
isequaldict(b1::Jacobi, b2::Jacobi) = length(b1)==length(b2) && iscompatible(b1,b2)
33-
34-
first_moment(b::Jacobi{T}) where {T} = (b.α+b.β+1≈0) ?
35-
T(2).^(b.α+b.β+1)*gamma(b.α+1)*gamma(b.β+1) :
36-
T(2).^(b.α+b.β+1)*gamma(b.α+1)*gamma(b.β+1)/(b.α+b.β+1)/gamma(b.α+b.β+1)
37-
# 2^(b.α+b.β) / (b.α+b.β+1) * gamma(b.α+1) * gamma(b.β+1) / gamma(b.α+b.β+1)
38-
39111
jacobi_α(b::Jacobi) = b.α
40112
jacobi_β(b::Jacobi) = b.β
41113

42114
measure(b::Jacobi) = JacobiWeight(b.α, b.β)
43115

44116

45-
isorthogonal(dict::Jacobi, measure::JacobiWeight) =
46-
dict.α == measure.α && dict.β == measure.β
47-
48-
gauss_rule(dict::Jacobi) = GaussJacobi(length(dict), dict.α, dict.β)
49117

50-
interpolation_grid(dict::Jacobi) = JacobiNodes(length(dict), dict.α, dict.β)
51-
iscompatiblegrid(dict::Jacobi, grid::JacobiNodes) = length(dict) == length(grid) && dict.α ≈ grid.α && dict.β ≈ grid.β
52-
isorthogonal(dict::Jacobi, measure::GaussJacobi) = jacobi_α(dict) ≈ jacobi_α(measure) && jacobi_β(dict) ≈ jacobi_β(measure) && opsorthogonal(dict, measure)
53-
54-
issymmetric(dict::Jacobi) = jacobi_α(dict)≈jacobi_β(dict)
118+
"""
119+
The basis of ultraspherical orthogonal polynomials on `[-1,1]`.
55120
56-
# See DLMF (18.9.2)
57-
# http://dlmf.nist.gov/18.9#i
58-
function rec_An(b::Jacobi{T}, n::Int) where {T}
59-
if (n == 0) && (b.α + b.β+1 == 0)
60-
one(T)/2*(b.α+b.β)+1
61-
else
62-
T(2*n + b.α + b.β + 1) * (2*n + b.α + b.β + 2) / T(2 * (n+1) * (n + b.α + b.β + 1))
63-
end
121+
They are orthogonal with respect to the weight function
122+
```
123+
w(x) = (1-x^2)^(λ-1/2).
124+
```
125+
"""
126+
struct Ultraspherical{T} <: AbstractJacobi{T}
127+
n :: Int
128+
λ :: T
64129
end
65130

66-
function rec_Bn(b::Jacobi{T}, n::Int) where {T}
67-
if (n == 0) && ((b.α + b.β + 1 == 0) || (b.α+b.β == 0))
68-
one(T)/2*(b.α-b.β)
69-
else
70-
T(b.α^2 - b.β^2) * (2*n + b.α + b.β + 1) / T(2 * (n+1) * (n + b.α + b.β + 1) * (2*n + b.α + b.β))
71-
end
72-
end
131+
const Gegenbauer = Ultraspherical
73132

74-
rec_Cn(b::Jacobi{T}, n::Int) where {T} =
75-
T(n + b.α) * (n + b.β) * (2*n + b.α + b.β + 2) / T((n+1) * (n + b.α + b.β + 1) * (2*n + b.α + b.β))
133+
jacobi_α(b::Ultraspherical{T}) where T = b.λ - one(T)/2
134+
jacobi_β(b::Ultraspherical{T}) where T = b.λ - one(T)/2
76135

77-
function dict_innerproduct_native(d1::Jacobi, i::PolynomialDegree, d2::Jacobi, j::PolynomialDegree, measure::JacobiWeight; options...)
78-
T = coefficienttype(d1)
79-
if iscompatible(d1, d2) && isorthogonal(d1, measure)
80-
if i == j
81-
a = d1.α
82-
b = d1.β
83-
n = value(i)
84-
2^(a+b+1)/(2n+a+b+1) * gamma(n+a+1)*gamma(n+b+1)/factorial(n)/gamma(n+a+b+1)
85-
else
86-
zero(T)
87-
end
88-
else
89-
dict_innerproduct1(d1, i, d2, j, measure; options...)
90-
end
91-
end
136+
measure(b::Ultraspherical{T}) where T = DomainIntegrals.UltrasphericalWeight(b.λ)
92137

93-
# function expansion_sum(src1::Jacobi, src2::Jacobi, coef1, coef2)
94-
# if iscompatible(src1, src2)
95-
# default_expansion_sum(src1, src2, coef1, coef2)
96-
# else
97-
# to_chebyshev_expansion_sum(src1, src2, coef1, coef2)
98-
# end
99-
# end
100-
101-
# function expansion_multiply(src1::Jacobi, src2::Jacobi, coef1, coef2)
102-
# if iscompatible(src1, src2)
103-
# default_poly_expansion_multiply(src1, src2, coef1, coef2)
104-
# else
105-
# to_chebyshev_expansion_multiply(src1, src2, coef1, coef2)
106-
# end
107-
# end
108138

109139
## Printing
110140
function show(io::IO, b::Jacobi{Float64})
@@ -123,14 +153,6 @@ function show(io::IO, b::Jacobi{T}) where T
123153
end
124154
end
125155

126-
# TODO: move to its own file and make more complete
127-
# Or better yet: implement in terms of Jacobi polynomials
128-
struct UltrasphericalBasis{T} <: OPS{T}
129-
n :: Int
130-
alpha :: T
131-
end
132-
133-
jacobi_α(b::UltrasphericalBasis) = b.α
134-
jacobi_β(b::UltrasphericalBasis) = b.α
135-
136-
weightfun(b::UltrasphericalBasis, x) = (1-x)^(b.α) * (1+x)^(b.α)
156+
show(io::IO, b::Ultraspherical{Float64}) = print(io, "Ultraspherical($(length(b)), $(b.λ))")
157+
show(io::IO, b::Ultraspherical{T}) where T =
158+
print(io, "Ultraspherical{$(T)}($(length(b)), $(b.λ))")

‎test/test_generic_OPS.jl‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -8,8 +8,8 @@ using Test
88
function test_half_range_chebyshev()
99
n = 60
1010
T = 2.
11-
b1 = HalfRangeChebyshevIkind(n,T)
12-
b2 = HalfRangeChebyshevIIkind(n,T)
11+
b1 = BasisFunctions.HalfRangeChebyshevIkind(n,T)
12+
b2 = BasisFunctions.HalfRangeChebyshevIIkind(n,T)
1313

1414
α = -.5
1515
nodes, weights = gaussjacobi(100, α, 0)
@@ -142,7 +142,7 @@ end
142142

143143
function test_roots_of_legendre_halfrangechebyshev()
144144
N = 10
145-
B = HalfRangeChebyshevIkind(N,2.)
145+
B = BasisFunctions.HalfRangeChebyshevIkind(N,2.)
146146
@test 1+maximum(abs.(B[N].(real(roots(resize(B,N-1))))))≈1
147147
B = Legendre(N)
148148
@test 1+maximum(abs.(B[N].(real(roots(resize(B,N-1))))))≈1

0 commit comments

Comments
 (0)