Function evaluates wrong when called from dblint

On 6/9/11, Neilen Marais <nmarais at gmail.com> wrote:

> (%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

I nominate this one for the Maxima "unintended consequences" hall of fame.

There is a variable z in dblint and the repeated evaluation and
substitution in f2 pulls in (courtesy of dynamic binding, i.e. non-lexical)
whatever is the current value of z. Try tracing subst -- you'll see
stuff like this:
 Enter subst [[1.529155181083518 = 1/2],
... that can't be good.

If you rename z in expr, sublis, and f2 to zz or something like that,
or rename z in dblint to something else, I think you'll get a different result.

I think the integrand function can be simplified substantially and that
should make it easier to see what it is supposed to do.
My general advice is to construct the integrand with the desired value
of z pasted into it (instead of trying to substitute the value at run time).

Some other minor items: do assignments outside of the list of local
variables in block, so the assignments are evaluated in sequence,
not in parallel (which makes a difference in this case because there is
an assignment depends on the previous one); and predicates in "if"
need not be enclosed in is(...).

I was tinkering with different ways to avoid the dynamic binding name
conflict problem. Here are two ideas. One is to use a macro (FOO)
to construct the integrand function with the z value pasted into it.
The other idea is to replace "block" in dblint with a lexical block.
I have a naive implementation of such a thing, called BLEX.
That makes the two variables named "z" distinct, so you don't need
to rename z in your code.

Here's the macro:

FOO (f, expr, zeq) ::= subst (zeq, buildq ([f, expr],
  f (x, y) := block ([r, theta, phi],
    r : sqrt (x^2 + y^2 + z^2),
    theta : acos (z / r),
    phi : if equal (x, 0) and equal (y, 0) then 0 else atan2 (y, x),
    float (expr))));

FOO (f2b, ''expr, ''sublis);
grind (f2b);
 =>
f2b(x,y):=block([r,theta,phi],r:sqrt(y^2+x^2+1/4),theta:acos(1/(2*r)),
    phi:if equal(x,0) and equal(y,0) then 0 else atan2(y,x),
    float(-200*%pi*cos(phi)*cos(1000000000*%pi*r/149896229-500000000*%pi/149896229)*cos(theta)
              *sin(theta)
           /r))$

dblint (f2b, lambda ([x], -1/2), lambda ([x], 1/2), -1/2, 1/2);
 => 1.32634644008552E-14


Here's the lexical block idea:

load ("/mnt/windows2/Users/robert/maxima/blex.mac");
load ("./dblint_blex.mac");
dblint_blex (f1, lambda ([x], -1/2), lambda ([x], 1/2), -1/2, 1/2);
 => 1.4210854715202E-14


I'm guessing that lexical scope would be a wonderful thing.
If someone wants to work on that with me, that would be terrific.

Robert Dodier

PS. I've attached the various scripts I was tinkering with,
in case anybody wants to take a look.
-------------- next part --------------
A non-text attachment was scrubbed...
Name: marais_dblint.mac
Type: application/octet-stream
Size: 1745 bytes
Desc: not available
URL: <http://www.math.utexas.edu/pipermail/maxima/attachments/20110612/c119a573/attachment-0004.obj>;
-------------- next part --------------
A non-text attachment was scrubbed...
Name: marais_dblint.out
Type: application/octet-stream
Size: 12342 bytes
Desc: not available
URL: <http://www.math.utexas.edu/pipermail/maxima/attachments/20110612/c119a573/attachment-0005.obj>;
-------------- next part --------------
A non-text attachment was scrubbed...
Name: dblint_blex.mac
Type: application/octet-stream
Size: 1057 bytes
Desc: not available
URL: <http://www.math.utexas.edu/pipermail/maxima/attachments/20110612/c119a573/attachment-0006.obj>;
-------------- next part --------------
A non-text attachment was scrubbed...
Name: blex.mac
Type: application/octet-stream
Size: 204 bytes
Desc: not available
URL: <http://www.math.utexas.edu/pipermail/maxima/attachments/20110612/c119a573/attachment-0007.obj>;