quadgk AbsTol/RelTol parameters combinations

2 次查看（过去 30 天）
Alejandro
Alejandro 2025-4-14
评论： Alejandro 2025-4-19
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:
  2 个评论
Torsten
Torsten 2025-4-15
What is "the desired result" ? Do you have integral values of high precision to compare with ?
Alejandro
Alejandro 2025-4-15
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.

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

Torsten
Torsten 2025-4-15
编辑：Torsten 2025-4-15
Ran in:
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
  3 个评论
显示 1更早的评论
Torsten
Torsten 2025-4-16
编辑：Torsten 2025-4-16
Ran in:
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
Alejandro 2025-4-19
Thanks for your useful feedback Torsten. Will generate the values and compere against my reference tables.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Special Functions 的更多信息

标签

产品


版本

R2024b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by