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.

Check my multi-root solver "scan range method" presentation at HHC2012 on YouTube. Click here.