A minimal version of nintegrate, called nint here,
uses defaults epsrel=1d-4,limit=600, key=3, and
uses either quad_qagi (infinite domain),
quad_qag (finite domain oscillation),
quad_qags (finite domain singular),
with the choice between the latter two functions
made on the basis of which method takes the
least number of function evaluations (as
suggested by Raymond Toy).
The slatec error messages are surpressed using
Ray's lisp code replacement j4save.lisp.
Two global lists are created for after the fact
diagnosis:
1. gargL includes the method used and
the list of input args.
2. goutL includes the method used and
the complete output of quadpack.
In addition, if the global flag details is set
to true, nint will print out the method used.
nint.mac sets details:false.
If the numerical integral cannot be found using
the defaults, appropriate error messages are
written to the screen (rather than error code
integers).
The two global lists gargL and goutL allow the
user to see the defaults being used, and this
can help the user go directly to a quadpack
method and play with the options.
Here are examples of an oscillating integrand,
an integrand with singular features,
an integral over an infinite domain,
and an example which the defaults cannot
cope with.
=============================
(%i1) load(nint);
(%o1) "c:/work2/nint.mac"
(%i2) details:true$
oscillating integrand
(%i3) g : (x-2)^2 * sin (4000*x)$
(%i4) tval : float (integrate (g,x,2,3));
(%o4) -1.5862464665268553E-4
(%i5) nint(g,x,2,3);
quad_qag
(%o5) -1.5862464665235301E-4
(%i6) relerr(%);
(%o6) 2.0962997448270408E-12
(%i7) gargL;
(%o7) [qag,(x-2)^2*sin(4000*x),x,2,3,3,epsrel = 1.0E-4,limit = 600]
(%i8) goutL;
(%o8) [qag,-1.5862464665235301E-4,1.5316233778733241E-8,11873,0]
singular integrand
(%i9) g:x^(1/2)*log(1/x);
(%o9) -sqrt(x)*log(x)
(%i10) tval : float (integrate(g,x,0,1));
(%o10) 0.44444444444444
(%i11) nint(g,x,0,1);
quad_qags
(%o11) 0.44444444444444
(%i12) relerr(%);
(%o12) 3.7470027081099033E-16
infinite domain
(%i13) g : x^2*exp (-4*x)$
(%i14) tval : float (integrate (g,x,0,inf));
(%o14) 0.03125
(%i15) nint (g,x,0,inf);
quad_qagi
(%o15) 0.031250000000007
(%i16) relerr(%);
(%o16) 2.1760371282653068E-13
fail behavior:
(%i17) g: cos (log (x)/x)/x;
(%o17) cos(log(x)/x)/x
(%i18) integrate (g,x,0,1);
(%o18) 'integrate(cos(log(x)/x)/x,x,0,1)
(%i19) nint(g,x,0,1);
keyval = 3 prec = 1.0E-4 limval = 600
quad_qag error code = too many subintervals done
quad_qags error code = integral is probably divergent or slowly convergent
(%o19) false
=================================
code file: (using print_file from Ch.2 mfiles.mac
which produces less space between lines when
copying Xmaxima output to a text file, as compared
with printfile)
(%i21) load(mfiles);
(%o21) "c:/work2/mfiles.mac"
(%i22) print_file("nint.mac")$
-------------------------------------------------------
/* nint.mac
oct. 12, 2011
nintegrate code */
load ("j4save.lisp")$
codelist :
[[1,"too many subintervals done"],
[2,"excessive roundoff error detected"],
[3,"extremely bad integrand behavior observed"],
[4, "integration failed to converge"],
[5, "integral is probably divergent or slowly convergent"],
[6, "the input is invalid"]]$
/* tval should hold "true" value */
relerr(u) := abs(u - tval)/abs(tval)$
reldiff(u,v) := 2.0*abs(u - v)/abs(u + v)$
/* oct. 12, 2011 nint is minimal version which sets
limit = 600, epsrel = 1d-4, chooses either
quad_qag or quad_qags for finite interval
based on min number of function evaluations
and calls quad_qagi for infinite domain.
two global lists are available for diagnosis
after invocation:
gargL is a list of method and args input,
goutL is a list of method and quadpack total output.
If details is set to true, nint will print out
the method used.
*/
nint (expr,_var%,_a%, _b%) :=
block ([prec : 1d-4,limval : 600, keyval:3, argL,rL, gL,gsL,ans ],
if _a% = minf or _b% = inf then (
argL : [expr,_var%,_a%, _b%, epsrel=prec,limit=limval],
gargL : cons (qagi, argL),
if details then print ("quad_qagi"),
rL : apply ('quad_qagi, argL),
goutL : cons (qagi,rL),
if rL[4] = 0 then (
ans : rL[1],
gargL : [qagi,expr,_var%,_a%, _b%, epsrel=prec,limit=limval],
goutL : cons (qagi,rL))
else (
ans : false,
print (" prec = ",prec," limval = ",limval),
print (" quad_qagi error code = ", part (codelist,rL[4],2))))
else (
/* finite interval case */
gL : apply ('quad_qag,
[expr,_var%,_a%,_b%,keyval,epsrel=prec,limit=limval]),
gsL : apply ('quad_qags,
[expr,_var%,_a%,_b%,epsrel=prec,limit=limval]),
/* select winner */
if gL[4] = 0 and gsL[4] = 0 then ( /* no errors case */
if gL[3] < gsL[3] then (
ans : gL[1],
if details then print ("quad_qag"),
gargL :
[qag,expr,_var%,_a%,_b%,keyval,epsrel=prec,limit=limval],
goutL : cons (qag,gL))
else (
ans : gsL[1],
if details then print ("quad_qags"),
gargL : [qags,expr,_var%,_a%,_b%,epsrel=prec,limit=limval],
goutL : cons (qags,gsL) ))
else if gL[4] = 0 then (
ans : gL[1],
if details then print ("quad_qag"),
gargL :
[qag,expr,_var%,_a%,_b%,keyval,epsrel=prec,limit=limval],
goutL : cons (qag,gL))
else if gsL[4] = 0 then (
ans : gsL[1],
if details then print ("quad_qags"),
gargL : [qags,expr,_var%,_a%,_b%,epsrel=prec,limit=limval],
goutL : cons (qags,gsL))
else (
ans : false,
print (" keyval = ",keyval," prec = ",prec," limval =
",limval),
print (" quad_qag error code = ", part (codelist,gL[4],2)),
print (" quad_qags error code = ", part (codelist,gsL[4],2)))),
ans)$
quad_j4save('control, 0)$
display2d:false$
details:false$
-------------------------------------------------------------------------
Ted Woollett