#650 integrate yields incorrect result on rational function

closed
nobody
5
2006-11-07
2004-11-25
Robert Dodier
No

"integrate" yields incorrect results on some rational
functions.

"Division by 0" is strange.

The definite integral below is certainly greater than 0
as the integrand is positive over [0, 1]. "integrate
(1/((x-3)^4+1/2), x)" returns the noun form, so maybe
(maybe) what happens is that the noun form is evaluated
at the limits of integration and it's the same, hence 0
is the result. (Just guessing there.)

Note that the difference between the two integrands is
that one is 1/(something + 1), while the other is
1/(same something + 1/2).

----------------------------------------
(%i1) integrate (1/((x-3)^4+1), x, 0, 1);

Division by 0
-- an error. Quitting. To debug this try
DEBUGMODE(TRUE);

(%i2) integrate (1/((x-3)^4+1/2), x, 0, 1);

(%o2) 0

(%i3) build_info ();
Maxima version: 5.9.1
Maxima build date: 21:24 9/23/2004
host type: i686-pc-linux-gnu
lisp-implementation-type: CMU Common Lisp
lisp-implementation-version: 19a

----------------------------------------
Same behavior observed in CVS build of 2004/11/24.

Discussion

  • Logged In: NO

    Two remarks:
    The method that Maxima uses to solve such integrals is
    integration by residues. Trace (residue) shows that this
    function is called with the correct arguments but answers
    zeroes for integrate (1/((x-3)^4+1/2), x, 0, 1);
    and runs into a division by zero for the other integral.

    The indefinite integrals can be solved with changevar:
    'Integrate(1/((x-3)^4+1/2), x);
    changevar (%, x - 3 - y ,y ,x);
    ev (%, Integrate);

     
  • Logged In: YES
    user_id=588346

    This is apparently another GCD problem.

    Fixed by setting the 'algebraic' flag:

    integrate (1/((x-3)^4+1), x, 0, 1),algebraic:true

    For the indefinite integral, you can factor over the
    Gaussians then use partfrac:

    integrate ( partfrac ( gfactor( 1/((x-3)^4+1) ), x ), x )

    You can simplify the result using ratsimp(rectform(...))

    -s

     
  • Robert Dodier
    Robert Dodier
    2006-04-10

    • labels: 460522 --> Lisp Core - Integration
     
  • Raymond Toy
    Raymond Toy
    2006-11-06

    Logged In: YES
    user_id=28849

    Recent changes in defint.lisp and residue.lisp has fixed the
    issue with the integrate(1/((x-3)^4+1),x,0,1).

    For integrate(1/((x-3)^4+1/2),x,0,1), it does seem to be
    either a gcd or residue problem.

    We can work around this issue by making 2 changes in
    residu.lisp. In the function RES, for the case of simple
    poles, we can use RES1 instead of $RESIDUE. In the function
    RES1, we need to call $RECTFORM for each pole to make is
    simpler.

    With these changes applied, this integral can be evaluated.
    However, the answer takes some time and is quite messy. We
    can get the simpler result by using
    factor(expand(sqrtdenest(<foo>))). The answer is:

    (2*atan((sqrt(2)-4*2^(1/4)+8)/(49*2^(3/4)+sqrt(2)-8))
    -log((2^(3/4)+12*sqrt(2)+73*2^(1/4)-2)/(33*2^(1/4)))
    +log((2^(3/4)-12*sqrt(2)+73*2^(1/4)+2)/(33*2^(1/4)))
    -2*atan((sqrt(2)+4*2^(1/4)+8)/(-49*2^(3/4)+sqrt(2)-8)))
    /(2*2^(3/4))

    Numerical evaluation of this compares favorably with the
    numerical result from quad_qags.

     
  • Raymond Toy
    Raymond Toy
    2006-11-07

    • status: open --> closed
     
  • Raymond Toy
    Raymond Toy
    2006-11-07

    Logged In: YES
    user_id=28849

    Fixed as suggested. This is really a workaround instead of
    the true fix.