How to monitor quad_qags error codes for the inner integral?
Subject: How to monitor quad_qags error codes for the inner integral?
From: Edwin Woollett
Date: Thu, 15 Dec 2011 12:12:20 -0800
I am so far not successful in trying to use quad_qags
for a double integral in which I monitor the error
code returned by the inner integral:
--------------------------------------
Maxima 5.25.1 http://maxima.sourceforge.net
using Lisp GNU Common Lisp (GCL) GCL 2.6.8 (a.k.a. GCL)
(%i1) float(integrate(integrate(x*y^2,y,0,1),x,0,1));
(%o1) 0.16666666666667
(%i2) quad_qags(quad_qags(x*y^2,y,0,1)[1],x,0,1);
(%o2) [0.16666666666667,1.8503717077085952E-15,21,0]
(%i3) qarea(ex,uL,vL) :=
block([%%u,u1,u2,%%v,v1,v2,%uu,verrL:[],
%ruL,%uans,%uerr,
exu : buildq([%%u,ex,v1,v2,%%v,verrL],lambda([u0],
block([%expr,%v1,%v2,%rL,%ans,%errcode],
%expr:subst(u0,%%u,ex),
%v1:subst(u0,%%u,v1),
%v2:subst(u0,%%u,v2),
%rL:quad_qags(%expr,%%v,%v1,%v2),
%ans : part(%rL,1),
%errcode : part(%rL,4),
if %errcode # 0 then verrL : cons(%errcode,verrL),
%ans)))],
[%%u,u1,u2] : uL,
[%%v,v1,v2] : vL,
display(%%u,u1,u2,%%v,v1,v2),
print (" exu(0.5) = ",exu(0.5)),
%ruL : quad_qags(exu(%uu),%uu,u1,u2,limit=700),
display(%ruL),
if not listp(%ruL) then
(print(" %ruL is not a list"),return(false)),
%uans : part(%ruL,1),
%uerr : part(%ruL,4),
if length(verrL) > 0 then
print("internal integral error code(s) = ", verrL),
if %uerr # 0 then
print(" external integral error code = ",%uerr),
%uans)$
(%i4) qarea (x*y^2, [x,0,1], [y,0,1] );
%%u = x
u1 = 0
u2 = 1
%%v = y
v1 = 0
v2 = 1
exu(0.5) = 0.16666666666667
%ruL = quad_qags(%uu*y^2,%uu,0,1,epsrel = 1.0E-8,epsabs = 0.0,limit = 700)
%ruL is not a list
(%o4) false
-----------------------------------------
Any suggestions about how I can make this work or use
another tactic?
Ted Woollett