+- HP Forums (https://www.hpmuseum.org/forum)
+-- Forum: HP Software Libraries (/forum-10.html)
+--- Forum: General Software Library (/forum-13.html)
+--- Thread: (71B) Bisection Plus for the HP-71B (/thread-419.html)
(71B) Bisection Plus for the HP-71B - Namir - 01-12-2014 05:44 PM
This program implements the new Bisection Plus algorithm for the HP-Prime. Click here to download a pdf file that discusses the new algorithm.
I have attached a ZIP file that contains the program for the Windows HP-71B emulator.
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| < toelrance
Return X2HP-71B Implementation
Code:
10 REM BISECTION PLUS
20 DESTROY ALL
30 DEF FNX(X)=EXP(X)-3*x^2
40 INPUT "A? ";A
50 INPUT "B? ";B
60 INPUT "TOLER? ";T
70 IF T<=0 OR T>=1 THEN T=1E-8
80 F7=FNX(A) @ F8=FNX(B)
90 IF F7*F8>0 THEN 380
95 X2=A @ I9=0
100 I9=I9+1 @ L=X2 @ X1=(A+B)/2
110 F1=FNX(X1)
120 IF F1*F7<0 THEN 160
130 S = (F7 – F1)/(A – X1)
140 I = F7 – S * A
150 GOTO 180
160 S = (F8 – F1)/(B – X1)
170 I = F8 – S * B
180 X2 = -I / S # DISP x2 @ WAIT .2
190 F2 = FNX(X2)
200 if F1*F2 > 0 THEN 260
210 A = X1
220 F7 = F1
230 B = X2
240 F8 = F2
250 GOTO 320
260 If F2*F7 < 0 THEN 300
270 A=X2
280 F7=F2
290 GOTO 320
300 B=X2
310 F8=F2
320 IF ABS(A - B) < T THEN 350
330 IF ABS(X2 - L) < T THEN 350
340 GOTO 100
350 DISP "ROOT = ";X2 @ PAUSE
360 DISP "ITERS = ";I9
370 GOTO 390
380 DISP "FA AND FB HAVE SAME SIGN"
390 ENDUsage
1. Press [RUN]. The program prompts you to enter the the value of A.
2. Key in the value of A and press [ENTER]. The program prompts you to enter the the value of B.
3. Key in the value of B and press [ENTER]. The program prompts you to enter the the value for the tolerance.
4. Key in the value for the tolerance and press [ENTER].
5. The program displays intermediate refinement for the guess. When it converges the program displays "ROOT=" followed by the root value.
6. Press [CONT]. The program displays the number of iterations.
Example
Find the root for f(x)=exp(x)-3*X^2 in the range [3, 4].
1. Press [RUN]. The program prompts you to enter the the value of A.
2. Key in 3 and press [ENTER]. The program prompts you to enter the the value of B.
3. Key in 4 and press [ENTER]. The program prompts you to enter the the value for the tolerance.
4. Key in 1E-8 and press [ENTER].
5. The program displays intermediate refinement for the guess. When it converges the program displays "ROOT=3.73308".
6. Press [CONT]. The program displays the number of iterations as "ITER=7".
Customization
Line 30 contains the definition of f(X). Edit this line to reflect your function.
RE: (71B) Bisection Plus for the HP-71B - C.Ret - 06-11-2021 06:37 PM
Hi,
I really appreciate this algorithm that stay in the inputted range.
Today, it serve me a bit seeking for the roots of a bunch of polynomials.
I just modify the original program to add two features :
- Function definition no more enter into the program listing (easier customization)
- A short cut exiting the loop when an exact root (Y=0) is fortunately found (this only happened with test polynomial with integer coefficients).
As useall, since I modified and enhance a few the original code, I post here to share for further user .
5 DIM F$[49]
10 INPUT "F(X)=",F$;F$ @ INPUT "[A,B]=";A,B @ INPUT "Tol=","1E-8";T @ INPUT "MaxIte=","30";M
12 X=B @ Y2=VAL(F$) @ X=A @ Y1=VAL(F$) @ IF Y1*Y2>0 THEN BEEP @ DISP "Fa & Fb = sign !" @ END
14 FOR N=1 TO M @ P=X @ X=(A+B)/2 @ Y=VAL(F$) @ IF ABS(B-A)<T THEN 24
16 IF Y*Y1<0 THEN S=(Y1-Y)/(A-X) @ I=S*A-Y1 ELSE S=(Y2-Y)/(B-X) @ I=S*B-Y2
18 L=X @ Z=Y @ X=I/S @ Y=VAL(F$) @ IF Y=0 OR ABS(P-X)<T THEN 24
20 IF Z*Y<0 THEN A=L @ Y1=Z @ B=X @ Y2=Y ELSE IF Y*Y2<0 THEN A=X @ Y1=Y ELSE B=X @ Y2=Y
22 DISP X @ NEXT N @ BEEP
24 DISP USING '3ZX,"F("K,")="SDE';N,T*(X DIV T),Y @ END
Usage:
- Press [RUN]. The program prompts you to enter the function definition. Be sure to enter the function without BASIC syntax error, using X (uppercase variable) as the independent variable. The previous function definition may be clear by pressing the [ON] key. Key in the new function definition or edit the proposed one and press [END LINE] to terminate the input.
- The program prompt you for the range [A,B]. Key in the values of A and the value of B separated by a comma. Press [END LINE] to validate.
- Key in the value for the tolerance or use the default one.
- Key in the maximum iteration count and press [END LINE].
The program displays intermediate refinement for the guess. Eventually, adjust display speed with DELAY command. This have to be tuned before starting the program.
When it converges or after the maximal iteration count, the program displays the root value in the format :
00n F( x.xxxxxxxx ) = ±#E-0##.
Where:
00n indicate iteration count.
x.xxxxxxxx is the root \( x_0 \) rounded to the tolerance.
±#E-0## is the sign and amplitude of the residual \( y=F(x_0)\approx 0 \)
Another root can be seek for the same function or for a new function, just press [ f ][CONT] or [RUN] to process. All the parameters can be kept or modified as wish... no edition of the program is needed.
Step 1: Root of \( f(x)=e^x-3x^2 \) in \( \left [\,-2\,,\,0\,\right ] \)
Code:
[RUN]
F(X)=_ [f][EXP(] X ) - 3 * X [g] ^ 2 [END LINE]
[A,B]=_ -2,0[END LINE]
Tol=1E-8 [→][→][→][4][END LINE] ; 1E-8 change to 1E-4
MaxIte=30 [END LINE]
-.275321257597 -.432885068285 -.457442411298 -.458918397673
005 F(-.4589)=+2E-006 [f][CONT]Code:
F(X)=EXP(X)-3*X^2 [END LINE]
[A,B]=_ 0,3[END LINE]
TOL=1E-8 [END LINE]
MaxIte=30 [END LINE]
.458952661568 .883433519492 .909669444836 .91000548197 .910004566002
006 F(.91000757)=+2E-011 [f][CONT]Code:
F(X)=EXP(X)-3*X^2 [END LINE]
[A,B]= 2,4[END LINE]
TOL=1E-8 [END LINE]
MaxIte=30 [END LINE]
3.67759480515 3.72781796368 3.73284101745 3.73307361257 3.73307896662 3.73307902828
007 F(3.73307902)=+1E-010 [OFF]Edited 25.06.2021 for typos and examples reformat