quadpack functions wrapper for numerical integration



I am beginning revisions of ch. 8, numerical integration,
(maxima by example), and want to include a simple
wrapper function for the quadpack functions for
one dimensional quadrature.

(area and volume integral wrappers are a separate project.)

If we ignore efficiency and principal value integrals,
we can get away with just using quad_qags (finite
domain) and quad_qagi.

A global variable 'details' (set to false inside nint.mac)
governs how much output nint produces.

---------------------------
(%i1) load(nint);
(%o1)                          c:/work2/nint.mac
(%i2) fpprintprec:8$
(%i3) nint (a*x,x,0,1);
                       undefined parameter in integrand

(%o3)                                false
(%i4) quad_qags(sqrt(x)*log(1/x),x,0,1);
(%o4)                 [0.444444, 4.93432455E-16, 315, 0]
(%i5) nint(sqrt(x)*log(1/x),x,0,1);
(%o5)                              0.444444
(%i6) details;
(%o6)                                false
(%i7) details:true$
(%i8) nint(sqrt(x)*log(1/x),x,0,1);
 case quad_qags
 est. abs. error =  4.93432455E-16
 no. of evaluations =  315
(%o8)                              0.444444
(%i9) details:false$
(%i10) quad_qagi(exp(-x^2),x,minf,inf);
(%o10)                [1.7724539, 1.42026368E-8, 270, 0]
(%i11) nint(exp(-x^2),x,minf,inf);
(%o11)                             1.7724539
(%i12) quad_qags(cos(50*x)*sin(3*x)*exp(-x),x,0,1);
(%o12) [-0.00191455,2.8536E-15,315,0]
(%i13) nint(cos(50*x)*sin(3*x)*exp(-x),x,0,1);
(%o13) -0.00191455
--------------------------------------------------------
The code so far is:

-------------------------------------
/* nint.mac
   nintegrate code */

 nint (expr,_var%,_a%, _b%, [_v%]) :=
   block ([argL,rL,nerr ],


     if not numberp (float (subst (_var%=1.0,expr))) then
       ( disp ("undefined parameter in integrand"),
         return() ),

     argL : flatten ( cons ([expr,_var%,_a%,_b%],[_v%])),


     if _a% = minf or _b% = inf then (
       if details then print (" case quad_qagi"),
         rL : apply ('quad_qagi, argL) )
     else  (
       if details then print (" case quad_qags "),
      rL : apply ('quad_qags,argL)),


      if not listp (rL) then (
       disp ("undefined parameter produces a noun form with default 
options"),
       return (rL)),


     nerr : part (rL,4),
     if nerr = 1 then disp ("too many subintervals done")
     else if nerr = 2 then disp ("excessive roundoff error detected")
     else if nerr = 3 then disp ("extremely bad integrand behavior 
observed")
     else if nerr = 4 then disp ("quadrature failed to converge")
     else if nerr = 5 then disp ("integral is probably divergent or slowly 
convergent")
     else if nerr = 6 then disp ("the input is invalid"),

     if details then (
        print (" est. abs. error = ",part (rL,2)),
        print (" no. of evaluations = ",part (rL,3))),
     part (rL,1))$


     details : false$
    -------------------------------------------

Ted Woollett