Devising an analytical approach to solving two infinite (improper) integrals of type $\int_0^\infty f(u) J_\mu (au) J_\nu(bu) \,\mathrm{d}u$

Question

While I was elaborating on a fluid mechanical problem, I came across the two following functions defined in terms of non-trivial infinite (improper) integrals over the wavenumber $u$ in the interval $(0,\infty)$ $$ \psi_\alpha (a,b) = \int_0^\infty \frac{u^\frac{3}{2}}{\left( u^2+\alpha^2 \right)^{\frac{1}{2}}} \, J_0 (au) J_{\frac{3}{2}}(bu) \,\mathrm{d}u $$ and $$ \phi_\alpha (a,b) = \int_0^\infty \left( \frac{u^2+\alpha^2}{u} \right)^{\frac{1}{2}} J_1 (au) J_{\frac{5}{2}}(bu) \,\mathrm{d}u $$ wherein $\alpha > 0$ is a small parameter, and $a>0$ and $b>0$ are the variables. Those two improper integrals arose while solving dual integral equations using the well-established analytical approaches devised by $\,$ S n e d d o n $\,$ and $\,$ C o p s o n.

It can easily be demonstrated that those two integrals are convergent. In particular, for $\alpha = 0$, the corresponding expression can readily be determined analytically. They are found to depend on whether $a<b$ or $a>b$.

In the general case of interest, however, I am unable to figure out how to proceed. I have tried to use the series or integral representations of the Bessel functions and try to make an analytical progress but unfortunately this did not seem to work. 

Any help is highly appreciated.

Thank you very much!

N O T E : Using the change of variables $v=u/\alpha$, $A=\alpha a$, and $B=\alpha b$, the integrands can me made independent of the parameter $\alpha$.

Poisson’s and Related Integrals (10.9.3 in DLMF) $$ J_\nu(z) = \frac{2 \left( \frac{z}{2} \right)^\nu}{\pi^\frac{1}{2} \Gamma \left( \nu + \frac{1}{2} \right)} \int_0^1 \left( 1-t^2 \right)^{\nu-\frac{1}{2}} \cos (zt) \, \mathrm{d}t \, . $$

Answer 1

I tried sketching on the start of a way forward with the Mellin transform for the first integral, i.e.

$$ \int_0^{\infty} \frac{t^{3/2}}{(t^2+\alpha^2)^{1/2}}J_0(at)J_{3/2}(bt)dt $$

Note that this is far from being an actual solution at the moment. It was rather too long to fit into the comments. I'm not sure all the terms I got are correct, I haven't double checked them, but what still remains is the inverse Mellin transform which really doesn't look like too much fun, even though I think it is possible to compute. Maybe someone can take up the calculations to compute the inverse Mellin transform?

I started with trying to write this as an integral on the form of

$$ \int_0^{\infty}h_1\left(\frac{x}{t}\right)h_2(t)\frac{dt}{t} $$

as the Mellin transform of this is $\tilde{h}_1(s)\tilde{h}_2(s)$ where the two are the Mellin transforms of $h_1(t)$ and $h_2(t)$ respectively.

With $h_1(t) = J_{\nu}\left(\frac{2}{\sqrt{t}}\right)$ and $h(t) = f(t)J_{\mu}(at), f(t)=\frac{t^{3/2}}{(t^2+\alpha^2)^{1/2}}$. the Mellin transform of $h_1(t)$ is

$$ \mathcal{M}[h_1(t);s] =\frac{\Gamma(\nu/2-s)}{\Gamma(1+\nu/2-s/2)} $$

so to get $J_\nu(bt) = h_1(2bt^{-2})$, one has

$$ \mathcal{M}[J_\nu(bt);s] = \mathcal{M}[h_1(2bt^{-2});s] = 2^{-s/2-1}b^{-s/2}\frac{\Gamma(\nu/2-s/2)}{\Gamma(s/2+\nu/2+1)} $$

