# Integration of products of spherical Hankel functions yield different results

Background

I am trying to solve the following integral $$J_{nm}^{(j)}(x) =\int_1^\infty y^{(-j)}h_n^{(1)}(x y)j_m(x^* y)\mathrm{d}y$$

where $$x$$ is complex valued, the asterisk denotes the complex conjugate, $$h_n^{(1)}$$ is the spherical Hankel function of the first kind, and $$j_m$$ is the spherical Bessel function of the first kind.

In the end, I will need numerical values for the integrals but in order to find an appropriate quadrature I was playing around with these integrals in Mathematica.

For certain values for $$n$$, $$m$$, $$j$$, and $$x$$ Mathematica will be able to evaluate the integrals analytically using the Mathematica's Integrate-function. The problem is that if I use the numerical quadrature NIntegrate the result will be very different. This might not be surprising, since Bessel functions are known to be notoriously difficult to numerically integrate. However, the integration also yields incorrect results on very small subintervals where integrand of $$H_{nm}^{(j)}$$ appears to be linear. So the reason does not seem to be the usual problems with the oscillations of the integrand.

An example:

The problem appears in particular if the imaginary part of $$x$$ is small.

Integrating on the subinterval $$[1, 1.1]$$ for $$x = (1 + i) \times 10^{-8}, n = 0, m = 1, j = 0$$

yields for the analytic integration $$\int_1^{1.1} y^{(-j)}h_n^{(1)}(x y)j_m(x^* y)\mathrm{d}y = 6.42699\times 10^7 + 6.42699\times10^7 i$$ while the numerical integration yields $$\int_1^{1.1} y^{(-j)}h_n^{(1)}(x y)j_m(x^* y)\mathrm{d}y = -0.0333333 - 3.5\times10^{-10} i.$$

Looking at the graph of the kernel, the result from the numerical quadrature looks much more reasonable.

Question

I looked the problem up in the literature and there are many articles discussing these types of integrals. However, the specifics never quite match my case. Now I have some questions:

• Where is the difference coming from and which result can I trust? Why?
• Is my assumption correct that integration of a complex integrand along the real axis in this case means that I do not have to deal with complex analysis (calculus of residues)?
• Why does the problem arise only for certain values of $$j$$, $$n$$, $$m$$, and $$x$$. In particular, why does it appear especially if $$\mathrm{Im(x)} \ll 1$$?

Minimal Working Example

Here is the code:

J[y_, x_, n_, m_, j_] := y^(-j) * SphericalHankelH1[n, x * y] * SphericalBesselJ[m, Conjugate[x] * y];
Integrate[J[z, (1+I) * 10^(-8), 0, 1, 0], {z, 1, 1.1}]
NIntegrate[J[z, (1+I) * 10^(-8), 0, 1, 0], {z, 1, 1.1}]


Thank you!

• Integrate[J[z, (1+I) * 10^(-8), 0, 1 0], {z, 1, 1.1}]; should be Integrate[J[z, (1+I) * 10^(-8), 0, 1, 0], {z, 1, 1.1}]; Commented Jan 3, 2022 at 18:16
• Thank you @bbgodfrey. I changed it. Also, I cannot use the  formatting to make the math nicer to view. Commented Jan 3, 2022 at 18:20
• It probably has to do with the inadequate precision due to extremely large/small terms in the integral. Take a look at the analytic integral (int = Integrate[J[z, (1 + I)*10^(-8), 0, 1, 0], z]) and its plot Plot[Evaluate@Re@int, {z, .9, 1.2}]. The plot looks much nices if you increase the precision Plot[Evaluate@Re@int, {z, .9, 1.2}, WorkingPrecision -> 30]. Commented Jan 3, 2022 at 18:22
• Therefore, increase the precision of the boundary in the integration: Integrate[J[z, (1 + I)*10^(-8), 0, 1, 0], {z, 150, 1.150}]. This yields -0.0333333329833333333333333 - 3.499999970577778*10^-10 I. Commented Jan 3, 2022 at 18:29
• @Domen: thank you very much; also to the others.Your help was very fast and solved/answered my question. Commented Jan 3, 2022 at 20:48

Integrate is a symbolic solver and sometimes messes up numerically when inexact coefficients are used.

Using a high-enough precision to avoid round-off error:

Integrate[jj[z, (1 + I)*10^(-8), 0, 1, 0], {z, 1, 1.150}]
(*  -0.0333333329833333333333333 - 3.499999970577778*10^-10 I  *)


Using N[] on an exact result; note that machine-precision is insufficient to evaluate the result accurately:

res = Integrate[jj[z, (1 + I)*10^(-8), 0, 1, 0], {z, 1, 11/10}];
N[res]
N[res, 16]
(*
1.32231*10^7 + 1.32231*10^7 I
-0.03333333298333333 - 3.5000000*10^-10 I
*)
`
• Oh, while I was playing, Domen gave a similar answer in the comments. Commented Jan 3, 2022 at 18:45