Re: Yet another incorrect integral
- To: mathgroup at smc.vnet.net
- Subject: [mg39344] Re: [mg39338] Yet another incorrect integral
- From: Andrzej Kozlowski <akoz at mimuw.edu.pl>
- Date: Wed, 12 Feb 2003 03:51:56 -0500 (EST)
- Sender: owner-wri-mathgroup at wolfram.com
As happens wiht quite a few such cases the problem seems to be with
Limit. It still sometimes helps to load in the Calculus`Limit`
package. Another useful trick is to evaluate the integral with a finite
bound, then use FullSimplify and then take take the limit as the bound
goes to infinity.
In[1]:=
<< "Calculus`Limit`"
In[2]:=
FullSimplify[Integrate[x^p/E^x^2, {x, z, t}]]
Out[2]=
(1/2)*((-t^(1 + p))*ExpIntegralE[1/2 - p/2, t^2] +
z^(1 + p)*ExpIntegralE[1/2 - p/2, z^2])
So much does not need the limit package and seems correct.
In[3]:=
Limit[%, t -> Infinity]
Out[3]=
(1/2)*z^(1 + p)*ExpIntegralE[1/2 - p/2, z^2]
The built-in Limit is not able to do evaluate this.
ALthough the answer returned by Mathematica does not look like yours I
believe they are in fact equal, at least in most of the complex plane
(I have not considered this carefully but quick numerical experiments
indicate this is almost certainly true). They are not equal for z==0;
however in this case Mathematica gives:
Integrate[x^p/E^x^2, {x, 0, Infinity}]
If[Re[p] > -1, (1/2)*Gamma[(1 + p)/2],
Integrate[x^p/E^x^2, {x, 0, Infinity}]]
which again agrees with your answer for values of p for which the
integral is convergent.
Andrzej Kozlowski
Yokohama, Japan
http://www.mimuw.edu.pl/~akoz/
http://platon.c.u-tokyo.ac.jp/andrzej/
On Tuesday, February 11, 2003, at 06:47 PM, Bob Stagat wrote:
> Consider the following integral...
>
> Integrate[x^p*E^(-x^2), {x, z, Infinity}]
>
> Using the substitution x = t^2 it is easy to show that the answer is
> half the incomplete gamma function, (1/2)*Gamma[(p + 1)/2, z^2].
>
> However, if I ask Mathematica to do this integral, here's what I get...
>
> Integrate[x^p E^(-x^2), {x, z, Infinity}]
> PowerExpand[%]
> % /. z -> 0
> %% /. z -> Infinity
>
> Out[146]=
> (1/2)*(z^(p + 1)*Gamma[(p + 1)/2, z^2]*
> (z^2)^((1/2)*(-p - 1)) + Gamma[(p + 1)/2])
>
> Out[147]=
> (1/2)*(Gamma[(p + 1)/2] + Gamma[(p + 1)/2, z^2])
>
> Out[148]=
> Gamma[(p + 1)/2]
>
> Out[149]=
> (1/2)*Gamma[(p + 1)/2]
>
> This is incorrect. The correct result should be:
>
> Out[147]=
> (1/2)*Gamma[(p + 1)/2, z^2]
>
>
> Out[148]=
> (1/2)*Gamma[(p + 1)/2]
>
> Out[149]=
> 0
>
> Does anyone understand why Mathematica screws up on such a simple
> integral?
>
> -Bob Stagat-
> --
>
>
>