+- HP Forums (https://www.hpmuseum.org/forum)
+-- Forum: HP Software Libraries (/forum-10.html)
+--- Forum: HP-65/67/97 Software Library (/forum-12.html)
+--- Thread: Bisection Plus program (/thread-459.html)
Bisection Plus program - Namir - 01-18-2014 04:43 AM
Here is an HP-67 implementation for the new Bisection Plus algorithm that calculators the root of a nonlinear function.
Algorithm
Code:
Given f(x)=0, the root-bracketing interval [A,B], and the tolerance for the root of f(x):
Calculate Fa = f(A) and Fb=f(B).
Exit if Fa*Fb > 0
X2 = A
Repeat
LastX2 = X2
X1=(A+B)/2
Fx1 = f(X1)
If Fx1*Fa > 0 then
Slope = (Fb – Fx1)/(B – X1)
Intercept = Fb – Slope * B
Else
Slope = (Fa – Fx1)/(A – X1)
Intercept = Fa – Slope * A
End If
X2=-Intercept / Slope
Fx2 = f(X2)
If Fx1*Fx2 < 0 then
A = X1
Fa = Fx1
B = X2
Fb = Fx2
Else
If Fx2*Fa > 0 then
A=X2
Fa=Fx2
Else
B=X2
Fb=Fx2
End If
End If
Until |A-B| < tolerance OR |X2 - LastX2| < tolerance
Return X2Memory Map
R0 = Tolerance
R1 = A
R2 = B
R3 = X1
R4 = X2
R5 = F(A)
R6 = F(B)
R7 = F(X1)
R8 = F(X2)
R9 = Slope
RA = LastX2
RB = Iterations
Listing
Code:
1 LBL A
2 STO 1 # Store A
3 RTN
4 LBL B
5 STO 2 # Store B
6 RTN
7 LBL 8 # Supply the defualt tolerance value
8 EEX
9 8
10 CHS
11 RTN
12 LBL C
13 X=0? # is the tolerance equal to zero?
14 GSB 8 # get the default value
15 STO 0 # Store tolerance
16 RCL 1
17 GSB E # Calculate f(A)
18 STO 5
19 RCL 2
20 GSB E # Calculate f(B)
21 STO 6
22 RCL 5
23 *
24 X>0? # Do the values of f(A) and f(B) have the same sign?
25 GTO 7 # jump to error message
26 RCL 1
27 STO 4 # X2=A
28 0
29 STO B # Iteration counter = 0
30 LBL 0 # ---------------- start iterations loop
31 1
32 RCL B
33 +
34 STO B # Increment iteration counter
35 RCL 4
36 STO A # LastX2 = X2
37 RCL 1
38 RCL 2
39 +
40 2
41 / # X1 = (A+B)/2
42 STO 3
43 GSB E # Calculate f(X1)
44 STO 7
45 RCL 5
46 *
47 X>0? # Do FX1 and FA have the same sign?
48 GTO 1
49 RCL 5
50 RCL 7
51 -
52 RCL 1
53 RCL 3
54 -
55 / # Calculate slope = (FA - FX)/(A - X)
56 STO 9
57 RCL 1
58 *
59 RCL 5
60 - # Calculate the negative value of the intercept
61 GTO 2 # Finish calculations of X2 and f(X2) in lbl 02
62 LBL 1
63 RCL 6
64 RCL 7
65 -
66 RCL 2
67 RCL 3
68 -
69 / # Calculate slope = (FB - FX)/(B - X)
70 STO 9
71 RCL 2
72 *
73 RCL 6
74 - # Calculate the negative value of the intercept
75 LBL 2
76 RCL 9
77 / # Calculate X2
78 STO 04
79 PAUSE # Pause to display X2
80 GSB E # Calculate f(X2)
81 STO 8
82 RCL 7
83 RCL 8
84 *
85 X>0? # DO f(X1) and f(X2) have the same sing?
86 GTO 3 # Replace A or B with X2
87 RCL 3
88 STO 1 # A = X1
89 RCL 4
90 STO 2 # B = X2
91 RCL 7
92 STO 5 # F(A) = F(X1)
93 RCL 8
94 STO 6 # F(B) = F(X2)
95 GTO 5 # jump to section that test for convergence
96 LBL 3
97 RCL 8
98 RCL 5
99 *
100 X>0? # Do f(X2) annd f(A) have the same signs?
101 GTO 4
102 RCL 4
103 STO 2 # B = X2
104 RCL 8
105 STO 6 # f(B) = f(X2)
106 GTO 5 # jump to section that test for convergence
107 LBL 4
108 RCL 4
109 STO 1 # A = X2
110 RCL 8
111 STO 5 # f(A) = f(X2)
112 LBL 5 # start testing for convergence
113 RCL 4
114 RCL A
115 -
116 ABS
117 RCL 0
118 X>Y? # |X2 - LastX2| < tolerance?
119 GTO 6 # exit loop and show the value for the root
120 RCL 1
121 RCL 2
122 -
123 ABS
124 X<>Y
125 X>Y? # |A - B| < tolerance?
126 GTO 6 # exit loop and show the value for the root
127 GTO 0 # ---------------- end of iterations loop
128 LBL 6
129 RCL 4
130 STOP
131 RCL B
132 RTN
133 LBL 7
134 0
135 /
136 RTN
137 LBL E # f(x) = GTO 3 ou can edit the code after LBL E and until RTN
138 EXP
139 LASTX
140 X^2
141 3
142 *
143 -
144 RTNUsage
1. Enter the value fo A and press the key A.
2. Enter the value for B and press the key B.
3. Enter the value for the tolerance (or enter 0 to use the default tolerance value of 1E-8) and press the key C.
4. The program displays intermediate refinements for the roots. It stops with the root value.
5. Press the R/S key to view the number of iterations.
6. (Optional) to evaluate f(X) at any value of X, enter that value and then press the key E.
Example
Find the root for f(x)=exp(x)-3*X^2 in the range [3,4].
1. Enter 3 and press the key A.
2. Enter 4 and press the key B.
3. Enter 0 (to use the default tolerance value of 1E-8) and press the key C.
4. The program displays intermediate refinements for the roots. It stops with the root value of 3.73308.
5. Press the R/S key to view the number of iterations. The program displays 7.
Customization
Label E contains the commands that implement f(X). Edit the code after label E to insert the code for your f(X). Registers C, D, E, and the secondary registers are available to store the values of X and any other intermediate values used in calculating f(X). To use the secondary registers, I recommend the following template for LBL E:
LBL E
P<>S
STO 0
commands to calculate f(x) using R0 to R9.
P<>S
RTN