9
$\begingroup$

I'm having some troubles with the following code I wrote in MMA 10 some time ago:

Clear[g, a, b, c, amp]

g[x_] := I (2 Sqrt[E])/(Pi*b) Exp[-(x - c/2)^2/(2 b^2)] Exp[I (2 Pi)/a (x)] (x - c/2);

amp[x_, k_] = Integrate[-I*g[xd] Exp[-I k xd ], {xd, 0, x}];

Plot3D[
 Abs[5 amp[x, k*50 + shift] /. {a -> 0.01, c -> 1, b -> 1/8, 
    shift -> 2 Pi/0.01}],
 {x, 0, 1}, {k, -1, 1},
 PlotRange -> {{0, 1}, {-1, 1}, {0, 1}},
 LabelStyle -> {24, Black},
 BoxStyle -> Black,
 TicksStyle -> Black,
 PlotPoints -> 50,
 ImageSize -> 500]

The output in MMA 10 looks like the following plot:

enter image description here

However, in MMA 11.3 I get this error message and this output plot:

General::munfl: Exp[-5704.68-289.16 I] is too small to represent as a normalized machine number; precision may be lost.

General::munfl: Exp[-5704.68-289.16 I] is too small to represent as a normalized machine number; precision may be lost.

enter image description here

It looks like a precision problem, but setting WorkingPrecision to different values doesn't help. I'm sure I've read something about it on stackexchange, but I can't remember where (it might be something related to the fact that MMA 11.3 doesn't work with arbitrary precision automatically, but I might be talking nonsense...)

$\endgroup$
  • 2
    $\begingroup$ Hm. I also observed that Mathematica 11.3 throws General::munflmuch more often than the version before... $\endgroup$ – Henrik Schumacher Aug 28 '18 at 13:07
  • $\begingroup$ I found something related in mathematica.stackexchange.com/questions/170416/… : the answers are a bit annoying though, it looks like there is no way of having old codes working again without rewriting parts of them $\endgroup$ – Fraccalo Aug 28 '18 at 13:58
  • 2
    $\begingroup$ WorkingPrecision will work if you replace the 0.01 with exact versions. $\endgroup$ – Brett Champion Aug 28 '18 at 14:38
  • 1
    $\begingroup$ @BrettChampion I tried WorkingPrecision->Infinity using 1/100 instead of 0.01 but it still doesn't work. It works with WorkingPrecision->50, 20, 10 and 5 (which confuses me a bit...) Also, this is a toy case, while in a realistic code I would have very weird numbers for a,b,c (e.g one value could be 0.00004622): how do I replace them with the corresponding exact version? $\endgroup$ – Fraccalo Aug 28 '18 at 14:47
  • 3
    $\begingroup$ You can't use Infinity for the working precision. For your weird numbers, just use Rationalize before plotting, e.g., Plot3D[Evaluate[Rationalize[.., 0]], ..] $\endgroup$ – Carl Woll Aug 28 '18 at 15:56
4
$\begingroup$

Turning a combination of great comments into an answer. The cause is Mathematica no longer checking for underflow since 11.3. Workaround is to turn your function description to a numerically exact form using Rationalize and making Plot3D to use arbitrary precision (which is no-longer the default behaviour) by explicitly setting WorkingPrecision -> 5.

rr = Rationalize[
  Abs[5 amp[x, k*50 + shift] /. {a -> 0.01, c -> 1, b -> 1/8, 
     shift -> 2 Pi/0.01}], 0];
Plot3D[rr, {x, 0, 1}, {k, -1, 1}, 
 PlotRange -> {{0, 1}, {-1, 1}, {0, 1}}, LabelStyle -> {24, Black}, 
 BoxStyle -> Black, TicksStyle -> Black, PlotPoints -> 50, 
 ImageSize -> 500, WorkingPrecision -> 5]

enter image description here

If I do not Rationalize before hand it gives a warning, but still calculates using WorkingPrecision. However, without Rationalize it seems to require larger WorkinPrecision for a reasonable result. I have an example, where even WorkingPrecision->1000 does not give the same plot, as I get with v11.1. If someone can explain it, I am happy to extend the answer.

I would also expect

rr = N[Abs[
   5 amp[x, k*50 + shift] /. {a -> 0.01, c -> 1, b -> 1/8, 
     shift -> 2 Pi/0.01}], 25]
Plot3D[rr, {x, 0, 1}, {k, -1, 1}, 
 PlotRange -> {{0, 1}, {-1, 1}, {0, 1}}, LabelStyle -> {24, Black}, 
 BoxStyle -> Black, TicksStyle -> Black, PlotPoints -> 50, 
 ImageSize -> 500]

to work, but it does not. It really seems, that WorkingPrecision and Rationalize are required.

$\endgroup$
  • 3
    $\begingroup$ N[1.0*Pi, 200] and N[1.0*Pi] give the same result, because the machine precision 1.0 "infects" the calculation. (Even though machine precision is ~16 digits, it's considered to be lower precision than any arbitrary precision.) So N by itself isn't enough to define rr. SetPrecision, on the other hand, can forcibly increase the precision. Even in that case, though, you still need to use WorkingPrecision for the plot, since it determines the precision of the specific values of x and k used to evaluate rr. Otherwise the default is to use machine precision and it underflows. $\endgroup$ – Brett Champion Aug 30 '18 at 3:55

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.