From 25503b0c82c9e8541b2ea0782e1f1628e0c1b082 Mon Sep 17 00:00:00 2001 From: Daan Huybrechs Date: Wed, 9 Oct 2024 21:19:14 +0200 Subject: [PATCH] add ultraspherical weight --- Project.toml | 6 ++- src/DomainIntegrals.jl | 1 + src/weights.jl | 103 +++++++++++++++++++++++++++++++---------- test/runtests.jl | 23 ++++++++- test/test_integrals.jl | 5 +- test/test_measures.jl | 22 +++++++++ test/test_suite.jl | 22 --------- 7 files changed, 131 insertions(+), 51 deletions(-) delete mode 100644 test/test_suite.jl diff --git a/Project.toml b/Project.toml index 5e79909..ad87ebf 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "DomainIntegrals" uuid = "cc6bae93-f070-4015-88fd-838f9505a86c" authors = ["Daan Huybrechs "] -version = "0.5" +version = "0.5.1" [deps] CompositeTypes = "b152e2b5-7a66-4b01-a709-34e65c35f657" @@ -12,10 +12,10 @@ HCubature = "19dc6840-f33b-545b-b366-655c7e3ffd49" IntervalSets = "8197267c-284f-5f27-9208-e0e47529a953" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc" +SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" [compat] -julia = "1.10" CompositeTypes = "0.1.3" DomainSets = "0.7.2" FastGaussQuadrature = "0.5" @@ -23,7 +23,9 @@ GaussQuadrature = "0.5.7" HCubature = "1.4.0" IntervalSets = "0.7" QuadGK = "2.3.1" +SpecialFunctions = "2" StaticArrays = "1.0" +julia = "1.10" [extras] Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" diff --git a/src/DomainIntegrals.jl b/src/DomainIntegrals.jl index 1b96757..87e7090 100644 --- a/src/DomainIntegrals.jl +++ b/src/DomainIntegrals.jl @@ -2,6 +2,7 @@ module DomainIntegrals using QuadGK, LinearAlgebra using GaussQuadrature, FastGaussQuadrature, HCubature +using SpecialFunctions using StaticArrays using IntervalSets, DomainSets diff --git a/src/weights.jl b/src/weights.jl index 26b06b6..8989829 100644 --- a/src/weights.jl +++ b/src/weights.jl @@ -7,6 +7,7 @@ export LebesgueMeasure, dx, LegendreWeight, JacobiWeight, + UltrasphericalWeight, ChebyshevWeight, ChebyshevTWeight, ChebyshevUWeight, @@ -103,6 +104,9 @@ struct LegendreWeight{T} <: LebesgueMeasure{T} end LegendreWeight() = LegendreWeight{Float64}() +jacobi_α(μ::LegendreWeight{T}) where {T} = zero(T) +jacobi_β(μ::LegendreWeight{T}) where {T} = zero(T) + similar(μ::LegendreWeight, ::Type{T}) where {T <: Real} = LegendreWeight{T}() support(μ::LegendreWeight{T}) where {T} = ChebyshevInterval{T}() @@ -116,12 +120,23 @@ abstract type AbstractJacobiWeight{T} <: Weight{T} end ==(μ1::AbstractJacobiWeight, μ2::AbstractJacobiWeight) = jacobi_α(μ1) == jacobi_α(μ2) && jacobi_β(μ1) == jacobi_β(μ2) +support(μ::AbstractJacobiWeight{T}) where {T} = ChebyshevInterval{T}() + +jacobi_moment(α, β) = 2^(α+β+1) * gamma(α+1)*gamma(β+1) / gamma(α+β+2) + +moment(μ::AbstractJacobiWeight) = jacobi_moment(jacobi_α(μ), jacobi_β(μ)) + +isnormalized(μ::AbstractJacobiWeight) = moment(μ) ≈ 1 + "The Jacobi weight on the interval `[-1,1]`." struct JacobiWeight{T} <: AbstractJacobiWeight{T} α :: T β :: T - JacobiWeight{T}(α = zero(T), β = zero(T)) where {T} = new(α, β) + function JacobiWeight{T}(α = zero(T), β = zero(T)) where {T} + @assert (α > -1) && (β > -1) + new(α, β) + end end JacobiWeight() = JacobiWeight{Float64}() JacobiWeight(α, β) = JacobiWeight(promote(α, β)...) @@ -131,7 +146,6 @@ JacobiWeight(α::N, β::N) where {N<:Number} = JacobiWeight(float(α), float(β) jacobi_weightfun(x, α, β) = (1-x)^α * (1+x)^β similar(μ::JacobiWeight, ::Type{T}) where {T <: Real} = JacobiWeight{T}(μ.α, μ.β) -support(μ::JacobiWeight{T}) where {T} = ChebyshevInterval{T}() unsafe_weightfun(μ::JacobiWeight, x) = jacobi_weightfun(x, μ.α, μ.β) jacobi_α(μ::JacobiWeight) = μ.α @@ -189,35 +203,74 @@ Base.show(io::IO, μ::ChebyshevUWeight) = print(io, "√(1+x)^2 dx (ChebyshevU) Display.object_parentheses(μ::ChebyshevUWeight) = true -convert(::Type{JacobiWeight}, μ::LegendreWeight{T}) where {T} = - JacobiWeight{T}(0, 0) -convert(::Type{JacobiWeight}, μ::ChebyshevTWeight{T}) where {T} = - JacobiWeight{T}(-one(T)/2, -one(T)/2) -convert(::Type{JacobiWeight}, μ::ChebyshevUWeight{T}) where {T} = - JacobiWeight{T}(one(T)/2, one(T)/2) +""" +The `Ultraspherical` (or `Gegenbauer`) weight is the measure on `[-1,1]` with +the weight function `w(x) = (1-x^2)^(λ - 1/2)`. +""" +struct UltrasphericalWeight{T} <: DomainIntegrals.AbstractJacobiWeight{T} + λ :: T + + function UltrasphericalWeight{T}(λ) where T + @assert λ > -one(T)/2 + new(λ) + end +end + +const GegenbauerWeight = UltrasphericalWeight + +UltrasphericalWeight(α::T) where {T<:AbstractFloat} = UltrasphericalWeight{T}(α) +UltrasphericalWeight(α::N) where {N<:Number} = UltrasphericalWeight(float(α), float(β)) + +ultraspherical_weightfun(x, λ::T) where T = (1-x^2)^(λ - one(T)/2) + +jacobi_α(μ::UltrasphericalWeight{T}) where {T} = μ.λ-one(T)/2 +jacobi_β(μ::UltrasphericalWeight{T}) where {T} = μ.λ-one(T)/2 + +similar(μ::UltrasphericalWeight, ::Type{T}) where {T <: Real} = + UltrasphericalWeight{T}(μ.λ) +unsafe_weightfun(μ::UltrasphericalWeight, x) = ultraspherical_weightfun(x, μ.λ) + +Base.show(io::IO, μ::UltrasphericalWeight) = print(io, "(1-x^2)^(λ - 1/2) dx (Ultraspherical)") +Display.object_parentheses(μ::UltrasphericalWeight) = true + + +convert(::Type{JacobiWeight}, μ::AbstractJacobiWeight{T}) where T = + JacobiWeight{T}(jacobi_α(μ), jacobi_β(μ)) +convert(::Type{JacobiWeight{T}}, μ::AbstractJacobiWeight) where T = + JacobiWeight{T}(jacobi_α(μ), jacobi_β(μ)) -function convert(::Type{ChebyshevTWeight}, μ::JacobiWeight{T}) where {T} +function convert(::Type{ChebyshevTWeight}, μ::AbstractJacobiWeight{T}) where T (jacobi_α(μ) ≈ -one(T)/2 && jacobi_β(μ) ≈ -one(T)/2) || throw(InexactError(:convert, ChebyshevTWeight, μ)) ChebyshevTWeight{T}() end - -function convert(::Type{ChebyshevUWeight}, μ::JacobiWeight{T}) where {T} +function convert(::Type{ChebyshevTWeight{T}}, μ::AbstractJacobiWeight) where T + (jacobi_α(μ) ≈ -one(T)/2 && jacobi_β(μ) ≈ -one(T)/2) || throw(InexactError(:convert, ChebyshevTWeight, μ)) + ChebyshevTWeight{T}() +end +function convert(::Type{ChebyshevUWeight}, μ::AbstractJacobiWeight{T}) where T (jacobi_α(μ) ≈ one(T)/2 && jacobi_β(μ) ≈ one(T)/2) || throw(InexactError(:convert, ChebyshevUWeight, μ)) ChebyshevUWeight{T}() end - -function convert(::Type{LegendreWeight}, μ::JacobiWeight{T}) where {T} - (μ.α ≈ 0 && μ.β ≈ 0) || throw(InexactError(:convert, LegendreWeight, μ)) +function convert(::Type{ChebyshevUWeight{T}}, μ::AbstractJacobiWeight) where T + (jacobi_α(μ) ≈ one(T)/2 && jacobi_β(μ) ≈ one(T)/2) || throw(InexactError(:convert, ChebyshevUWeight, μ)) + ChebyshevUWeight{T}() +end +function convert(::Type{LegendreWeight}, μ::AbstractJacobiWeight{T}) where {T} + (jacobi_α(μ) == 0 && jacobi_β(μ) == 0) || throw(InexactError(:convert, LegendreWeight, μ)) LegendreWeight{T}() end - -jacobi_α(μ::LegendreWeight{T}) where {T} = zero(T) -jacobi_β(μ::LegendreWeight{T}) where {T} = zero(T) - -==(μ1::AbstractJacobiWeight, μ2::LegendreWeight) = - jacobi_α(μ1) == jacobi_α(μ2) && jacobi_β(μ1) == jacobi_β(μ2) -==(μ1::LegendreWeight, μ2::AbstractJacobiWeight) = - jacobi_α(μ1) == jacobi_α(μ2) && jacobi_β(μ1) == jacobi_β(μ2) +function convert(::Type{LegendreWeight{T}}, μ::AbstractJacobiWeight) where T + (jacobi_α(μ) == 0 && jacobi_β(μ) == 0) || throw(InexactError(:convert, LegendreWeight, μ)) + LegendreWeight{T}() +end +function convert(::Type{UltrasphericalWeight}, μ::AbstractJacobiWeight{T}) where T + (jacobi_α(μ) == jacobi_β(μ)) || throw(InexactError(:convert, UltrasphericalWeight, μ)) + UltrasphericalWeight{T}(jacobi_α(μ)+one(T)/2) +end +function convert(::Type{UltrasphericalWeight{T}}, μ::AbstractJacobiWeight) where T + (jacobi_α(μ) == jacobi_β(μ)) || throw(InexactError(:convert, UltrasphericalWeight, μ)) + UltrasphericalWeight{T}(jacobi_α(μ)+one(T)/2) +end "The generalised Laguerre measure on the halfline `[0,∞)`." @@ -240,6 +293,8 @@ unsafe_weightfun(μ::LaguerreWeight, x) = laguerre_weightfun(x, μ.α) laguerre_α(μ::LaguerreWeight) = μ.α +==(μ1::LaguerreWeight, μ2::LaguerreWeight) = laguerre_α(μ1) == laguerre_α(μ2) + Base.show(io::IO, μ::LaguerreWeight) = laguerre_α(μ) == 0 ? print(io, "exp(-x)dx (Laguerre)") : print(io, "x^$(laguerre_α(μ))exp(-x)dx (Laguerre)") Display.object_parentheses(μ::LaguerreWeight) = true @@ -262,13 +317,13 @@ Base.show(io::IO, μ::HermiteWeight) = print(io, "exp(-x^2)dx (Hermite)") Display.object_parentheses(μ::HermiteWeight) = true -"The Gaussian measure with weight exp(-|x|^2/2)." +"The Gaussian measure with weight 1/sqrt(2π)^(n/2)*exp(-|x|^2/2)." struct GaussianWeight{T} <: Weight{T} end GaussianWeight() = GaussianWeight{Float64}() gaussian_weightfun(x, ::Type{T} = prectype(x)) where {T} = - 1/(2*convert(T, pi))^(length(x)/2) * exp(-norm(x)^2) + 1/(2*T(pi))^(length(x)/2) * exp(-norm(x)^2) similar(μ::GaussianWeight, ::Type{T}) where {T} = GaussianWeight{T}() isnormalized(μ::GaussianWeight) = true diff --git a/test/runtests.jl b/test/runtests.jl index c194712..22aeb3f 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1 +1,22 @@ -include("test_suite.jl") + +using Test +using DomainSets, HCubature, LinearAlgebra, StaticArrays + +using DomainIntegrals + +include("test_measures.jl") +include("test_integrals.jl") +include("test_gauss.jl") + +@testset "Measures" begin + test_measures() + test_discrete_measures() +end + +@testset "Quadrature rules" begin + test_gauss() +end + +@testset "Selected integrals" begin + test_integrals() +end diff --git a/test/test_integrals.jl b/test/test_integrals.jl index 2faeaf0..aa3aa59 100644 --- a/test/test_integrals.jl +++ b/test/test_integrals.jl @@ -71,12 +71,13 @@ function test_some_integrals() @test abs(integral(strategy, x->1, UnitDisk()) - pi) < 1e-5 end -function test_discrete_measures() +function test_integrals_with_discrete_measures() m1 = DomainIntegrals.GenericDiscreteWeight([0.2], [0.5]) - @test integral(exp, m1) == 0.5*exp(0.2) + @test integral(exp, m1) ≈ 0.5*exp(0.2) end function test_integrals() test_some_integrals() test_readme_examples() + test_integrals_with_discrete_measures() end diff --git a/test/test_measures.jl b/test/test_measures.jl index b7c129b..00dc579 100644 --- a/test/test_measures.jl +++ b/test/test_measures.jl @@ -84,6 +84,28 @@ function test_measures() @test jacobi_α(m8) == 1/2 @test jacobi_β(m8) == 1/2 @test m8 == JacobiWeight(1/2, 1/2) + + m9 = JacobiWeight(0.4, 0.7) + @test support(m9) == (-1..1) + @test !isdiscrete(m9) + @test iscontinuous(m9) + @test !isnormalized(m9) + @test domaintype(m8) == Float64 + @test jacobi_α(m9) == 0.4 + @test jacobi_β(m9) == 0.7 + @test weightfun(m9, 0.22) ≈ (1-0.22)^0.4 * (1+0.22)^0.7 + @test_throws AssertionError JacobiWeight(-2, 0) + @test_throws AssertionError JacobiWeight(0, -1) + + m10 = UltrasphericalWeight(0.3) + @test support(m10) == (-1..1) + @test !isdiscrete(m10) + @test iscontinuous(m10) + @test !isnormalized(m10) + @test domaintype(m10) == Float64 + @test jacobi_α(m10) == 0.3-0.5 + @test jacobi_β(m10) == 0.3-0.5 + @test weightfun(m10, 0.22) ≈ (1-0.22^2)^(0.3-0.5) end function test_discrete_measures() diff --git a/test/test_suite.jl b/test/test_suite.jl deleted file mode 100644 index dc36559..0000000 --- a/test/test_suite.jl +++ /dev/null @@ -1,22 +0,0 @@ - -using Test -using DomainSets, HCubature, StaticArrays - -using DomainIntegrals - -include("test_measures.jl") -include("test_integrals.jl") -include("test_gauss.jl") - -@testset "Measures" begin - test_measures() - test_discrete_measures() -end - -@testset "Quadrature rules" begin - test_gauss() -end - -@testset "Selected integrals" begin - test_integrals() -end