From bc46bccbddd75a82bf9f4a44621924b6e48aeedc Mon Sep 17 00:00:00 2001 From: Federico Cerisola Date: Sat, 28 Oct 2023 15:56:07 +0100 Subject: [PATCH] Refactored singular integrals into their own sub-module --- src/SingularIntegrals.jl | 93 +++++++++++++++++++++++++++++++++++++++ src/SpectralDensities.jl | 14 +++++- src/WeakCoupling.jl | 94 ---------------------------------------- test/runtests.jl | 8 ++-- 4 files changed, 109 insertions(+), 100 deletions(-) create mode 100644 src/SingularIntegrals.jl diff --git a/src/SingularIntegrals.jl b/src/SingularIntegrals.jl new file mode 100644 index 0000000..92e069f --- /dev/null +++ b/src/SingularIntegrals.jl @@ -0,0 +1,93 @@ +""" + cauchy_quadgk(g, a, b, x0=0; kws...) + +Computes the Cauchy principal value of the integral of a function `g` +with singularity at `x0` over the interval `[a, b]`, that is +```math +\\mathrm{P.V.}\\int_a^b\\frac{g(x)}{x-x_0}\\mathrm{d}x +``` +Note that `x0` must be contained in the interval `[a, b]`. +The actual integration is performed by the `quadgk` method of the +`QuadGK.jl` package and the keyword arguments `kws` are passed +directly onto `quadgk`. + +Note: If the function `g` contains additional integrable singularities, +the user should manually split the integration interval around them, +since currently there is no way of passing integration break points +onto `quadgk`. + +## Arguments +- `g`: The function to integrate. +- `a`: The lower bound of the interval. +- `b`: The upper bound of the interval. +- `x0`: The location of the singularity. Default value is the origin. +- `kws...`: Additional keyword arguments accepted by `quadgk`. + +## Returns +A tuple `(I, E)` containing the approximated integral `I` and +an estimated upper bound on the absolute error `E`. + +## Throws +- `ArgumentError`: If the interval `[a, b]` does not include the singularity `x0`. + +## Examples +```julia +julia> cauchy_quadgk(x -> 1/(x+2), -1.0, 1.0) +(-0.549306144334055, 9.969608472104596e-12) + +julia> cauchy_quadgk(x -> x^2, 0.0, 2.0, 1.0) +(4.0, 2.220446049250313e-16) +``` +""" +function cauchy_quadgk(g, a, b, x0=zero(promote_type(typeof(a),typeof(b))); kws...) + a < x0 < b || throw(ArgumentError("domain must include the singularity")) + g₀ = g(x0) + g₀int = (b-x0) == -(a-x0) ? zero(g₀) : g₀ * log(abs((b-x0)/(a-x0))) / (b - a) + return quadgk(x -> (g(x)-g₀)/(x-x0) + g₀int, a, x0, b; kws...) +end + +""" + hadamard_quadgk(g, g′, a, b, x0=0; kws...) + +Computes the Hadamard finite part of the integral of a function `g` +with a singularity at `x0` over the interval `[a, b]`, that is +```math +\\mathcal{H}\\int_a^b\\frac{g(x)}{(x-x_0)^2}\\mathrm{d}x +``` +Note that `x0` must be contained in the interval `[a, b]`. +The actual integration is performed by the `quadgk` method of the +`QuadGK.jl` package and the keyword arguments `kws` are passed +directly onto `quadgk`. + +Note: If the function `g′` contains additional integrable singularities, +the user should manually split the integration interval around them, +since currently there is no way of passing integration break points +onto `quadgk`. + +## Arguments +- `g`: The function to integrate. +- `g′`: The derivative of the function `g`. +- `a`: The lower bound of the interval. +- `b`: The upper bound of the interval. +- `x0`: The location of the singularity. Defaults value is the origin. +- `kws...`: Additional keyword arguments accepted by `quadgk`. + +## Returns +A tuple `(I, E...)` containing the approximated Hadamard finite part integral `I` and +an additional error estimate `E` from the quadrature method. + +## Throws +- `ArgumentError`: If the interval `[a, b]` does not include the singularity `x0`. + +## Examples +```julia +julia> hadamard_quadgk(x -> log(x+1), x -> 1/(x+1), 0.0, 2.0, 1.0) +(-1.6479184330021648, 9.969608472104596e-12) +``` +""" +function hadamard_quadgk(g, g′, a, b, x0=zero(promote_type(typeof(a),typeof(b))); kws...) + a < x0 < b || throw(ArgumentError("domain must include the singularity")) + I0 = g(a)/(a-x0) - g(b)/(b-x0) + I1 = cauchy_quadgk(g′, a, b, x0; kws...) + return (I0 + I1[1], I1[2:end]...) +end \ No newline at end of file diff --git a/src/SpectralDensities.jl b/src/SpectralDensities.jl index 8d957f1..f611681 100644 --- a/src/SpectralDensities.jl +++ b/src/SpectralDensities.jl @@ -39,16 +39,26 @@ include("InversePolyKernelSD.jl") export InversePolyKernelSD +module SingularIntegrals + using QuadGK + + include("SingularIntegrals.jl") + + export cauchy_quadgk, hadamard_quadgk +end + +export SingularIntegrals + module WeakCoupling using ForwardDiff using QuadGK using ..SpectralDensities + using ..SingularIntegrals include("WeakCoupling.jl") export weak_coupling_Δ, weak_coupling_Δprime, - weak_coupling_Σ, weak_coupling_Σprime, - cauchy_quadgk, hadamard_quadgk + weak_coupling_Σ, weak_coupling_Σprime end export WeakCoupling diff --git a/src/WeakCoupling.jl b/src/WeakCoupling.jl index da542cd..660c66d 100644 --- a/src/WeakCoupling.jl +++ b/src/WeakCoupling.jl @@ -107,98 +107,4 @@ function weak_coupling_Σprime(J::AbstractSD, ωB) g(ω) = 4*J(ω)*abs(ωB)*ω/(ω + abs(ωB))^2 g′(ω) = ForwardDiff.derivative(g,ω) return hadamard_quadgk(g, g′, zero(ωB), Inf, abs(ωB))[1] -end - -""" - cauchy_quadgk(g, a, b, x0=0; kws...) - -Computes the Cauchy principal value of the integral of a function `g` -with singularity at `x0` over the interval `[a, b]`, that is -```math -\\mathrm{P.V.}\\int_a^b\\frac{g(x)}{x-x_0}\\mathrm{d}x -``` -Note that `x0` must be contained in the interval `[a, b]`. -The actual integration is performed by the `quadgk` method of the -`QuadGK.jl` package and the keyword arguments `kws` are passed -directly onto `quadgk`. - -Note: If the function `g` contains additional integrable singularities, -the user should manually split the integration interval around them, -since currently there is no way of passing integration break points -onto `quadgk`. - -## Arguments -- `g`: The function to integrate. -- `a`: The lower bound of the interval. -- `b`: The upper bound of the interval. -- `x0`: The location of the singularity. Default value is the origin. -- `kws...`: Additional keyword arguments accepted by `quadgk`. - -## Returns -A tuple `(I, E)` containing the approximated integral `I` and -an estimated upper bound on the absolute error `E`. - -## Throws -- `ArgumentError`: If the interval `[a, b]` does not include the singularity `x0`. - -## Examples -```julia -julia> cauchy_quadgk(x -> 1/(x+2), -1.0, 1.0) -(-0.549306144334055, 9.969608472104596e-12) - -julia> cauchy_quadgk(x -> x^2, 0.0, 2.0, 1.0) -(4.0, 2.220446049250313e-16) -``` -""" -function cauchy_quadgk(g, a, b, x0=zero(promote_type(typeof(a),typeof(b))); kws...) - a < x0 < b || throw(ArgumentError("domain must include the singularity")) - g₀ = g(x0) - g₀int = (b-x0) == -(a-x0) ? zero(g₀) : g₀ * log(abs((b-x0)/(a-x0))) / (b - a) - return quadgk(x -> (g(x)-g₀)/(x-x0) + g₀int, a, x0, b; kws...) -end - -""" - hadamard_quadgk(g, g′, a, b, x0=0; kws...) - -Computes the Hadamard finite part of the integral of a function `g` -with a singularity at `x0` over the interval `[a, b]`, that is -```math -\\mathcal{H}\\int_a^b\\frac{g(x)}{(x-x_0)^2}\\mathrm{d}x -``` -Note that `x0` must be contained in the interval `[a, b]`. -The actual integration is performed by the `quadgk` method of the -`QuadGK.jl` package and the keyword arguments `kws` are passed -directly onto `quadgk`. - -Note: If the function `g′` contains additional integrable singularities, -the user should manually split the integration interval around them, -since currently there is no way of passing integration break points -onto `quadgk`. - -## Arguments -- `g`: The function to integrate. -- `g′`: The derivative of the function `g`. -- `a`: The lower bound of the interval. -- `b`: The upper bound of the interval. -- `x0`: The location of the singularity. Defaults value is the origin. -- `kws...`: Additional keyword arguments accepted by `quadgk`. - -## Returns -A tuple `(I, E...)` containing the approximated Hadamard finite part integral `I` and -an additional error estimate `E` from the quadrature method. - -## Throws -- `ArgumentError`: If the interval `[a, b]` does not include the singularity `x0`. - -## Examples -```julia -julia> hadamard_quadgk(x -> log(x+1), x -> 1/(x+1), 0.0, 2.0, 1.0) -(-1.6479184330021648, 9.969608472104596e-12) -``` -""" -function hadamard_quadgk(g, g′, a, b, x0=zero(promote_type(typeof(a),typeof(b))); kws...) - a < x0 < b || throw(ArgumentError("domain must include the singularity")) - I0 = g(a)/(a-x0) - g(b)/(b-x0) - I1 = cauchy_quadgk(g′, a, b, x0; kws...) - return (I0 + I1[1], I1[2:end]...) end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 06c6199..631efee 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -32,10 +32,10 @@ using Test @test Q(Jpoly_gauss) ≈ reorganisation_energy(Jpoly_gauss) end - @testset "Weak coupling" begin - @test WeakCoupling.cauchy_quadgk(x -> 1/(x+2), -1.0, 1.0)[1] ≈ -log(3)/2 - @test WeakCoupling.cauchy_quadgk(x -> x^2, 0.0, 2.0, 1.0)[1] ≈ 4.0 - @test WeakCoupling.hadamard_quadgk(x -> log(x+1), x -> 1/(x+1), 0.0, 2.0, 1.0)[1] ≈ -3*log(9)/4 + @testset "Singular Integrals" begin + @test SingularIntegrals.cauchy_quadgk(x -> 1/(x+2), -1.0, 1.0)[1] ≈ -log(3)/2 + @test SingularIntegrals.cauchy_quadgk(x -> x^2, 0.0, 2.0, 1.0)[1] ≈ 4.0 + @test SingularIntegrals.hadamard_quadgk(x -> log(x+1), x -> 1/(x+1), 0.0, 2.0, 1.0)[1] ≈ -3*log(9)/4 end @testset "Memory kernels" begin