sagemath / sage

Main repository of SageMath
https://www.sagemath.org
Other
1.45k stars 481 forks source link

erroneous integral #26138

Open 97c6705f-be81-4d4a-a395-b6442e123d13 opened 6 years ago

97c6705f-be81-4d4a-a395-b6442e123d13 commented 6 years ago

The following computation is wrong (Sage 7.5.1):

sage: integral(arccos(cos(x)/(1+2*cos(x))),(x,0,pi/2))
1/4*pi^2

Indeed, the actual value is 5/24*\pi^2 .

Reference: Paul J. Nahin, "Inside interesting integrals", Springer, confirmed by a numerical computation with Maple.

Upstream: Reported upstream. No feedback yet.

Component: symbolics

Keywords: integral

Issue created by migration from https://trac.sagemath.org/ticket/26138

seblabbe commented 6 years ago
comment:1

Which is also confirmed by numerical_integral:

sage: integral(arccos(cos(x)/(1+2*cos(x))),(x,0,pi/2))
1/4*pi^2
sage: n(_)
2.46740110027234
sage: numerical_integral(arccos(cos(x)/(1+2*cos(x))),0,pi/2)
(2.056167583560283, 2.2828045937564967e-14)
sage: n(5/24*pi^2)
2.05616758356028
fchapoton commented 5 years ago
comment:2

Apparently, maxima returns a wrong primitive.

sage: f = arccos(cos(x)/(1+2*cos(x)))
sage: g = f.integrate(x); g
2*arccos((sin(x)^2/(cos(x) + 1)^2 - 1)/((2*(sin(x)^2/(cos(x) + 1)^2 - 1)/(sin(x)^2/(cos(x) + 1)^2 + 1) - 1)*(sin(x)^2/(cos(x) + 1)^2 + 1)))*arctan(sin(x)/(cos(x) + 1))
sage: g.simplify_full()
2*arccos(cos(x)/(2*cos(x) + 1))*arctan(sin(x)/(cos(x) + 1))

EDIT: But maxima itself does not return this answer:

maxima.eval('integrate(arccos(cos(x)/(2*cos(x)+1)),x)')
kcrisman commented 5 years ago
comment:3

Maxima uses acos, as it turns out. Anyway, the problem here somehow involves domain:complex, found in sage/maxima_lib.py:

Maxima 5.42.2 http://maxima.sourceforge.net
using Lisp ECL 16.1.2
Distributed under the GNU Public License. See the file COPYING.
Dedicated to the memory of William Schelter.
The function bug_report() provides bug reporting information.
(%i1) display2d:false;

(%o1) false
(%i2) domain:complex;

(%o2) complex
(%i3) integrate(acos(cos(x)/(2*cos(x)+1)),x);

(%o3) 2*atan(sin(x)/(cos(x)+1))
       *acos((1-sin(x)^2/(cos(x)+1)^2)/((sin(x)^2/(cos(x)+1)^2+1)
                                       *((2*(1-sin(x)^2/(cos(x)+1)^2))
                                        /(sin(x)^2/(cos(x)+1)^2+1)
                                        +1)))

Interestingly, without that flag, the much-abused abs_integrate gave something different (no idea if it was right, of course):

Maxima 5.42.2 http://maxima.sourceforge.net
using Lisp ECL 16.1.2
Distributed under the GNU Public License. See the file COPYING.
Dedicated to the memory of William Schelter.
The function bug_report() provides bug reporting information.
(%i1) display2d:false;

(%o1) false
(%i2) integrate(acos(cos(x)/(2*cos(x)+1)),x);   

(%o2) 'integrate(acos(cos(x)/(2*cos(x)+1)),x)
(%i3) load(abs_integrate);

(%o3) "/Users/.../sage/local/share/maxima/5.42.2/share/contrib/integration/abs_integrate.mac"
(%i4) integrate(acos(cos(x)/(2*cos(x)+1)),x);

(%o4) -(2*%pi*atan(sin(x)/(cos(x)+1))-%pi*x)/2

See Maxima bug report 3566.

kcrisman commented 5 years ago

Upstream: Reported upstream. No feedback yet.