Natural Logarithm of Bessel Function, Overflow

Question

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?

Plump

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? – Zellers

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. – Zellers

>> besselk(750,1) is bigger then 1.7977e+308, however logarithm of besselk(750,1) is small number. I am calculating log(besselk(nu, Z)). – Plump

That is the most ambiguous comment ever! yeah, that gives Inf, we know, you put that in your question. But that is the result! A number bigger than 1.7977e+308. ill ask again: do you need to play with numerical values that big? Why would you? – Zellers

Yes actually I have such a big numbers (Simulation staff...). Have u any idea to deal with this problem? – Plump

I do not know how to deal with it, but I know that numerical stuff using this big numbers will bite your ass. The numerical errors you will get due to floating point errors are going to be HUGE. I highly disencourage you to use numbers this big in simulations. Read carefully: docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html . – Zellers

If bessel is bigger than 1.8e+308, then its ln would be bigger than 700. That is not small; you have a mistake somewhere. Are you mixing nu and Z? – Occident

You are right, log(besselk(nu, Z)) should be small. However, I have no idea how to solve your problem.... Isn't there a way you can know the order of besselk(nu, Z) without computing it? – Zellers

@Occident U are right the expression could be 800, 1000 or even 5000 but small compare to inf. – Plump

@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. – Dogcart

How accurate does your result need to be? Have you looked at en.wikipedia.org/wiki/Bessel_function#Asymptotic_forms – Occident

@Occident That is good approximation only for a small nu. In my case I have very big nu. – Plump

Guys I am trying to write something like that: f = double(sym('ln(besselk(double(nu), double(Z)))'));. But 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); Is there any way to avoid this error? – Plump

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. – Howzell

The answer to this question on Mathoverflow may also be useful. – Howzell

Answer 1

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.

Alteration

Thanks. I will try to use this approximation, however I am not sure whether the the precision will be enough or not... – Plump

@JamesGreen: I just tried it for nu=750 and z=1. See my edit. – Alteration

Answer 2

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×10²⁰⁵⁵, 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.

Heavyduty

Answer 3

Have you tried the integral representation?

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

Molech

Thanks for suggestion. Integration is very slow, I have to do lots of iterations.... – Plump

Answer 4

Here is a recursive implementation in python of the logarithm of the Bessel Function of the second kind. More information are in the technical report.

Bream