Difference between NIntegrate and Integrate results

Question

I want to integrate a function of two variables that is defined via a piecewise function. I can compute the integral analytically as follows:

Integrate[Exp[-y x^α1 ], {x, 0, x0}, Assumptions -> {1 < α1, x0 > 0, y > y0, y0 > 0}] + Integrate[Exp[-x0^(α1 - α2) y x^α2], {x, x0, Infinity}, Assumptions -> {1 < α1, 0 < α2 < 1, x0 > 0, y > y0, y0 > 0, x0 > 0}]

that gives the solution

(x0 ExpIntegralE[(-1 + α2)/α2, x0^α1 y])/α2 +(y^(-1/α1) (Gamma[1/α1] -Gamma[1/α1, x0^α1 y]))/α1

Then, I integrate over y

Integrate[(x0 ExpIntegralE[(-1 + α2)/α2, x0^α1 y])/α2 + (y^(-1/α1) (Gamma[1/α1] - Gamma[1/α1, x0^α1 y]))/α1, {y, y0, Infinity}, Assumptions -> {α1 > 1, 1 > α2 > 0, x0 > 1, y0 > 0}]

and obtain

(y0^(1 - 1/α1) Gamma[1/α1])/(1 - α1) + (y0^((-1 + α1)/α1) Gamma[1/α1, x0^α1 y0])/(-1 + α1) + ( x0 ((E^(-x0^α1 y0)x0^-α1 (α1 - α2))/(-1 + α1) - x0^(-(α1/α2)) y0^((-1 + α2)/α2) Gamma[1/α2, x0^α1 y0]))/(-1 + α2)

I can evaluate it with the parameters that interest me and obtain a result:

With[{α1 = 1.4, α2 = 0.8, x0 = 8, y0 = 1}, ( y0^(1 - 1/α1) Gamma[1/α1])/(1 - α1) + ( y0^((-1 + α1)/α1)   Gamma[1/α1, x0^α1 y0])/(-1 + α1) + (  x0 ((E^(-x0^α1 y0)       x0^-α1 (α1 - α2))/(-1 + α1) -      x0^(-(α1/α2)) y0^((-1 + α2)/α2)      Gamma[1/α2, x0^α1 y0]))/(-1 + α2)]

which evaluates as -3.18998. This can't be true, because the original function was positive in all its domain.

If I try to integrate the function numerically, it warns me about possible singularities.

With[{α1 = 1.4, α2 = 0.8, x0 = 8, R = 10, y0 = 1},  NIntegrate[  Piecewise[{{Exp[-y x^α1 ],      x <= x0}, {Exp[-x0^(α1 - α2) y x^α2],      x > x0}}], {x, 0, Infinity}, {y, y0, Infinity}]]

And gives me a completely different result: 3.21302*10^10

My questions are: - Why does Mathematica produce an analytical solution that appears to be wrong? - How can I obtain a good estimate of the numerical integral? I've tried all Methods and played around with the options, but couldn't get an estimate that I could trust.

Thank you

Answer 1

In Mathematica 10 I can evaluate this directly, without any intermediate steps.

With[{α1 = 1.4, α2 = 0.8, x0 = 8, y0 = 1}, 
 Integrate[
  Piecewise[{{Exp[-y x^α1], 
     x <= x0}, {Exp[-x0^(α1 - α2) y x^α2], 
     x > x0}}], {x, 0, Infinity}, {y, y0, Infinity}, 
  PrincipalValue -> True]]
(* ∞ *)