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: Fri, 16 Dec 2011 12:53:15 -0800
On Dec 15, 2011, Raymond Toy wrote:
----------------------------
>I think you need to do something like this:
>
>quad_qags(lambda([x], block([r : quad_qags(x*y^2,y,0,1)], print(r[4]),
> r[1])), x, 0,
> 1);
>
>0
>0
>0
>...
>0
> [.166666, 1.85e-15, 21, 0]
----------------------------
Thanks for the suggested approach.
I need to encase this in a function (eventually
part of nint) which accepts a general expression
and lists which include the variables of integration
and the limits.
I want to pass the general expression and limit args to
quad_qags for both the inner and outer integrals.
I try to do this in the following.
The lambda version you suggest works ok with
strictly numerical limits, but not for variable
limits on 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) qarea(ex,uL,vL) :=
block([u,u1,u2,v,v1,v2,u0,expr,v11,v22],
[u,u1,u2] : uL,
[v,v1,v2] : vL,
display(u,u1,u2,v,v1,v2),
expr:subst(u0,u,ex),
/* allow for variable limits on the inner integral */
v11:subst(u0,u,v1),
v22:subst(u0,u,v2),
display(expr,v11,v22),
quad_qags(lambda([u0], block([r : quad_qags(expr,v,v11,v22)],
if listp(r) then print(r[4])
else print(r),
r[1])), u0, u1, u2))$
/* numerical limits work: */
(%i2) qarea(x*y^2,[x,0,1],[y,0,1]);
u = x
u1 = 0
u2 = 1
v = y
v1 = 0
v2 = 1
expr = u0*y^2
v11 = 0
v22 = 1
0
0
...
0
(%o2) [0.16666666666667,1.8503717077085952E-15,21,0]
--------------------------------
/* variable limits on the inner integral
doesn't work: */
(%i3) qarea(x*y^2,[x,0,1],[y,0,x]);
u = x
u1 = 0
u2 = 1
v = y
v1 = 0
v2 = x
expr = u0*y^2
v11 = 0
v22 = u0
quad_qags(u0*y^2,y,0,u0,epsrel = 1.0E-8,epsabs = 0.0,limit = 200)
(%o3) quad_qags(lambda([u0],
block([r:quad_qags(expr,v,v11,v22)],
if listp(r) then print(r[4]) else print(r),
r[1])),u0,0,1,epsrel = 1.0E-8,epsabs = 0.0,limit = 200)
----------------------------------
Ted