Integration error
in [Mathematica]
|
From: Jason Alexander on 6 Apr 2010 07:22 Hello all, I'm getting a strange result when calculating what I thought was a relatively straightforward integral. If I evaluate the following expression, I receive an imaginary result when it should, in fact, equal 1. (This happens in Mathematica 7.0.1 on Mac OS X): In[1]:= Integrate[( 4 (8/(2320 + 3 x (-96 + 3 x)) + 8/(1460 + 3 x (-76 + 3 x)) + 6/( 800 + 3 x (-56 + 3 x)) + 4/(340 + 3 x (-36 + 3 x)) + 1/( 80 + 3 x (-16 + 3 x))))/(9 \[Pi]), {x, -Infinity, Infinity}] Out[1]:= 1/3 + (52 I)/9 However, if I break the integral into two parts, as below, then I get an answer of 1. In[2]:= Integrate[( 4 (8/(2320 + 3 x (-96 + 3 x)) + 8/(1460 + 3 x (-76 + 3 x)) + 6/( 800 + 3 x (-56 + 3 x)) + 4/(340 + 3 x (-36 + 3 x)) + 1/( 80 + 3 x (-16 + 3 x))))/(9 \[Pi]), {x, -Infinity, 0}] + Integrate[( 4 (8/(2320 + 3 x (-96 + 3 x)) + 8/(1460 + 3 x (-76 + 3 x)) + 6/( 800 + 3 x (-56 + 3 x)) + 4/(340 + 3 x (-36 + 3 x)) + 1/( 80 + 3 x (-16 + 3 x))))/(9 \[Pi]), {x, 0, Infinity}] // Simplify Out[2]:= 1 Furthermore, if I compute the antiderivative, I get the following: (1/(9 \[Pi]))4 (-(1/24) ArcTan[1/4 (8 - 3 x)] - 1/4 ArcTan[1/4 (28 - 3 x)] - 1/3 ArcTan[1/4 (38 - 3 x)] + 2/3 ArcTan[3/4 (-16 + x)] + 1/3 ArcTan[3/4 (-6 + x)] + 1/3 ArcTan[1/4 (-38 + 3 x)] + 1/4 ArcTan[1/4 (-28 + 3 x)] + 1/24 ArcTan[1/4 (-8 + 3 x)]) And if I define f[x] to be the above, and then ask Mathematica to calculate Limit[f[n] - f[-n], n -> Infinity] I get 1, as expected. Is this a bug, or am I overlooking some subtle point in how Mathematica computes improper integrals? Cheers, Jason -- Dr. J. McKenzie Alexander Department of Philosophy, Logic and Scientific Method London School of Economics Houghton Street, London WC2A 2AE
From: dh on 6 Apr 2010 08:21 On 06.04.2010 13:22, Jason Alexander wrote: > Hello all, > > I'm getting a strange result when calculating what I thought was a relatively straightforward integral. If I evaluate the following expression, I receive an imaginary result when it should, in fact, equal 1. (This happens in Mathematica 7.0.1 on Mac OS X): > > In[1]:= > Integrate[( > 4 (8/(2320 + 3 x (-96 + 3 x)) + 8/(1460 + 3 x (-76 + 3 x)) + 6/( > 800 + 3 x (-56 + 3 x)) + 4/(340 + 3 x (-36 + 3 x)) + 1/( > 80 + 3 x (-16 + 3 x))))/(9 \[Pi]), {x, -Infinity, Infinity}] > > Out[1]:= > 1/3 + (52 I)/9 > > However, if I break the integral into two parts, as below, then I get an answer of 1. > > In[2]:= > Integrate[( > 4 (8/(2320 + 3 x (-96 + 3 x)) + 8/(1460 + 3 x (-76 + 3 x)) + 6/( > 800 + 3 x (-56 + 3 x)) + 4/(340 + 3 x (-36 + 3 x)) + 1/( > 80 + 3 x (-16 + 3 x))))/(9 \[Pi]), {x, -Infinity, 0}] + > Integrate[( > 4 (8/(2320 + 3 x (-96 + 3 x)) + 8/(1460 + 3 x (-76 + 3 x)) + 6/( > 800 + 3 x (-56 + 3 x)) + 4/(340 + 3 x (-36 + 3 x)) + 1/( > 80 + 3 x (-16 + 3 x))))/(9 \[Pi]), {x, 0, Infinity}] // Simplify > > Out[2]:= > 1 > > Furthermore, if I compute the antiderivative, I get the following: > > (1/(9 \[Pi]))4 (-(1/24) ArcTan[1/4 (8 - 3 x)] - > 1/4 ArcTan[1/4 (28 - 3 x)] - 1/3 ArcTan[1/4 (38 - 3 x)] + > 2/3 ArcTan[3/4 (-16 + x)] + 1/3 ArcTan[3/4 (-6 + x)] + > 1/3 ArcTan[1/4 (-38 + 3 x)] + 1/4 ArcTan[1/4 (-28 + 3 x)] + > 1/24 ArcTan[1/4 (-8 + 3 x)]) > > And if I define f[x] to be the above, and then ask Mathematica to calculate > > Limit[f[n] - f[-n], n -> Infinity] > > I get 1, as expected. > > Is this a bug, or am I overlooking some subtle point in how Mathematica computes improper integrals? > > > Cheers, > > Jason > > -- > Dr. J. McKenzie Alexander > Department of Philosophy, Logic and Scientific Method > London School of Economics > Houghton Street, London WC2A 2AE > > Hi Jason, looks like a bug to me. You have a rational real integrand without poles in R, therefore, the definite integral must be real. Please report this bug to Wolfram. Daniel -- Daniel Huber Metrohm Ltd. Oberdorfstr. 68 CH-9100 Herisau Tel. +41 71 353 8585, Fax +41 71 353 8907 E-Mail:<mailto:dh(a)metrohm.com> Internet:<http://www.metrohm.com>
From: J. McKenzie Alexander on 7 Apr 2010 03:20 Just a follow-up... Here's another problem with Mathematica's integration engine: take the PDF for the (0,1) - L=E9vy distribution, shift it 1/4 to the right, then calculate the area under the curve between -1/4 and 0. I get an imaginary answer... In[41]:== PDF[LevyDistribution[0, 1]][x] /. {x -> (x - 1/4)} Out[41]== E^(-(1/(2 (-(1/4) + x))))/(Sqrt[2 \[Pi]] (-(1/4) + x)^(3/2)) and In[42]:== Integrate[ E^(-(1/(2 (-(1/4) + x))))/( Sqrt[2 \[Pi]] (-(1/4) + x)^(3/2)), {x, -1/4, 0} ] Out[42]== -I (Erfi[1] - Erfi[Sqrt[2]]) Cheers, Jason > Hi Jason, > looks like a bug to me. You have a rational real integrand without > poles in R, therefore, the definite integral must be real. > Please report this bug to Wolfram. > Daniel > > > -- > > Daniel Huber > Metrohm Ltd. > Oberdorfstr. 68 > CH-9100 Herisau > Tel. +41 71 353 8585, Fax +41 71 353 8907 > E-Mail:<mailto:dh(a)metrohm.com> > Internet:<http://www.metrohm.com> > > -- Dr J. McKenzie Alexander Department of Philosophy, Logic and Scientific Method London School of Economics and Political Science Houghton Street, London WC2A 2AE
From: dh on 7 Apr 2010 03:20 Hi Jason, here the integrand is a complex function In your case, the trouble maker is the factor: (-(1/4) + x)^(-3/2); this is clearly imaginary for x<1/4 Daniel On 06.04.2010 15:47, J. McKenzie Alexander wrote: > Just a follow-up... > > Here's another problem with Mathematica's integration engine: take the PDF for the (0,1) - L=E9vy distribution, shift it 1/4 to the right, then calculate the area under the curve between -1/4 and 0. I get an imaginary answer... > > In[41]:== PDF[LevyDistribution[0, 1]][x] /. {x -> (x - 1/4)} > > Out[41]== E^(-(1/(2 (-(1/4) + x))))/(Sqrt[2 \[Pi]] (-(1/4) + x)^(3/2)) > > and > > In[42]:== Integrate[ > E^(-(1/(2 (-(1/4) + x))))/( > Sqrt[2 \[Pi]] (-(1/4) + x)^(3/2)), {x, -1/4, 0} > ] > > Out[42]== -I (Erfi[1] - Erfi[Sqrt[2]]) > > Cheers, > > Jason > >> Hi Jason, >> looks like a bug to me. You have a rational real integrand without >> poles in R, therefore, the definite integral must be real. >> Please report this bug to Wolfram. >> Daniel >> >> >> -- >> >> Daniel Huber >> Metrohm Ltd. >> Oberdorfstr. 68 >> CH-9100 Herisau >> Tel. +41 71 353 8585, Fax +41 71 353 8907 >> E-Mail:<mailto:dh(a)metrohm.com> >> Internet:<http://www.metrohm.com> >> >> > > -- > Dr J. McKenzie Alexander > Department of Philosophy, Logic and Scientific Method > London School of Economics and Political Science > Houghton Street, London WC2A 2AE > > > > > > -- Daniel Huber Metrohm AG International Headquarters Oberdorfstr. 68, CH-9101 Herisau / Switzerland Phone +41 71 353 8606, Fax +41 71 353 89 01 Mail <mailto:dh(a)metrohm.com> Web <http://www.metrohm.com
From: DrMajorBob on 7 Apr 2010 03:20
You (and all of us) are overlooking MANY subtle points buried in Integrate. (We have no choice.) Suffice it to say, always check your results another way. Break things down into smaller intervals, as you have done. And don't trust ANYBODY. Bobby On Tue, 06 Apr 2010 06:23:03 -0500, Jason Alexander <jalex(a)lse.ac.uk> wrote: > Hello all, > > I'm getting a strange result when calculating what I thought was a > relatively straightforward integral. If I evaluate the following > expression, I receive an imaginary result when it should, in fact, equal > 1. (This happens in Mathematica 7.0.1 on Mac OS X): > > In[1]:= > Integrate[( > 4 (8/(2320 + 3 x (-96 + 3 x)) + 8/(1460 + 3 x (-76 + 3 x)) + 6/( > 800 + 3 x (-56 + 3 x)) + 4/(340 + 3 x (-36 + 3 x)) + 1/( > 80 + 3 x (-16 + 3 x))))/(9 \[Pi]), {x, -Infinity, Infinity}] > > Out[1]:= > 1/3 + (52 I)/9 > > However, if I break the integral into two parts, as below, then I get an > answer of 1. > > In[2]:= > Integrate[( > 4 (8/(2320 + 3 x (-96 + 3 x)) + 8/(1460 + 3 x (-76 + 3 x)) + 6/( > 800 + 3 x (-56 + 3 x)) + 4/(340 + 3 x (-36 + 3 x)) + 1/( > 80 + 3 x (-16 + 3 x))))/(9 \[Pi]), {x, -Infinity, 0}] + > Integrate[( > 4 (8/(2320 + 3 x (-96 + 3 x)) + 8/(1460 + 3 x (-76 + 3 x)) + 6/( > 800 + 3 x (-56 + 3 x)) + 4/(340 + 3 x (-36 + 3 x)) + 1/( > 80 + 3 x (-16 + 3 x))))/(9 \[Pi]), {x, 0, Infinity}] // Simplify > > Out[2]:= > 1 > > Furthermore, if I compute the antiderivative, I get the following: > > (1/(9 \[Pi]))4 (-(1/24) ArcTan[1/4 (8 - 3 x)] - > 1/4 ArcTan[1/4 (28 - 3 x)] - 1/3 ArcTan[1/4 (38 - 3 x)] + > 2/3 ArcTan[3/4 (-16 + x)] + 1/3 ArcTan[3/4 (-6 + x)] + > 1/3 ArcTan[1/4 (-38 + 3 x)] + 1/4 ArcTan[1/4 (-28 + 3 x)] + > 1/24 ArcTan[1/4 (-8 + 3 x)]) > > And if I define f[x] to be the above, and then ask Mathematica to > calculate > > Limit[f[n] - f[-n], n -> Infinity] > > I get 1, as expected. > > Is this a bug, or am I overlooking some subtle point in how Mathematica > computes improper integrals? > > > Cheers, > > Jason > > -- > Dr. J. McKenzie Alexander > Department of Philosophy, Logic and Scientific Method > London School of Economics > Houghton Street, London WC2A 2AE > > -- DrMajorBob(a)yahoo.com |