Bisection Plus Program

[Namir]

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 X2

Memory 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 RTN

Usage

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).