Wrong integration results of trigonometric functions

[smeunier-amd]

I have a lot of problems using trigonometric functions, Piecewise, Abs and Integral.

Here's the simplest case :

from sympy import *

theta = Symbol("theta", real=True)
theta_range = (theta, -pi / 2, pi / 2)

f = Abs(cos(theta))

I = Integral(f, theta_range)
I.doit()

It outputs :

0

I use sympy 1.10 and python 3.9.7.

[ThePauliPrinciple]

It seems to me that Piecewise doesn't take into account periodicity. If we redefine the problem in the range 0-2pi we get the right result:

>>> (Integral(f, (theta, 0, 2*pi)) - Integral(f, (theta, pi/2, 3*pi/2))).doit()
2

Not sure if this is the systematic problem you are encountering "I have a lot of problems".

>>> Integral(abs(cos(theta)), theta).doit()
Piecewise(
  (-sin(theta), theta <= 0), 
  (sin(theta), theta <= pi/2), 
  (2 - sin(theta), theta <= 3*pi/2), 
  (sin(theta) + 4, theta <= 2*pi), 
  (4 - sin(theta), True)
)

Note that wolfram alpha simply gives sin(x)*sign(cos(x)) (which seems to be correct) and it should be possible to provide this answer in SymPy too.

[smeunier-amd]

Indeed, a lot of my problems are about the periodicity. I was able to workaround a few of them like you did ! But in some cases, with equations using different shifted and scaled angles, splitting integrals rapidly become an hassle.

Here's another :

from sympy import *

theta = Symbol("theta", real=True)

f = Abs(sin(theta))

I = Integral(f, (theta, 0, 1 * pi))
display(I.doit())

I = Integral(f, (theta, 0, 2 * pi))
display(I.doit())

I = Integral(f, (theta, 0, 3 * pi))
display(I.doit())

I = Integral(f, (theta, 0, 4 * pi))
display(I.doit())

It outputs :

2
4
2
4

[smeunier-amd]

Another one leading to an exception :

from sympy import *

theta = Symbol("theta", real=True)
phi = Symbol("phi", real=True)

f = cos(theta - phi)

I = Integral(f, (theta, 0, pi))
display(I.doit())

I = Integral(Abs(f), (theta, 0, pi))
display(I.doit())

It outputs :

2 * sin(phi)

KeyError
...

I don't expect any simplification from sympy.

KeyError log.

[oscarbenjamin]

Internally the integrand is rewritten as

p = Piecewise((cos(theta), cos(theta) >= 0), (-cos(theta), True))

sympy/sympy/integrals/integrals.py

Lines 532 to 539 in d09723f

	if function.has(Abs, sign) and (
	(len(xab) < 3 and all(x.is_extended_real for x in xab)) or
	(len(xab) == 3 and all(x.is_extended_real and not x.is_infinite for
	x in xab[1:]))):
	# some improper integrals are better off with Abs
	xr = Dummy("xr", real=True)
	function = (function.xreplace({xab[0]: xr})
	.rewrite(Piecewise).xreplace({xr: xab[0]}))

The condition cos(x)>=0 is awkward.

That then gives a valid antiderivative (with the same awkward condition):

In [1]: x = Symbol('x', real=True)

In [2]: f = Abs(cos(x))

In [3]: f
Out[3]: │cos(x)│

In [4]: f.rewrite(Piecewise)
Out[4]: 
⎧cos(x)   for cos(x) ≥ 0
⎨                       
⎩-cos(x)    otherwise   

In [5]: f.rewrite(Piecewise).piecewise_integrate(x)
Out[5]: 
⎧sin(x)   for cos(x) ≥ 0
⎨                       
⎩-sin(x)    otherwise 

That then goes wrong here:

sympy/sympy/integrals/integrals.py

Line 711 in d09723f

evalued_pw = piecewise_fold(Add(*piecewises))._eval_interval(x, a, b)

In [4]: f.rewrite(Piecewise).piecewise_integrate(x)._eval_interval(x, -pi/2, pi/2)
Out[4]: 0

At this point we get some weird output:

sympy/sympy/functions/elementary/piecewise.py

Lines 611 to 615 in d09723f

	# handle a Piecewise with lo <= hi and no x-independent relationals
	# -----------------------------------------------------------------
	ok, abei = self._intervals(x)
	if not ok:
	from sympy.integrals.integrals import Integral

