1
|
7888 Views
|
24 Replies
|
4 Total Likes
View groups...
Share
GROUPS:

# Problem with Integrate

Posted 10 years ago
 Hello. I was trying to do the following integral: Integrate[1/Sqrt[2 (1 + x)^3 + 1], {x, 0, 1}] // N  The output in Mathematica 9 is a complex number that is not correct -2.96267 - 1.92762 I  When I try to solve by using NIntegrate, the output is the correct one NIntegrate[1/Sqrt[2 (1 + x)^3 + 1], {x, 0, 1}] 0.376065  Is this a bug in Integrate or I am missing something? Can someone test in Mathematica 10? Thanks
24 Replies
Sort By:
Posted 9 years ago
 Let me add two remarks1) the exact value of the jump position seems to be xj = 2^(2/3) - 1 N[2^(2/3) - 1, 20] (* 0.58740105196819947475 *) 2) an antiderivative without a jump can be calculated as follows ad2F1[x_] = Integrate[1/Sqrt[1 + 2 x^3], {x, 0, t}, Assumptions -> t > 0] /. t -> 1 + x (* Out= (1 + x) Hypergeometric2F1[1/3, 1/2, 4/3, -2 (1 + x)^3] *) Hence we find a hypergeometric function instead of the "bizarre" elliptic integral.Proof D[ad2F1[x], x] // Simplify (* Out= 1/Sqrt[1 + 2 (1 + x)^3] *) 
Posted 10 years ago
 Did somebody notice already that the simple substitution x+1->y, dx->dy, {x,0,1} -> {y,1,2} leads to the correct result?In fact: sol2 = Integrate[1/Sqrt[2 y^3 + 1], {y, 1, 2}] (* Out= Sqrt Hypergeometric2F1[1/6, 1/2, 7/6, -(1/2)] - Hypergeometric2F1[1/6, 1/2, 7/6, -(1/16)] *) sol2 // N (* Out= 0.376065 *) solN = NIntegrate[1/Sqrt[2 (1 + x)^3 + 1], {x, 0, 1}] (* Out= 0.376065 *) $Version (* Out= "8.0 for Microsoft Windows (64-bit) (October 7, 2011)" *) Best regards, Wolfgang Posted 10 years ago  Yes, for the indefinite integral you're right; it is still in Version 10.0.0 the case. The difference of valued indefinite integrals does not convert into the Hypergeometric2F1[] form and the bug reappears In:= N[(ReplaceAll[#, y -> 2] - ReplaceAll[#, y -> 1]) &[Integrate[1/Sqrt[2 y^3 + 1], y]]] Out= -2.96267 - 1.92762 I In:= N[FullSimplify[(ReplaceAll[#, y -> 2] - ReplaceAll[#, y -> 1])] &[Integrate[1/Sqrt[2 y^3 + 1], y]]] Out= -2.96267 - 1.92762 I In:= (Integrate[1/Sqrt[2 y^3 + 1], y] /. y -> x + 1) === Integrate[1/Sqrt[2 (1 + x)^3 + 1], x] Out= True that means this is a bug not about a path singularity, because the indefinite integral does not go along a path. It's something about turning back to the right leaf. For example, look at this, the derivative does not turn back to the original integrand, but presents an equation of six-th degree instead: Posted 10 years ago  One gets the intended result with a simple shift of the integration variable of the definite integral In:= NIntegrate[1/Sqrt[2 (1 + x)^3 + 1], {x, 0, 1}] Out= 0.376065 In:=$Version Out= "10.0 for Microsoft Windows (64-bit) (June 29, 2014)" Let y = x + 1, dy = dx and In:= NIntegrate[1/Sqrt[2 y^3 + 1], {y, 1, 2}] Out= 0.376065 In:= Integrate[1/Sqrt[2 y^3 + 1], {y, 1, 2}] Out= Sqrt*Hypergeometric2F1[1/6, 1/2, 7/6, -1/2] - Hypergeometric2F1[1/6, 1/2, 7/6, -1/16] In:= % // N Out= 0.376065 With other words, from 1 to 2 one does not encounter a singularity like the one Daniel mentioned, because In:= (Integrate[1/Sqrt[2 y^3 + 1], y] /. y -> x + 1) === Integrate[1/Sqrt[2 (1 + x)^3 + 1], x] Out= True of course.
