+- HP Forums (https://www.hpmuseum.org/forum)
+-- Forum: HP Software Libraries (/forum-10.html)
+--- Forum: HP-41C Software Library (/forum-11.html)
+--- Thread: Bisection Plus Program (/thread-398.html)
Bisection Plus Program - Namir - 01-11-2014 11:29 AM
This program implements the new Bisection Plus algorithm for the HP-41C. The attachment contains a ZIP file. This file has a .raw file that you can load with HP-41C emulators that read .raw files. Click here to download a pdf file that discusses the new algorithm.
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
R00 = Tolerance
R01 = A
R02 = B
R03 = X1
R04 = X2
R05 = F(A)
R06 = F(B)
R07 = F(X1)
R08 = F(X2)
R09 = Slope
R10 = LastX2
R11 = Iterations
Listing
Code:
1 LBL "BISPLS"
2 LBL A
3 "A^B?"
4 PROMPT # Promot for [A, B]
5 STO 02
6 X<>Y
7 STO 01
8 "TOLER"
9 CF 22
10 PROMPT # Prompt for tolerance value
11 FC? 22
12 1E-8 # default tolerance is 1E-08
13 STO 00
14 RCL 01
15 XEQ E # Calculate f(A)
16 STO 05
17 RCL 02
18 XEQ E # Calculate f(B)
19 STO 06
20 RCL 05
21 *
22 X>0? # Do the values of f(A) and f(B) have the same sign?
23 GTO 07 # jump to error message
24 RCL 01
25 STO 04 # X2=A
26 0
27 STO 11 # Iteration counter = 0
28 LBL 00 # ---------------- start iterations loop
29 1
30 STO+ 11 # Increment iteration counter
31 RCL 04
32 STO 10 # LastX2 = X2
33 RCL 01
34 RCL 02
35 +
36 2
37 / # X1 = (A+B)/2
38 STO 03
39 XEQ E # Calculate f(X1)
40 STO 07
41 RCL 05
42 *
43 X>0? # Do FX1 and FA have the same sign?
44 GTO 01
45 RCL 05
46 RCL 07
47 -
48 RCL 01
49 RCL 03
50 -
51 / # Calculate slope = (FA - FX)/(A - X)
52 STO 09
53 RCL 01
54 *
55 RCL 05
56 - # Calculate the negative value of the intercept
57 GTO 02 # Finish calculations of X2 and f(X2) in lbl 02
58 LBL 01
59 RCL 06
60 RCL 07
61 -
62 RCL 02
63 RCL 03
64 -
65 / # Calculate slope = (FB - FX)/(B - X)
66 STO 09
67 RCL 02
68 *
69 RCL 06
70 - # Calculate the negative value of the intercept
71 LBL 02
72 RCL 09
73 / # Calculate X2
74 STO 04
75 PSE # Pause to display X2
76 XEQ E # Calculate f(X2)
77 STO 08
78 RCL 07
79 RCL 08
80 *
81 X>0? # DO f(X1) and f(X2) have the same sing?
82 GTO 03 # Replace A or B with X2
83 RCL 03
84 STO 01 # A = X1
85 RCL 04
86 STO 02 # B = X2
87 RCL 07
88 STO 05 # F(A) = F(X1)
89 RCL 08
90 STO 06 # F(B) = F(X2)
91 GTO 05 # jump to section that test for convergence
92 LBL 03
93 RCL 08
94 RCL 05
95 *
96 X>0? # Do f(X2) annd f(A) have the same signs?
97 GTO 04
98 LBL 04
99 RCL 04
100 STO 02 # B = X2
101 RCL 08
102 STO 06 # f(B) = f(X2)
103 GTO 05 # jump to section that test for convergence
104 LBL 04
105 RCL 04
106 STO 01 # A = X2
107 RCL 08
108 STO 05 # f(A) = f(X2)
109 LBL 05 # start testing for convergence
110 RCL 04
111 RCL 10
112 -
113 ABS
114 RCL 00
115 X>Y? # |X2 - LastX2| < tolerance?
116 GTO 06 # exit loop and show the value for the root
117 RCL 01
118 RCL 02
119 -
120 ABS
121 X<>Y
122 X>Y? # |A - B| < tolerance?
123 GTO 06 # exit loop and show the value for the root
124 GTO 00 # ---------------- end of iterations loop
125 LBL 06
126 "RT="
127 RCL 04
128 ARCL X
129 PROMPT # Display the refined root value
130 "ITRS="
131 ARCL 11
132 PROMPT
133 RTN
134 LBL 07
135 "FA*FB>0"
136 PROMPT
137 RTN
138 LBL E # f(x) = 0. ou can edit the code after LBL E and until RTN
139 EXP
140 LASTX
141 X^2
142 3
143 *
144 -
145 RTNUsage
1. Press [A}. The program prompts you to enter the root-bracketing interval [A, B].
2. Key in the value of A, press [ENTER}, key in the value of B, and then press [R/S].
3. The program prompts you to enter the tolerance value.
4. Enter the tolerance value and then press [R/S]. To accept the default value of 1E-8, just press the [R/S].
5. The program displays intermediate refinment for the guess. When it converges the program displays "RT=" followed by the root value. Th eprogram also displays the number of iterations.
8. (Optional) to evaluate f(X) at any value of X, enter that value and then press [E].
Example
Find the root for f(x)=exp(x)-3*X^2 in teh range [3,4].
1. Set the display mode to FIX 5.
2. Press [A]. The program displays the prompt "A^B?".
3 Key in 3, press {ENTER], key in 4, and then press [R/S].
4. The program prompts you to enter the tolerance value. Accept the defualt value of 1E-8 by simply pressing the [R/S] key.
5. The program pauses to display intermediate refinement for the root. When done, the program displays "RT=3.73308".
6. Press [R/S] to view the number of iterations. The program displays "ITRS=7.00000".
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 12 and up are available to store the values of X and any other intermediate values used in calculating f(X).