I am trying to compute the following integral (numerically, since I couldn't find an analytical expression) for very large $b \gtrapprox 10^9$:
$$ \int_0^1 J_0\left(\exp(b\cdot x)\right)dx $$
($J_0$ is the Bessel function of first kind.)
The problem is that the argument to the Bessel function quickly overflows and so I have to truncate the integrand to something like (on a 64-bit architecture):
$$ \begin{cases} J_0\left(\exp(b\cdot x)\right) &, b\cdot x \leq 709.78\\ 0 & , \textrm{otherwise} \end{cases} $$
However at these cutoff values the value of the integral is not stable within desired accuracy ($\Delta \lessapprox 10^{-16}$; especially when using scipy.integrate the value would be different by a factor of up to $10^2$ between the different methods; the plot was computed with mpmath):
I tried employing mpmath for using arbitrary floating point accuracy however when I set the cutoff larger by only a factor of 100 the computation of the integral becomes prohibitively long and only yields an accuracy of $\approx 10^{-12}$.
So does anyone know a way to compute these kind of integrals? Many thanks in advance!