#2520 abs_integrate fails on abs(sin(x)) and abs(cos(x))

None
open
nobody
5
2014-12-15
2012-12-17
No

Summary: integrate(abs(sin(x)),x) and integrate(abs(cos(x)),x) are computed incorrectly. integrate(abs(sin(x)),x,0,2*%pi) causes crashes, and integrate(abs(cos(x)),x,0,2*%pi) returns nonsense.

In more detail: Maxima 5.28.0 (SBCL/OSX 10.6.8 - downloaded from SF version), Maxima 5.26.0 (ECL 11.1 / Sage), Maxima 5.29.1 (ECL 12.12 and CLISP) all compute integrate(abs(sin(x)),x) and integrate(abs(cos(x)),x) incorrectly.

On Maxima 5.26.0 integrate(abs(sin(x)),x,0,2*%pi) returns the unevaluated integral, on Maxima 5.29.1 it causes segfaults, on Maxima 5.28.0 one gets "Control stack exhausted (no more space for function call frames)."
On Maxima 5.28.0 integrate(abs(cos(x)),x,0,2*%pi) returns 0, after a long wait and many messages "log: encountered log(0)."

See more details on http://trac.sagemath.org/sage_trac/ticket/13364 (this is where we stumbled upon this problem).

Related

Bugs: #2520