-> if not ok:
(Pdb) l
611  	                return rv
612  	
613  	        # handle a Piecewise with lo <= hi and no x-independent relationals
614  	        # -----------------------------------------------------------------
615  	        ok, abei = self._intervals(x)
616  ->	        if not ok:
617  	            from sympy.integrals.integrals import Integral
618  	            # not being able to do the interval of f(x) can
619  	            # be stated as not being able to do the integral
620  	            # of f'(x) over the same range
621  	            return Integral(self.diff(x), (x, lo, hi))  # unevaluated
(Pdb) p ok
True
(Pdb) p abei
[(0, pi/2, sin(x), 0), (3*pi/2, 2*pi, sin(x), 0), (-oo, oo, -sin(x), 1)]

I think that the intention is that the _intervals function would return a list of intervals over which the piecewise takes different forms e.g.:

In [4]: p = Piecewise((y, x>2), (z, x>1), (t, x>0))

In [5]: p
Out[5]: 
⎧y  for x > 2
⎪            
⎨z  for x > 1
⎪            
⎩t  for x > 0

In [6]: p._intervals(x)
Out[6]: (True, [(2, oo, y, 0), (1, oo, z, 1), (0, oo, t, 2)])

The result [(0, pi/2, sin(x), 0), (3*pi/2, 2*pi, sin(x), 0), (-oo, oo, -sin(x), 1)] is only valid for 0 <= x <= 2*pi though. I think that wrong result comes from here:

In [12]: from sympy.solvers.inequalities import _solve_inequality

In [13]: _solve_inequality(cos(x)>0, x)
Out[13]: 
⎛            π⎞   ⎛3⋅π              ⎞
⎜0 ≤ x ∧ x < ─⎟ ∨ ⎜─── < x ∧ x < 2⋅π⎟
⎝            2⎠   ⎝ 2               ⎠

That is called here:

sympy/sympy/functions/elementary/piecewise.py

Lines 681 to 684 in d09723f

	try:
	rv = _solve_inequality(r, sym)
	except NotImplementedError:
	return False, 'Unable to solve relational %s for %s.' % (r, sym)

Clearly if solve_inequality can't properly handle trigonometric functions then its output should not be trusted here without some checking.

[smeunier-amd]

I have some problems using solveset too. Is that related ? I was trying to cleanup an example or two.

[ThePauliPrinciple]

sympy/sympy/integrals/integrals.py

The comment in this code says "# some improper integrals are better off with Abs". Perhaps some proper integrals are also better off with Abs? (considering the sin(x)*sign(cos(x)) solution of wolfram)

[smeunier-amd]

Trying to integrate various equations, I tried to find intervals using solveset.

A simple example :

from sympy import *

theta = Symbol("theta", real=True)
phi = Symbol("phi", real=True)

solveset(cos(theta - phi) >= 0, theta, domain=Reals)

It outputs :

TypeError
...

TypeError log.

[jksuom]

sin(x)*sign(cos(x)) is only a "piecewise antiderivative" because it is discontinuous at each point where cos(x) is zero. It can be used to compute definite integrals only over intervals in which cos(x) does not change sign.

[ThePauliPrinciple]

I see, so that would have the same/a similar problem as the current piecewise function that SymPy gives?

[oscarbenjamin]

The integral of a piecewise integrand should be split into parts. The current code does that but fails for periodic integrals because the old solvers don't handle periodic equations properly.

Given an integral like Integral(Abs(f(x)), (x, a, b)) the obvious approach would be to solve for f(x)=0 in the interval (a, b) and then split the integral over several intervals like (a, x1), (x2, x3), ..., (xn, b). Doing that would potentially be better than rewriting Abs to Piecewise.

For periodic integrals the obvious approach would be to transform the integral so it is periodic with period 2*pi. Then the integral itself is over a first interval (a % (2*pi), 2*pi) a last interval (0, b % (2*pi)) and then some number of periods giving a multiple of the integral over (0, 2*pi). Then all integrals are in the range (0, 2*pi) which is currently handled correctly (at least in this case).

The OP example would transform to

In [1]: integrate(abs(cos(x)), (x, 3*pi/2, 2*pi)) + integrate(abs(cos(x)), (x, 0, pi/2))
Out[1]: 2

I think the inequality solvers should have an option to raise an error instead of returning incorrect results for trigonometric functions. The domain argument could also be made to handle these cases properly:

In [22]: solve_univariate_inequality(cos(x)>0, x, domain=Interval(-pi/2, pi/2))
Out[22]: 
            π
0 ≤ x ∧ x < ─
            2

[oscarbenjamin]

This should be changed:

sympy/sympy/solvers/inequalities.py

Lines 515 to 518 in d09723f

	inf, sup = domain.inf, domain.sup
	if sup - inf is S.Infinity:
	domain = Interval(0, period, False, True).intersect(_domain)
	_domain = domain

[oscargus]

The type error seems to relate to a symbolic Interval.