Bisection, Secant and Newton Raphson Solver - J E Patterson

J E Patterson - jepspectro.com - 20231119

Description


Introduction

The three solvers in this program have their own particular strengths. With the additional choice of the built-in solver, in the hp-15c Simulator by Torsten Manz, solutions to most root finding problems can be found.

Bisection Solver


GSB A runs a bisection solver which acts on f(x) = 0 under label E. This program was written to test the hp-15c Simulator by Torsten Manz.

Label E can hold any equation of the form f(x) = 0. Note that some of the hp-15C Owner's Handbook examples may require some extra ENTER statements at the beginning as the stack is expected to be filled with x. The open box example on page 189 requires three ENTER statements at the beginning.

To compare with the built-in solver press 1 ENTER 2 f SOLVE E.

A problem from the hp-15c Owner's Handbook - page 189:

A 4 decimetre by 8 decimetre metal sheet is available, i.e. 400 mm by 800 mm.
The box should hold a volume V of 7.5 cubic decimetres, i.e. 7.5 litres.
How should the metal be folded for the tallest box in decimetres.
We are using decimetres rather than mm because equation entry is simplified.
Let x be the height.
Volume V = (8 - 2x)(4-2x)x. There are two sides and two ends of height x.
Rearrange to f(x) = 4((x - 6)x + 8)x - V = 0 and solve for the height x.

Instructions:

There are 3 solutions depending on the guesses.

0 ENTER 1 GSB A gives x = 0.2974 decimetres or 29.74 mm - a flat box.
1 ENTER 2 GSB A gives x = 1.5 decimetres or 150 mm - a reasonable height box.
2 ENTER 3 GSB A can't find a root as there are none in this interval. Press any key to stop.
3 ENTER 4 GSB A can't find a root as there are none in this interval. Press any key to stop.
4 ENTER 5 GSB A gives x = 4.206 decimetres or 420.26 mm - an impossible box height.

Bisection Solver Procedure:

Choose two guesses x1 and x0.
Do
Set flag 1
Calculate f(x0)
Let x2 = (x1 - x0)/2
Calculate f(x2)
If f(x2)*f(x0) < 0 let x1 = x2, clear flag 1
If flag 1 is set let x0 = x2
Loop until f(x2) = 0 within an error tolerance determined by the FIX, SCI or ENG significant digits setting.
Display root x2.

Note:

Choose guesses which bracket the root of interest. Adjust for other roots.

Secant Solver

GSB B runs a secant solver which acts on f(x) = 0 under label E.

RAN# is substituted for f(x0)-f(x1) if a divide by zero error would occur. If the equation to be solved has a square root term an error will occur for negative inputs. Use g TEST 2, 0 as the first two equation program statements. This is a test for a negative input and sets it to zero as the lowest valid input. Negative inputs may occur during the iteration.

Instructions:

The Secant solver uses the box problem to display two roots outside the bracketed guesses.
There are 3 roots depending on the guesses.

0 ENTER 1 GSB B gives x = 0.2974 decimetrers or 29.74 mm - a flat box.
1 ENTER 2 GSB B gives x = 1.5 decimetres or 150 mm - a reasonable height box.
2 ENTER 3 GSB B gives x = 1.5 decimetres or 150 mm - a reasonable height box.
3 ENTER 4 GSB B gives x = 4.2026 decimetres or 420.26 mm - an impossible box.
4 ENTER 2 GSB B points to the correct answer of x = 1.5 decimetres in spite of a potential divide by zero error.

Secant Solver Procedure:

Choose two starting points x0 and x1.
If x0 = x1 add 1E-7 to x1
Set f(x0) = RAN# a non-zero, non-integer starting value to avoid possible divide by zero. f(x0) is updated by f(x1) after one iteration.
Do
Calculate f(x1)
Let x2 = x1 - f(x1)*(x0 - x1)/(f(x0) - f(x1))
if f(x0)-f(x1) = 0 then substitute RAN#
If f(x1) = 0 let root = x1, display x1 and exit
else
x0 = x1
f(x0) = f(x1)
x1 = x2
Loop until f(x1) = 0 within an error tolerance determined by the FIX, SCI or ENG significant digits setting.
Display root x1.

