Raymond Toy <rtoy@earthlink.net> writes:
> I'm open to suggestions on what arguments these functions should take
> and what they should return.
>
> Also, since these are in Lisp, it opens up the possibility of doing
> multiple integrals adaptively
Below you find a quick prototype hack which implements a QUANC8 (from
share/numeric/qq.lisp) or ROMBERG like interface for the new QAGS
function. Here's an iterated integral from qq.dem (the other examples
in that demo also work)
(C2) qags(qags(exp(x)*sin(y),x,0,2),y,0,1);
(D2) 2.937034356357518
(C3) qags_info;
(D3) [2.937034356357518, 3.26076316654324E-14, 21, 0]
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ cut ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
(in-package "MAXIMA")
(defmvar $qags_epsrel 1d-8)
(defmvar $qags_limit 200)
(defmvar $qags_info)
(defmspec $qags (form)
(if (cdddr (setq form (cdr form)))
(apply #'call-qags
(meval `((lambda) ((mlist) ,(cadr form)) ,(car form)))
(mapcar #'meval (cddr form)))
(apply #'call-qags (mapcar #'meval form))))
(def%tr $qags (form)
`($float
,@(cdr (tr-lisp-function-call
(if (cdddr (setq form (cdr form)))
`((call-qags)
((lambda) ((mlist) ,(cadr form)) ,(car form)) ,@(cddr form))
`((call-qags) ,@form)) nil))))
(defun call-qags (fun a b)
(bind-tramp1$
fun fun
(mqags fun (mto-float a) (mto-float b))))
;; Adapted from quadpack.lisp
(defun mqags (fun a b)
(let* ((lenw (* 4 $qags_limit))
(work (make-array lenw :element-type 'double-float))
(iwork (make-array $qags_limit :element-type 'f2cl-lib:integer4)))
(multiple-value-bind (junk z-a z-b z-epsabs z-epsrel result abserr neval ier z-limit z-lenw last)
(slatec::dqags fun a b 0d0 (float $qags_epsrel)
0d0 0d0 0 0 $qags_limit lenw 0 iwork work)
(declare (ignore junk z-a z-b z-epsabs z-epsrel z-limit z-lenw))
(setq $qags_info (list '(mlist) result abserr neval ier))
result)))
;;; mqags.lisp ends here
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ cut ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Wolfgang