Skip to content

Commit acc682e

Browse files
MikaelSlevinskyMikaelSlevinsky
authored
Merge pull request #87 from JuliaApproximation/dunkl-xu-disk
2 parents da04724 + 6145dcb commit acc682e

File tree

4 files changed

+286
-6
lines changed

4 files changed

+286
-6
lines changed

‎src/MultivariateOrthogonalPolynomials.jl

Lines changed: 11 additions & 6 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

@@ -20,14 +20,19 @@ import InfiniteArrays: InfiniteCardinal
2020

2121
import ClassicalOrthogonalPolynomials: jacobimatrix, Weighted, orthogonalityweight, HalfWeighted
2222
import HarmonicOrthogonalPolynomials: BivariateOrthogonalPolynomial, MultivariateOrthogonalPolynomial, Plan,
23-
PartialDerivative, BlockOneTo, BlockRange1, interlace
23+
PartialDerivative, AngularMomentum, 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, AngularMomentum,
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: 189 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,189 @@
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+
@simplify function *(Dx::PartialDerivative{1}, w_P::WeightedDunklXuDisk)
61+
wP, P = w_P.args
62+
@assert P.β == wP.β
63+
β = P.β
64+
n = mortar(Fill.(oneto(∞),oneto(∞)))
65+
k = mortar(Base.OneTo.(oneto(∞)))
66+
dat = BlockBroadcastArray(hcat,
67+
(-4 .* (k .+ (β - 1)).*(n .- k .+ 1)./(2k .+ (2β - 1))),
68+
0 .* n,
69+
(-k .* (k .+ 1) ./ ((2k .+ (2β - 1)) .* (k .+ β)) .* (n .+ k .+ 2β))
70+
)
71+
WeightedDunklXuDisk(β-1) * _BandedBlockBandedMatrix(dat', axes(k,1), (1,-1), (2,0))
72+
end
73+
74+
@simplify function *(Dy::PartialDerivative{2}, w_P::WeightedDunklXuDisk)
75+
wP, P = w_P.args
76+
@assert P.β == wP.β
77+
k = mortar(Base.OneTo.(oneto(∞)))
78+
T = promote_type(eltype(Dy), eltype(P)) # avoid bug in convert
79+
WeightedDunklXuDisk(P.β-1) * _BandedBlockBandedMatrix((T(-2).*k)', axes(k,1), (1,-1), (1,-1))
80+
end
81+
82+
# P^{(β)} ↗ P^{(β+1)}
83+
function dunklxu_raising(β)
84+
n = mortar(Fill.(oneto(∞),oneto(∞)))
85+
k = mortar(Base.OneTo.(oneto(∞)))
86+
dat = PseudoBlockArray(Vcat(
87+
(-(k .+ (β - 1)) .* (2n .+ (2β - 1)) ./ ((2k .+ (2β - 1)) .* (4n .+ 4β)))', # n-2, k-2
88+
(0 .* n)', # n-2, k-1
89+
(-(k .+ 2β) .* (k .+ (2β + 1)) ./ ((2k .+ (2β - 1)) .* (2k .+ 2β)) .* (2n .+ (2β-1)) ./ (8n .+ 8β))', # n-2, k
90+
(0 .* n)', # n-1, k-2
91+
(0 .* n)', # n-1, k-1
92+
(0 .* n)', # n-1, k
93+
(2 .* (k .+ (β - 1)) .* (n .- k .+ 1) .* (n .- k .+ 2) ./ ((2k .+ (2β-1)) .* (2n .+ 2β) .* (2n .+ (2β + 1))))', # n, k-2
94+
(0 .* n)', # n, k-1
95+
((k .+ 2β) .* (k .+ (2β + 1)) ./ ((2k .+ (2β - 1)) .* (2k .+ 2β)) .* (n .+ k .+ 2β) .* (n .+ k .+ (2β + 1)) ./ ((2n .+ 2β) .* (2n .+ (2β + 1))))' # n, k
96+
), (blockedrange(Fill(3,3)), axes(n,1)))
97+
_BandedBlockBandedMatrix(dat, axes(k,1), (0,2), (0,2))
98+
end
99+
100+
function \(A::DunklXuDisk, B::DunklXuDisk)
101+
if A.β == B.β
102+
Eye((axes(B,2),))
103+
elseif A.β == B.β + 1
104+
dunklxu_raising(B.β)
105+
else
106+
error("not implemented for $A and $B")
107+
end
108+
end
109+
110+
# (1-x^2-y^2) P^{(β)} ↘ P^{(β-1)}
111+
function dunklxu_lowering(β)
112+
n = mortar(Fill.(oneto(∞),oneto(∞)))
113+
k = mortar(Base.OneTo.(oneto(∞)))
114+
dat = PseudoBlockArray(Vcat(
115+
((2n .+ (2β - 1)) ./ (n .+ β) .* (k .+ (β - 1)) ./ (2k .+ (2β - 1)))', # n, k
116+
(0 .* n)', # n, k+1
117+
((n .+ (β - 1/2)) ./ (n .+ β) .* k .* (k .+ 1) ./ ((2k .+ (2β - 1)) .* (2k .+ 2β)))', # n, k+2
118+
(0 .* n)', # n+1, k
119+
(0 .* n)', # n+1, k+1
120+
(0 .* n)', # n+1, k+2
121+
(-8 .* (n .- k .+ 1) .* (n .- k .+ 2) ./ ((2n .+ 2β) .* (2n .+ (2β + 1))) .* (k .+ (β - 1)) ./(2k .+ (2β - 1)))', # n+2, k
122+
(0 .* n)', # n+2, k+1
123+
(-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
124+
), (blockedrange(Fill(3,3)), axes(n,1)))
125+
_BandedBlockBandedMatrix(dat, axes(k,1), (2,0), (2,0))
126+
end
127+
128+
function \(w_A::WeightedDunklXuDisk, w_B::WeightedDunklXuDisk)
129+
wA,A = w_A.args
130+
wB,B = w_B.args
131+
132+
@assert wA.β == A.β
133+
@assert wB.β == B.β
134+
135+
if A.β == B.β
136+
Eye((axes(B,2),))
137+
elseif A.β + 1 == B.β
138+
dunklxu_lowering(B.β)
139+
else
140+
error("not implemented for $A and $B")
141+
end
142+
end
143+
144+
\(w_A::DunklXuDisk, w_B::WeightedDunklXuDisk) =
145+
(DunklXuDiskWeight(0) .* w_A) \ w_B
146+
147+
# Actually Jxᵀ
148+
function jacobimatrix(::Val{1}, P::DunklXuDisk)
149+
β = P.β
150+
n = mortar(Fill.(oneto(∞),oneto(∞)))
151+
k = mortar(Base.OneTo.(oneto(∞)))
152+
dat = PseudoBlockArray(Vcat(
153+
((2n .+ (2β - 1)) ./ (4n .+ 4β))', # n-1, k
154+
(0 .* n)', # n, k
155+
((n .- k .+ 1) .* (n .+ k .+ 2β) ./ ((n .+ β) .* (2n .+ (2β + 1))))', # n+1, k
156+
), (blockedrange(Fill(1, 3)), axes(n,1)))
157+
_BandedBlockBandedMatrix(dat, axes(k,1), (1,1), (0,0))
158+
end
159+
160+
# Actually Jyᵀ
161+
function jacobimatrix(::Val{2}, P::DunklXuDisk)
162+
β = P.β
163+
n = mortar(Fill.(oneto(∞),oneto(∞)))
164+
k = mortar(Base.OneTo.(oneto(∞)))
165+
dat = PseudoBlockArray(Vcat(
166+
((k .+ (β - 1)) .* (2n .+ (2β - 1)) ./ ((2k .+ (2β - 1)) .* (2n .+ 2β)))', # n-1, k-1
167+
(0 .* n)', # n-1, k
168+
(-k .* (k .+ 2β) .* (2n .+ (2β - 1)) ./ ((2k .+ (2β - 1)) .* (2k .+ 2β) .* (4n .+ 4β)))', # n-1, k+1
169+
(0 .* n)', # n, k-1
170+
(0 .* n)', # n, k
171+
(0 .* n)', # n, k+1
172+
(-(2k .+ (2β - 2)) .* (n .- k .+ 1) .* (n .- k .+ 2) ./ ((2k .+ (2β - 1)) .* (n .+ β) .* (2n .+ (2β + 1))))', # n+1, k-1
173+
(0 .* n)', # n+1, k
174+
(k .* (k .+ 2β) .* (n .+ k .+ 2β) .* (n .+ k .+ (2β + 1)) ./ ((2k .+ (2β - 1)) .* (2k .+ 2β) .* (n .+ β) .* (2n .+ (2β + 1))))', # n+1, k+1
175+
), (blockedrange(Fill(3, 3)), axes(n,1)))
176+
_BandedBlockBandedMatrix(dat, axes(k,1), (1,1), (1,1))
177+
end
178+
179+
@simplify function *(A::AngularMomentum, P::DunklXuDisk)
180+
β = P.β
181+
n = mortar(Fill.(oneto(∞),oneto(∞)))
182+
k = mortar(Base.OneTo.(oneto(∞)))
183+
dat = PseudoBlockArray(Vcat(
184+
(2 .* (k .+ (β - 1)) .* (n .- k .+ 1) ./ (2k .+ (2β - 1)))', # n, k-1
185+
(0 .* n)', # n, k
186+
(-k .* (k .+ 2β) .* (n .+ k .+ 2β) ./ ((2k .+ (2β - 1)) .* (2k .+ 2β)))', # n, k+1
187+
), (blockedrange(Fill(3, 1)), axes(n,1)))
188+
DunklXuDisk(β) * _BandedBlockBandedMatrix(dat, axes(k,1), (0,0), (1,1))
189+
end

‎test/runtests.jl

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -2,5 +2,6 @@ using MultivariateOrthogonalPolynomials, Test
22

33
include("test_modalinterlace.jl")
44
include("test_disk.jl")
5+
include("test_rectdisk.jl")
56
include("test_triangle.jl")
67
# include("test_dirichlettriangle.jl")

‎test/test_rectdisk.jl

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

0 commit comments

Comments
 (0)