Just a lazy solver algortihm
Message #1 Posted by PGILLET on 28 June 2013, 7:01 p.m.
I've had much fun implementing this easy but lazy & fat imperfect solver algorithm:
The user gives an interval and the program tries to find all the roots of a function in this interval, doing all the guesses automatically.
Each time the program calls the built-in solver (with Guess in the middle of an interval) and finds a root or extremum, a new interval is submitted (with new Lower and Upper bounds) and, if necessary, another one is put on a stack for later use (PUSH).
If an interval is too small (the user defines tolerance xtol) or if the built-in solvers does not returns a root or extremum, the next interval is taken from the stack (POP), until the stack is empty.
If the built-in solver returns a root (or extremum), the program displays it (if root) and explores these cases:
Case 1: Lower+xtol < Root/Extremum < Guess < Upper:
"Assume" that Interval(Root/Extremum,Guess) is empty (what a shame:)
PUSH Interval(Lower,Root/Extremum)
SOLVE with guess in the middle of Interval(Guess,Upper)
Case 2: Lower < Guess < Root/Extremum < Upper-xtol
"Assume" that Interval(Guess,Root/extremum) is empty (what a shame:)
PUSH Interval(Lower,Guess)
SOLVE with guess in the middle of Interval(Root/Extremum,Upper)
Case 3: Lower < Guess < Upper-xtol <= Root/extremum
"Assume" Interval(Guess,Upper) is empty (what a shame:)
SOLVE with guess in the middle of Interval(Lower,Guess)
Case 4: Root/Extremum <= Lower+xtol < Guess < Upper
"Assume" Interval(Lower,Guess) is empty (what a shame:)
SOLVE with guess in the middle of Interval(Guess,Upper)
IN the subroutine called by the built-in solver for evaluating the function to solve, I added code to more quickly round the function to 0 or detect extremum (f(x) very close to previous one) (the user defines tolerance ytol). A root or extremum will be then more easily found by the built-in solver.
The stack uses 2 registers per interval (Lower and Upper), 1 register as stack pointer and indirect adresssing.
This lazy algorithm seems to be working and useful for a quick tour of the roots of a "well-behaved" function.
For finding the roots of sin(x) in [-10,10], with xtol=1E-3 and ytol=1E-15, on Free42, this algorithm quickly finds the 7 roots and needed at most 3 intervals on stack and ... hum ... 95 calls to solver (what a shame :).
For HP42S:
00 { 531-Byte Prgm }
01>LBL "FXSLV"
SAVE CALCULATOR STACK
02 STO "FXX"
03 Rv
04 STO "FXY"
05 Rv
06 STO "FXZ"
07 Rv
08 STO "FXT"
09 Rv
INITIALIZATIONS
10 EXITALL
11 "LOWER ?"
12 -100
13 PROMPT
14 STO "L"
15 "UPPER ?"
16 100
17 PROMPT
18 STO "U"
19 4
20 STO 01 (Stack pointer:=Stack Top)
21 0
22 STO 02 (Highest Stack pointer)
23 STO 03 (Counts calls to built-in solver)
24 "DIFF X ?" (xtol)
25 1E-3
26 PROMPT
27 STO "P"
28 "DIFF Y ?" (ytol)
29 1E-15
30 PROMPT
31 STO "P2"
32 PGMSLV "FX" (Function is defined in label FX)
33 GTO "SCAN"
GET NEXT INTERVAL
34>LBL "NEXT"
35 2
36 RCL 01
37 X<=Y?
38 GTO "STOP"
39 XEQ "POP"
SOLVE WITH GUESS IN THE MIDDLE OF INTERVAL
40>LBL "SCAN"
41 RCL "U"
42 RCL "L"
43 -
44 RCL "P"
45 X>=Y? (Too small interval ?=> Get next interval)
46 GTO "NEXT"
47 1
48 STO+ 03 (One more call to built-in solver)
49 0
50 STO "D" (Previous f(x):=0)
51 RCL "L"
52 RCL "U"
53 +
54 2
55 ÷
56 STO "G" (Guess=(Lower+Upper)/2
57 STO "X"
58 ENTER (Guess2=Guess1)
59 SOLVE "X"
60 R^ (X=Return code (normally 0 if Root (or extremum in my case (see code at end of FX label))
61 X!=0? (Not root or extremum => Get next interval)
62 GTO "NEXT"
63 RCL "U"
64 RCL "P"
65 -
66 RCL "X"
67 X<Y? (Root>=Upper-xtol ? => See Case 3)
68 GTO 01
69 RCL "G"
70 STO "U"
71 GTO "SCAN" (Submit Interval(Lower,Guess)
72>LBL 01
73 RCL "L"
74 RCL "P"
75 +
76 RCL "X"
77 X>Y? (Root<=Lower+xtol ? => See Case 4)
78 GTO "FOUND" (Root or extremum found in interval)
79 RCL "G"
80 STO "L"
81 GTO "SCAN" (Submit Interval(Guess,Upper)
ROOT OR EXTREMUM FOUND
82>LBL "FOUND"
83 "ROOT="
84 FC? 01
85 PROMPT (Display if root)
86 RCL "G"
87 RCL "X"
88 X>Y? (Root/Extremum<=Guess ? => See Case 1)
89 GTO 02
90 RCL "X"
91 XEQ "PUSH" (PUSH Interval(Lower,Root/Extremum))
92 RCL "G"
93 STO "L"
94 GTO "SCAN" (Submit Interval(Guess,Upper)
95>LBL 02 (Root/Extremum>Guess ? => See Case 2)
96 RCL "G"
97 XEQ "PUSH" (PUSH Interval(Lower,Guess))
98 RCL "X"
99 STO "L"
100 GTO "SCAN" (Submit Interval(Root/Extremum,Upper)
STOP AND EXIT
101>LBL "STOP"
102 "STOP"
103 AVIEW
104 RCL "FXT" (Restore calculator stack)
105 RCL "FXZ"
106 RCL "FXY"
107 RCL "FXX"
108 CLV "FXX" (Clear variables)
109 CLV "FXY"
110 CLV "FXZ"
111 CLV "FXT"
112 CLV "L"
113 CLV "U"
114 CLV "G"
115 CLV "P"
116 CLV "P2"
117 CLV "D"
118 CF 01
119 RTN
FUNCTION DEFINITION + CODE
120>LBL "FX"
121 MVAR "X"
122 RCL "X"
123 X^2
124 3
125 -
126 CF 01 (Clear Flag1, will be set if extremum found)
127 0
128 RCL ST Y
129 ABS
130 RCL "P2" (TZYX=(f(x),0,f(x),ytol))
131 X>=Y? (f(x)<=ytol ? => Round to 0)
132 R^ (TZYX=(0,f(x),ytol,f(x)))
133 R^ (TZYX=(f(x),ytol,f(x),0) or (0,f(x),ytol,f(x))
134 RCL "D"
135 X<>Y
136 STO "D" (Previous f(x):=f(x) or 0)
137 X=0? (f(x)=0 or rounded to 0 ? => Return)
138 GTO 04
139 -
140 ABS
141 RCL "P2"
142 X>=Y? (ABS(fx)-previous(fx))<=ytol ? => Extremum => Set flag 1
143 SF 01
144 RCL "D" (f(x))
145 FS? 01 (flag 1 set ? (extremum) => Force value 0 for function)
146 0
147>LBL 04
148 RTN
PUSH
149>LBL "PUSH"
150 1
151 STO+ 01
152 RCL "L" (Lower unchanged)
153 STO IND 01
154 1
155 STO+ 01
156 R^ (Future Upper=Calculator stack register X)
157 STO IND 01
158 RCL 01
159 RCL 02
160 X>=Y? (New highest Stack pointer ? => Put it in R02
161 GTO 05
162 X<>Y
163 STO 02 (New highest Stack pointer)
164>LBL 05
165 RTN
POP
166>LBL "POP"
167 RCL IND 01
168 STO "U" (New Upper)
169 1
170 STO- 01
171 RCL IND 01
172 STO "L" (New Lower)
173 1
174 STO- 01
175 RTN
176 .END.
| 
Go back to the main exhibit hall