Fix for Larry Bates's integration bug



Camm Maguire <camm@enhanced.com> writes:

> (C1) f(a):=integrate(sin(a)/(1-cos(a)^2*sin(t)^2), t,0,2*%Pi);
> 
>                                       SIN(a)
> (D1)           f(a) := INTEGRATE(-------------------, t, 0, 2 %PI)
>                                         2       2
>                                  1 - COS (a) SIN (t)
> (C2) f(%Pi/2);
> 
> (D2)                      2 %PI
> (C3) f(3.14159/2);

[...]

> 
> (D3)                        0
> 
> The correct answer is 2Pi, and in fact does not depend on the choice
> of the variable a if 0<a<Pi.
> 
> Larry Bates
> 

The origin of the bug is revealed by tracing the function PERIOD (in
defint), which tests if some expression is periodic with respect to a
given period and variable. I trace some more functions to tell the
whole story.

(C1) trace(?intsc1,?period,?sinint);
(D1) 			   [INTSC1, PERIOD, SININT]
(C2) f(a):=integrate(sin(a)/(1-cos(a)^2*sin(t)^2), t,0,2*%pi);
				      SIN(a)
(D2) 	      f(a) := INTEGRATE(-------------------, t, 0, 2 %PI)
				       2       2
				1 - COS (a) SIN (t)
(C3) f(1.5);

RAT replaced -.005003751699777 by -1107//221234 = -.005003751683738

RAT replaced -.005003751699777 by -1107//221234 = -.005003751683738
					 1
1 Enter ?INTSC1 [0, 2 %PI, -----------------------------]
						    2
			   1 - .0050037516997773 SIN (t)
				       1
 1 Enter ?PERIOD [2 %PI, -----------------------------, t]
						  2
			 1 - .0050037516997773 SIN (t)

RAT replaced -.005003751699777 by -1107//221234 = -.005003751683738
 1 Exit  ?PERIOD FALSE
1 Exit  ?INTSC1 FALSE

RAT replaced -.005003751699777 by -1107//221234 = -.005003751683738
			       1
1 Enter ?SININT [-----------------------------, t]
					  2
		 1 - .0050037516997773 SIN (t)

RAT replaced -.005003751699777 by -1107//221234 = -.005003751683738
			      220127 TAN(t)
		221234 ATAN(-----------------)
			    SQRT(48699576718)
1 Exit  ?SININT ------------------------------
		      SQRT(48699576718)
(D3) 				       0
(C4) untrace(?intsc1,?period,?sinint);
(D4) 			   [INTSC1, PERIOD, SININT]

The bug in PERIOD (which results from comparing expressions which have
been simplified in different ways) has the consequence that defint
thinks that INTSC1 can't be used, so SININT is tried to find an
indefinite integral; the resulting expression, however, is valid only
in (open) intervals where tan has no pole...

After applying the patch, I get the result

(D7) 		  2.235435669225693E-6 SQRT(800452634953) %PI
(C8) %/%pi,numer;
(D8) 			       1.999999999999074

BTW, I replaced your 3.14159/2 by 1.5 in order to get slightly less
ugly numbers and also since otherwise apparently some numbers become
too small (in absolute value) and that vile `LOG(0) has been
generated' error turns up somehow.

This seems to completely fix the case where a is a float.  Sadly, in
other cases we run into some unrelated bug.

(C9) f(3/2);
quotient is not exact

I haven't fixed this one.  Note, however, that other methods like
doing a partial fraction decomposition wrt sin(t) beforehand will
work.

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ cut ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Index: defint.lisp
===================================================================
RCS file: /cvsroot/maxima/maxima/src/defint.lisp,v
retrieving revision 1.3
diff -C2 -r1.3 defint.lisp
*** defint.lisp	11 Jul 2001 05:06:22 -0000	1.3
--- defint.lisp	8 Jul 2002 15:35:43 -0000
***************
*** 1492,1496 ****
  
  (DEFUN PERIOD (P E VAR)
!        (ALIKE1 E (NO-ERR-SUB (m+ P VAR) E))) 
  
  (DEFUN INFR (A)
--- 1492,1498 ----
  
  (DEFUN PERIOD (P E VAR)
!   (and (ALIKE1 (no-err-sub var E) (setq e (NO-ERR-SUB (m+ P VAR) E)))
!        ;; means there was no error
!        (not (eq e t))))
  
  (DEFUN INFR (A)
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ cut ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Wolfgang
-- 
wjenkner@inode.at