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