I have been trying to evaluate the following integral involving modified Bessel functions:

\begin{align} \int r I_1 (r) K_1(r) dr \end{align}

This integral has an explicit expression given in the Wolfram documentation in terms of other Bessel functions. Granted, this expression has divergent components for integer-order Bessel functions but still has a well-defined limit at least for the above case. After performing these manually, one gets that the integral evaluates to

1/2 (-1 + r BesselI[0, r] (r BesselK[0, r] + BesselK[1, r]) + 
   r BesselI[1, r] (-BesselK[0, r] + r BesselK[1, r]))

However, when evaluating in Mathematica, this returns

In[1]:= Integrate[r BesselI[1, r] BesselK[1, r], r]

Out[1]= MeijerG[{{1, 3/2}, {}}, {{1, 2}, {0, 0}}, r, 1/2]/(4 Sqrt[\[Pi]])

Furthermore, this Meijer-G function refuses to simplify with FunctionExpand.

Is there a way get Mathematica to simplify the output of the integration? As far as I can see, the Mathematica internals should already know how to perform this simplification (or at least return the indefinite integral not in Meijer-G form).


Furthermore, if one tries to evaluate

Integrate[z BesselI[\[Nu], a z] BesselK[\[Nu], a z], z]

as in the docs, Mathematica doesn't actually apply the rule discussed above, instead expressing in terms of hypergeometricPFQ functions. Why isn't this integration rule ever being used?

  • 1
    $\begingroup$ Moreover, Series[MeijerG[{{1, 3/2}, {}}, {{1, 2}, {0, 0}}, r, 1/2]/( 4 Sqrt[\[Pi]]), {r, 0, 2}] produces two "Power::infy: Infinite expression 1/0 encountered." and $$ \frac{G_{2,4}^{2,2}\left(0,\frac{1}{2}| \begin{array}{c} 1,\frac{3}{2} \\ 1,2,0,0 \\ \end{array} \right)}{4 \sqrt{\pi }}+\text{ComplexInfinity} r+\text{ComplexInfinity} r^2+O\left(r^3\right).$$ $\endgroup$
    – user64494
    Aug 23, 2023 at 5:45

1 Answer 1


You could try

 Integrate[r BesselI[n,a r] BesselK[n,a r],r] //FunctionExpand, {n,1,0}] 
(* r BesselI[1, a r] BesselK[1, a r]  *)


By trying different values of n the general case seems to be

  • $\begingroup$ Perhaps you can shed some light on how you obtained the bottom expression? At least, that's not what I get from naive FunctionExpand and FullSimplify routines on the generic form of the integral. $\endgroup$
    – Aaron
    Aug 24, 2023 at 3:31
  • $\begingroup$ @Aaron if the limit n-> integer evaluates to indeterminate one can try the series approach as I did above. Then I tried for different n, integer and non integer and checked the derivative symbolically and numerically. $\endgroup$
    – Andreas
    Aug 24, 2023 at 8:44
  • $\begingroup$ the term -(n/(2*a^2)) may be omitted as it is just an integration constant $\endgroup$
    – Andreas
    Aug 24, 2023 at 12:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.