#1040 find_root evaluates its arguments strangely

closed
nobody
5
2008-07-30
2006-12-11
Anonymous
No

Someone asked a question on the mailing list about using the output of solve to do numerical root finding.

we found that you could not use the output of solve directly in find_root because it evaluates its first argument strangely. Instead it required explicit evaluation with the '' operator.

This behavior is very confusing even for relatively advanced maxima users, but especially for new users.

example of the problem below:

-----Transcript------

Maxima 5.10.0 http://maxima.sourceforge.net
Using Lisp GNU Common Lisp (GCL) GCL 2.6.7 (aka GCL)

(%i1) h1(s):=1/(1+tau[1]*s);
h2(s):=1/(1+tau[2]*s);
x(s):=1/s;
tau[1]:33e-6;
tau[2]:33e-6;
eqn:first(solve(ilt(h1(s)*h2(s)*x(s),s,t)=0.1,t));
1
(%o1) h1(s) := ----------
1 + tau s
1
(%i2) 1
(%o2) h2(s) := ----------
1 + tau s
2
(%i3) 1
(%o3) x(s) := -
s
(%i4) (%o4) 3.2999999999999996E-5
(%i5) (%o5) 3.2999999999999996E-5
(%i6)
`rat' replaced 3.2999999999999996E-5 by 33//1000000 = 3.3000000000000003E-5

`rat' replaced 0.9 by 9//10 = 0.9
1000000 t
---------
33
297 %e - 330
(%o6) t = ---------------------
10000000
(%i7) find_root(%,t,0,2);

Maxima encountered a Lisp error:

Error in MACSYMA-TOP-LEVEL [or a callee]: ((MEQUAL SIMP) $T
((MTIMES SIMP)
((RAT SIMP) 1 10000000)
((MPLUS SIMP) -330
((MTIMES SIMP) 297
((MEXPT SIMP) $%E
((MTIMES SIMP)
((RAT SIMP) 1000000
33)
$T)))))) is not of type (OR
RATIONAL
LISP:FLOAT).

Discussion

  • Robert Dodier

    Robert Dodier - 2006-12-27
    • labels: 460523 --> Lisp Core - Solving equations
     
  • Robert Dodier

    Robert Dodier - 2008-07-30

    Logged In: YES
    user_id=501686
    Originator: NO

    Since a revision or two ago, find_root evaluates its arguments as an ordinary function.

    This particular example can't be solved by find_root(..., 0, 2) because the expression overflows at 2. However find_root(..., 0, 1/1000) succeeds and yields 1.7549783076857198E-5.

    Closing this report as fixed.

     
  • Robert Dodier

    Robert Dodier - 2008-07-30
    • status: open --> closed
     

Log in to post a comment.