6
$\begingroup$

For the following integral,

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

Mathematica gives the following analytic result:

-(Exp[a^2/(8b)] π 
    (a^2 BesselI[-1/4, a^2/(8b)] - 
      (a^2 + 4b) BesselI[1/4, a^2/(8b)] + 
        a^2 (BesselI[3/4, a^2/(8b)] - BesselI[5/4, a^2/(8b)]))) /
    (8 Sqrt[2] Sqrt[a b^3])

Plotting this as a function of $a$ (with $b = 1.0$, say) along with the numerical integration result

NIntegrate[x^2 Exp[-a x^2 - b x^4] /. b -> 1.0, {x, -∞, ∞}]

we obtain the following:

plot

My question is: why does this occur with the analytic function? The exponential $\exp(a^2/8b)$ grows very large as $a^2/8b$ increases, so my guess is that the analytic result given by Mathematica is not valid everywhere -- but there is no indication of this given. Is there some way that I can see what is going on here so that I can obtain the correct analytic result for all values of $a^2/8b$?

Edit

While the solution may be similar to that given in the other question that this has been marked as a potential duplicate of, this question is not a duplicate, as it is not obvious that this expression involving Bessel functions would require the same sort of "finessing" as high-degree polynomials.

Improve this question
$\endgroup$
1
  • $\begingroup$ In retrospect, anybody familiar with the asymptotic behavior of modified Bessel functions would quickly find the returned result to be numerically suspect, since $I$ is large for large arguments when the final result is supposed to be small. This is akin to computing $\exp(-x)$ from $\cosh$ and $\sinh$ for large $x$. $\endgroup$
    – J. M.'s missing motivation ♦
    Commented May 21, 2016 at 4:43

2 Answers 2

Reset to default
10
$\begingroup$

It seems that the analytic result is correct, but the precision is lost when converting it to a number. For example, if we use a higher precision, we get consistent results between numerical and analytical integration:

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

ListPlot[{
  Table[{a, f[a, 1`50]}, {a, 1/10, 14, 1/10}],
  Table[{a, g[a, 1.]}, {a, 1/10, 14, 1/10}]
  }, Joined -> True]

enter image description here

You can also set the precision in Plot using WorkingPrecision:

Plot[{f[a, 1], g[a, 1]}, {a, 1, 14}, WorkingPrecision -> 50]
Improve this answer
$\endgroup$
3
  • $\begingroup$ Thank you! It seems as though the precision is needed for the numerator of the analytic function -- with $a=15.0, b=1.0$, the numerator "equals" $2.85959 \times 10^{10}$, whereas with $a = 15.0\unicode{x0060} 50$ and $b=1.0\unicode{x0060} 50$, the numerator equals $0.65766425...$. Is there a way that I can add this precision into the analytic function / its input parameters, $a$ and $b$, so that I can plot it without resorting to ListPlot? $\endgroup$
    – AnInquiringMind
    Commented May 21, 2016 at 3:40
  • $\begingroup$ @PhysicsCodingEnthusiast You can use WorkingPrecision in Plot to do that. See the update. $\endgroup$
    – xslittlegrass
    Commented May 21, 2016 at 3:45
  • $\begingroup$ That wasn't working for me before (I tried WorkingPrecision in Plot even before posting this question), but I think it was because I was setting $b=1.0$ instead of with the extra precision. Thank you for your help! $\endgroup$
    – AnInquiringMind
    Commented May 21, 2016 at 4:02
9
$\begingroup$

The culprit, as suspected by xslittlegrass, is indeed numerical instability; in particular, this is because of the perverse combination of modified Bessel functions exhibited in the result returned by Mathematica.

Using a recurrence identity satisfied by the modified Bessel function of the first kind, we can simplify the expression returned, like so:

Integrate[x^2 Exp[-a x^2 - b x^4], {x, -∞, ∞}, Assumptions -> {a > 0, b > 0}] /.
          BesselI[5/4, a^2/(8 b)] ->
          BesselI[-3/4, a^2/(8 b)] - 4 b/a^2 BesselI[1/4, a^2/(8 b)] // FullSimplify
   (a^2 E^(a^2/(8 b)) (- BesselK[1/4, a^2/(8 b)] + BesselK[3/4, a^2/(8 b)]))/(8 Sqrt[a b^3])

where the result is now in terms of the function of the second kind.

(As to why Mathematica is unable to return this result automatically, I've no idea.)

Thus,

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

Plot[pcefunction[a, 1], {a, 0, 15}]

nice and stable

Improve this answer
$\endgroup$
5
  • $\begingroup$ Thank you for reminding me about that recurrence identity -- I had been finding ones valid only for integer or half-integer modified Bessel functions. It seems as though that identity still needs to be coded into Mathematica! $\endgroup$
    – AnInquiringMind
    Commented May 21, 2016 at 4:35
  • 2
    $\begingroup$ In this very case, your observed instability boils down to "trying to compute a very tiny number by subtracting a bunch of big numbers together", so I initially tried to re-express everything in terms of $K$, which is tiny for large arguments. The cancellations I then observed prompted me to try the recurrence, and thus… $\endgroup$
    – J. M.'s missing motivation ♦
    Commented May 21, 2016 at 4:41
  • $\begingroup$ A follow-up for you: I am graphing a log-log plot of pcefunction with $b = 10^{-3}$ and $10^{-10} \leq a \leq 0.1$, and I obtain a graph of pcefunction that has values below ~110 for the entire domain. If I try to generate the same plot with WorkingPrecision -> 50, I obtain a result for pcefunction with values on the order of $10^{46}$. This persists even when I try to explicitly set the values of $b$ and the bounds of $a$ to have a precision of 50. Why should changing the WorkingPrecision give me such a different answer, and how can I get a consistent answer? Thanks so much! $\endgroup$
    – AnInquiringMind
    Commented May 23, 2016 at 2:37
  • $\begingroup$ That seems like it should be a separate question. Please include the plots if you do ask it. $\endgroup$
    – J. M.'s missing motivation ♦
    Commented May 23, 2016 at 3:41
  • $\begingroup$ Please see the follow-up question $\endgroup$
    – AnInquiringMind
    Commented May 24, 2016 at 2:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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