The natural logarithm of the Bessel function, overflow

Question

I am trying to calculate the logarithm of a modified Bessel function of the second type in MATLAB, i.e. something like that:

log(besselk(nu, Z)) 

where for example

 nu = 750; Z = 1; 

I have a problem because the value of log(besselk(nu, Z)) goes to infinity, because besselk(nu, Z) is infinity. However, log(besselk(nu, Z)) should be really small.

I'm trying to write something like

 f = double(sym('ln(besselk(double(nu), double(Z)))')); 

However, I get the following error:

Error using mupadmex Error in MuPAD: DOUBLE command cannot convert input expression to double array. If the input expression contains a symbolic variable, use the VPA function instead.

Error in sym / double (line 514) Xstr = mupadmex ('symobj :: double', Ss, 0) `;

How to avoid this error?

Answer 1

As Newuff pointed out, DLMF gives asymptotic expansions of K_nu (z) for large nu. From 10.41.2 we find for real positive arguments z:

 besselk(nu,z) ~ sqrt(pi/(2nu)) (ez/(2nu))^-nu 

which gives after some simplification

 log( besselk(nu,z) ) ~ 1/2*log(pi) + (nu-1/2)*log(2nu) - nu(1 + log(z)) 

So this is O (nu log (nu)). Not surprisingly, direct calculation is not performed for nu> 750.

I do not know how accurate this approximation is. Perhaps you can compare it for values ​​where besselk is less than numeric infinity to see if this suits your purpose?

EDIT: I just tried for nu = 750 and z = 1: the above approximation gives 4.7318e + 03, and with the result of horchler we get log (1.02 * 10 ^ 2055) = 2055 * log (10) + log (1.02) = 4.7318 e + 03. Therefore, this is correct for at least 5 significant digits, for nu> = 750 and z = 1! If this is enough for you, it will be much faster than symbolic math.

Answer 2

You are doing a few things wrong. It makes no sense to use double for your two arguments in besselk and convert the output to symbolic. You should also avoid typing the old line into sym . Instead, you want to conditionally evaluate besselk (which will return about 1.02 × 10 ²⁰⁵⁵ much more than realmax ), take the result log symbolically, and then convert it back to double precision.

The following is enough: if one or more input arguments are symbolic, the symbolic version of besselk will be used:

 f = double(log(besselk(sym(750), sym(1)))) 

or in the old string form:

 f = double(sym('log(besselk(750, 1))')) 

If you want to keep your parameters symbolic and evaluate at a later time:

 syms nu Z; f = log(besselk(nu, Z)) double(subs(f, {nu, Z}, {750, 1})) 

Make sure you do not flip the nu and Z values ​​in your math, since large orders ( nu ) are not very common.

Answer 3

Have you tried the integral representation?

Journal [Integration [Cosh [Nu t] / E ^ (Z Cosh [t]), {t, 0, infinity}]]