Subject: Function evaluates wrong when called from dblint
From: Neilen Marais
Date: Thu, 09 Jun 2011 11:37:27 +0200
I'm trying to integrate a function using dblint. The original integrand
is specified in spherical coordinates, but I need to evaluate it on
rectangular domains in order to compare it to the output of a numerical
program I'm writing. In the example below the two functions f1() and
f2() are identical apart from the debugging print statement. Note how in
both cases f1(0.0, 0.0) and f2(0.0, 0.0) both evaluate to zero when
called directly; when called from dblint the same parameter values
result in something completely different as shown by the print
statement!
Am I doing something wrong here, or is there a bug going on? I'm using
Maxima version 5.21.1 as packaged in 64-bit Ubuntu 10.10.
Thanks
Neilen
(%i1) ratprint:false$
(%i2) sublis:[z=1/2]$
(%i3) f(x,y):=block([
r:sqrt(x^2+y^2+z^2), theta:acos(z/r),
phi:(if is(equal(x, 0)) and is(equal(y, 0)) then 0 else atan2(y,x)),
tmp],
r:float(subst(sublis, r)),
theta:float(subst(sublis,theta)),
phi:float(subst(sublis,phi)),
tmp:float(ev(subst(sublis, expr), eval)),
/*(if tmp > 10000 then print([x, y, tmp, r, theta, phi]) else 0), */
print([x, y, tmp, r, theta, phi]),
tmp)$
(%i4) expr:-(200*%pi*cos(phi)*cos(theta)*sin(theta)*cos((1000000000*%pi*z)/149896229-(1000000000*%pi*r)/149896229))/r$
(%i5) f(float(0.0), float(0.0));
[0.0,0.0,0.0,0.5,acos(0.5/r),0.0]
(%o5) 0.0
(%i6) load (dblint);
(%o6) "/usr/share/maxima/5.21.1/share/numeric/dblint.mac"
(%i7) f2(x,y):=block([
r:sqrt(x^2+y^2+z^2), theta:acos(z/r),
phi:(if is(equal(x, 0)) and is(equal(y, 0)) then 0 else atan2(y,x)),
tmp],
r:float(subst(sublis, r)),
theta:float(subst(sublis,theta)),
phi:float(subst(sublis,phi)),
tmp:float(ev(subst(sublis, expr), eval)),
(if tmp > 10000 then print([x, y, tmp, r, theta, phi]) else 0),
/*print([x, y, tmp, r, theta, phi])*/
tmp)$
(%i8) f2(float(0), float(0));
(%o8) 0.0
(%i9) dblint(f2, lambda([x],-1/2), lambda([x], 1/2), -1/2, 1/2);
[0.0,0.0,6.304713613808152*10^58,1.4850279263678455*10^-64,3.141592653589793-1.0*acos((1.4850279263678455*10^-64)/r),0.0]
(%o9) 7.0052373486757241*10^55