quadgk AbsTol/RelTol parameters combinations

Dear network.

I am having trouble getting the desired result of an integral involving Bessel functions Jo and Yo.

Need your help with a powerful set of combinations of the AbsTol/RelTol parameters that will help me get a low-error result

This is the equation I am trying to solve, with t as a parameter:

[Torsten]

What is "the desired result" ? Do you have integral values of high precision to compare with ?

[Alejandro]

Hi Torsten, yes.

I have figures from various papers and books to compare with.

The current results I am obtaining in MATLAB using either the quadgk or integral commands are off by +- 10%, which requires an optimization of the AbsTol/RelTol parameters.

Answer 1

[Torsten]

umin = 1e-16;
f = @(t,u) exp(-t*u.^2)./(u.*(besselj(0,u).^2+bessely(0,u).^2));
g = @(u) pi/2 * atan((2*double(eulergamma)-log(4)+2*log(u))/pi);
qD = @(t) 1 + 4/pi^2*( g(umin) + quadgk(@(u)f(t,u),umin,Inf) );
format long
t = 0.1:0.1:10;
plot(t,arrayfun(@(t)qD(t),t))
xlabel('t')
ylabel('qD')
grid on

[Alejandro]

Thanks for your help. I see that you split the integral into two terms.

I will review in detail to understand how you came up with the g-term integral for short values of u.

[Torsten]

Consider

syms u
f = u*(bessely(0,u)^2+1);

f = 

series(f)

ans = 

g = u*(4*(eulergamma-log(sym('2'))+log(u))^2/sym(pi)^2+1)

g = 

●

int(1/g)

ans = 

Limit for int(1/g) as u -> 0+ is pi/2 * atan(-Inf) = -pi^2/4.

Thus for f(t,u) = exp(-t*u.^2)./(u.*(besselj(0,u).^2+bessely(0,u).^2)) I computed

int(f,0,Inf) = int(f,0,umin) + int(f,umin,Inf) ~ int(1/g,0,umin) + int(f,umin,Inf) = pi/2*atan((2*eulergamma-log(4)+2*log(umin))/pi) + pi^2/4 + int(f,umin,Inf)

Now multiply by 4/pi^2 to get qD.

[Alejandro]

Thanks for your useful feedback Torsten. Will generate the values and compere against my reference tables.