Natural Logarithm of Bessel Function, Overflow

8
votes

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

log(besselk(nu, Z))

where e.g.

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 small indeed.

I am 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 command: DOUBLE cannot convert the input expression into a 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', S.s, 0)`;

How can I avoid this error?

3
algorithmmatlabmathbessel-functions
Is it overflow? or the result is actually Infinite? You value is bigger than 1.7977e+308, do you need to play with numerical values that big? - Ander Biguri
Side note: it doesn't overflow, but gives Inf. Technically overflowing would be that the number starts from negative again, but that only happens with integers, not floating points. - Ander Biguri
@AnderBiguri Technically, overflowing would be that the number starts from negative again, but that only happens with integers, not floating points. What you're referring to is called "integer wrap around". Overflow is a general term for "hitting the limit", and applies to floating-point numbers as well. - jub0bs
The DLMF gives asymptotic expansions for large arguments in section 10.40 and for large order in section 10.41 which may be suitable for your use case. - njuffa
The answer to this question on Mathoverflow may also be useful. - njuffa

3 Answers

4
votes

As njuffa 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)) (e z/(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 it is O(nu log(nu)). No surprise the direct calculation fails for nu > 750.

I dont know how accurate this approximation is. Perhaps you can compare it for the values where besselk is smaller than the numerical infinity, to see if it fits your purpose?

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

5
votes

You're doing a few things incorrectly. It makes no sense to use double for your two arguments to to besselk and the convert the output to symbolic. You should also avoid the old string based input to sym. Instead, you want to evaluate besselk symbolically (which will return about 1.02×102055, much greater than realmax), take the log of the result symbolically, and then convert back to double precision.

The following is sufficient – when one or more of the input arguments is 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 that you haven't flipped the nu and Z values in your math as large orders (nu) aren't very common.

0
votes

Have you tried the integral representation?

Log[Integrate[Cosh[Nu t]/E^(Z Cosh[t]), {t, 0, Infinity}]]