3
$\begingroup$

I'm looking for a well-thought quadrature rule for this measure $d\mu(t)=\frac{dt}{t^s}$ for $s\in(0,1)$, the underlying motivation is to compute this integral $$ \lambda^{s-1}=\frac{1}{\Gamma(1-s)}\int_{0}^{\infty} e^{-t\lambda}\frac{dt}{t^s} $$ with $\lambda > 0$, see this reference, Page 5, equation (5), https://arxiv.org/abs/1808.05159. So I'm looking for a way of dealing with the singularity at $t=0$ as well as the limit to $\infty$. Observe that $$ \int_{T}^\infty \frac{dt}{t^s} = \infty $$ I was using the QuadGK.jl package and it is working just fine

λ = 2; s = 0.8
λs = λ^(s-1)

f(t) = exp(-λ*t)/t^s/gamma(1-s)
res, err = quadgk(f,0,Inf)

The predicted error and the real error is the same order of magnitude. However, when using the package for another formula $$ \lambda^s = \frac{1}{\Gamma(-s)}\int_0^\infty (e^{-t\lambda}-1)\frac{dt}{t^{1+s}} $$ same reference as before. I again do the quadrature

λs = λ^s
g(t) = (exp(-λ*t)-1)/t^(1+s)/gamma(-s)
res, err = quadgk(g,0,Inf)

it works just fine, but this time the real error is of order 1e-3 while the predicted error is around 1e-8. And if I put a positive matrix instead of a scalar

A = # some positive definite matrix
h(t) = (exp(-A.*t)-I)/t^(1+s)/gamma(-s)
res, err = quadgk(h,0,Inf)

I get an ERROR: DomainError with 0.5: integrand produced NaN in the interval (0.0, 1.0), and, to my understanding, quadgk does not evaluate the function at $0$ as per the documentation https://juliamath.github.io/QuadGK.jl/stable/. And if I write

res, err = quadgk(h,1e-50,Inf)

the integral blows up. If I augment the lower bound to 1e-20 then I get a norm2 error of 0.24 w.r.t the fractional matrix, if I keep augmenting the lower bound (1e-10, 1e-5), the norm2 error keeps augmenting. I computed the error using the exact eigenvalues and eigenvectors. Has somebody got an intuition to where the instability behind the scenes might be? I would also like to know whether the Gauss-Kronrod quadrature is optimal for this integral, and if I'm doing fine, as I'm relatively new to this :)

EDIT (following a comment): unfortunately does not make a lot of difference, I leave the fully reproducible code here

using QuadGK
using SpecialFunctions 

s = 0.8
nx = 3; Δx = 1/(nx+1) # change nx to whatever
A = -[2. -1 0; -1 2 -1; 0 -1 2]/Δx^2
Χ = zeros(Float64,nx,nx)
for p=1:nx
    Χ[:,p] = sin.(p*pi*(1:nx)*Δx)
 end
Λ = diagm(2*(cos.((1:nx)*pi*Δx).-1)/Δx^2)
As = -Χ*(diagm(diag(-Λ).^s))/Χ
h(t) = (exp(A.*t)-I)/t^(1+s)/gamma(-s)
res1, err1 = quadgk(h,1e-10,Inf)
res21, err21 = quadgk(h,1e-10,1e-5)
res22, err22 = quadgk(h,1e-5,Inf)
res2 = res21 + res22; err2 = err21 + err22
@show res2
@show norm(res2+As)/norm(As)
@show norm(res1+As)/norm(As)
@show norm(res1-res2)
$\endgroup$

1 Answer 1

4
$\begingroup$

You can split the integral up such that you don't get the singularity at $t=0$. Something like this

using QuadGK
using SpecialFunctions 

λ = 2.0; s = 0.8
λs = λ^(s-1)

f(t) = exp(-λ*t)/t^s/gamma(1-s)
@show res, err = quadgk(f,0,Inf)

λs = λ^s
g(t) = (exp(-λ*t)-1)/t^(1+s)/gamma(-s)
@show res, err = quadgk(g,0,Inf)

A = [2 -1 0; -1 2 -1; 0 -1 2] # some positive definite matrix
h(t) = (exp(-A.*t)-I)/t^(1+s)/gamma(-s)

singularity_integral, singularity_error = quadgk(h, 1e-10, 1e-5)
regular_integral, regular_error = quadgk(h, 1e-5, Inf)
@show result = singularity_integral + regular_integral
@show error = singularity_error + regular_error

which will print

(res, err) = quadgk(f, 0, Inf) = (0.8705505365052336, 1.2164780075652479e-8)
(res, err) = quadgk(g, 0, Inf) = (1.7402976255082268, 2.515036776166618e-8)
result = singularity_integral + regular_integral = [1.6837899364901885 -0.7050512093960599 -0.03988519645256732; -0.7050512093960599 1.6439047399317526 -0.7050512098966073; -0.03988519645256732 -0.7050512098966073 1.683789936502625]
error = singularity_error + regular_error = 3.822238569033589e-8

Alternatively, you can also look into something like Cauchy principal values and this example

--Edit: update to new information--

I'm afraid I can't reproduce the same error. When I run your code (I've added using LinearAlgebra for diagm) I get

res2 = [15.354655632950173 -6.419790221101635 -0.3665284746983265; -6.419790221101635 14.988127158828647 -6.4197902215468154; -0.3665284746983265 -6.4197902215468154 15.354655632887033]
norm(res2 + As) / norm(As) = 0.018647374162762107
norm(res1 + As) / norm(As) = 0.01864737561229984
norm(res1 - res2) = 4.9872522079524005e-8
4.9872522079524005e-8
$\endgroup$
1

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.