Newton Raphson solver

GSB C runs a Newton Raphson solver which acts on f(x) = 0 under label E. Only one guess is required.

Instructions

For the box problem there are 3 solutions depending on the guesses.
0 GSB C gives x = 0.2974 decimetres or 29.74 mm - a flat box
1 GSB C gives x = 1.5 decimetres or 150 mm - a reasonable height box.
2 GSB C gives x = 1.5 decimetres or 150 mm - a reasonable height box.
3 GSB C gives x = 0.2974 decimetres or 29.74 mm - a flat box. Here the secant line now intersects the x axis below the smallest root.
4 GSB C gives x = 4.2026 decimetres or 420.26 mm - an impossible box.

Newton-Raphson Solver Procedure:

Choose a starting point x1
Do
Calculate f(x1)
If x1 = 0 use 1 instead to avoid a divide by zero from the derivative f'(x1)
define h = 1/10000
calculate f'(x1) ≈ (f(x1+h) - f(x1))/h
calculate f(x1)/f'(x1)
subtract from x1
Loop until f(x1) = 0 within an error tolerance determined by the FIX, SCI or ENG significant digits setting.
Display root x1.

Notes:

A TEST=0 ex statement in the Newton-Raphson solver program is a way of entering 1 without it attaching to the following EEX statement. TEST=0 10x or TEST=0 COS could also be used. Also the obvious TEST=0 1 ENTER would also do.

There are useful discussions on Newton's solver for the hp-12C here and here.

SOLVEkey.pdf has a good explanation of the additional tricks used to solve difficult equations using the built-in Solver.

In SCI and ENG modes on the DM15 and a real hp-15C some roots are not always found. This is because RND is implemented slightly differently in the hp-15C Simulator. Just stop the iteration by pressing any key and examine the register where x is held. This can be done with GSB D. Alternatively change to FIX mode - e.g. FIX 4 before solving.

Guesses can be tested with GSB E. For the Bisection solver look for a sign change. For the Secant solver with multiple roots try choosing both guesses to be on the same side of a root. This can alter the direction of the iteration. For the Newton-Raphson solver try moving the guess closer to the root or to the other side.

Equations are entered under label E. For example f(x) = 3x2 - 5x + 1 = 0 is entered using GTO E, changing to program mode with g P/R and deleting any entered function. X is placed on the stack twice as it is needed in two terms.

ENTER ENTER g x2 3 X X<>Y 5 X - 1 +

The roots are 0.2324 and 1.4343, using guesses 0,1 and 1,2.

An alternative form of this equation is: f(x) = x(3x-5)+1 = 0

The program is shorter: ENTER ENTER 3 X 5 - X 1 +.

Program Resources

Labels

Name Description Name Description
 A Bisection solve routine - needs two guesses in x and y  2 Bisection - Store x0 to x1 and x2 to x2
 B Secant solve routine - needs two guesses in x and y  3 Bisection - Store x1 to x0 and x2 to x1
 C Newton - Raphson Solve routine - needs one guess in x  4 Secant - Iteration loop
 D Recover root after program stop  5 Newton - iteration loop
 E Formula to be Solved, f(x) = 0  6 If two guesses are the same, add 1E-7 to R.1
 1 Bisection loop updates x2 and either x0 to x1 or x1 to x0

Storage Registers

Name Description
.0 x0
.1 x1
.2 Bisection - x2 = (x0+x1)/2, Secant x0-x1, Newton - h
.3 Bisection - f(x0), Secant - f(x0)
.4 Bisection f(x2), Secant - f(x1), Newton - f(x1)

Flags

