Playing with a tricky integral in Maple and Mathematica
A colleague recently sent me this issue. Consider the following integral
Attempting to evaluate this with Mathematica 9 gives 0:
f[a_, b_] := Exp[I*(a*x^3 + b*x^2)];
Integrate[f[a, b], {x, -Infinity, Infinity},
Assumptions -> {a > 0, b \[Element] Reals}]
Out[2]= 0
So, for all valid values of a and b, Mathematica tells us that the resulting integral will be 0.
However, Let’s set a=1/3 and b=0 in the function f and ask Mathematica to evaluate that integral
In[3]:= Integrate[f[1/3, 0], {x, -Infinity, Infinity}]
Out[3]= (2*Pi)/(3^(2/3)*Gamma[2/3])
This is definitely not zero as we can see by numerically evaluating the result:
In[5]:= (2*Pi)/(3^(2/3)*Gamma[2/3])//N Out[5]= 2.23071
Similarly if we put a=1/3 and b=1 we get another non-zero result
In[7]:= Integrate[f[1/3, 1], {x, -Infinity, Infinity}] // InputForm
Out[7]=2*E^((2*I)/3)*Pi*AiryAi[-1]
In[8]:=
2*E^((2*I)/3)*Pi*AiryAi[-1]//N
Out[8]= 2.64453 + 2.08083 I
We are faced with a situation where Mathematica is disagreeing with itself. On one hand, it says that the integral is 0 for all a and b but on the other it gives non-zero results for certain combinations of a and b. Which result do we trust?
One way to proceed would be to use the NIntegrate[] function for the two cases where a and b are explicitly defined. NIntegrate[] uses completely different algorithms from Integrate. In particular it uses numerical rather than symbolic methods (apart from some symbolic pre-processing).
NIntegrate[f[1/3, 0], {x, -Infinity, Infinity}]
gives 2.23071 + 0. I and
NIntegrate[f[1/3, 1], {x, -Infinity, Infinity}]
gives 2.64453 + 2.08083 I
What we’ve shown here is that evaluating these integrals using numerical methods gives the same result as evaluating using symbolic methods and numericalizing the result. This gives me some confidence that its the general, symbolic evaluation that’s incorrect and hence I can file a bug report with Wolfram Research.
Maple 17.01 on the general problem
Since my University has just got a site license for Maple, I wondered what Maple 17.01 would make of the general integral. Using Maple’s 1D input/output we get:
> myint := int(exp(I*(a*x^3+b*x^2)), x = -infinity .. infinity); myint := (1/12)*3^(1/2)*2^(1/3)*((4/3)*Pi^(5/2)*(((4/27)*I)*b^3/a^2)^(1/3)*exp(((2/27)*I)*b^3/a^2) *BesselI(-1/3, ((2/27)*I)*b^3/a^2)+((2/3)*I)*Pi^(3/2)*3^(1/2)*b*2^(2/3) *hypergeom([1/2, 1], [2/3, 4/3],((4/27)*I)*b^3/a^2)/(-I*a)^(2/3)-(8/27)*2^(1/3)*Pi^(5/2)*b^2* exp(((2/27)*I)*b^3/a^2)*BesselI(1/3, ((2/27)*I)*b^3/a^2)/((-I*a)^(4/3)*(((4/27)*I)*b^3/a^2)^(1/3))) /(Pi^(3/2)*(-I*a)^(1/3))+(1/12)*3^(1/2)*2^(1/3)*((4/3)*Pi^(5/2)*(((4/27)*I)*b^3/a^2) ^(1/3)*exp(((2/27)*I)*b^3/a^2)*BesselI(-1/3, ((2/27)*I)*b^3/a^2)+((2/3)*I)*Pi^(3/2)*3^(1/2)*b*2^(2/3) *hypergeom([1/2, 1], [2/3, 4/3], ((4/27)*I)*b^3/a^2)/(I*a)^(2/3)-(8/27)*2^(1/3)*Pi^(5/2)*b^2* exp(((2/27)*I)*b^3/a^2)*BesselI(1/3, ((2/27)*I)*b^3/a^2)/((I*a)^(4/3)*(((4/27)*I)*b^3/a^2)^(1/3))) /(Pi^(3/2)*(I*a)^(1/3))
That looks scary! To try to determine if it’s possibly a correct general result, let’s turn this expression into a function and evaluate it for values of a and b we already know the answer to.
>f := unapply(%, a, b): >result1:=simplify(f(1/3,1)); result1 := (2/27)*3^(1/2)*Pi*exp((2/3)*I)*((-(1/3)*I)^(1/3)*3^(2/3)*(-1)^(1/6)*BesselI(-1/3, (2/3)*I) +3*(-(1/3)*I)^(2/3)*BesselI(-1/3, (2/3)*I)+3*(-(1/3)*I)^(2/3)*BesselI(1/3, (2/3)*I)-BesselI(1/3, (2/3)*I) *3^(1/3)*(-1)^(1/3))/(-(1/3)*I)^(2/3) evalf(result1); 2.644532853+2.080831872*I
Recall that Mathematica returned 2*E^((2*I)/3)*Pi*AiryAi[-1]=2.64453 + 2.08083 I for this case. The numerical results agree to the default precision reported by both packages so I am reasonably confident that Maple’s general solution is correct.
Not so simple simplification?
I am also confident that Maple’s expression involving Bessel functions is equivalent to Mathematica’s expression involving the AiryAi function. I haven’t figured out, however, how to get Maple to automatically reduce the Bessel functions down to AiryAi. I can attempt to go the other way though. In Maple:
>convert(2*exp((2*I)/3)*Pi*AiryAi(-1),Bessel); (2/3)*exp((2/3)*I)*Pi*(-1)^(1/6)*(-BesselI(1/3, (2/3)*I)*(-1)^(2/3)+BesselI(-1/3, (2/3)*I))
This is more compact than the Bessel function result I get from Maple’s simplify so I guess that Maple’s simplify function could have done a little better here.
Not so general general solution?
Maple’s solution of the general problem should be good for all a and b right? Let’s try it with a=1/3, b=0
f(1/3,0); Error, (in BesselI) numeric exception: division by zero
Uh-Oh! So it’s not valid for b=0 then! We know from Mathematica, however, that the solution for this particular case is (2*Pi)/(3^(2/3)*Gamma[2/3])=2.23071. Indeed, if we solve this integral directly in Maple, it agrees with Mathematica
>myint := int(exp(I*(1/3*x^3+0*x^2)), x = -infinity .. infinity); myint := (2/9)*3^(5/6)*(-1)^(1/6)*Pi/GAMMA(2/3)-(2/9)*3^(5/6)*(-1)^(5/6)*Pi/GAMMA(2/3) >simplify(myint); (2/3)*3^(1/3)*Pi/GAMMA(2/3) >evalf(myint); 2.230707052+0.*I
Going to back to the general result that Maple returned. Although we can’t calculate it directly, We can use limits to see what its value would be at a=1/3, b=0
>simplify(limit(f(1/3, x), x = 0)); (2/9)*Pi*3^(2/3)/((-(1/3)*I)^(1/3)*GAMMA(2/3)*((1/3)*I)^(1/3)) >evalf(%) evalf(%); 2.230707053+0.1348486379e-9*I
The symbolic expression looks different and we’ve picked up an imaginary part that’s possibly numerical noise. Let’s use more numerical precision to see if we can make the imaginary part go away. 100 digits should do it
>evalf[100](%%); 2.230707051824495741427486519543450239771293806030489125938415383976032081571780278667202004477941904 -0.855678513686459467847075286617333072231119978694352241387724335424279116026690601018858453303153701e-100*I
Well, more precision has made the imaginary part smaller but it’s still there. No matter how much precision I use, it’s always there…getting smaller and smaller as I ramp up the level of precision.
What’s going on?
All I’m doing here is playing around with the same problem in two packages. Does anyone have any further insights into some of the issues raised?
Weirdly, in Mathematica 9.0.1 I get the answer:
ConditionalExpression[0, a \[Element] Reals && Im[b] > 0]
for the general integral. That has presumably just skipped out both of your test cases, saying that they are invalid.
Testing with a=1, b=2i in Mathematica with NIntegrate gives 1.16379 + 0. I, while Integrate directly gives (4 E^(16/27) BesselK[1/3, 16/27])/(3 Sqrt[3]) in this case, which comes to the same thing.
>Weirdly, in Mathematica 9.0.1 I get the answer:
>ConditionalExpression[0, a \[Element] Reals && Im[b] > 0]
I’m definitely just getting 0 but then I am using 9.0
In[3]:= $Version
Out[3]= “9.0 for Microsoft Windows (64-bit) (January 25, 2013)”
“My” 9.0.1 gives 0 and WolframAlpha times out (I don’t have the pro version so that could be the issue here). It is always worrying to see how many errors are still present in these systems…
My $Version is “9.0 for Mac OS X x86 (64-bit) (January 24, 2013)” – how bizarre.
In Maple: try simplify(simplify(myint), size) and the answer goes from length 891 to length 378. In the presence of a and b symbolic, try combine(denom(myint), symbolic) and you see that this answer cannot be used right away for a = 0 or b = 0. In the case b=0, as you say, you can still use limits to evaluate the answer, but that doesn’t help for a = 0. So going back to the integrand, substitute a = 0 and recompute the integral: you get a piecewise solution completing the picture, I think. Maple could have returned a piecewise in the first place including this case a = 0.
Forgot to mention: to express Bessels and others in terms of Airy functions, you can use convert(myint, Airy).
Forgot to mention: to express myint in terms of Airy functions, you can use convert(myint, Airy).
I posted two questions on that:
http://mathematica.stackexchange.com/questions/29563/bug-in-integrate-command
http://math.stackexchange.com/questions/456424/hard-integral-that-standard-cas-get-totally-wrong
If you intend that `a` and `b` are purely real (or positive) then it can help to supply such information to the system. That may be apply to both the integration and the simplification steps. And there may be other avenues for getting the “best” simplification.
For example, in Maple 17.02,
> restart:
> myint:=int(exp(I*(a*x^3+b*x^2)),x=-infinity..infinity) assuming real:
> simp_myint:=simplify(simplify(myint),size) assuming real:
> lprint(simp_myint);
-2/27/(I*a)^(2/3)/(-I*a)^(2/3)/(1/a^2)^(1/3)*(-(I*b^3)^(2/3)*a*(1/a^2)^(2/3)*((I
*a)^(1/3)*(-I*a)^(2/3)+(I*a)^(2/3)*(-I*a)^(1/3))*BesselI(-1/3,2/27*I*b^3/a^2)+I*
((I*a)^(2/3)-(-I*a)^(2/3))*b^2*BesselI(1/3,2/27*I*b^3/a^2))*3^(1/2)*exp(2/27*I*b
^3/a^2)*Pi/(I*b^3)^(1/3)/a
> evalc(eval(limit(myint,b=0),a=1/3));
1/3 1/3 1/3
Pi 3 2 4
1/3 —————–
GAMMA(2/3)
> evalf(%);
2.230707052
I reported the bug to Wolfram Research and they suggested using ExpToTrig to get around the problem. I’m using a different version of Mathematica right now because I’m at home
In[1]:= $Version
Out[1]= “8.0 for Linux x86 (64-bit) (October 10, 2011)”
In[2]:= ExpToTrig[Exp[I*(a*x^3 + b*x^2)]]
Out[2]= Cos[b x^2 + a x^3] + I Sin[b x^2 + a x^3]
Mathematica can use this to get a very nice looking result.
In[3]:= Integrate[ExpToTrig[Exp[I*(a*x^3 + b*x^2)]], {x, -Infinity, Infinity},
Assumptions -> {a > 0, b \[Element] Reals}] // InputForm
Out[3]//InputForm=
(2*E^((((2*I)/27)*b^3)/a^2)*Pi*AiryAi[-b^2/(3*3^(1/3)*a^(4/3))])/
(3^(1/3)*a^(1/3))
If we put a=1/3,b=1 or a=1/3,b=0 into this expression, we recover the specific cases shown in the main blog post.
In response to my question in the Mathematica Stackexchange Daniel Lichtblau from Wolframresearch has acknowledged the bug: “[…] the complex exponential route has trouble in Slater convolution when it hits things like Exp[2*Log[I*a]], which it cannot discern from Exp[2*Log[-I*a]] (both show up because we integrate from -infinity to zero by negating and going 0 to infinity). An effort is made to unravel such things, but apparently this is not always getting it correct.”
and: “Have put in place a provisional fix. Now hoping nothing breaks as a result of that– would be nice for the fix to survive until the next release.”
Source: http://mathematica.stackexchange.com/questions/29563/have-i-found-a-bug-in-integrate
Very nice! Thanks for doing that vonjd.