Posted 10 years ago
 Dear Udo, thanks for your comments.But if you try (I am using Mathematica 9)In:= sol = Integrate[1/Sqrt[2 y^3 + 1], y]Out= (1/(3^(1/4) Sqrt[1 + 2 y^3]))(-1)^(1/6) 2^( 2/3) Sqrt[(-1)^(5/6) (-1 + (-2)^(1/3) y)] Sqrt[ 1 + (-2)^(1/3) y + (-2)^(2/3) y^2] EllipticF[ArcSin[Sqrt[-(-1)^(5/6) - I (-2)^(1/3) y]/3^(1/4)], (-1)^( 1/3)]In:= (sol /. y -> 2) - (sol /. y -> 1) // NOut= -2.96267 - 1.92762 IYou still have a wrong answer. So, in my opnion, there is a bug in Integrate function.
Posted 10 years ago
 Although possible (but you commented that you think not) it's unlikely that they are related--there are a variety of crashing bugs that one will come across in Mathematica (not frequently, but they happen) and I am guessing that they happen because of things deep in the works.
Posted 10 years ago
 Regarding the crash associated with plotting a long expression, I wonder if the bug is related to this post: http://community.wolfram.com/groups/-/m/t/314527?p_p_auth=3uJCatzrPostscript: I tried the plot in http://community.wolfram.com/groups/-/m/t/314527?p_p_auth=3uJCatzr, but using a shorter expression for the first argument to Plot. The plot in that post still crashes (randomly, it seems). So, perhaps the bugs are not relatives.
Posted 10 years ago
 No crash in Mathematica 9
Posted 10 years ago
 I filed a bug report.
Posted 10 years ago
 Yes, that's what happens on my computer also. Weird.
Posted 10 years ago
 Plotting theDerivative using Plot[theDerivative,{x,0,1}] does not crash my Kernel, but plotting the full expression for theDerivative as written in my post does crash it reproducibly.
