Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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?

share|improve this question
1  
In case it helps, Maxima gives 2*%pi/b * (2*sqrt(1-b^2)+b^2-2)/(sqrt(1-b^2)+b^2-1). When computing with +i instead of -i in the exponential, it gives instead 2*%pi/b * (sqrt(1-b^2)-1)/sqrt(1-b^2), but it's actually the same. –  Jean-Claude Arbaut May 19 at 10:36
    
Have a play with this Manipulate to see the dependence on b of the singularities of the integrand in the complex Theta plane: Manipulate[ContourPlot[Abs[Exp[-I Theta]/(1 + b Cos[Theta]) /. {Theta -> ThetaR + I ThetaI}], {ThetaR, -1, 2 Pi + 1}, {ThetaI, -Pi - 1, Pi + 1}, PlotRange -> {Automatic, 30}, Contours -> 50, Epilog -> {Red, Thick, Line[{{0, 0}, {2 Pi, 0}}]}], {{b, 0.5}, -1, 1}]. You need to ensure that the path of integration goes appropriately around the singularities in order to get the intended result. –  Stephen Luttrell May 19 at 11:07
    
@StephenLuttrell Is there a way to Tell Mathematica that theta is real, and integration is done on the real interval (0,2pi). I don't have Mma at hand right now to do the test, but there seem to be an option Assumption->Element[x,Reals] –  Jean-Claude Arbaut May 19 at 11:44
    
The Manipulate in my previous comment shows that the position of the singularities of the integrand depends strongly on the value of b, so I can't see how any automatic integration scheme - e.g. integration along real Theta - will necessarily give you the "right" answer. In general, the answer depends on how the integration path winds around the singularities. However, I concede that in this particular case I haven't looked closely at exactly what is going on, other than to observe the "wandering singularities", which alerts me to the need to take care. –  Stephen Luttrell May 19 at 11:59
    
I found a (dubious, but it works) series expansion in b approach that works. Compute the coefficient of b^n using coeff[n_] = Integrate[Exp[-I Theta]/(1 + b Cos[Theta]) // SeriesCoefficient[#, {b, 0, n}] & // Evaluate, {Theta, 0, 2 Pi}, Assumptions -> (n >= 0 && n \[Element] Integers)] // FullSimplify - this takes a while to compute. Sum the series using integral[b_] = (Sum[b^n coeff[n], {n, 0, Infinity}]) // FullSimplify, which gives ((-2 + 2 b^4 + (-1 + b^2) Sqrt[1 + b^2] + Sqrt[1 + b^2] Sqrt[1 - b^4] + Sqrt[2 - 2 b^4] Sqrt[1 + Sqrt[1 - b^4]]) \[Pi])/(b(-1 + b^4)). –  Stephen Luttrell May 19 at 12:28

2 Answers 2

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.

share|improve this answer
1  
(+1) "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". I am totally agree. I can also say that the NIntegrate framework is implemented through the step-by-step approach in rather transparent and controllable way for a user (although insufficiently documented as always in Mathematica). –  Alexey Popkov May 19 at 16:21

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
*)
share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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