4
$\begingroup$

I have a function that, while the maths itself is unimportant, at certain values it results in a very large number multiplying a very small number. E.g. 10^450000 * 10^-449998. As you can see, this should output the more-sensible number 10^2. However Mathematica is rounding the small number to zero thus the whole calculation breaks down. How can I prevent this? I've played with MinNumber and MachinePrecision but neither seem to fix the issue.

Thanks in advance!

Edit: Including the equation as requested:

$\exp[\cfrac{\omega^2}{\sigma^2}] BesselK[1,\cfrac{\sqrt(\cfrac{\omega}{\sigma^2})}{\sqrt(\cfrac{\sigma^2}{\omega^3})}]$

Equation breaks down for $\sigma<0.01*\omega$ and generally unreliable below $\sigma<0.04*\omega$

Edit2: And the Mathematica code!

bessktot[ω0_, σ_] := BesselK[1, Sqrt[ω0/σ^2]/Sqrt[σ^2/ω0^3]]
expcalctot[ω0_, σ_] := E^(ω0^2/σ^2);
ω0 = 2 \[Pi] 10^12;
σ[BWpc_] := BWpc/100 ω0;
σt = σ[1]
bessktot[ω0, σt]
expcalctot[ω0, σt]
expcalctot[ω0, σt]*bessktot[ω0, σt]

Edit3: Thanks Bob Hanlon, (I can't comment back on your answer yet). Your answer works for certain values input to the bessel. The problem appears to be even more fundamental than I realised. It appears Mathematica can't calculate non-integer $BesselK[1,x]$ functions when $x>741$. Is there a way around this?

$\endgroup$
  • $\begingroup$ On my machine, even with machine precision exponents, 10^45000000000000.*10^-44999999999998. outputs 1.0 * 10^2 from a fresh startup. It should only start rounding to 0 if it runs out of precision during the calculation, one of the intermediate functions does not handle arbitrary precision arithmetic (not common, but some don't), if it's explicitly told to start rounding somewhere, or if something is actually multiplied by 0. We will need a bit more code to diagnose what's actually going on, so please consider posting your function. $\endgroup$ – eyorble Apr 16 at 13:55
  • $\begingroup$ @eyorblade 0^45000000000000. is certainly no a machine precision number: Precision[110^45000000000000.]. $\endgroup$ – Henrik Schumacher Apr 16 at 13:59
  • $\begingroup$ @HenrikSchumacher In a_^b_, the b is machine precision, is it not? The resulting value is not, obviously, because there's any precision tracking at all. I suspect Mathematica knows what it's doing here, in that with exact a the precision isn't lost very quickly at all, but the exponent isn't the source of the uncertainty in this case. $\endgroup$ – eyorble Apr 16 at 14:03
  • $\begingroup$ If the issue is with Mathematica code, please post the Mathematica code rather than a LaTeX representation. $\endgroup$ – JimB Apr 16 at 14:14
  • 2
    $\begingroup$ @Andrew When I approved your edit, I notice that you have two different accounts, both with the same name. I do not know how that happened, but you should be able to combine the two. Then, you will not need the approval of others to edit your own question. $\endgroup$ – bbgodfrey Apr 16 at 17:39
9
$\begingroup$

There was a change in handling underflow in later versions.

$Version

(* "11.3.0 for Mac OS X x86 (64-bit) (March 7, 2018)" *)

Clear["Global`*"]

bessktot[ω0_, σ_] := 
  BesselK[1, Sqrt[ω0/σ^2]/Sqrt[σ^2/ω0^3]];
expcalctot[ω0_, σ_] := E^(ω0^2/σ^2);
ω0 = 2 π 10^12;
σ[BWpc_] := BWpc/100 ω0;
σt = σ[1];

MachinePrecision is insufficient,

(expr = expcalctot[ω0, σt]*bessktot[ω0, σt]) // N

(* 0. *)

Use arbitrary precision

expr // N[#, $MachinePrecision] &

(* 0.01253361135127051 *)

expr // N[#, 20] &

(* 0.012533611351270505734 *)

EDIT: For large input to BesselK you need to control the precision of the input.

BesselK[1, #] & /@ {801., 801.`20, SetPrecision[801., 20]}

(* {0., 5.9781508629496523*10^-350, 5.9781508629496523*10^-350} *)
$\endgroup$
  • 1
    $\begingroup$ Note that, although 801.`20 and SetPrecision[801., 20] are equivalent, in general these approaches can give different results for numeric literals. E.g. 1.1`20 != SetPrecision[1.1, 20]. That's because in the latter case 1.1 is first converted to a machine number then extended the equivalent arbitrary-precision value, while 1.1`20 becomes arbitrary-precision directly at the parsing stage. $\endgroup$ – Ruslan Jul 25 at 17:08
3
$\begingroup$

Another way to go about this is to use the representation of your function in terms of Meijer $G$:

With[{ω = N[2 π 1*^12], σ = N[2 π 1*^12]/100},
     -MeijerG[{{1/2}, {}}, {{-1, 1}, {}}, 2 (ω/σ)^2]/Sqrt[π]]
   0.012533611351270506

Other computing environments, however, usually provide an exponentially-scaled modified Bessel function (e.g. MATLAB) because of numerical problems like these; since Mathematica does not implement these convenient functions, we have to either use arbitrary precision, or make do.

$\endgroup$
  • 1
    $\begingroup$ (That's it for me for now, Wolfram One tends to not work right if the Internet connection is spotty.) $\endgroup$ – J. M. will be back soon Jul 25 at 14:19
  • $\begingroup$ great to see you around again $\endgroup$ – b3m2a1 Jul 26 at 3:29
  • 1
    $\begingroup$ Unfortunately, I can only come sporadically depending on whether I can go to an Internet café, but I hope I can be regular again soon. $\endgroup$ – J. M. will be back soon Jul 28 at 2:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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