Skip to content

Commit a677010

Browse files
MikaelSlevinskyMikaelSlevinsky
committed
start adding Dunkl-Xu recurrences
1 parent 04704a8 commit a677010

File tree

3 files changed

+234
-5
lines changed

3 files changed

+234
-5
lines changed

‎src/MultivariateOrthogonalPolynomials.jl

Lines changed: 10 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
module MultivariateOrthogonalPolynomials
2-
using ClassicalOrthogonalPolynomials, FastTransforms, BlockBandedMatrices, BlockArrays, DomainSets,
3-
QuasiArrays, StaticArrays, ContinuumArrays, InfiniteArrays, InfiniteLinearAlgebra,
2+
using ClassicalOrthogonalPolynomials, FastTransforms, BlockBandedMatrices, BlockArrays, DomainSets,
3+
QuasiArrays, StaticArrays, ContinuumArrays, InfiniteArrays, InfiniteLinearAlgebra,
44
LazyArrays, SpecialFunctions, LinearAlgebra, BandedMatrices, LazyBandedMatrices, ArrayLayouts,
55
HarmonicOrthogonalPolynomials
66

@@ -22,12 +22,17 @@ import ClassicalOrthogonalPolynomials: jacobimatrix, Weighted, orthogonalityweig
2222
import HarmonicOrthogonalPolynomials: BivariateOrthogonalPolynomial, MultivariateOrthogonalPolynomial, Plan,
2323
PartialDerivative, BlockOneTo, BlockRange1, interlace
2424

25-
export UnitTriangle, UnitDisk, JacobiTriangle, TriangleWeight, WeightedTriangle, PartialDerivative, Laplacian,
26-
MultivariateOrthogonalPolynomial, BivariateOrthogonalPolynomial, Zernike, RadialCoordinate,
27-
zerniker, zernikez, Weighted, Block, ZernikeWeight, AbsLaplacianPower
25+
export MultivariateOrthogonalPolynomial, BivariateOrthogonalPolynomial,
26+
UnitTriangle, UnitDisk,
27+
JacobiTriangle, TriangleWeight, WeightedTriangle,
28+
DunklXuDisk, DunklXuDiskWeight, WeightedDunklXuDisk,
29+
Zernike, ZernikeWeight, zerniker, zernikez,
30+
PartialDerivative, Laplacian, AbsLaplacianPower,
31+
RadialCoordinate, Weighted, Block
2832

2933
include("ModalInterlace.jl")
3034
include("disk.jl")
35+
include("rectdisk.jl")
3136
include("triangle.jl")
3237

3338

‎src/rectdisk.jl

