11
$\begingroup$

Bug introduced in 7.0 or earlier and fixed in 11.0

This is a follow-up to this question regarding numerical instabilities occurring with modified Bessel functions. In trying to explore J.M.'s answer to that question in greater depth by setting WorkingPrecision -> 50 (as advised in the other answer by xslittlegrass), I obtained vastly different (and incorrect!) results than when not explicitly setting a non-default WorkingPrecision (see below for plots). I am trying to figure out why this is occurring and was advised to make a separate question regarding this.

The function that I am attempting to plot (for fixed $b$, as a function of $a$) is

f[a_, b_] := (a^2 Exp[a^2/(8 b)] (BesselK[3/4, a^2/(8 b)] - 
  BesselK[1/4, a^2/(8 b)]))/(8 Sqrt[a b^3])

which is the analytic solution to

Integrate[x^2 Exp[-a x^2 - b x^4], {x, -∞, ∞}, Assumptions -> {a > 0, b > 0}]

I also plot the numerical solution to this,

g[a_,b_] := NIntegrate[x^2 Exp[-a x^2 - b x^4], {x, -∞, ∞}]

Without any regards to WorkingPrecision, I obtain the following plot:

amin = 10^-10; amax = 0.1; bval = 10^-3;
LogLogPlot[{f[a, bval], g[a, bval]}, {a, amin, amax}, PlotRange -> All]

enter image description here

Even with explicitly giving the variable/bounds higher precision,

amin = 1`50 10^-10; amax = 0.1`50; bval = 1`50 10^-3;

the plot remains the same as above. If, however, I add in WorkingPrecision -> 50, a plot with function values many orders of magnitude higher is obtained:

amin = 1`50 10^-10; amax = 0.1`50; bval = 1`50 10^-3;
LogLogPlot[{f[a, bval], g[a, bval]}, {a, amin, amax}, WorkingPrecision -> 50,
             PlotRange -> All]

enter image description here

The fact that this is occurring is worrisome, as I would generally expect the results obtained with higher WorkingPrecision to be more trustworthy; in this case, however, they are completely wrong -- here's a representative plot:

bval = 10^-3;
Plot[x^2 Exp[-a x^2 - bval x^4] /. a -> 10^-5, {x, -15, 15}]

enter image description here

Why is this happening? Thank you very much!

Improve this question
$\endgroup$
6
  • $\begingroup$ Do you mean Integrate instead of NIntegrate in your definition of g? $\endgroup$
    – bbgodfrey
    May 24, 2016 at 4:19
  • $\begingroup$ I believe he did mean NIntegrate[], @bb; he was comparing the solution I derived in the other thread with the results of numerical integration. $\endgroup$
    – J. M.'s slightly less busy ♦
    May 24, 2016 at 4:27
  • $\begingroup$ @J.M. I understand your point. However, Assumptions is not a legal option for NIntegrate. $\endgroup$
    – bbgodfrey
    May 24, 2016 at 4:29
  • $\begingroup$ Yes, prolly an error of copying, @bb... $\endgroup$
    – J. M.'s slightly less busy ♦
    May 24, 2016 at 4:32
  • 5
    $\begingroup$ Possible duplicate of LogPlot axes labels destroyed when working in high precision $\endgroup$
    – Alexey Popkov
    May 24, 2016 at 5:52

1 Answer 1

Reset to default
12
$\begingroup$

There appears to be a bug, not in Integrate or in BesselK, but in the vertical-axis Ticks of LogLogPlot. Consider the simple case,

LogLogPlot[Exp[x], {x, 10^-10, 1}, PlotRange -> All]

enter image description here

as it should be. However,

LogLogPlot[Exp[x], {x, 10^-10, 1}, WorkingPrecision -> 50, PlotRange -> All]

enter image description here

In fact, any value of WorkingPrecision except MachinePrecision triggers this error. Even WorkingPrecision -> $MachinePrecision produces the error.

Improve this answer
$\endgroup$
6
  • $\begingroup$ Oh my. That would be quite unfortunate. I will have to see if this issue arises for various functions when I'm at my computer again tomorrow morning – this is something that would somehow need to be reported to the appropriate Wolfram team, right? $\endgroup$
    – AnInquiringMind
    May 24, 2016 at 4:40
  • $\begingroup$ @PhysicsCodingEnthusiast Absolutely. Will you do so, or do you wish me to? Your call. $\endgroup$
    – bbgodfrey
    May 24, 2016 at 4:41
  • $\begingroup$ I can do so – is there a good way to do this? Should I provide a link to this question? I'll wait to do it until I'm at my computer again (currently on my phone) – in the meantime, perhaps others will chime in with their own findings. Thanks! $\endgroup$
    – AnInquiringMind
    May 24, 2016 at 4:44
  • $\begingroup$ @Physics, yes, please do link to this question as evidence if and when you make a report. $\endgroup$
    – J. M.'s slightly less busy ♦
    May 24, 2016 at 4:50
  • $\begingroup$ @PhysicsCodingEnthusiast Report the bug here. Do not be surprised, if the response is slow. $\endgroup$
    – bbgodfrey
    May 24, 2016 at 4:54

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.