#1906 Definite Integration of 1/(a-b*cos(x)) wrong

closed
nobody
5
2010-03-25
2010-02-17
Anonymous
No

Maxima 5.20.1 with wxMaxima.

integrate(1/(a-b*cos(x)),x,0,%pi); where a>0, 0<b<a yields 0.

Discussion

  • Barton Willis
    Barton Willis
    2010-02-26

    Notice how a float enters into the asksign:

    (%i4) integrate(1/(1-a*cos(x)),x);
    Is a^2-1.0 positive or negative?neg;
    (%o4) (2*atan(((2*a+2)*sin(x))/(2*sqrt(1-a^2)*(cos(x)+1))))/sqrt(1-a^2)

     
  • Raymond Toy
    Raymond Toy
    2010-03-23

    The incorrect result comes from polelist failing to identify the locations of the poles.

    Since the integrand is even, we can integrate from -%pi to %pi (or 0 to 2*%pi) and take half of the result. This integral is converted to the contour integral of -2/(b*yy^2-2*a*y+b) around the unit circle. This is evaluated by residues. We want to find the poles inside the unit circle and polelist is supposed to do that. The poles are correctly determined, but unfortunately for these poles, polelist cannot find the one pole that is in the circle. Therefore the function res thinks there are no poles in the unit circle and returns 0.

    When a and b are numbers, polelist does a better job and normally determines the pole that is within the unit circle.

    Perhaps the function that determines whether the pole is in the unit circle needs to be enhanced?

     
  • Dan Gildea
    Dan Gildea
    2010-03-24

    possible solution:

    (defun unitcir (grand var)
    (numden grand)
    (let ((result (princip (res nn* dn*
    #'(lambda (pt)
    (eq (let ((limitp nil))
    ($asksign (m+ -1 (cabs pt))))
    '$neg))
    #'(lambda (pt)
    (alike1 1 (cabs pt)))))))
    (cond (result (m* '$%pi result))
    (t nil))))

     
  • Raymond Toy
    Raymond Toy
    2010-03-25

    Yes, that works nicely. Need to make the second lambda also call asksign. And since the question is the same, it's nice to cache the answer from the first lambda.

     
  • Raymond Toy
    Raymond Toy
    2010-03-25

    dgildea's basic solution checked in.

    Closing bug.

     
  • Raymond Toy
    Raymond Toy
    2010-03-25

    • status: open --> closed