complex number from the integration of a real integrand

Question

Consider this integral: $$ \int_0^\pi d\theta \sin\theta\, {1 \over \gamma^q}\left((1-\gamma^2)^2\sqrt{1 \over 1-\left(1-(1-\gamma^2)^4\right)\cos^2\theta} - 1\right)^q $$ for $0\leq \gamma < 1$. The integrand is real, but the result of the integration is a complex number.

Here is the code:

q = 14;
A1 = Assuming[
   a \[Element] Reals && 0 <= a < 1 && b \[Element] Reals && 
    0 <= b < 1 && theta >= 0 && theta \[Element] Reals,
   Integrate[Sin[theta] (b Sqrt[1/(1 - a Cos[theta]^2)] - 1)^q, 
    theta]];
B1 = Assuming[
   a \[Element] Reals && 0 <= a < 1 && b \[Element] Reals && 
    0 <= b < 1, 
   Limit[A1, theta -> Pi, Direction -> "FromBelow"] - 
    Limit[A1, theta -> 0, Direction -> "FromAbove"]
   ];
C1 = 1/gamma^q B1 /. {a -> (1 - (1 - gamma^2)^4), 
    b -> (1 - gamma^2)^2};

Print["q = ", q, " --> C1/.gamma\[Rule] 0.1 = ", C1 /. gamma -> 0.1];
Print["q = ", q, " --> C1/.gamma\[Rule] 0.5 = ", C1 /. gamma -> 0.5];
Print["q = ", q, " --> C1/.gamma\[Rule] 0.9 = ", C1 /. gamma -> 0.9];

and the results are

q = 14 --> C1/.gamma-> 0.1 = 454.747 +0.355271 I
q = 14 --> C1/.gamma-> 0.5 = 0.101885 +7.27596*10^-12 I
q = 14 --> C1/.gamma-> 0.9 = -25.9344+4.85305*10^-16 I

As we see the results are complex numbers and especially for q=14, and gamma = 0.1 the imaginary part is not negligible. So, the question is when the integrand is real, why the result of the integration is complex and has an imaginary part?

Any idea what is the problem and how to fix it?

Answer 1

Mathematica's antiderivative has poles and or not continuous, so FTOC can not be used. Either that (I have not figured which it is), or the anti-derivative is simply wrong.

But using Rubi, it produces anti-derivative which allows using FTOC on it and does not produce complex values. And now the result of integration is real.

<< Rubi`
q = 14;
Arubi = Int[Sin[theta] (b Sqrt[1/(1 - a Cos[theta]^2)] - 1)^q, theta];
B1 = Assuming[
  a \[Element] Reals && 0 <= a < 1 && b \[Element] Reals && 
   0 <= b < 1, 
  Limit[Arubi, theta -> Pi, Direction -> "FromBelow"] - 
   Limit[Arubi, theta -> 0, Direction -> "FromAbove"]];

C1 = 1/gamma^q B1 /. {a -> (1 - (1 - gamma^2)^4),  b -> (1 - gamma^2)^2}
Print["q = ", q, " --> C1/.gamma\[Rule] 0.1 = ", C1 /. gamma -> 0.1]
Print["q = ", q, " --> C1/.gamma\[Rule] 0.5 = ", C1 /. gamma -> 0.5]
Print["q = ", q, " --> C1/.gamma\[Rule] 0.9 = ", C1 /. gamma -> 0.9]

For an example, for gamma -> 0.5, here is Rubi's anitderivative

  Arubi0 = Arubi /. {a -> (1 - (1 - gamma^2)^4), b -> (1 - gamma^2)^2};
  Plot[Arubi0 /. gamma -> 0.5, {theta, -4 Pi, 4 Pi}]

how to fix it

I do not know now how to fix it in Mathematica if any. But you can try Rubi to integrate this.

Update to answer comment

Is there any way to test which one is correct?

One way is to check the derivative of the antiderivative gives back the integrand, or something free from the variable of inegration.

This below shows Mathematic's result does not do this. While Rubi's does.

To make it easier to do, I used this substitution {a -> (1 - (1 - gamma^2)^4), b -> (1 - gamma^2)^2} /. gamma -> 1/2 on both, otherwise, it will take very long time. This function below returns 1 if it verifies, else 0.

verifyAntiDerivative[anti_, integrand_, x_] := Module[{tmp},
  tmp = D[anti, x] - integrand;
  If[FreeQ[Simplify[tmp], x], Return[1, Module]];
  tmp = D[Simplify[anti], x] - integrand;
  If[FreeQ[Simplify[tmp], x], Return[1, Module]];
  tmp = Simplify[D[anti, x]] - integrand;
  If[FreeQ[Simplify[tmp], x], Return[1, Module]];
  0];

Now do the following

<< Rubi`
q = 14;
integrand = Sin[theta] (b Sqrt[1/(1 - a Cos[theta]^2)] - 1)^q;
anti = Int[integrand, theta];
anti0 = anti /. {a -> (1 - (1 - gamma^2)^4), b -> (1 - gamma^2)^2}/. gamma -> 1/2
integrand0 =integrand/.{a ->(1 - (1 -gamma^2)^4),b ->(1 - gamma^2)^2} /. gamma -> 1/2;
verifyAntiDerivative[anti0, integrand0, theta]

(*1*)

Now do the same for Mathematica

q = 14;
integrand = Sin[theta] (b Sqrt[1/(1 - a Cos[theta]^2)] - 1)^q;
anti = Integrate[integrand, theta];
anti0 = anti /. {a -> (1 - (1 - gamma^2)^4), b -> (1 - gamma^2)^2} /.gamma -> 1/2
integrand0=integrand/.{a -> (1 -(1 - gamma^2)^4),b -> (1 - gamma^2)^2}/. gamma -> 1/2;
verifyAntiDerivative[anti0, integrand0, theta]

  (* 0 *)