#1797 definite integral - bad answer


This bug has been reported at Sage forum http://groups.google.cz/group/sage-devel/t/c9c147c4dd4ea840

integrate(cos(x)^2 * (1 + sin(x)^2)^-3,x,0,%pi/2) returns a value wich is numericaly about 1.4.... but the correct value is close to 0.49

This is my bug_report(), but according to the discussion, the same problem is in the CVS version.

Maxima version: 5.18.1
Maxima build date: 13:5 5/3/2009
host type: i686-pc-linux-gnu
lisp-implementation-type: CMU Common Lisp
lisp-implementation-version: CVS 19d 19d-release (19D)


  • Andrej Vodopivec

    I think the error comes from same-sheet-subs (defint.lisp). The function substitutes the limits of integration into expressions which look like atan(something). It substitutes the limits plus it adds contributions from poles of something between the limits. The error is that it also adds contributions from poles on the limits, which I think should be ignored.

    The patch below fixes this behavior and with the patch applied we get

    (%i3) integrate(cos(x)^2 * (1 + sin(x)^2)^-3,x,0,%pi/2);
    (%o3) 7*%pi/2^(11/2)
    (%i4) %, numer;
    (%o4) 0.485940321361071

    Index: defint.lisp

    RCS file: /cvsroot/maxima/maxima/src/defint.lisp,v
    retrieving revision 1.66
    diff -u -r1.66 defint.lisp
    --- defint.lisp 6 Sep 2009 13:50:43 -0000 1.66
    +++ defint.lisp 18 Oct 2009 12:19:35 -0000
    @@ -938,7 +938,7 @@
    ((not (equal (sratsimp ipart) 0)) ())
    (t (let ((pole (poles-in-interval (let (($algebraic t))
    (sratsimp (cadr exp)))
    - var ll ul)))
    + var ll ul :on-boundary nil)))
    (cond ((and pole (not (or (eq pole '$unknown)
    (eq pole '$no))))
    (do ((l pole (cdr l)) (llist ()))
    @@ -3259,7 +3259,7 @@
    ;;; Returns $no $unknown or a list of poles in the interval (ll ul)
    ;;; for exp w.r.t. var.
    ;;; Form of list ((pole . multiplicity) (pole1 . multiplicity) ....)
    -(defun poles-in-interval (exp var ll ul)
    +(defun poles-in-interval (exp var ll ul &key (on-boundary t))
    (let* ((denom (cond ((mplusp exp)
    ($denom (sratsimp exp)))
    ((and (mexptp exp)
    @@ -3290,7 +3290,7 @@
    (t '$no)))
    (let* ((soltn (caar dummy))
    ;; (multiplicity (cdar dummy)) (not used? -- cwh)
    - (root-in-ll-ul (in-interval soltn ll ul)))
    + (root-in-ll-ul (in-interval soltn ll ul :on-boundary on-boundary)))
    (cond ((eq root-in-ll-ul '$no) '$no)
    ((eq root-in-ll-ul '$yes)
    (let ((lim-ans (is-a-pole exp soltn)))
    @@ -3365,7 +3365,7 @@
    'epsilon 0 '$plus))

    -(defun in-interval (place ll ul)
    +(defun in-interval (place ll ul &key (on-boundary t))
    ;; real values for ll and ul; place can be imaginary.
    (let ((order (ask-greateq ul ll)))
    (cond ((eq order '$yes))
    @@ -3374,8 +3374,8 @@
    (list '(mlist simp) ll ul)))))
    (if (not (equal ($imagpart place) 0))
    - (let ((lesseq-ul (ask-greateq ul place))
    - (greateq-ll (ask-greateq place ll)))
    + (let ((lesseq-ul (ask-greateq ul place :on-boundary on-boundary))
    + (greateq-ll (ask-greateq place ll :on-boundary on-boundary)))
    (if (and (eq lesseq-ul '$yes) (eq greateq-ll '$yes)) '$yes '$no))))

    (defun real-roots (exp var)
    @@ -3401,7 +3401,7 @@
    (cadr dummy))

    -(defun ask-greateq (x y)
    +(defun ask-greateq (x y &key (on-boundary t))
    ;;; Is x > y. X or Y can be $MINF or $INF, zeroA or zeroB.
    (let ((x (cond ((among 'zeroa x)
    (subst 0 'zeroa x))
    @@ -3432,8 +3432,10 @@
    ((eq y '$minf)
    (t (let ((ans ($asksign (m+ x (m- y)))))
    - (cond ((member ans '($zero $pos) :test #'eq)
    + (cond ((eq ans '$pos)
    + ((eq ans '$zero)
    + (if on-boundary '$yes '$no))
    ((eq ans '$neg)
    (t '$unknown)))))))

  • Dieter Kaiser

    Dieter Kaiser - 2009-11-12

    For the record:

    This bug is related to the bugs:

    ID: 2890315 - integrate(cot(x)^4,x,%pi/6,%pi/2)
    ID: 2880797 - bad answer in integrate(sqrt(sin(t)^2+cos(t)^2),t,0,2*%pi)

    Again we get wrong poles for an atan expression. The suggested patch in the thread to the bug report ID: 2890315 - integrate(cot(x)^4,x,%pi/6,%pi/2) does not help for this examle.

    Dieter Kaiser

  • Dan Gildea

    Dan Gildea - 2009-11-13
    • status: open --> closed
  • Dan Gildea

    Dan Gildea - 2009-11-13

    Fixed in defint.lisp rev 1.67 with a simpler version of the patch below.


Log in to post a comment.

Get latest updates about Open Source Projects, Conferences and News.

Sign up for the SourceForge newsletter:

No, thanks