A colleague recently sent me this issue. Consider the following integral<\/p>\n
<\/p>\n
Attempting to evaluate this with Mathematica<\/a> 9 gives 0:<\/p>\n So, for all valid values of a and b, Mathematica tells us that the resulting integral will be 0.<\/p>\n However, Let’s set a=1\/3 and b=0 in the function f and ask Mathematica to evaluate that integral<\/p>\n This is definitely not zero as we can see by numerically evaluating the result:<\/p>\n Similarly if we put a=1\/3 and b=1 we get another non-zero result<\/p>\n We are faced with a situation where Mathematica is disagreeing with itself.\u00a0 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.\u00a0 Which result do we trust?<\/p>\n One way to proceed would be to use the NIntegrate[]<\/a> function for the two cases where a and b are explicitly defined.\u00a0 NIntegrate[] uses completely different algorithms from Integrate.\u00a0 In particular it uses numerical rather than symbolic methods (apart from some symbolic pre-processing).<\/p>\n gives 2.23071 + 0. I and<\/p>\n gives 2.64453 + 2.08083 I<\/p>\n 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.\u00a0 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.<\/p>\n Maple 17.01 on the general problem Since my University has just got a site license for Maple<\/a>, I wondered what Maple 17.01 would make of the general integral. Using Maple’s 1D input\/output<\/a> we get:<\/p>\n 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.<\/p>\n 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.<\/p>\n Not so simple simplification?<\/strong><\/p>\n 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:<\/p>\n 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.<\/p>\n Not so general general solution?<\/strong><\/p>\n Maple’s solution of the general problem should be good for all a and b right?\u00a0 Let’s try it with a=1\/3, b=0<\/p>\n 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<\/p>\n 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<\/p>\n 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<\/p>\n 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.<\/p>\n What’s going on?<\/strong><\/p>\n All I’m doing here is playing around with the same problem in two packages.\u00a0 Does anyone have any further insights into some of the issues raised?<\/p>\n","protected":false},"excerpt":{"rendered":" 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 […]<\/p>\n","protected":false},"author":3,"featured_media":0,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"standard","meta":{"jetpack_post_was_ever_published":false,"footnotes":"","jetpack_publicize_message":"","jetpack_publicize_feature_enabled":true,"jetpack_social_post_already_shared":true,"jetpack_social_options":{"image_generator_settings":{"template":"highway","default_image_id":0,"font":"","enabled":false},"version":2}},"categories":[25,4,8],"tags":[],"class_list":["post-5031","post","type-post","status-publish","format-standard","hentry","category-maple","category-math-software","category-mathematica"],"jetpack_publicize_connections":[],"jetpack_featured_media_url":"","jetpack_shortlink":"https:\/\/wp.me\/p3swhs-1j9","jetpack_sharing_enabled":true,"_links":{"self":[{"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=\/wp\/v2\/posts\/5031","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=\/wp\/v2\/users\/3"}],"replies":[{"embeddable":true,"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=%2Fwp%2Fv2%2Fcomments&post=5031"}],"version-history":[{"count":20,"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=\/wp\/v2\/posts\/5031\/revisions"}],"predecessor-version":[{"id":5062,"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=\/wp\/v2\/posts\/5031\/revisions\/5062"}],"wp:attachment":[{"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=%2Fwp%2Fv2%2Fmedia&parent=5031"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=%2Fwp%2Fv2%2Fcategories&post=5031"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=%2Fwp%2Fv2%2Ftags&post=5031"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}} f[a_, b_] := Exp[I*(a*x^3 + b*x^2)];\r\nIntegrate[f[a, b], {x, -Infinity, Infinity}, \r\n Assumptions -> {a > 0, b \\[Element] Reals}]\r\n\r\nOut[2]= 0<\/pre>\nIn[3]:= Integrate[f[1\/3, 0], {x, -Infinity, Infinity}]\r\n\r\nOut[3]= (2*Pi)\/(3^(2\/3)*Gamma[2\/3])<\/pre>\nIn[5]:=\r\n(2*Pi)\/(3^(2\/3)*Gamma[2\/3])\/\/N\r\n\r\nOut[5]= 2.23071<\/pre>\n
In[7]:= Integrate[f[1\/3, 1], {x, -Infinity, Infinity}] \/\/ InputForm\r\n\r\nOut[7]=2*E^((2*I)\/3)*Pi*AiryAi[-1]\r\n\r\nIn[8]:=\r\n2*E^((2*I)\/3)*Pi*AiryAi[-1]\/\/N\r\n\r\nOut[8]= 2.64453 + 2.08083 I<\/pre>\nNIntegrate[f[1\/3, 0], {x, -Infinity, Infinity}]<\/pre>\nNIntegrate[f[1\/3, 1], {x, -Infinity, Infinity}]<\/pre>\n
\n<\/strong><\/p>\n> myint := int(exp(I*(a*x^3+b*x^2)), x = -infinity .. infinity);\r\n\r\nmyint := (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)\r\n*BesselI(-1\/3, ((2\/27)*I)*b^3\/a^2)+((2\/3)*I)*Pi^(3\/2)*3^(1\/2)*b*2^(2\/3)\r\n*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*\r\nexp(((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)))\r\n\/(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)\r\n^(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)\r\n*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*\r\nexp(((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)))\r\n\/(Pi^(3\/2)*(I*a)^(1\/3))<\/pre>\n
>f := unapply(%, a, b):\r\n>result1:=simplify(f(1\/3,1));\r\n\r\nresult1 := (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)\r\n+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)\r\n*3^(1\/3)*(-1)^(1\/3))\/(-(1\/3)*I)^(2\/3)\r\n\r\nevalf(result1);\r\n2.644532853+2.080831872*I<\/pre>\n
>convert(2*exp((2*I)\/3)*Pi*AiryAi(-1),Bessel);\r\n\r\n(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))<\/pre>\n
f(1\/3,0);\r\nError, (in BesselI) numeric exception: division by zero<\/pre>\n
>myint := int(exp(I*(1\/3*x^3+0*x^2)), x = -infinity .. infinity);\r\n\r\nmyint := (2\/9)*3^(5\/6)*(-1)^(1\/6)*Pi\/GAMMA(2\/3)-(2\/9)*3^(5\/6)*(-1)^(5\/6)*Pi\/GAMMA(2\/3)\r\n\r\n>simplify(myint);\r\n(2\/3)*3^(1\/3)*Pi\/GAMMA(2\/3)\r\n\r\n>evalf(myint);\r\n2.230707052+0.*I<\/pre>\n
>simplify(limit(f(1\/3, x), x = 0));\r\n\r\n(2\/9)*Pi*3^(2\/3)\/((-(1\/3)*I)^(1\/3)*GAMMA(2\/3)*((1\/3)*I)^(1\/3))\r\n\r\n>evalf(%)\r\n\r\nevalf(%);\r\n2.230707053+0.1348486379e-9*I<\/pre>\n
>evalf[100](%%);\r\n\r\n2.230707051824495741427486519543450239771293806030489125938415383976032081571780278667202004477941904\r\n-0.855678513686459467847075286617333072231119978694352241387724335424279116026690601018858453303153701e-100*I<\/pre>\n