I'd like to submit you the program attached (Hp42_11_EV.raw), which
attempts to find all real eigenvalues of a real (square) matrix, using
SOLVE & doing automatically all the necessary guesses.
It's based on a personal algorithm which could also be used for finding
roots of a function: it "scans" intervals (Guess=middle of interval,
then SOLVE) and uses a stack (feeded with intervals to scan later
(PUSH) or discarded if to scan now (POP))
1) Input the initial interval [L,U] (which should contain all roots)
(For the matrix eigenvalues, I used [-FNRM,FNRM] automatically )
2) Input the minimum difference in X ("Diff X") (in order to reject a
proposed root (not enough different) or to consider that an interval is
too little to explore, if difference is below or equal)
(e.g. 1E-03)
3) Input the maximum difference in f(X) ("Diff Y") (to detect an
extremum OR to round f(X) to zero, more quickly)
(e.g. 1E-15)
if ABS(f(X))<=DiffY => return f(X)=0
if ABS(f(X)-old f(X))<=DiffY => Set Flag 1 and return f(X)=0
(Extremum detected)
Set Old f(X) = f(X)
MAIN LOOP:
4) If Interval length (U-L) <= DiffX => too little interval =>
POP other interval =>Goto 4) or STOP if Stack empty
5) Solve, i.e.: Scan Interval using Guess G=(L+U)/2 (Result X OK if
return code 0, (if <>0=> POP other interval
=> Goto 4 or STOP if Stack empty))
6) if Result X >=(U-DiffX) => Upper half [G,U] is empty =>
Scan Lower Half [L,G] => Goto 4)
if Result X<=(L+DiffX) => Lower half [L,G] is empty =>
Scan Upper Half [G,U] => Goto 4)
if Result X>(L+DiffX) AND X<(U-DiffX)
=> Root found ! (Display value if not extremum; if extremum
(Flag 1 set) => treat it as a root, without displaying it)
if Result X>G => Upper part [G,X] is empty, PUSH
Lower half [L,G] and Scan now Upper part [X,U] => Goto 4)
if Result X<=G => Lower part [X,G] is empty, PUSH
Lower part [L,X] and Scan now Upper Half [G,U] => Goto 4)
- The stack begins at Register R03 and uses R01 as stack pointer:
e.g. R01=8 => there are 3 intervals in the stack: [L1,U1], [L2,U2],
[L3,U3] => R03=L1 R04=U1 R05=L2 R06=U2 R07=L3 R08=U3
if POP interval => L=R07 (L3), U=R08 (U3), Set R01=6 (8-2)
if PUSH a 4th interval [L4, U4] => R09=L4 R10=U4, Set R01=10 (8+2)
- The automatic extremum detection furnished by SOLVE (return code 2)
caused problems sometimes (Stuck => I did my own (see above, using
Flag 1))
- I had to round f(X) to 0 more quickly than SOLVE (otherwise was stuck
sometimes)
Pierre GILLET (pierregillet AT swing DOT be)