Posted 10 years ago
 If I calculatesln = DSolve[{y'[x] == 1/Sqrt[2 (1 + x)^3 + 1], y == 0}, y, x]yp[x_] = D[y[x] /. sln, x]Plot[yp[x] - 1/Sqrt[2 (1 + x)^3 + 1], {x, 0, 1}]works giving me a plot that stays with 1.5 x 10^-15Also, plotting "theDerivative" does not crash my computer.
Posted 10 years ago
 A side observation...If one computes theIntegral = Integrate[1/Sqrt[2 (1 + x)^3 + 1], x] and then computes theDerivative = D[theIntegral, x] // FullSimplify which yields -(((-1)^(5/6) Sqrt Sqrt[ 1 + 2 (1 + x)^3])/(Sqrt[-(-1)^(5/6) (1 + 2^(1/3) + 2^(1/3) x)] Sqrt[(-1)^(5/6) (-1 + (-2)^(1/3) (1 + x))] Sqrt[ 1 + (-2)^(1/3) (1 + x) + (-2)^(2/3) (1 + x)^2] Sqrt[ 3 - (-1)^(1/6) Sqrt + (1 + x) Root[108 + #1^6 &, 1]] Sqrt[ 3 + (-1)^(5/6) Sqrt + (1 + x) Root[108 + #1^6 &, 2]])) and then attempts to Plot this using Plot[-(((-1)^(5/6) Sqrt Sqrt[ 1 + 2 (1 + x)^3])/(Sqrt[-(-1)^(5/6) (1 + 2^(1/3) + 2^(1/3) x)] Sqrt[(-1)^(5/6) (-1 + (-2)^(1/3) (1 + x))] Sqrt[ 1 + (-2)^(1/3) (1 + x) + (-2)^(2/3) (1 + x)^2] Sqrt[ 3 - (-1)^(1/6) Sqrt + (1 + x) Root[108 + #1^6 &, 1]] Sqrt[ 3 + (-1)^(5/6) Sqrt + (1 + x) Root[108 + #1^6 &, 2]])) , {x, 0, 1}] the result is that it crashes the Mathematica Kernel. (I will report this.) In:= SystemInformation["Small"] Out= {"Kernel" -> {"SystemID" -> "MacOSX-x86-64", "ReleaseID" -> "10.0.0.0 (5098698, 5098537)", "CreationDate" -> DateObject[{2014, 6, 29}, TimeObject[{20, 38, 32}]]}, "FrontEnd" -> {"OperatingSystem" -> "MacOSX", "ReleaseID" -> "10.0.0.0 (5098698, 2014062704)", "CreationDate" -> DateObject[{2014, 6, 27}, TimeObject[{23, 17, 40.}]]}} 
Posted 10 years ago
 Slipping in Evaluate prevents the crash. fyi. Plot[ Evaluate[-(((-1)^(5/6) Sqrt Sqrt[ 1 + 2 (1 + x)^3])/(Sqrt[-(-1)^(5/6) (1 + 2^(1/3) + 2^(1/3) x)] Sqrt[(-1)^(5/ 6) (-1 + (-2)^(1/3) (1 + x))] Sqrt[ 1 + (-2)^(1/3) (1 + x) + (-2)^(2/3) (1 + x)^2] Sqrt[ 3 - (-1)^(1/6) Sqrt + (1 + x) Root[108 + #1^6 &, 1]] Sqrt[ 3 + (-1)^(5/6) Sqrt + (1 + x) Root[108 + #1^6 &, 2]]))], {x, 0, 1}] Ditto pre-processing the expression with FullSimplify. expr2 = FullSimplify[ expr , x >= 0] I don't want to discourage your reporting the crash.
Posted 10 years ago
 Interesting. I've reported the crash and it's in the bug system. And I just added your example to the report.
Posted 10 years ago
 Much simpler example that causes the crash: Plot[Root[108 + #1^6 &, 2], {x, 0, 1}] 
Posted 10 years ago
 @otsocol, Yes, I was looking at that antiderivative.As to whether the antiderivative itself is correct, About all I can say is it gives a numerical result consistent with NIntegrate if one only integrates from 0 to 1/2 (which excludes that singularity point). This might or might not answer your question.
Posted 10 years ago
 Hi Daniel. Thanks for your answer.When you say that the antiderivative has a jump discontinuity near x=0.59, I think you are looking at the Mathematica answerIn:= Integrate[1/Sqrt[2 (1 + x)^3 + 1], x]Out= ((-1)^(1/6) 2^(2/3) Sqrt[(-1)^(5/6) (-1 + (-2)^(1/3) (1 + x))] Sqrt[1 + (-2)^(1/3) (1 + x) + (-2)^(2/3) (1 + x)^2] EllipticF[ ArcSin[Sqrt[-(-1)^(5/6) - I (-2)^(1/3) (1 + x)]/3^(1/4)], (-1)^( 1/3)])/(3^(1/4) Sqrt[1 + 2 (1 + x)^3])Is that correct? But is it possible to trust in this result in this case?
Posted 10 years ago
 I see what you meansln = DSolve[{y'[x] == 1/Sqrt[2 (1 + x)^3 + 1], y == 0}, y]Plot [Im[y[x] /. sln], {x, 0, 1}, PlotRange -> All ] shows that y[x] develops an imaginary part at 0.59
Posted 10 years ago
 sln = DSolve[{y'[x] == 1/Sqrt[2 (1 + x)^3 + 1], y == 0}, y] sln represents a syntax error as such. If one is concerned about newbies and rookies not sticking to the syntax, one should stick to it at least oneself.
Posted 10 years ago
 Frank: In brief: Mathematica is using the Newton-Leibniz method to do the definite integral. It computes the antiderivative. It looks for points on the path of integration where the antiderivative might have discontinuities. It evaluates by taking limits at endpoints and all these points.The antiderivative has a jump discontinuity near x=0.59, which of course is between 0 and 1. The code that looks for discontinuities is failing to find it. Not entirely a surprise (finding them involves some hard-core equation solving, often transcendental). But it's worth further investigation.
Posted 10 years ago
 Daniel, I don't understand your post.
Posted 10 years ago
 There is a path singularity around .59 that is not found. Will investigate.
Posted 10 years ago
 Using Mathematica 10In:= Integrate[1/Sqrt[2 (1 + x)^3 + 1], {x, 0, 1}]Out= (1/(3^(1/4)))(-1)^(1/12) 2^( 2/3) (EllipticF[ ArcSin[((-1)^(11/12) Sqrt[1 + 2^(1/3)])/3^(1/4)], (-1)^(1/3)] - EllipticF[ ArcSin[((-1)^(11/12) Sqrt[1 + 2 2^(1/3)])/3^(1/4)], (-1)^(1/3)])In:= N[%]Out= -2.96267 - 1.92762 ISeeing -1 raised to fractional powers makes a branch cut problem a possibility.