quadpack functions wrapper for numerical integration
Subject: quadpack functions wrapper for numerical integration
From: Edwin Woollett
Date: Wed, 5 Oct 2011 13:09:12 -0700
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