Matlab numerical integration involving Bessel functions returns NaN

Question

I need to numerically compute integrals such as this (some parameters omitted for simplicity):

$$ \int_{0}^{\infty} e^{-x^2} I_{0}(x) K_{0}(x) \mathrm{d}x $$ where $I_{0}$ and $K_{0}$ denote the modified Bessel functions of the first and second kind, respectively.

As $x \to \infty$, we have $K_{0} \to 0$ and $I_{0} \to \infty$, but the integrand goes to zero rather fast. Yet, instead of smartly truncating the integration interval, Matlab complains about NaN, apparently due to $I_{0}$ diverging to infinity.

>> integral(@(x) exp(-x.^2) .* besseli(0,x) .* besselk(0,x) , 0 , inf)
Warning: Infinite or Not-a-Number value encountered. 
> In funfun/private/integralCalc>iterateScalarValued at 349
  In funfun/private/integralCalc>vadapt at 133
  In funfun/private/integralCalc at 84
  In integral at 89 

I have tried quad, quadgk and quadl too, with similar results. So what is a robust approach to evaluate these kinds of integrals? Or how do you teach Matlab to automatically cut off the tail if it's expected to be very thin?

Note: I am not very familiar with the inner workings of quadrature methods, so I'd prefer a simple generic fix over a super-customized one, if there is.

Answer 1

Try using the optional scale arguments for besseli (documentation) and besselk (documentation), which change how these functions are calculated internally:

integral(@(x)exp(-x.^2).*besseli(0,x,1).*besselk(0,x,1),0,Inf)

which returns 1.161822658007091. The scaling factors for the two Bessel functions cancel each other in this case so you don't need to do anything else. You can increase the accuracy by adjusting the absolute and relative tolerances, e.g.:

f = @(x)exp(-x.^2).*besseli(0,x,1).*besselk(0,x,1);
integral(f,0,Inf,'AbsTol',eps,'RelTol',100*eps)

which returns 1.161822643395565, and agrees more closely with the variable precision result from Mathematica. You can use f = @(x)exp(x.*(1-x)).*besseli(0,x,1).*besselk(0,x); as well.

I'm not sure why integral has a problem with your original formulation, but it may be fundamental to do with the growth and decay rates of your sub-functions and the quadrature routines used. Even Matlab's symbolic and variable precision int did not return a result (and I get errors when trying to use MuPAD's numeric::quadrature directly). You might consider reporting this case to The MathWorks by filing a bug report.