Strange phenomenon occurring in analytic integration result involving Bessel functions

Question

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:

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.

Answer 1

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]

You can also set the precision in Plot using WorkingPrecision:

Plot[{f[a, 1], g[a, 1]}, {a, 1, 14}, WorkingPrecision -> 50]

Answer 2

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}]