Skip to content

Fix ultraspherical Integral indexing #212

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 4 commits into from
Mar 9, 2023
Merged

Fix ultraspherical Integral indexing #212

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
a396f3f
ultraspherical integral indexing
jishnub Mar 9, 2023
855774d
Add Legendre tests
jishnub Mar 9, 2023
1d3529b
type inference in integral
jishnub Mar 9, 2023
0642035
Check for inference only on v1.8+
jishnub Mar 9, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
24 changes: 18 additions & 6 deletions src/Spaces/Ultraspherical/UltrasphericalOperators.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -65,7 +65,9 @@ function Integral(sp::Ultraspherical{<:Any,<:IntervalOrSegment}, m::Number)
ConcreteIntegral(sp,m)
else # Convert up
nsp = Ultraspherical(m,domain(sp))
IntegralWrapper(ConcreteIntegral(nsp,m)*Conversion(sp,nsp),m)
Integralop = ConcreteIntegral(nsp,m)
C = Conversion(sp,nsp)
IntegralWrapper(Integralop * C, m, sp, rangespace(Integralop))
end
end

Expand Down Expand Up @@ -100,8 +102,10 @@ linesum(f::Fun{<:Ultraspherical{LT,DD}}) where {LT,DD<:IntervalOrSegment} =
sum(setcanonicaldomain(f))*arclength(d)/2


rangespace(D::ConcreteIntegral{<:Ultraspherical{LT,DD}}) where {LT,DD<:IntervalOrSegment} =
order(domainspace(D)) == 1 ? Chebyshev(domain(D)) : Ultraspherical(order(domainspace(D))-D.order,domain(D))
function rangespace(D::ConcreteIntegral{<:Ultraspherical{LT,DD}}) where {LT,DD<:IntervalOrSegment}
k = order(domainspace(D))-D.order
k == 0 ? Chebyshev(domain(D)) : Ultraspherical(k, domain(D))
end

function getindex(Q::ConcreteIntegral{<:Ultraspherical{LT,DD}},k::Integer,j::Integer) where {LT,DD<:IntervalOrSegment}
T=eltype(Q)
Expand All @@ -112,7 +116,12 @@ function getindex(Q::ConcreteIntegral{<:Ultraspherical{LT,DD}},k::Integer,j::Int

if λ == 1 && k==j+1
C = complexlength(d)/2
strictconvert(T,C./(k-1))
strictconvert(T, C/(k-1))
elseif λ == m && k == j + m
C = complexlength(d)/2
U1toC = C/(k-1)
UmtoU1 = pochhammer(one(T)*λ,-(m-1))*(complexlength(d)/4)^(m-1)
strictconvert(T, U1toC * UmtoU1)
elseif λ > 1 && k==j+m
strictconvert(T,pochhammer(one(T)*λ,-m)*(complexlength(d)/4)^m)
else
Expand Down Expand Up @@ -300,8 +309,11 @@ end

# TODO: include in getindex to speed up
function Integral(sp::Chebyshev{DD},m::Integer) where {DD<:IntervalOrSegment}
IntegralWrapper(TimesOperator([Integral(Ultraspherical(m,domain(sp)),m),
Conversion(sp,Ultraspherical(m,domain(sp)))]),m)
usp = Ultraspherical(m,domain(sp))
I = Integral(usp,m)
C = Conversion(sp,usp)
T = TimesOperator(I, C)
IntegralWrapper(T, m, sp, sp)
end


Expand Down
22 changes: 22 additions & 0 deletions test/UltrasphericalTest.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -234,5 +234,27 @@ using Static
gexp = gexp - coefficients(gexp)[1]
@test g ≈ gexp
end

@testset for n in 3:6, d in ((), (0..1,)),
sp in (Chebyshev(d...), Ultraspherical(1, d...), Ultraspherical(0.5, d...))
f = Fun(sp, Float64[zeros(n); 2])
@test Integral(1) * (Derivative(1) * f) ≈ f
@test Integral(2) * (Derivative(2) * f) ≈ f
@test Integral(3) * (Derivative(3) * f) ≈ f
@test Derivative(1) * (Integral(1) * f) ≈ f
@test Derivative(2) * (Integral(2) * f) ≈ f
@test Derivative(3) * (Integral(3) * f) ≈ f
end

if VERSION >= v"1.8"
@inferred Integral(Chebyshev(), 3)
@inferred (() -> Integral(Ultraspherical(1), 3))()
@inferred (() -> Integral(Ultraspherical(3), 1))()
# without constant propagation, this should be a small union
TA = typeof(Integral(Ultraspherical(3), 1))
TB = typeof(Integral(Ultraspherical(1), 3))
@inferred Union{TA, TB} Integral(Ultraspherical(3), 1)
@inferred Union{TA, TB} Integral(Ultraspherical(1), 3)
end
end
end