Real integral giving imaginary answer

Question

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?

Answer 1

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"]

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

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

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"]

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));

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"]

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);

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

Answer 2

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:

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.