[Maxima-commits] CVS: maxima/archive/share/trash qq.dem, NONE, 1.1 qq.lisp, NONE, 1.2 qq.usg, NONE, 1.2

[Robert D. <rob...@us...>]

Update of /cvsroot/maxima/maxima/archive/share/trash
In directory sc8-pr-cvs7.sourceforge.net:/tmp/cvs-serv19120/archive/share/trash

Added Files:
	qq.dem qq.lisp qq.usg 
Log Message:
Cut out package qq (containing function quanc8) and expunge all mention
of the qq package and quanc8. Move qq files (qq.dem, qq.lisp, qq.usg) to archive/share/trash/. 
quanc8 is superseded by QUADPACK functions. Topic was discussed on the mailing list
circa 2007/01/02.

--- NEW FILE: qq.dem ---

/* Numerical Integration Demo for QUANC8*/

/* SETUP_AUTOLOAD(QQ,QUANC8); */

/* show the time and the garbage collection time, and use the 
   GC-demon to make better use of flonum space */

showtime:'all;
/* IF STATUS(FEATURE,ITS)=TRUE THEN LOAD("GCDEMN"); */

/* Sample highly oscillatory problem */

g(x):=(mode_declare(x,float),sin(1.0/x));

/* inefficient way to use QUANC8 */

quanc8(g(x),x,0.01,0.1);

