# Numerical integration of modified bessel function

I need to compute the following integral:

NIntegrate[ BesselI[-nu, k x]/x ,{x, r1, r}]


in which nu=-(2m-1)/2 and I have to increase m, k is a complex number,r1=1, r=1.1 as an example(r>r1).

When I run my code it generates errors and unreasonable answers for some values of m. I know that there's a closed form answer but I have to compute it numerically. How can I improve my results?

• You are reposting this and I doubt you did not notice the code formatting in the thread where you originally asked this question, so please refer to the help centre and format your post accordingly. – Sektor Nov 4 '14 at 20:36

You can integrate as follows.

Integrate[BesselI[-nu, k*x]/x, {x, r1, r},
Assumptions :> {k \[Element] Complexes, r1 \[Element] Reals,
r \[Element] Reals, nu \[Element] Reals, r1 > 0, r > r1}


The result is a complicated expression in terms of Gamma and HypergeometricPFQRegularized functions. Nevertheless, it can be evaluated at large values of m. For example,

Block[{k=2.0536*10^(-6) + 2.3066*I, nu, r1=1.0, r=1.1, m=30},
nu = -(2m - 1) / 2.0;
(* long expression from NIntegrate command above *)
]

(* -4.96714*10^-31 + 4.9674*10^-31 I *)


For numeric integration you need to assign numeric values for all parameters (e.g., k) and the product of k and x must include a space (k x) or an asterisk (k*x).

With[{k = 1 + I, r1 = 1, r = 1.1},
Table[
nu = -(2 m - 1)/2.;
NIntegrate[BesselI[-nu, k x]/x, {x, r1, r}],
{m, 6}]]


{0.0691516 + 0.0652726 I, 0.00791225 + 0.0454624 I, -0.00714297 + 0.0116356 I, -0.00279089 + 0.000771435 I, -0.000421279 - 0.000226613 I, -0.0000195765 - 0.0000616956 I}

EDIT:

The integral can be done analytically

f[k_, r1_, r_, m_] =
Assuming[{r > r1 > 0}, Integrate[BesselI[(2 m - 1)/2, k x]/x, {x, r1, r}]]


For relatively large values of m the values of the integral are essentially zero.

With[{r1 = 1, r = 1.1},
Prepend[
Table[{m,
f[1 + I, r1, r, m] // N,
f[2.0536 10^-6 + 2.3066 I, r1, r, m] // N},
{m, 30, 40, 2}],
{"m", "k = 1 + I", "k = 2.0536*10^-6 + 2.3066I"}]] //
Grid[#, Frame -> All] &


• Thank you so much but when I increase m to 30,40,... ,I can't solve the errors. – shahrzad mirhosseini Nov 5 '14 at 18:55
• I guess it's highly oscillatory and even with LevineRule the errors don't vanish. If k is real there is no problem but when I have complex numbers... – shahrzad mirhosseini Nov 5 '14 at 18:57
• [Kappa] = 2.0536 10^-6 + 2.3066 I – shahrzad mirhosseini Nov 5 '14 at 18:58