Lines changed: 155 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,155 @@
1+
"""
2+
DunklXuDisk(β)
3+
4+
are the orthogonal polynomials on the unit disk with respect to (1-x^2-y^2)^β, defined as `P_{n-k}^{k+β+1/2,k+β+1/2}(x)*(1-x^2)^{k/2}*P_k^{β,β}(y/sqrt{1-x^2})`.
5+
"""
6+
struct DunklXuDisk{T, V} <: BivariateOrthogonalPolynomial{T}
7+
β::V
8+
end
9+
10+
DunklXuDisk(β::T) where T = DunklXuDisk{float(T), T}(β)
11+
DunklXuDisk() = DunklXuDisk(0)
12+
13+
==(D1::DunklXuDisk, D2::DunklXuDisk) = D1.β == D2.β
14+
15+
axes(P::DunklXuDisk{T}) where T = (Inclusion(UnitDisk{T}()),blockedrange(oneto(∞)))
16+
17+
copy(A::DunklXuDisk) = A
18+
19+
Base.summary(io::IO, P::DunklXuDisk) = print(io, "DunklXuDisk($(P.β))")
20+
21+
22+
"""
23+
DunklXuDiskWeight(β)
24+
25+
is a quasi-vector representing `(1-x^2-y^2)^β`.
26+
"""
27+
struct DunklXuDiskWeight{T, V} <: Weight{T}
28+
β::V
29+
end
30+
31+
DunklXuDiskWeight(β::T) where T = DunklXuDiskWeight{float(T),T}(β)
32+
33+
axes(P::DunklXuDiskWeight{T}) where T = (Inclusion(UnitDisk{T}()),)
34+
35+
Base.summary(io::IO, P::DunklXuDiskWeight) = print(io, "(1-x^2-y^2)^$(P.β) on the unit disk")
36+
37+
const WeightedDunklXuDisk{T} = WeightedBasis{T,<:DunklXuDiskWeight,<:DunklXuDisk}
38+
39+
WeightedDunklXuDisk(β) = DunklXuDiskWeight(β) .* DunklXuDisk(β)
40+
41+
@simplify function *(Dx::PartialDerivative{1}, P::DunklXuDisk)
42+
β = P.β
43+
n = mortar(Fill.(oneto(∞),oneto(∞)))
44+
k = mortar(Base.OneTo.(oneto(∞)))
45+
dat = BlockBroadcastArray(hcat,
46+
((k .+ (β - 1)) ./ (2k .+ (2β - 1)) .* (n .- k .+ 1)),
47+
0 .* n,
48+
(((k .+ 2β) .* (k .+ (2β + 1))) ./ ((2k .+ (2β - 1)) .* (2k .+ 2β)) .* (n .+ k .+ 2β) ./ 2)
49+
)
50+
DunklXuDisk(β+1) * _BandedBlockBandedMatrix(dat', axes(k,1), (-1,1), (0,2))
51+
end
52+
53+
@simplify function *(Dy::PartialDerivative{2}, P::DunklXuDisk)
54+
β = P.β
55+
k = mortar(Base.OneTo.(oneto(∞)))
56+
T = promote_type(eltype(Dy), eltype(P)) # avoid bug in convert
57+
DunklXuDisk(β+1) * _BandedBlockBandedMatrix(((k .+ T(2β)) ./ 2)', axes(k,1), (-1,1), (-1,1))
58+
end
59+
60+
# P^{(β)} ↗ P^{(β+1)}
61+
function dunklxu_raising(β)
62+
n = mortar(Fill.(oneto(∞),oneto(∞)))
63+
k = mortar(Base.OneTo.(oneto(∞)))
64+
dat = PseudoBlockArray(Vcat(
65+
(-(k .+ (β - 1)) .* (2n .+ (2β - 1)) ./ ((2k .+ (2β - 1)) .* (4n .+ 4β)))', # n-2, k-2
66+
(0 .* n)', # n-2, k-1
67+
(-(k .+ 2β) .* (k .+ (2β + 1)) ./ ((2k .+ (2β - 1)) .* (2k .+ 2β)) .* (2n .+ (2β-1)) ./ (8n .+ 8β))', # n-2, k
68+
(0 .* n)', # n-1, k-2
69+
(0 .* n)', # n-1, k-1
70+
(0 .* n)', # n-1, k
71+
(2 .* (k .+ (β - 1)) .* (n .- k .+ 1) .* (n .- k .+ 2) ./ ((2k .+ (2β-1)) .* (2n .+ 2β) .* (2n .+ (2β + 1))))', # n, k-2
72+
(0 .* n)', # n, k-1
73+
((k .+ 2β) .* (k .+ (2β + 1)) ./ ((2k .+ (2β - 1)) .* (2k .+ 2β)) .* (n .+ k .+ 2β) .* (n .+ k .+ (2β + 1)) ./ ((2n .+ 2β) .* (2n .+ (2β + 1))))' # n, k
74+
), (blockedrange(Fill(3,3)), axes(n,1)))
75+
_BandedBlockBandedMatrix(dat, axes(k,1), (0,2), (0,2))
76+
end
77+
78+
function \(A::DunklXuDisk, B::DunklXuDisk)
79+
if A.β == B.β
80+
Eye((axes(B,2),))
81+
elseif A.β == B.β + 1
82+
dunklxu_raising(B.β)
83+
else
84+
error("not implemented for $A and $B")
85+
end
86+
end
87+
88+
# (1-x^2-y^2) P^{(β)} ↘ P^{(β-1)}
89+
function dunklxu_lowering(β)
90+
n = mortar(Fill.(oneto(∞),oneto(∞)))
91+
k = mortar(Base.OneTo.(oneto(∞)))
92+
dat = PseudoBlockArray(Vcat(
93+
((2n .+ (2β - 1)) ./ (n .+ β) .* (k .+ (β - 1)) ./ (2k .+ (2β - 1)))', # n, k
94+
(0 .* n)', # n, k+1
95+
((n .+ (β - 1/2)) ./ (n .+ β) .* k .* (k .+ 1) ./ ((2k .+ (2β - 1)) .* (2k .+ 2β)))', # n, k+2
96+
(0 .* n)', # n+1, k
97+
(0 .* n)', # n+1, k+1
98+
(0 .* n)', # n+1, k+2
99+
(-8 .* (n .- k .+ 1) .* (n .- k .+ 2) ./ ((2n .+ 2β) .* (2n .+ (2β + 1))) .* (k .+ (β - 1)) ./(2k .+ (2β - 1)))', # n+2, k
100+
(0 .* n)', # n+2, k+1
101+
(-4 .* (n .+ k .+ 2β) .* (n .+ k .+ (2β + 1)) ./ ((2n .+ 2β) .* (2n .+ (2β + 1))) .* k .* (k .+ 1) ./ ((2k .+ (2β - 1)) .* (2k .+ 2β)))' # n+2, k+2
102+
), (blockedrange(Fill(3,3)), axes(n,1)))
103+
_BandedBlockBandedMatrix(dat, axes(k,1), (2,0), (2,0))
104+
end
105+
106+
function \(w_A::WeightedDunklXuDisk, w_B::WeightedDunklXuDisk)
107+
wA,A = w_A.args
108+
wB,B = w_B.args
109+
110+
@assert wA.β == A.β
111+
@assert wB.β == B.β
112+
113+
if A.β == B.β
114+
Eye((axes(B,2),))
115+
elseif A.β + 1 == B.β
116+
dunklxu_lowering(B.β)
117+
else
118+
error("not implemented for $A and $B")
119+
end
120+
end
121+
122+
\(w_A::DunklXuDisk, w_B::WeightedDunklXuDisk) =
123+
(DunklXuDiskWeight(0) .* w_A) \ w_B
124+
125+
# Actually Jxᵀ
126+
function jacobimatrix(::Val{1}, P::DunklXuDisk)
127+
β = P.β
128+
n = mortar(Fill.(oneto(∞),oneto(∞)))
129+
k = mortar(Base.OneTo.(oneto(∞)))
130+
dat = PseudoBlockArray(Vcat(
131+
((2n .+ (2β - 1)) ./ (4n .+ 4β))', # n-1, k
132+
(0 .* n)', # n, k
133+
((n .- k .+ 1) .* (n .+ k .+ 2β) ./ ((n .+ β) .* (2n .+ (2β + 1))))', # n+1, k
134+
), (blockedrange(Fill(1, 3)), axes(n,1)))
135+
_BandedBlockBandedMatrix(dat, axes(k,1), (1,1), (0,0))
136+
end
137+
138+
# Actually Jyᵀ
139+
function jacobimatrix(::Val{2}, P::DunklXuDisk)
140+
β = P.β
141+
n = mortar(Fill.(oneto(∞),oneto(∞)))
142+
k = mortar(Base.OneTo.(oneto(∞)))
143+
dat = PseudoBlockArray(Vcat(
144+
((k .+ (β - 1)) .* (2n .+ (2β - 1)) ./ ((2k .+ (2β - 1)) .* (2n .+ 2β)))', # n-1, k-1
145+
(0 .* n)', # n-1, k
146+
(-k .* (k .+ 2β) .* (2n .+ (2β - 1)) ./ ((2k .+ (2β - 1)) .* (2k .+ 2β) .* (4n .+ 4β)))', # n-1, k+1
147+
(0 .* n)', # n, k-1
148+
(0 .* n)', # n, k
149+
(0 .* n)', # n, k+1
150+
(-(2k .+ (2β - 2)) .* (n .- k .+ 1) .* (n .- k .+ 2) ./ ((2k .+ (2β - 1)) .* (n .+ β) .* (2n .+ (2β + 1))))', # n+1, k-1
151+
(0 .* n)', # n+1, k
152+
(k .* (k .+ 2β) .* (n .+ k .+ 2β) .* (n .+ k .+ (2β + 1)) ./ ((2k .+ (2β - 1)) .* (2k .+ 2β) .* (n .+ β) .* (2n .+ (2β + 1))))', # n+1, k+1
153+
), (blockedrange(Fill(3, 3)), axes(n,1)))
154+
_BandedBlockBandedMatrix(dat, axes(k,1), (1,1), (1,1))
155+
end

‎test/test_rectdisk.jl

Lines changed: 69 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,69 @@
1+
using MultivariateOrthogonalPolynomials, StaticArrays, BlockArrays, BlockBandedMatrices, ArrayLayouts,
2+
QuasiArrays, Test, ClassicalOrthogonalPolynomials, BandedMatrices, FastTransforms, LinearAlgebra
3+
import MultivariateOrthogonalPolynomials: dunklxu_raising, dunklxu_lowering
4+
5+
@testset "Dunkl-Xu disk" begin
6+
@testset "basics" begin
7+
P = DunklXuDisk()
8+
@test copy(P) ≡ P
9+
10+
xy = axes(P,1)
11+
x,y = first.(xy),last.(xy)
12+
@test xy[SVector(0.1,0.2)] == SVector(0.1,0.2)
13+
@test x[SVector(0.1,0.2)] == 0.1
14+
@test y[SVector(0.1,0.2)] == 0.2
15+
end
16+
17+
@testset "operators" begin
18+
N = 5
19+
20+
β = 0.123
21+
22+
P = DunklXuDisk(β)
23+
Q = DunklXuDisk(β+1)
24+
WP = WeightedDunklXuDisk(β)
25+
WQ = WeightedDunklXuDisk(β+1)
26+
27+
x, y = first.(axes(P, 1)), last.(axes(P, 1))
28+
29+
L = WP \ WQ
30+
R = Q \ P
31+
32+
X = P \ (x .* P)
33+
Y = P \ (y .* P)
34+
35+
@test (L*R)[Block.(1:N), Block.(1:N)] ≈ (I - X^2 - Y^2)[Block.(1:N), Block.(1:N)]
36+
37+
∂x = PartialDerivative{1}(P)
38+
∂y = PartialDerivative{2}(P)
39+
40+
Dx = Q \ (∂x * P)
41+
Dy = Q \ (∂y * P)
42+
43+
Mx = Q \ (x .* Q)
44+
My = Q \ (y .* Q)
45+
46+
A = Mx[Block.(1:N), Block.(1:N+1)]*Dy[Block.(1:N+1), Block.(1:N)] - My[Block.(1:N), Block.(1:N+1)]*Dx[Block.(1:N+1), Block.(1:N)]
47+
48+
B = (Q \ P)[Block.(1:N), Block.(1:N)]
49+
50+
C = B \ A
51+
52+
@test C ≈ Tridiagonal(C)
53+
54+
λ = eigvals(Matrix(C))
55+
56+
@test λ ≈ im*imag(λ)
57+
58+
@testset "truncations" begin
59+
KR,JR = Block.(1:N),Block.(1:N)
60+
61+
@test L[KR,JR] isa BandedBlockBandedMatrix
62+
@test R[KR,JR] isa BandedBlockBandedMatrix
63+
@test X[KR,JR] isa BandedBlockBandedMatrix
64+
@test Y[KR,JR] isa BandedBlockBandedMatrix
65+
@test Dx[KR,JR] isa BandedBlockBandedMatrix
66+
@test Dy[KR,JR] isa BandedBlockBandedMatrix
67+
end
68+
end
69+
end

0 commit comments

Comments
 (0)