Number Description
1 Bisection - If flag 1 is set x0 = x2 else x1 = x2

Program

Line Display Key Sequence Line Display Key Sequence Line Display Key Sequence
000 042 44 .1 STO . 1 084 44 .2 STO . 2
001 42,21,11 f LBL A 043 42 36 f RAN # 085 45,40, .1 RCL + . 1
002 44 .0 STO . 0 044 44 .3 STO . 3 086 32 15 GSB E
003 34 x↔y 045 42,21, 4 f LBL 4 087 45,30, .4 RCL − . 4
004 44 .1 STO . 1 046 45 .0 RCL . 0 088 45,10, .2 RCL ÷ . 2
005 42,21, 1 f LBL 1 047 45,30, .1 RCL − . 1 089 45 .4 RCL . 4
006 43, 4, 1 g SF 1 048 44 .2 STO . 2 090 34 x↔y
007 45 .0 RCL . 0 049 45 .1 RCL . 1 091 10 ÷
008 32 15 GSB E 050 32 15 GSB E 092 44,30, .1 STO − . 1
009 44 .3 STO . 3 051 44 .4 STO . 4 093 45 .4 RCL . 4
010 45 .1 RCL . 1 052 45 .1 RCL . 1 094 43 34 g RND
011 45,40, .0 RCL + . 0 053 34 x↔y 095 43,30, 0 g TEST x≠0
012 2 2 054 45,20, .2 RCL × . 2 096 22 5 GTO 5
013 10 ÷ 055 45 .3 RCL . 3 097 45 .1 RCL . 1
014 44 .2 STO . 2 056 45,30, .4 RCL − . 4 098 43 32 g RTN
015 32 15 GSB E 057 43 20 g x=0 099 42,21,14 f LBL D
016 44 .4 STO . 4 058 42 36 f RAN # 100 45 .1 RCL . 1
017 45,20, .3 RCL × . 3 059 10 ÷ 101 43 32 g RTN
018 43,30, 2 g TEST x<0 060 30 − 102 42,21, 6 f LBL 6
019 32 3 GSB 3 061 45 .1 RCL . 1 103 26 EEX
020 43, 6, 1 g F? 1 062 44 .0 STO . 0 104 16 CHS
021 32 2 GSB 2 063 34 x↔y 105 7 7
022 45 .4 RCL . 4 064 44 .1 STO . 1 106 40 +
023 43 34 g RND 065 45 .4 RCL . 4 107 43 32 g RTN
024 43,30, 0 g TEST x≠0 066 44 .3 STO . 3 108 42,21,15 f LBL E
025 22 1 GTO 1 067 43 34 g RND 109 36 ENTER
026 45 .2 RCL . 2 068 43,30, 0 g TEST x≠0 110 36 ENTER
027 43 32 g RTN 069 22 4 GTO 4 111 36 ENTER
028 42,21, 2 f LBL 2 070 45 .1 RCL . 1 112 6 6
029 45 .2 RCL . 2 071 43 32 g RTN 113 30 −
030 44 .0 STO . 0 072 42,21,13 f LBL C 114 20 ×
031 43 32 g RTN 073 44 .1 STO . 1 115 8 8
032 42,21, 3 f LBL 3 074 42,21, 5 f LBL 5 116 40 +
033 45 .2 RCL . 2 075 45 .1 RCL . 1 117 20 ×
034 44 .1 STO . 1 076 32 15 GSB E 118 4 4
035 43, 5, 1 g CF 1 077 44 .4 STO . 4 119 20 ×
036 43 32 g RTN 078 45 .1 RCL . 1 120 7 7
037 42,21,12 f LBL B 079 43 20 g x=0 121 48 .
038 44 .0 STO . 0 080 12 eˣ 122 5 5
039 34 x↔y 081 26 EEX 123 30 −
040 43,30, 5 g TEST x=y 082 4 4 124 43 32 g RTN
041 32 6 GSB 6 083 10 ÷