8

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?

jub0bs
  • 60,866
  • 25
  • 183
  • 186
James Green
  • 79
  • 5
  • 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 Sep 09 '15 at 16:28
  • 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 Sep 09 '15 at 16:32
  • >> 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)). – James Green Sep 09 '15 at 16:32
  • 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? – Ander Biguri Sep 09 '15 at 16:34
  • Yes actually I have such a big numbers (Simulation staff...). Have u any idea to deal with this problem? – James Green Sep 09 '15 at 16:41
  • 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: https://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html . – Ander Biguri Sep 09 '15 at 16:43
  • 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`? – Teepeemm Sep 09 '15 at 16:43
  • 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? – Ander Biguri Sep 09 '15 at 16:45
  • @Teepeemm U are right the expression could be `800`, `1000` or even `5000` but small compare to `inf`. – James Green Sep 09 '15 at 16:52
  • 1
    @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](https://en.wikipedia.org/wiki/Integer_overflow#Origin)". Overflow is a general term for "hitting the limit", and [applies to floating-point numbers as well](https://en.wikipedia.org/wiki/IEEE_floating_point#Exception_handling). – jub0bs Sep 09 '15 at 17:06
  • How accurate does your result need to be? Have you looked at https://en.wikipedia.org/wiki/Bessel_function#Asymptotic_forms – Teepeemm Sep 09 '15 at 17:12
  • @Teepeemm That is good approximation only for a small `nu`. In my case I have very big `nu`. – James Green Sep 09 '15 at 17:15
  • 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? – James Green Sep 09 '15 at 17:22
  • 2
    The DLMF gives asymptotic expansions for large arguments in [section 10.40](http://dlmf.nist.gov/10.40) and for large order in [section 10.41](http://dlmf.nist.gov/10.41) which may be suitable for your use case. – njuffa Sep 09 '15 at 17:39
  • 1
    The answer to [this question on Mathoverflow](http://mathoverflow.net/questions/158598/asymptotic-expansion-of-modified-bessel-function-k-alpha) may also be useful. – njuffa Sep 09 '15 at 17:49

3 Answers3

5

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.

horchler
  • 18,384
  • 4
  • 37
  • 73
4

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.

Andreas H.
  • 5,557
  • 23
  • 32
0

Have you tried the integral representation?

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

Louis Ricci
  • 20,804
  • 5
  • 48
  • 62