Edwin,
There are three issues here. One is a simple matter of quoting; one is a
variable scope problem; and the third is quad_qags's poor error handling.
When you write find_root(myg(y),y,1,3), all the arguments are evaluated
normally, and as it happens myg(y) evaluates to -5.0 (I'll discuss that
below). So you are calculating find_root(-5.0,y,1,3) which of course has
the same sign at the endpoints.
Maxima works this way so that you can write
expr: x^3-1$
find_root(expr,x,-2,2);
or
find_root(diff((x^3-1)/(x-2),x),x,-10,10);
Otherwise, you could only use explicit expressions as the argument to
find_root.
I see you also tried ''(myg(y)) and ev(myg(y)), but those are going in the
wrong direction. You don't want *more* evaluation of myg(y), but *less*:
find_root( '(myg(y)), y, 1, 3) /* This is preferable */
or
find_root('myg(y)), y, 1, 3)
Alas, this doesn't work, either. Why? Because find_root binds the value of
'y' when evaluating myg, and myg internally uses a free instance of the
variable 'y' -- remember that Maxima has dynamic scope. So you need to
explicitly make y local:
myg(x) := block([val, qlist, numer,y], /* make y local */
numer : true, qlist : quad_qags(2 *y, y, sqrt(5.0), x), val :
qlist[1] , val);
or quote it:
myg(x) := block([val, qlist, numer],
numer : true, qlist : quad_qags( 2*'y, 'y, sqrt(5.0), x), val :
qlist[1] , val);
or simply avoid the variable 'y' in the outer call (this works, but is not a
clean, general solution):
find_root( '(myg(q)),q,1,3 )
Finally, you might wonder why find_root(myg(y),y,1,3) didn't cause an error
in quad_qags in the original case where you called myg('y). It turns out
that quad_qags is sloppily written, and treats non-numerical endpoints as
zero rather than throwing an error:
quad_qags(x,x,a,b) => [0,0,21,0] /* !!! should give an error */
It also doesn't throw errors in the standard Maxima way, but returns an
error code as the 4th element of the list, which makes it easy to ignore
(your code, for example, doesn't check it). So even if it *had* checked the
non-numerical endpoint case, your code wouldn't have noticed.
I consider this a bug in quad_qags, or at least poor practice.
-s
-------------------------
On Nov 5, 2007 7:28 PM, Edwin Woollett <woollett at charter.net> wrote:
> myg(x):=block([val,qlist,numer],numer:true,qlist:quad_qags(2*y,y,sqrt(5.0
> ),x),val:qlist[1],
> val)$
> (%o3) done
> (%i4) myg(1.0);
> (%o4) - 4.000000000000001
> (%i5) myg(3.0);
> (%o5) 3.999999999999999
> (%i6) find_root(myg(y),y,1,3);
> function has same sign at endpoints
> [f(1.0) = - 5.000000000000002, f(3.0) = - 5.000000000000002]
> -- an error. To debug this try debugmode(true);
> (%i7) find_root(''myg(y),y,1,3);
> function has same sign at endpoints
> [f(1.0) = - 5.000000000000002, f(3.0) = - 5.000000000000002]
> -- an error. To debug this try debugmode(true);
> (%i8) find_root(ev(myg(y)),y,1,3);
> function has same sign at endpoints
> [f(1.0) = - 5.000000000000002, f(3.0) = - 5.000000000000002]
> -- an error. To debug this try debugmode(true);
>
>