/* give G as the thing to compile, then a semi-colon to read it in */
compile('g);

/* see how much faster now, using fast version of QUANC8 */

quanc8('g,0.01,0.1);

/* note that romberg doesn't easily manage oscillatory behaviors */

errcatch(romberg('g,0.01,0.1));

/* a function with two large and narrow peaks,
   so hard to integrate with a fixed step-size */

h(x):=(mode_declare(x,float),1.0/(.0001+(x-.3)^2)+1.0/(.0025+(x-.9)^2));
/* give H and a semi-colon in response */
compile('h);

/* the exact answer for integrate(h(x),x) is
   100.0*atan(100.0*x-30.0)+20.0*atan(20.0*x-18.0) */
/* the exact answer for integrate(h(x),x,0.0,1.0) is 361.847622.
   note that the romberg result is more accurate in this case */
romberg('h,0.0,1.0);
/* but at a large time cost compared to adaptive */
quanc8('h,0.0,1.0);
/* reduce QUANC8's relative error tolerance by a factor of 10 */
quanc8_relerr:1.0E-5;
/* does not take much longer, and gives much better result now */
quanc8('h,0.0,1.0);
/* put QUANC8_RELERR back to 1.0e-4 */
quanc8_relerr:1.0E-4$
/* the exact answer for integrate(h(x),x,0.0,10.0) is 372.33606.
   while too tough for fixed step size method */
errcatch(romberg('h,0.0,10.0));
/*QUANC8 is still super speedy even now */
quanc8('h,0.0,10.0);

/* double integral of sample function exp(x)*sin(y) 
   from x=0.0 to 2.0, y=0.0 to 1.0 */
/* you must be sure that the quanc8 source is loaded
   (load("qq");) before you  translate a call to quanc8 in a function, 
   else an error in macro expansion will occur when it can't figure 
   out if the 3 or 4 arg version is used */

/* the exact answer for 
   integrate(integrate(exp(x)*sin(y),x,0.0,2.0),y,0.0,1.0) is 2.93703434.
   integrate over y to get the final answer */
quanc8(quanc8(exp(x)*sin(y),x,0.0,2.0),y,0.0,1.0);

p():=quanc8(lambda([y],mode_declare(y,float),
                   quanc8(lambda([x],mode_declare(x,float),
                                 exp(x)*sin(y)),
                          0.0,2.0)),0.0,1.0);

p();
translate(p);
p();
/* give p and a semi-colon */
compile('p);
p();

--- NEW FILE: qq.lisp ---
;;;;;;;;;;;;;;;;;;; -*- Mode: Lisp; Package: Macsyma -*- ;;;;;;;;;;;;;;;;;;;
;;;     (c) Copyright 1982 Massachusetts Institute of Technology         ;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(in-package :maxima)

(macsyma-module qq)

;; (load-macsyma-macros Numerm)

;;; 10:55 pm Sept 25, 1981 - LPH & GJC
;;; added $numer and $float set to T so that quanc8(x,x,0,1/2); works. this
;;; is done inside the prog for $quanc8 so that they are put back when done.
;;; 3:30 pm Feb 12, 1982 - JPG, LPH, & GJC
;;; removed the $numer:$float:true and replaced with mto-float which is
;;; defined in maxsrc;numer > .
;;; 3:45 am April 11, 1982 - LPH
;;; fixed an error in estimating TOLERR due to neglect of sub-interval size
;;; when abserr check is made.

(defmvar $quanc8_flag 0.0
  "integer part-how many regions failed to converge. fractional part-
  how much of the interval was left when it failed")
(defmvar $quanc8_errest 0.0
  "an estimated error based on the difference between one- and two-panel
  calculation of an interval, summed over the whole region")
(defmvar $quanc8_abserr 1.0e-8
  "absolute error tolerance for reasonable single-precision answers")
(defmvar $quanc8_relerr 1.0e-4
  "relative error condition for reasonable run-times")	

;; (DEFUN ($QUANC8 TRANSLATED-MMACRO) (F A B &OPTIONAL (C NIL C-P))
;;   (IF (NOT C-P)
;;       `((CALL-QUANC8) ,F ,A ,B)
;;       `((CALL-QUANC8) ((LAMBDA) ((MLIST) ,A) ,F) ,B ,C)))

(defmspec $quanc8 (form)
  (if (cdddr (setq form (cdr form))) 
      (apply #'call-quanc8
	     (meval `((lambda) ((mlist) ,(cadr form)) ,(car form)))
	     (mapcar #'meval (cddr form)))
      (apply #'call-quanc8 (mapcar #'meval form))))

(def%tr $quanc8 (form)
  `($float
    ,@(cdr (tr-lisp-function-call
	    (if (cdddr (setq form (cdr form)))
		`((call-quanc8)
		  ((lambda) ((mlist) ,(cadr form)) ,(car form)) ,@(cddr form))
		`((call-quanc8) ,@form)) nil))))

(defvar quanc8-free-list ()
  "For efficient calls to quanc8 keep arrays we need on a free list.")

(defvar quanc8-|^]| ())
(defun quanc8-|^]| () (setq quanc8-|^]| t))

(defun call-quanc8 (fun a b)
  (bind-tramp1$
   fun fun
   (let ((vals (if quanc8-free-list
		   (pop quanc8-free-list)
		   (list (*array nil 'flonum 17.)
			 (*array nil 'flonum 17.)
			 (*array nil 'flonum 9. 31.)
			 (*array nil 'flonum 9. 31.)
			 (*array nil 'flonum 32.))))
	 (user-timesofar (cons #'quanc8-|^]| user-timesofar))) 
     (prog1
      (quanc8 fun
	      (mto-float a)
	      (mto-float b)
	      (car vals)
	      (car (cdr vals))
	      (car (cddr vals))
	      (car (cdddr vals))
	      (car (cddddr vals)))
      (push vals quanc8-free-list)))))

(defun quanc8 (fun a b x-arr f-arr xsave-arr fsave-arr qright-arr)
  (declare (type (simple-array cl:float)
		 x-arr f-arr xsave-arr fsave-arr qright-arr))
  ;; local macros for typing convenience.
  (macrolet ((x (j) `(arraycall flonum x-arr ,j))
	     (f (j) `(arraycall flonum f-arr ,j))
	     (xsave (j k) `(arraycall flonum xsave-arr ,j ,k))
	     (fsave (j k) `(arraycall flonum fsave-arr ,j ,k))
	     (qright (j) `(arraycall flonum qright-arr ,j))) 
    ;; Rudimentary (non-ansi GCL compatible) error handling.
    (let (errset)
      (or (car (errset
    (prog ((levmin 1.) (levmax 30.) (levout 6.) (nomax 5000.) (nofin 0)
	 (w0 (//$ 3956.0 14175.0)) (w1 (//$ 23552.0 14175.0))
	 (w2 (//$ -3712.0 14175.0)) (w3 (//$ 41984.0 14175.0))
	 (w4 (//$ -18160.0 14175.0))
	 (result 0.0) (cor11 0.0) (area 0.0)
	 (nofun 0.) (lev 0.) (nim 1.) (qprev 0.0)
	 (stone (//$ (-$ b a) 16.0))
	 (i 0.)
	 (step 0.0) (qleft 0.0) (qnow 0.0) (qdiff 0.0)
	 (esterr 0.0) (tolerr 0.0) (temp 0.0)
	 ($numer t) ($float t))
	
	(declare (cl:float w0 w1 w2 w3 w4 result cor11 area qprev stone step
			 qleft qnow qdiff esterr tolerr temp)
		 (fixnum i levmin levmax levout nomax nofin nofun lev nim)
		 (boolean $numer $float))

	(setq
	 $quanc8_flag 0.0
	 $quanc8_errest 0.0
	 nofin (- nomax (* 8 (+ levmax (* -1. levout) (expt 2. (1+ levout))))))
	      
	(cond ((= a b)
	       (return 0.0)))

	(setf (x 0.)
	      a)

	(setf (x 16.) 
	      b)

	(setf (x 8.)
	      (*$ 0.5 (+$ (x 0.)
			  (x 16.))))

	(setf (x 4.)
	      (*$ 0.5 (+$ (x 0.)
			  (x 8.))))

	(setf (x 12.)
	      (*$ 0.5 (+$ (x 16.)
			  (x 8.))))

	(setf (x 2.)
	      (*$ 0.5 (+$ (x 0.)
			  (x 4.))))

	(setf (x 6.)
	      (*$ 0.5 (+$ (x 4.)
			  (x 8.))))

	(setf (x 10.)
	      (*$ 0.5 (+$ (x 12.)
			  (x 8.))))

	(setf (x 14.)
	      (*$ 0.5 (+$ (x 12.)
			  (x 16.))))


  do-25
        (when quanc8-|^]|
	      (mtell "QUANC8 calculating at X= ~S" (x i))
	      (setq quanc8-|^]| nil))
	      
	(setf (f i)
	      (fcall$ fun (x i)))
	
	(setq i (+ 2. i))

	(if (> i 16.)
	    (go do-25-done))

	(go do-25)

do-25-done
	(setq nofun 9.)

 tag-30
	(setq i 1.)

  do-30
	(setf (x i)
	      (*$ 0.5 (+$ (x (1- i))
			  (x (1+ i)))))
        (when quanc8-|^]|
	      (mtell "QUANC8 calculating at X= ~S" (x i))
	      (setq quanc8-|^]| nil))

	(setf (f i)
	      (fcall$ fun (x i)))

	(setq i (+ 2. i))

	(if (> i 15.)
	    (go do-30-done))

	(go do-30)

do-30-done
	(setq nofun (+ 8. nofun))
	
	(setq step (//$ (-$ (x 16.)
			    (x 0.))
			16.0))

	(setq qleft (*$ step (+$ (*$ w0 (+$ (f 0.)
					    (f 8.)))
				 (*$ w1 (+$ (f 1.)
					    (f 7.)))
				 (*$ w2 (+$ (f 2.)
					    (f 6.)))
				 (*$ w3 (+$ (f 3.)
					    (f 5.)))
				 (*$ w4 (f 4.)))))

	(setf (qright (1+ lev))
	      (*$ step (+$ (*$ w0 (+$ (f 8.)
				      (f 16.)))
			   (*$ w1 (+$ (f 9.)
				      (f 15.)))
			   (*$ w2 (+$ (f 10.)
				      (f 14.)))
			   (*$ w3 (+$ (f 11.)
				      (f 13.)))
			   (*$ w4 (f 12.)))))
	
	(setq qnow (+$ qleft (qright (1+ lev))))

	(setq qdiff (-$ qnow qprev))

	(setq area (+$ area qdiff))

	(setq esterr (//$ (abs qdiff)
			  1023.0))

	(setq tolerr (*$
		      (//$ step stone)
		      (max $quanc8_abserr
			   (*$ (abs area)
			       $quanc8_relerr))))

	(if (< lev levmin)
	    (go tag-50))

	(if (or (> lev levmax)
		(= lev levmax))
	    (go tag-62))

	(if (> nofun nofin)
	    (go tag-60))

	(if (< esterr tolerr)
	    (go tag-70))

 tag-50
	(setq nim (* 2. nim)
	      lev (1+ lev))

	(setq i 1.)

  do-52
	(setf (fsave i lev)
	      (f (+ 8. i)))

	(setf (xsave i lev)
	      (x (+ 8. i)))

	(setq i (1+ i))

	(if (> i 8.)
	    (go do-52-done))

	(go do-52)

do-52-done
	(setq qprev qleft)

	(setq i 1.)

  do-55
	(setf (f (+ 18. (* -2. i)))
	      (f (+ 9. (* -1. i))))

	(setf (x (+ 18. (* -2. i)))
	      (x (+ 9. (* -1. i))))

	(setq i (1+ i))

	(if (> i 8.)
	    (go tag-30))

	(go do-55)

 tag-60
	(setq nofin (* 2. nofin))

	(setq levmax levout)

	(setq $quanc8_flag (+$ $quanc8_flag (//$ (-$ b (x 0.))
				 (-$ b a))))

	(go tag-70)

 tag-62
	(setq $quanc8_flag (+$ $quanc8_flag 1.0))

 tag-70
	(setq result (+$ result qnow)
	      $quanc8_errest (+$ $quanc8_errest esterr)
	      cor11 (+$ cor11 (//$ qdiff 1023.0)))

 tag-72
	(if (= nim (* 2. (// nim 2.)))
	    (go tag-75))

	(setq nim (// nim 2.))

	(setq lev (1- lev))

	(go tag-72)

 tag-75
	(setq nim (1+ nim))

	(if (or (< lev 0.)
		(= lev 0.))
	    (go tag-80))

	(setq qprev (qright lev))

	(setq i 1.)

	(setf (x 0.)
	      (x 16.))

	(setf (f 0.)
	      (f 16.))

  do-78
	(setf (f (* 2. i))
	      (fsave i lev))

	(setf (x (* 2. i))
	      (xsave i lev))
	
	(setq i (1+ i))

	(if (> i 8.)
	    (go tag-30))
	
	(go do-78)

 tag-80
	(setq result (+$ result cor11))

	(if (= $quanc8_errest 0.0)
	    (return result))
	
 tag-82
    (setq temp (+$ $quanc8_errest (abs result)))
	
	(if (not (= temp (abs result)))
	    (return result))

	(setq $quanc8_errest (*$ 2.0 $quanc8_errest))
	
	(go tag-82))))

	  ;; For whatever reason and with a suitable
	  ;; discrete meaning of convergence...
	  (merror "QUANC8 failed to converge.")))))

--- NEW FILE: qq.usg ---

The file SHARE1;QQ FASL contains a function QUANC8 which can take either 
3 or 4 arguments. The 3 arg version computes the integral of the function 
specified as the first argument over the interval from lo to hi as in
QUANC8('function name,lo,hi).  The function name should be quoted.
The 4 arg version will compute the integral of the function or expression
(first arg) with respect to the variable (second arg) over the interval 
from lo to hi as in QUANC8(<f(x) or expression in x>,x,lo,hi).

The method used is the Newton-Cotes 8th order polynomial quadrature, and 
the routine is adaptive. It will thus spend time dividing the interval only 
when necessary to achieve the error conditions specified by the global 
variables QUANC8_RELERR (default value=1.0e-4) and QUANC8_ABSERR
(default value=1.0e-8) which give
the relative error test 
|integral(function)-computed value|< quanc8_relerr*|integral(function)|
and the absolute error test 
|integral(function)-computed value|<quanc8_abserr. 

The error from each subinterval is estimated and the contribution from a 
subinterval is accepted only when the integral over the subinterval satisfies
the error test over the subinterval. The total estimated error of the integral
is contained in the global variable QUANC8_ERREST (default value=0.0). 

The global variable QUANC8_FLAG (default value=0.0) will contain valuable 
information if the computation fails to satisfy the error conditions. The 
integer part will tell you how many subintervals failed to converge and the 
fractional part will tell you where the singular behavior is, as follows: 
singular point=lo+(1.-frac part)*(hi-lo).
Thus quanc8(tan(x),x,1.57,1.6); gives frac=.97 so trouble is at 1.57+.03*.03=
1.5709 (=half pi). If QUANC8_FLAG is not 0.0, you should be cautious in using 
the return value, and should try ROMBERG or a Simpson method and see if the
result checks. Analysis of possible singular behavior might be advisable. You
may get QUANC8_FLAG=<integer>.0 and an error message (such as division by 0)
when a singular point is hit in the interval. You will have to find the
singularity and eliminate it before QUANC8 will get an answer. Functions 
which have very large derivatives may throw the error estimate way off and
cause the wrong points to be used, and a wrong answer returned. Try
romberg(exp(-.002*x^2)*cos(x)^2,x,0.,100.); with the default tolerance, and
quanc8(exp(-.002*x^2)*cos(x)^2,x,0.,100.); with quanc8_relerr=1.e-7 and 1.e-8.
The last result is consistent with romberg while the previous one is off by
a factor of 2! This is due to the bad behavior of the derivatives near
x=10.0 which cause the adaptive routine to have trouble. If you use 
quanc8('f,a,c)+quanc8('f,c,b) where a<c<b, you will do better in such cases.

Typing a control-right-bracket (^]) will give a printout of where the 
computation is at the moment, and how much time has been used.

You can do batch("qq demo"); for some comparisons with the ROMBERG numerical
integrator (which is not adaptive). Note that ROMBERG usually gives more 
accurate answers for comparable tolerances, while QUANC8 will get the same 
answer faster even with a smaller tolerance, because ROMBERG subdivides the
whole interval if the total result is not within error tolerance, while QUANC8
improves only where needed, thus saving many function calls. ROMBERG will also
fail to converge when oscillatory behavior is overwhelming, while QUANC8 will 
adapt in the regions as it sees fit. (The global variable ROMBERGMIN is 
designed to allow you a minimum number of function calls in such cases, so 
that exp(-x)*sin(12*x) can be integrated from 0 to 4*%pi without erroneously 
giving 0. from the first few function calls.)

To make your macsyma user function callable in the fastest way, you must use
MODE_DECLARE and then translate and compile the function. Read TRANSL;TRANSL
NEWS for info on this, or ask me for more info. The speed of the computation
may be increased by well over an order of magnitude when compilation is used.
If you do multiple integrals, it is really necessary to compile the function
in order to avoid the time spent on function calls. A sample use of QUANC8
for a double integral is in the batch("qq demo"); and compilation is nearly 
a hundred times faster in doing the work!

Comments, bugs, etc. should be sent to LPH@MIT-MC.