Just a lazy solver algortihm « Next Oldest | Next Newest »

 ▼ PGILLET Junior Member Posts: 16 Threads: 6 Joined: Dec 2011 06-28-2013, 07:01 PM 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=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. ▼ Namir Posting Freak Posts: 2,247 Threads: 200 Joined: Jun 2005 06-28-2013, 11:47 PM Check my multi-root solver "scan range method" presentation at HHC2012 on YouTube. Click here.

 Possibly Related Threads... Thread Author Replies Views Last Post hp-prime solver and variable name fabrice48 22 6,880 12-10-2013, 03:25 AM Last Post: fabrice48 HP Prime Triangle solver BruceH 29 7,034 11-28-2013, 12:03 AM Last Post: Dale Reed Using units in Numeric Solver Harold A Climer 1 990 10-13-2013, 10:44 AM Last Post: Tim Wessman Does Prime Have a Multiple Equation Solver? Norman Dziedzic 2 1,086 09-20-2013, 09:43 AM Last Post: Norman Dziedzic [43s] : How the solver will be implemented Miguel Toro 3 1,307 03-14-2013, 06:09 PM Last Post: Walter B TVM-Solver for the PC fhub 14 3,280 12-26-2012, 03:24 PM Last Post: fhub [WP34s] New TVM-solver version fhub 43 8,578 12-26-2012, 06:12 AM Last Post: fhub HP-Solver Mike (Stgt) 2 926 10-10-2012, 02:44 AM Last Post: Mike (Stgt) WP34S solver question Reth 22 5,007 07-13-2012, 06:55 PM Last Post: Paul Dale Iterative vs. direct solver Don Shepherd 12 2,923 05-18-2012, 06:28 PM Last Post: Paul Dale

Forum Jump: