Raymond Toy <rtoy@earthlink.net> writes:
> I wonder about that estimated error.
Yes I know, the $QAGS_INFOs from the inner integral are overwritten...
> Still looking for the appropriate interface. I don't like globals
> that much, but perhaps that's appropriate for epsrel, limit, and
> other parameters. I prefer functional interfaces.
Actually, mutatis mutandis I just copied the defmspec and def%tr from
what I checked in for qq.lisp a few days ago. Strangely, I realized
only after that commit that I could have adapted the ROMBERG interface
instead, which does a similar thing in a different fashion. This
seems to be mostly a matter of taste. However, the ROMBERG interface
also binds some variables (for optimizing numerical frobbing) which I
didn't care about since the original (broken) $QUANC8 didn't care
about them. I still have to figure out if they are really necessary
or if they are implied by the bind-tramp stuff.
I think it would be best to write an interface common to the old
(ROMBERG, QUANC8) and new numerical integration functions. For both
DEFMSPEC and DEF%TR we have the complete calling form at our disposal,
so we could write something like
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ cut ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
(defun maxima-numerical-integration (form)
(let ((integrator (get (car (pop form)) 'numerical-integrator)))
(if (cdddr form)
(apply #'call-numerical-integrator integrator
(meval `((lambda) ((mlist) ,(cadr form)) ,(car form)))
(mapcar #'meval (cddr form)))
(apply #'call-numerical-integrator integrator (mapcar #'meval form)))))
(defun call-numerical-integrator (integrator fun a b)
(bind-tramp1$
fun fun
(funcall integrator fun (mto-float a) (mto-float b))))
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ cut ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
and the definition for the user-level function would be just
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ cut ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
(setf (get '$qags 'mfexpr*) #'maxima-numerical-integration)
(setf (get '$qags 'numerical-integrator) 'mqags)
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ cut ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
where MQAGS is the function adapted from quadpack.lisp (see my
previous message). We could also passe optional arguments to the
low-level integrator.
> We should also do something other than force the user to look at the
> error code to figure out if the integration was ok. Conditions
> would be ideal, but it seems gcl doesn't have conditions and/or
> handler-case/bind, etc.
The low-level integrator could signal conditions and
CALL-NUMERICAL-INTEGRATOR could handle them.
Wolfgang