Bug in mathematica analytic integration?

Question

Bug introduced in 9.0 or earlier and fixed in 13.0 or earlier

---

I found Mathematica provides me a wrong answer for a relatively simple analytically solvable integral:

Integrate[Exp[-I θ]/(1 + b Cos[θ]), {θ, 0, 2 π}, 
     Assumptions -> b < 1 && b > -1]

Provides as a result (2π)/b, which is incorrect. The integral can be computed easily writing the integrand as a geometric series.

Additionally, you can compute it numerically

NIntegrate[Exp[-I θ]/(1 + b Cos[θ]), {θ, 0, 2 π}]

which produces the right result, different from (2π)/b.

Is this a very bad bug? or am I doing something wrong? What could be the reason for this error?

Answer 1

This is indeed a serious and problematic issue.
We know many similar problems with symbolic integration which provides Integrate. There were some improvments in newer versions of the system but also some issues become worse, see e.g. Mathematica 9 can't integrate this function but earlier versions could.

). One can find more problems looking for tags calculus-and-analysis and bugs.
The bigger problem is that Wolfram Research Inc. seems not to be sure how to remove completely such bugs from the system.
There are a few things to say:

• possibly most general symbolic integration is a really difficult issue which needs many supplementing rules.
• basically there is no system free of any bugs (see e.g. The Art of Computer Programming by Donald Knuth).
• a natural way of development of such a sophisticated system as Mathematica assumes a step by step process rather than finding a universal key opening all doors.

The problem we encounter here most likely comes from the fact that underlying complex functions behind the symbolic result are not defined in the whole complex plane, see closely related post: Why does Integrate declare a convergent integral divergent?.
We can see that we can compute the indefinite integral:

f[θ_] = FullSimplify[ Integrate[ Exp[-I θ]/(1 + b Cos[θ]), θ], 
                                   0 <= θ <= 2 π && -1 < b < 1]

-((I ((2 ArcTan[(1 + b E^(-I θ))/Sqrt[-1 + b^2]])/Sqrt[-1 + b^2] 
- Log[ 2 E^(-I θ) (1 + b Cos[θ])]))/b)

In fact we have

D[ f[θ], θ] // FullSimplify

E^(-I θ)/(1 + b Cos[θ])

However as b is close to -1 or 1 then this term 1 + b Cos[θ] in the denominator becomes large when θ is close to 0 and then ArcTan and Log functions should be supplemented with an appropriate phase rule to get reliable results.

Interestingly this definite integral could be computed in Mathematica 7

Integrate[ Exp[-I θ]/(1 + b Cos[θ]), {θ, 0, 2 π}, Assumptions -> -1 < b < 1]

(2 (1 - 1/Sqrt[1 - b^2]) π)/b

The last result seems to be numerically compatibile with that of NIntegrate however it could be computed with certain phase rule which in general would be unnecessarily too restrictive yielding incorrect results in another cases. Taking into account the above remarks we would rather say that these incorrect results are rather due to imperfectness of symbolic integration (expected in future versions as well, in different forms though) than just a simple bugs.

Answer 2

Artes gives a good overview of the problems inherent in symbolic integration. Changing the form of the integral can yield different results, but then the problem is determining the correct one.

Using ExpToTrig on the integrand yields a different result:

ans1 = Integrate[Exp[-I θ]/(1 + b Cos[θ]),
  {θ, 0, 2 π}, Assumptions -> b < 1 && b > -1]
ans2 = Integrate[Exp[-I θ]/(1 + b Cos[θ]) // ExpToTrig,
  {θ, 0, 2 π}, Assumptions -> b < 1 && b > -1]
(*
  (2 π)/b
  ConditionalExpression[(2 (1 - 1/Sqrt[1 - b^2]) π)/b, b != 0]
*)

We can compare the two answers. Turning off symbolic preprocessing ensures that the NIntegrate will use the formula just as we pass it. (One can achieve a similar effect by defining f[θ_?NumericQ] := Exp[-I θ]/(1 + b Cos[θ]) and using f[θ] for the integrand.

Block[{b = 0.5},
 {NIntegrate[Exp[-I θ]/(1 + b Cos[θ]),
   {θ, 0, 2 π}, Method -> {Automatic, "SymbolicProcessing" -> 0}],
  ans1, ans2}
 ]
(*
  {-1.94402 + 3.65803*10^-16 I, 12.5664, -1.94402}
*)

So the second method seems to give a correct result. More values of b could be checked for further verification.

Another way to modify the form of the integral is to change the interval:

Integrate[Exp[-I θ]/(1 + b Cos[θ]),
 {θ, -π, π}, Assumptions -> b < 1 && b > -1]
(*
  ConditionalExpression[(
   2 (1 - 1/Sqrt[1 - b^2]) π)/b, (b != 0 && Re[1/Sqrt[b]] != 0) || b < 0]
*)

Or one can break down the real and imaginary parts (the third one, with the parts combine is basically what ExpToTrig does):

Integrate[Cos[θ]/(1 + b Cos[θ]),
 {θ, 0, 2 π}, Assumptions -> b < 1 && b > -1]
Integrate[Sin[θ]/(1 + b Cos[θ]),
 {θ, 0, 2 π}, Assumptions -> b < 1 && b > -1]

Integrate[(Cos[θ] - I Sin[θ])/(1 + b Cos[θ]),
 {θ, 0, 2 π}, Assumptions -> b < 1 && b > -1]

(*
  ConditionalExpression[(2 (1 - 1/Sqrt[1 - b^2]) π)/b, b != 0]
  ConditionalExpression[0, b != 0]

  ConditionalExpression[(2 (1 - 1/Sqrt[1 - b^2]) π)/b, b != 0]
*)

One thing worth noting is that in this form, the integrands are real-valued functions (ignoring the complex coefficient in the third form) of real variables. In fact, they are rational functions of sine and cosine, which can be converted to rational function via substitution, although I do not know how Mathematica handles them. The Exp[I θ] probably causes Mathematica to invoke a different algorithm to deal with functions of a complex variable.

If one wants, may do the substitution θ -> 2 ArcTan[t], although this will be equivalent to the integral from -π to π:

Integrate[TrigExpand@FullSimplify[
    Exp[-I θ]/(1 + b Cos[θ]) Dt[θ] /.
      θ -> 2 ArcTan[t], t ∈ Reals] /. Dt[t] -> 1,
 {t, -Infinity, Infinity}, Assumptions -> -1 < b < 1]
(*
  (2 (π - π/Sqrt[1 - b^2]))/b
*)