Error in trigonometric definite integration

[StanczakDominik]

I think this could be another case of #20360.

Consider this snippet:

>>> from sympy import *
>>> from sympy.abc import x
>>> a = symbols('a', positive=True)
>>> expr = (1 + a * cos(x))**-1

note that this expression, for a < 1, is always positive. With a concrete value this results in a positive result:

>>> integrate(expr.subs(a, Rational(1, 2)), (x, 0, 2 * pi))
4*sqrt(3)*pi/3

but with the general a, we get a result that doesn't actually depend on a:

>>> integrate(expr, (x, 0, 2 * pi))
0

I'd like to help solve this however I can!

Note: I've tested this on a bunch of versions between 1.6.2 and 1.0. The result is the same in each of them.

[ghost]

@StanczakDominik integral will be zero as you can write

cos(x) = (e^(ix) + e^(-ix))/2 using eulers formula

substitute it in the integral , you will get some function of e^(ix) inside the integral

now when you substitute e^(ix) = t

as limits of x are from 0 to 2pi

limits of t will be from e^(i * 0) to e^(i * 2pi) i.e from 1 to 1

integral with the same limits is eventually zero

[StanczakDominik]

There must be an error in your argument, as - for a < 1 - the function is everywhere positive and we are integrating it over a non zero interval. So that can't be right. It doesn't make sense graphically.

Also, if that was right, you'd expect the a = 1/2 result to be zero. It clearly isn't.

I've also done this numerically, instead. I seem to get, off the top of my head as i just got up, something like 2pi/sqrt(1-a^2), with a blow up at a = 1.

[gunvirranu]

@akshatsood2249 Unfortunately that doesn't work because you're moving your integration into the complex plane and converting it into a contour integral along the unit circle. Even though the path is closed as you hinted, the integrand is not analytic inside the unit circle, so the line integral doesn't necessarily go to zero. If you're interested, you can look into some complex analysis, very interesting stuff.

[ghost]

@StanczakDominik you are right brother, my approach is wrong, the integral boils down to pi/sqrt(1-a^2). @gunvirranu will surely have a look at complex stuff

[StanczakDominik]

I've added a test for this in #20890.

I also found a derivation (using the tan(y/2) substitution) for a related case that confirms my 2pi/sqrt(1-a^2) numerical result.

Also, I suspect this needs a "Wrong Result" label :(

[ghost]

Here first notice that a should lie in (-1,1) else the integral blows up, since then we end up having 1/0 .
0 is actually the principal cauchy value. This is some text on the same.

Just try for yourself, what happens when a become greater than 1. The computed integral in the pr will be correct only when a lies in (-1,1) , since the integral can only be calculated when a lies in (-1,1)

The integral in the link is actually different than the original integral, infact very different. The original integral has a as the coefficient while the the one in the link has 1/a as the coefficient. In the latter case a ends up having values in the domain
(-oo,-1) U (1,oo)

[StanczakDominik]

I might not have said so but yes, I completely agree that for a outside of (-1, 1), the integral blows up. But Sympy should still return the correct value for a in (-1, 1), no?

[ghost]

You mean the function should return an integral with the the cases being handled seperately.
somehwat like this ->

integrate( ( 1+a*cos)**-1,x) 
>>>2*a/(sqrt(1-a**2))    , -1<x<1
    0       ,  otherwise

[StanczakDominik]

Yeah, I suspect a Piecewise result would make the most sense!

[oscargus]

In the current master this gives an antiderivative of

This evaluates to the same value for both 0 and 2*pi.

However, in the current version this doesn't happen, but an exception is thrown:

Traceback (most recent call last):

  File "C:\Users\Oscar\sympy\sympy\solvers\solveset.py", line 597, in _solve_trig
    sol = _solve_trig1(f, symbol, domain)

  File "C:\Users\Oscar\sympy\sympy\solvers\solveset.py", line 705, in _solve_trig1
    raise _SolveTrig1Error("polynomial solutions must form FiniteSet")

_SolveTrig1Error: polynomial solutions must form FiniteSet


During handling of the above exception, another exception occurred:

Traceback (most recent call last):

  File "C:\Users\Oscar\AppData\Local\Temp/ipykernel_13956/1429305692.py", line 1, in <module>
    integrate(expr, (x, 0, 2 * pi))

  File "C:\Users\Oscar\sympy\sympy\integrals\integrals.py", line 1553, in integrate
    return integral.doit(**doit_flags)

  File "C:\Users\Oscar\sympy\sympy\integrals\integrals.py", line 486, in doit
    did = self.xreplace(reps).doit(**hints)

  File "C:\Users\Oscar\sympy\sympy\integrals\integrals.py", line 697, in doit
    evalued_pw = piecewise_fold(Add(*piecewises))._eval_interval(x, a, b)

  File "C:\Users\Oscar\sympy\sympy\functions\elementary\piecewise.py", line 551, in _eval_interval
    irv = self._handle_irel(x, handler)

  File "C:\Users\Oscar\sympy\sympy\functions\elementary\piecewise.py", line 419, in _handle_irel
    expr = handler(self.xreplace(reps))

  File "C:\Users\Oscar\sympy\sympy\functions\elementary\piecewise.py", line 550, in handler
    return ipw._eval_interval(x, lo, hi)

  File "C:\Users\Oscar\sympy\sympy\core\expr.py", line 980, in _eval_interval
    singularities = singularities | solveset(logterm.args[0], x,

  File "C:\Users\Oscar\sympy\sympy\solvers\solveset.py", line 2217, in solveset
    rv = solveset(f.xreplace({symbol: x}), x, domain)

  File "C:\Users\Oscar\sympy\sympy\solvers\solveset.py", line 2241, in solveset
    return _solveset(f, symbol, domain, _check=True)

  File "C:\Users\Oscar\sympy\sympy\solvers\solveset.py", line 1026, in _solveset
    result = _solve_trig(f, symbol, domain)

  File "C:\Users\Oscar\sympy\sympy\solvers\solveset.py", line 600, in _solve_trig
    sol = _solve_trig2(f, symbol, domain)

  File "C:\Users\Oscar\sympy\sympy\solvers\solveset.py", line 765, in _solve_trig2
    solns = solveset(g, y, S.Reals) - solveset(h, y, S.Reals)

  File "C:\Users\Oscar\sympy\sympy\solvers\solveset.py", line 2217, in solveset
    rv = solveset(f.xreplace({symbol: x}), x, domain)

  File "C:\Users\Oscar\sympy\sympy\solvers\solveset.py", line 2241, in solveset
    return _solveset(f, symbol, domain, _check=True)

  File "C:\Users\Oscar\sympy\sympy\solvers\solveset.py", line 1072, in _solveset
    result += _solve_radical(equation, u,

  File "C:\Users\Oscar\sympy\sympy\solvers\solveset.py", line 870, in _solve_radical
    for s in result:

  File "C:\Users\Oscar\sympy\sympy\sets\sets.py", line 1428, in __iter__
    raise TypeError(

TypeError: The computation had not completed because of the undecidable set membership is found in every candidates.

[oscarbenjamin]

That's a separate issue in solveset. I've opened #22058 for that.