with $\nu=3/2$, this simplifies to $\mathcal{M}[J_\nu(bt);s] = 2^{-s/2-1}b^{-s/2}\frac{\Gamma(3/4-s/2)}{\Gamma(s/2+5/2)}$ with the fundamental strip $-5<\Re(s)<3/2$.

For $h_2(t)$, the Mellin transform is in turn

$$ \tilde{h}_2(s) = \int_0^{\infty} t^{s-1}f(t)J_0(at)dt = \int_0^{\infty} \frac{t^{s+1/2}}{(t^2+\alpha^2)^{1/2}} J_0(at)dt $$

So this is also on the form of a Mellin transform. But this integral is luckily given in 6.567.8 (p. 679) in Table of integrals, series and products, Gradshteyn, Ryzhik, 7th ed. in the form of

$$ \int_0^{\infty} \frac{x^{\rho-1}J_{\nu}(ax)}{(x^2+k^2)^{\mu+1}}dx $$ which is the sum of two hypergeometric functions of ${}_1F_2$. For the parameters in this case, one has that

$$ \tilde{h}_2(s) = \int_0^{\infty} \frac{t^{s+1/2}}{(t^2+\alpha^2)^{1/2}} J_0(at)dt = \frac{\alpha^{s-3/2}\Gamma(s/2+3/4)\Gamma(-s/2-1/4)}{2\sqrt{\pi}}{}_1F_2\left(\frac{s}{2}+\frac{3}{4};\frac{s}{2}+\frac{5}{4},1;\left(\frac{a\alpha}{2}\right)^2 \right) + \frac{a^{-s-1/2}\Gamma(s/2+1/4)}{2^{-s+1/2}\Gamma(3/4-s/2)}{}_1F_2\left(\frac{1}{2};\frac{3}{4}-\frac{s}{2},\frac{5}{4}-\frac{s}{2};\left(\frac{a\alpha}{2} \right)^2 \right) $$

with $-3/2<\Re(s)<-1/2$ if I've managed to put in all terms correctly. Now, as $\tilde{h}_1(s)\tilde{h}_2(s)$ is known, it is "only" a matter of inverting this. And well, if the above is correct, one should then get whatever the original integral computes to. It looks like one would get infinite sums of the ${}_1F_2$ evaluated at each residual coming from the Gamma functions.

Answer 2

You are correct that to get back to the original integral, you would need to take the inverse Mellin transform of h~1(s)h~2(s) h ~ 1(s) h ~ 2(s).

The Gamma function poles will result in exponential terms in the inverse transform. The Mellin-transform approach can be powerful, but the resulting expressions can indeed become quite complicated. Some things to keep in mind:

Make sure you have the correct fundamental strip for convergence of the Mellin transform, and correspondingly the correct inversion contour for the inverse transform Double-check all the details in the transformations, as there are many terms and factors involved. Small errors can creep in easily. Consider also other approaches, like the series expansions or integral representations mentioned earlier. These may lead to more tractable results, especially for small α. If possible, check any analytical results you get against numerical evaluation of the original integral, to ensure the method is working correctly.

So the Mellin transform approach seems like a reasonable thing to try, but do be careful with the details, and consider other methods as well.

Some good options for numerically evaluating integrals like yours include:

Quadrature methods, like Gaussian quadrature. These work well for integrands that are well-behaved and unimodal (have only one peak). Adaptive quadrature, which uses quadrature but refines the subintervals to get higher accuracy where needed. Monte Carlo methods, which use random sampling to estimate the integral. These are very general but can be slow to converge. Romberg integration, which uses Richardson extrapolation to accelerate the convergence of trapezoid rule approximations.

For your integrals, I would suggest trying Gaussian quadrature or adaptive quadrature first. The integrands are well-behaved, but have exponential decay at infinity, so adaptive methods may be efficient. You can find implementations of these numerical integration methods in many numerical libraries, like SciPy or NumPy. Just be sure to evaluate the non-elementary functions (Bessel functions) numerically as well, to high precision. Comparing numerical results to your analytical results will be a good check to ensure the analysis is accurate. Even some digits of agreement would be a good sign, though more digits are better if possible.