HP Forums
Bisection Plus program - Printable Version

+- 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 X2

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

Usage

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