Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

So I have some variables $\,x_{1},\, x_{2},\, \nu\, =\, 12.654,\, 13.487,\, 0\,$ and the following function:

$\dfrac{(x_{1}\cdot(-BesselK(\nu,x_{1}\cdot125))\cdot BesselI(\nu,x_{2}\cdot125))-(x_{2}\cdot BesselK(\nu+1,x_{1}\cdot 125)\cdot BesselI(\nu+1,x_{2}\cdot 125))}{\dfrac{-1}{125\cdot x_{1}}}$

What's my problem? Well even though the final output is a reasonable number, the output of the Bessel functions are either super-huge or super-small and the calculation requires quad-precision (which is very sluggish). I need to scale down the value going into the Bessel function but I need to be able to get the end result to be the same as if I used quad-precision.

Thank you.

share|cite|improve this question
Is it possible to simply divide by ten within the bessel functions? – Rapid Dec 12 '12 at 9:42
I just checked and a simple division does not work :( I had a feeling it wouldn't be that simple. – Rapid Dec 12 '12 at 9:50
up vote 2 down vote accepted

So it turns out I can use the exponentially scaled modified bessel functions, namely:

$BesselI_{exp}(\nu,x) = \exp(-x)BesselI(\nu,x)\\BesselK_{exp}(\nu,x) = \exp(x)BesselK(\nu,x)$

and so,

$BesselI(\nu,x) = \exp(\log(BesselI_{exp}(\nu,x))+x)\\BesselK(\nu,x) = \exp(\log(BesselK_{exp}(\nu,x))-x)$

These scaled variants must be somehow magically coded into Pythons scipy package to avoid using huge numbers. Very impressive and quite fortunate.

This allows me to modify my equations to use the scaled functions which returns normal sized numbers and then modify the end result to account for the scale.

share|cite|improve this answer

If $125 x_1$ and $125 x_2$ are fairly large compared to $\nu^2$, you might want to use the asymptotic forms of these Bessel functions:

$$ \eqalign{K_\nu(x) &\sim \sqrt{\frac{\pi}{2x}} e^{-x} \left(1+{\frac {4\,{\nu}^{2}-1^2}{1!\; 8x}}+{\frac {(4 \nu^2-1^2)(4 \nu^2-3^2)}{2!\; (8x)^{2}}}\right.\cr &\left.+{\frac {(4 \nu^2 - 1^2)(4 \nu^2 - 3^2)(4 \nu^2 - 5^2)}{3! \;(8x)^{3}}}+{\frac {(4 \nu^2 - 1^2)(4 \nu^2 - 3^2)(4 \nu^2 - 5^2)(4 \nu^2 - 7^2)}{4! \;(8x)^{4}}}+\ldots \right)\cr I_\nu(x) &\sim \frac{e^x}{\sqrt{2\pi x}}\left(1-{\frac {4\,{\nu}^{2}-1^2}{1!\; 8x}}+{\frac {(4 \nu^2-1^2)(4 \nu^2-3^2)}{2!\; (8x)^{2}}}\right.\cr &\left.-{\frac {(4 \nu^2 - 1^2)(4 \nu^2 - 3^2)(4 \nu^2 - 5^2)}{3! \;(8x)^{3}}}+{\frac {(4 \nu^2 - 1^2)(4 \nu^2 - 3^2)(4 \nu^2 - 5^2)(4 \nu^2 - 7^2)}{4! \;(8x)^{4}}}-\ldots \right)\cr}$$

share|cite|improve this answer
Hi, thank you for this, the only problem is that it's the value the Bessel function produces that is too big. e.g. $BesselI(0,1000) = 2.486\dots\times 10^{432}$ unfortunately this number is larger than double precision floats can hold, hence the need to scale down the values inputted to the Bessel function so the output is not so large – Rapid Dec 12 '12 at 11:29

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.