4
$\begingroup$

Hello I am trying to evaluate the following integral

Integrate[4279/Sqrt[0.6817 + 0.3183*(1 + x)^3], {x, 0, 20}]

my mathematica 9 gives me -17605.1 - 16224.2 I which is impossible, as the integrand is real in the range 0->20.

while Wolfram gives me the right result 10496.

What did I do wrong?

$\endgroup$
  • 2
    $\begingroup$ That is the result of Nintegrate. $\endgroup$ – b.gates.you.know.what Mar 4 '15 at 8:12
  • $\begingroup$ How do I tell mathematica to not use NIntegrate? Also I would like to define the integral as a function of x (replacing the upperlimit as x) and then plot, how do I do that? Thanks $\endgroup$ – HuShu Mar 4 '15 at 8:19
  • $\begingroup$ I meant NIntegrate gives the correct result so you should consider using it. As to the second question lookup NumericQ and PatternTest in the documentation. $\endgroup$ – b.gates.you.know.what Mar 4 '15 at 8:32
  • $\begingroup$ This effect (bug(?) in Integrate[]) is due to discontinous antiderivatives. The Topic has been discussed several times here. See e.g. mathematica.stackexchange.com/questions/56957/… $\endgroup$ – Dr. Wolfgang Hintze Mar 4 '15 at 9:50
  • 1
    $\begingroup$ FWIW, translation by 1 produces the correct result: Integrate[4279/Sqrt[0.6817 + 0.3183*(x)^3], {x, 1, 21}]. $\endgroup$ – Michael E2 Mar 4 '15 at 11:19
6
$\begingroup$

The antiderivative is correct, in the sense its derivative gives back the original integrand

Clear[x]
integrand = 4279/Sqrt[6817/10000 + 3183/10000*(1 + x)^3];
mmaResult = Integrate[integrand, x];
integrandBack = Simplify[D[mmaResult, x]];
Plot[{integrandBack, integrand}, {x, 0, 20}, PlotTheme -> "Detailed"]

Mathematica graphics

But the definite integral does not give the correct result. Using limits

Limit[mmaResult, x -> 20] - Limit[mmaResult, x -> 0] // N

Mathematica graphics

This leads to one suspecting that it is the evaluation of the antiderivative which is the problem.

Plot[{Re[mmaResult], Im[mmaResult]}, {x, 0, 20}, PlotRange -> All, PlotTheme -> "Detailed"]

Mathematica graphics

Notice that Mathematica is saying there is discontinuity around x=1. Let us compare this to Maple, which does give the correct result for the definite integral

restart; 
integrand:=4279/sqrt(6817/10000 + 3183/10000*(1 + x)^3):
rMaple:=simplify(int(integrand,x));

Mathematica graphics

And

 evalf(int(integrand,x=0..20));
      (*10495.99273*)

I Copied Maple's antiderivative to Mathematica, to plot and compare with Mathematica's

mapleResult = -((427900/69066360513)*I)*69066360513^(1/3)*
 (-I*(I*3^(1/2)*69066360513^(1/3) + 69066360513^(1/3) - 6366*x - 6366))^
   (1/2)*2^(1/2)*6817^(1/2)*((3183*x + 69066360513^(1/3) + 3183)/
    (I*3^(1/2) + 3))^(1/2)*(-I*(I*3^(1/2)*69066360513^(1/3) - 
   69066360513^(1/3) + 6366*x + 6366))^(1/2)*
   EllipticF[(1/130191066)*6817^(1/2)*3^(3/4)*2^(1/2)*69066360513^
    (1/3)*(-I*(I*3^(1/2)*69066360513^(1/3) + 69066360513^(1/3) - 
    6366*x - 6366))^(1/2), (1/2)*3^(1/2) + (1/2)*I]/(3183*x^3 + 
       9549*x^2 + 9549*x + 10000)^(1/2);

 Plot[{Re[mapleResult], Im[mapleResult]}, {x, 0, 20}, PlotRange -> All, 
       PlotTheme -> "Detailed"]

Mathematica graphics

So, even though both Maple and Mathematica antiderivatives when differentiated, give back the integrand, Mathematica's plot shows there is a discontinuity in the antiderivative between x=1 and x=2, while Maple's version does not have this. But is this really true?

To make sure this is not just the way Mathematica evaluates this, due to the EllipticF function, I copied Mathematica antiderivative to maple to see how Maple evaluates it

restart; 
with(MmaTranslator); #load the package
mmaResult:=FromMma(`(855800 (-1)^(1/6) (6817/1061)^(1/3) Sqrt[(-1)^(5/6) (-1 + 
(-(3183/6817))^(1/3) (1 + x))] Sqrt[ 1 + (-(3183/6817))^(1/3) (1 + x) + 
 (-(3183/6817))^(2/3) (1 + x)^2]  EllipticF[ArcSin[Sqrt[-(-1)^(5/6) - I 
 (-(3183/6817))^(1/3) (1 + x)]/3^(1/4)], (-1)^(1/3)])/(3^(7/12) Sqrt[6817 + 3183 
    (1 + x)^3])`);

plot(Re(mmaResult),x=0..20);

Mathematica graphics

Which is different from Mathematica plot of the same function! No discontinuity. Notice it starts the same, at around 15000, but at near x=1.5 it does not break, like with Mathematica case. The same function! So this is the difference.

So this mean that the problem is in Mathematica evaluation of the EllipticF part of the antiderivative. Let see what Maple thinks the definite integral should be, using Mathematica antiderivative:

evalf(limit(mmaResult,x=20)-limit(mmaResult,x=0));
      (* 10495.99220-0.7825782301e-5*I  *)

Which is the correct result!

So, the antiderivative generated by Mathematica is correct. It is the evaluation phase, where the issue is, i.e. the evaluation of the ellipticF function. This results in the complex result you saw. It looks like a branch cut issue. Since the anitderivative needs to be evaluated to give the definite integral, hence the problem you saw.

Maple 18.2, and Mathematica 10.02

$\endgroup$
  • 1
    $\begingroup$ Nice analysis! Perhaps the final statement should be that the OP should consider to report this to support@wolfram.com. $\endgroup$ – Sjoerd C. de Vries Mar 4 '15 at 12:18
1
$\begingroup$

It might be Errr... a bug or some singular points(not on the real axis) might be involved and Integrate chose a strange Integrate route.

define

f[x_] := Integrate[4279/Sqrt[0.6817 + 0.3183*(1 + x)^3], x]

then

fdat = Transpose@{Range[0, 20, 0.1], f[x] /. x -> Range[0, 20, 0.1]}

plot the f generated by Integrate:

ListLinePlot[{{#[[1]], Im@#[[2]]} & /@ fdat, {#[[1]], Re@#[[2]]} & /@ 
fdat}, PlotRange -> All, PlotStyle -> {Red, Green}, 
PlotLegends -> {Im, Re}]

get: enter image description here

we found that at around x=1.5~1.6,there is a jump in Re part and Large Im part occured. If you didn't pass through this discontinuity point, Integrate will give right answer.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.