Discussion

  • Rupert Swarbrick

    Confusingly, on my machine, both the current git version and 2.28.0 give noun forms:

    rupert@skate /s/n/maxima> ./maxima-local
    Maxima 5.28.0_130_g925e3f2_dirty http://maxima.sourceforge.net
    using Lisp SBCL 1.0.57.0.debian
    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) integrate(abs(sin(x)),x,0,2*%pi);
                                  2 %pi
                                 /
                                 [
    (%o1)                        I      abs(sin(x)) dx
                                 ]
                                 /
                                  0
    (%i2) integrate(abs(cos(x)),x,0,2*%pi);
                                  2 %pi
                                 /
                                 [
    (%o2)                        I      abs(cos(x)) dx
                                 ]
                                 /
                                  0
    (%i3) rupert@skate /s/n/maxima> maxima
    
    Maxima 5.28.0 http://maxima.sourceforge.net
    using Lisp GNU Common Lisp (GCL) GCL 2.6.7 (a.k.a. GCL)
    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) integrate(abs(cos(x)),x,0,2*%pi);
                                  2 %pi
                                 /
                                 [
    (%o1)                        I      abs(cos(x)) dx
                                 ]
                                 /
                                  0
    (%i2) integrate(abs(sin(x)),x,0,2*%pi);
                                  2 %pi
                                 /
                                 [
    (%o2)                        I      abs(sin(x)) dx
                                 ]
                                 /
                                  0
    
     
  • Dima Pasechnik

    Dima Pasechnik - 2012-12-17

    one needs to explicitly
    load(abs_integrate);

     
  • Rupert Swarbrick

    Some more poking shows that, with abs(sin(x)), everything dies because Maxima is trying to find the limit of some truly hideous expression, returned by

    (ANTIDERIV '((MABS SIMP) ((%SIN SIMP) $X)))
    

    namely:

    ~~~~
    4((signum(1/(cos(x)+1))signum(sin(x))log(abs(sin(x)/(cos(x)+1)+%i))
    +signum(1/(cos(x)+1))
    signum(sin(x))log(abs(sin(x)/(cos(x)+1)-%i))
    -signum(1/(cos(x)+1))
    signum(sin(x))log(2sin(x)^2/(cos(x)+1)^2+2)
    +signum(1/(cos(x)+1))signum(sin(x))

    (log(2)-2sin(x)atan(sin(x)/(cos(x)+1))/(cos(x)+1)))
    /4
    +abs(sin(x))(atan(sin(x)/(cos(x)+1))/2
    +sin(x)/((cos(x)+1)
    (2*sin(x)^2/(cos(x)+1)^2+2)))
    /abs(cos(x)+1))
    ~~~~~

    It looks plausibly correct when you check it between, say, 0.1 and
    0.5, but then when you plot it you see that it's not even continuous!
    (I think that at every point it is continuous, the function is smooth
    and its derivative is probably equal to abs(sin(x)), but that's not
    particularly helpful!)

    A picture: [-img src=antideriv.png alt=graph of antiderivative: missing =-]

     
    Attachments
  • Dima Pasechnik

    Dima Pasechnik - 2012-12-17

    It is slightly better with abs(cos(x))
    (%i7) integrate(abs(cos(x)),x);

    log: encountered log(0).
    [..]
    (%o7) 2((1/2-sin(x)/((cos(x)+1)(sin(x)^2/(cos(x)+1)^2+1)))
    signum(sin(x)/(cos(x)+1)-1)
    +1)
    signum(sin(x)/(cos(x)+1)+1)

    Plotting this function shows a jump discontinuity at %pi, and maxima computes the limit at %pi wrongly (it claims both left and right limits to be 0), while the left must be 3, and the right must be -1, according to the plot. So it seems to be a bug in the limit computation.

     
  • Rupert Swarbrick

    In the abs(sin(x)) case, I hunted a bit further. Basically, Maxima is failing to compute some slightly unpleasant looking limit, and seemingly hanging somewhere in the gruntz1 routine. I managed to trim the testcase a little and

    u: log(abs(%i*cos(x)+%i)/abs(cos(x)+1))$
    v: abs(sin(x))*atan(sin(x)/(cos(x)+1))$
    
    limit([u,v], x, 2*%pi);
    => [0, 0]
    
    limit(u+v, x, 2*%pi);
    => Hangs...
    
     
  • Rupert Swarbrick

    And a bit further... We're dying in a call to $taylor. A related call that blows the stack is

    taylor (abs(sin(x))*atan(sin(x)/(cos(x)+1)), x, 0, 0);
    

    Removing abs(sin(x)) stops things going haywire.

     
  • Barton Willis

    Barton Willis - 2012-12-17

    Using GCL and Clozure CL, I'm unable to reproduce segfaults or crashes. Maybe these errors have something to do with ECL.

    In neighborhoods of zero, integrate(abs(sin(x)),x) and integrate(abs(cos(x)),x) are both correct. Unfortunately these antiderivaties are not valid on the reals. Maxima incorrectly assumes that these antiderivaties are valid on larger sets, and thus it incorrectly computes some definite integrals.

    Most all of this happens outside of the abs_integrate code. Maxima makes some fairly dubious substitutions, and the abs_integrate code (as far as I can tell) correctly integrates an expression that is mostly rational, but involves some absolute values.

     
  • Rupert Swarbrick

    Barton: That's surprising. I only have git master compiled with sbcl, but on ECL with 5.28.0, I get:

    rupert@skate ~> maxima
    Maxima 5.28.0 http://maxima.sourceforge.net
    using Lisp GNU Common Lisp (GCL) GCL 2.6.7 (a.k.a. GCL)
    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) load("abs_integrate");
    (%o1) /usr/share/maxima/5.28.0/share/contrib/integration/abs_integrate.mac
    (%i2) integrate (abs(sin(x)), x, 0, 2*%pi);
    
    Maxima encountered a Lisp error:
    
     Error in MACSYMA-TOP-LEVEL [or a callee]: Bind stack overflow.
    
    Automatically continuing.
    To enable the Lisp debugger set *debugger-hook* to nil.
    (%i3) 
    

    which seems to be what the bug reporter described. As far as I can tell, it's to do with Taylor going off the deep end on the expression in my previous note. I plan to trace and debug it further this evening, but I'd be very surprised if the behaviour varied between lisps.

     
  • Robert Dodier

    Robert Dodier - 2012-12-17
    • Description has changed:

    Diff:

    --- old
    +++ new
    @@ -1,10 +1,8 @@
    -Summary: integrate(abs(sin(x)),x) and integrate(abs(cos(x)),x) are computed incorrectly. integrate(abs(sin(x)),x,0,2*%pi) causes crashes, and integrate(abs(cos(x)),x,0,2*%pi) returns nonsense.
    +Summary: integrate(abs(sin(x)),x) and integrate(abs(cos(x)),x) are computed incorrectly. integrate(abs(sin(x)),x,0,2\*%pi) causes crashes, and integrate(abs(cos(x)),x,0,2\*%pi) returns nonsense.
    
     In more detail: Maxima 5.28.0 (SBCL/OSX 10.6.8 - downloaded from SF version), Maxima 5.26.0 (ECL 11.1 / Sage), Maxima 5.29.1 (ECL 12.12 and CLISP) all compute integrate(abs(sin(x)),x) and integrate(abs(cos(x)),x) incorrectly.
    
    -On Maxima 5.26.0 integrate(abs(sin(x)),x,0,2*%pi) returns the unevaluated integral, on Maxima 5.29.1 it causes segfaults, on Maxima 5.28.0 one gets "Control stack exhausted (no more space for function call frames)."
    -On Maxima 5.28.0 integrate(abs(cos(x)),x,0,2*%pi) returns 0, after a long wait and many messages "log: encountered log(0)."
    +On Maxima 5.26.0 integrate(abs(sin(x)),x,0,2\*%pi) returns the unevaluated integral, on Maxima 5.29.1 it causes segfaults, on Maxima 5.28.0 one gets "Control stack exhausted (no more space for function call frames)."
    +On Maxima 5.28.0 integrate(abs(cos(x)),x,0,2\*%pi) returns 0, after a long wait and many messages "log: encountered log(0)."
    
     See more details on http://trac.sagemath.org/sage_trac/ticket/13364 (this is where we stumbled upon this problem).
    -
    -
    
     
  • Robert Dodier

    Robert Dodier - 2012-12-17

    Fix up formatting: throw in backslashes to make literal asterisks.

     
  • Barton Willis

    Barton Willis - 2012-12-17

    $ ./maxima-local
    Maxima 5.28post http://maxima.sourceforge.net
    using Lisp Clozure Common Lisp Version 1.8-r15286M (WindowsX8664)
    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) load(abs_integrate);
    (%o1) c:/maxima-clean/share/contrib/integration/abs_integrate.mac
    (%i2) integrate (abs(sin(x)), x, 0, 2*%pi);

    --> hostile abort. maybe I didn't try this case with Clozure. With GCL, the result is a noun form.

    --Barton


    From: Rupert Swarbrick [rswarbrick@users.sf.net]
    Sent: Monday, December 17, 2012 15:16
    To: [maxima:bugs]
    Subject: [maxima:bugs] #2520 abs_integrate fails on abs(sin(x)) and abs(cos(x))

    Barton: That's surprising. I only have git master compiled with sbcl, but on ECL with 5.28.0, I get:

    rupert@skate ~> maxima
    Maxima 5.28.0 http://maxima.sourceforge.net
    using Lisp GNU Common Lisp (GCL) GCL 2.6.7 (a.k.a. GCL)
    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) load("abs_integrate");
    (%o1) /usr/share/maxima/5.28.0/share/contrib/integration/abs_integrate.mac
    (%i2) integrate (abs(sin(x)), x, 0, 2*%pi);

    Maxima encountered a Lisp error:

    Error in MACSYMA-TOP-LEVEL [or a callee]: Bind stack overflow.

    Automatically continuing.
    To enable the Lisp debugger set debugger-hook to nil.
    (%i3)

    which seems to be what the bug reporter described. As far as I can tell, it's to do with Taylor going off the deep end on the expression in my previous note. I plan to trace and debug it further this evening, but I'd be very surprised if the behaviour varied between lisps.


    [bugs:#2520]http://sourceforge.net/p/maxima/bugs/2520/ abs_integrate fails on abs(sin(x)) and abs(cos(x))

    Status: open
    Labels: abs_integrate
    Created: Mon Dec 17, 2012 10:22 AM UTC by Dima Pasechnik
    Last Updated: Mon Dec 17, 2012 08:50 PM UTC
    Owner: nobody

    Summary: integrate(abs(sin(x)),x) and integrate(abs(cos(x)),x) are computed incorrectly. integrate(abs(sin(x)),x,0,2%pi) causes crashes, and integrate(abs(cos(x)),x,0,2%pi) returns nonsense.

    In more detail: Maxima 5.28.0 (SBCL/OSX 10.6.8 - downloaded from SF version), Maxima 5.26.0 (ECL 11.1 / Sage), Maxima 5.29.1 (ECL 12.12 and CLISP) all compute integrate(abs(sin(x)),x) and integrate(abs(cos(x)),x) incorrectly.

    On Maxima 5.26.0 integrate(abs(sin(x)),x,0,2%pi) returns the unevaluated integral, on Maxima 5.29.1 it causes segfaults, on Maxima 5.28.0 one gets "Control stack exhausted (no more space for function call frames)."
    On Maxima 5.28.0 integrate(abs(cos(x)),x,0,2%pi) returns 0, after a long wait and many messages "log: encountered log(0)."

    See more details on http://trac.sagemath.org/sage_trac/ticket/13364 (this is where we stumbled upon this problem).


    Sent from sourceforge.net because you indicated interest in https://sourceforge.net/p/maxima/bugs/2520/

    To unsubscribe from further messages, please visit https://sourceforge.net/auth/prefs/

     

    Related

    Bugs: #2520

  • Rupert Swarbrick

    Assuming that the Taylor example I found is the root of this error, I carried on hunting. A simpler example again that causes the same error is

    taylor (cos(x)*sin(x)/abs(sin(x)), x, 0, 0);
    

    If you trace taylor2 beforehand, you'll see an infinite recursion. What causes it is that you get calls of the form

            5: (TAYLOR2
                ((MTIMES) ((RAT SIMP) 1 1)
                 ((MTIMES SIMP) ((%COS SIMP) $X) ((%SIN SIMP) $X)
                  ((MEXPT SIMP) ((MABS SIMP) ((%SIN SIMP) $X)) -1))
                 ((MEXPT) $X 1)))
    

    to evaluate this, you end up calling

                      10: (TAYLOR2 ((MABS SIMP) ((%SIN SIMP) $X)))
    

    and, unfortunately, this falls through to the bottom of the cond in taylor2, which looks like

    (t (let ((exact-poly () ))  ; Taylor series aren't exact
          (taylor2 (diff-expand e tlist))))
    

    DIFF-EXPAND, when called on this abs form, returns

                             ((MPLUS)
                              ((MTIMES) ((RAT SIMP) 1 1)
                               ((MTIMES SIMP) ((%COS SIMP) $X) ((%SIN SIMP) $X)
                                ((MEXPT SIMP) ((MABS SIMP) ((%SIN SIMP) $X))
                                 -1))
                               ((MEXPT) $X 1))
                              0)
    

    and we recurse. Bleugh.

     
  • Rupert Swarbrick

    Continuing the detective work, write f(x) for abs(sin(x)). The unhelpful result from DIFF-EXPAND is because it is trying to evaluate "f(0) + f'(0) + ..." (up to order one in this case) and calls EVAL-DERIV to actually evaluate the derivative it gets. This is of the form

    at(cos(x)*sin(x)/abs(sin(x)),x = 0)
    

    and EVAL-DERIV returns cos(x)*sin(x)/abs(sin(x)) rather than throwing an error because of the division by zero. As such, the result from DIFF-EXPAND actually has f'(x) rather than f'(0).

     
  • Rupert Swarbrick

    The attached patch fixes the stack overflow bug. The test for %derivative seems a bit hacky, but I'm not sure how to do better. Ideas, anyone?

    The only unexpected failure in the resulting testsuite run is a numerical error in rtest_gamma, but I'm pretty certain that's not my fault :-)

    With this patch, the sin integral returns a noun form. Not ideal, but better than crashing Maxima! The cos integral still returns the wrong answer, but that was tickling a different bug.

    Comments? (Which I'll read tomorrow. It's midnight and I've still got to cycle home!)

     
    • Dima Pasechnik

      Dima Pasechnik - 2012-12-18

      deleted a meaningless comment.

       
      Last edit: Dima Pasechnik 2012-12-18
  • Robert Dodier

    Robert Dodier - 2012-12-18

    Rupert, thanks for working on this problem. Looks like your analysis is right on the mark. Please keep up the good work.

    There has to be a way to express the fact that there is no Taylor expansion for some functions (e.g. absolute value). The most useful behavior for EVAL-DERIV is to throw something to indicate the derivative doesn't exist at the specified point. Returning a noun expression suggests that the derivative exists but EVAL-DERIV wasn't able to calculate it.

    Likewise the Maxima function 'taylor' would be more useful if it could indicate that the Taylor expansion doesn't exist. As it stands, there is no documented way for 'taylor' to express that. Maybe that is the first thing to work on: what should be the behavior of 'taylor' when there is no such Taylor expansion?

     
    • Rupert Swarbrick

      I'm not sure what you mean here. At the moment, $taylor throws an error (see tay-err), which is caught by $limit (yielding a noun form for the limit). Since there's not really a way to throw and catch errors in the maxima language, what else could happen?

       

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

Sign up for the SourceForge newsletter:





No, thanks