(12C) Solve f(x)=0 with modified Regula Falsi method

[Dieter]

Yesterday Gamo posted a Newton solver for the 12C, and I wondered if a Regula Falsi method was also possible. Since the 12C does not feature subroutine calls – which are usually required to evaluate f(x) – such a program is a bit tricky.

The following program implements a modified Regula Falsi method, here the Illinois variant. This is a well known algorithm that has been used since the days of the HP65. As the 12C does not support subroutines the following approach was used:

The evaluation of f(x) is required for the two initial guesses a and b, i.e. for f(a) and f(b), and for the function value of the new approximation c, i.e. f(c). After the function value is calculated, the program in all three cases jumps back to a common return point (line 09). Here a flag decides how to continue: the flag is a number in R0 which can be negative, positive or zero. The first case encodes the f(a) call, so that the program continues with storing f(a) and calculates f(b) next. The second case is set for f(b) so that the program stores this value and proceeds with the start of the iteration loop. Finally the flag is set to zero, in this case the program knows that it returns from an f(c) call, so it continues with the respective part of the Regula Falsi algorithm.

Here is the program. Please note: this is experimental code that has not seen much testing, so use it at your own risk, and please report any errors you find.

Code:

01 STO 2   store b
02 x<>y
03 STO 1   store a
04 1
05 CHS
06 STO 0   flag := -1
07 RCL 1
08 GTO 74  calculate f(a)
09 RCL 0   all f(x) calls return here
10 x=0?    flag = 0?
11 GTO 50  this was an f(c) call, continue with iteration loop
12 0
13 x<=y?   flag > 0?
14 GTO 22  this was an f(b) call, continue with storing f(b)
15 Rv
16 Rv      so flag is < 0
17 STO 4   this was an f(a) call, continue with storing f(a)
18 1
19 STO 0   flag := +1
20 RCL 2
21 GTO 74  calculate f(b)
22 STO 0   flag := 0
23 Rv
24 Rv
25 STO 5   store f(b)
26 RCL 4
27 x       f(a) x f(b)
28 CHS     if no sign change between a and b
29 SQRT    stop with error message
30 RCL 2
31 RCL 2
32 RCL 1   step 30 - 40
33 -       calculate the
34 RCL 5   new approximation c
35 RCL 4
36 -
37 /
38 RCL 5
39 x
40 -
41 STO 3   store c
42 RND     round c to display precision
43 RCL 2
44 RND     round b to display precision
45 -
46 x=0?    do both agree?
47 GTO 71  then exit
48 RCL 3
49 GTO 74  calculate f(c)
50 Rv      return from f(c) call, f(c) is in Y
51 x=0?    f(c)=0?
52 GTO 71  then exit
53 RCL 5
54 x<>y
55 x       f(b) x f(c), here f(c) is saved in LstX
56 0
57 x<=y?   no sign change between b and c?
58 GTO 64  then apply Illinois correction
59 RCL 2
60 STO 1   else a := b
61 RCL 5   and f(a) := f(b)
62 STO 4
63 GTO 66
64 2       Illinois adjustment:
65 STO/4   f(a) := 1/2 f(a)
66 RCL 3
67 STO 2   b := c
68 LstX
69 STO 5   f(b) := f(c)
70 GTO 30  and do another iteration
71 RCL 2   exit routine: return previous
72 RCL 3   and final approximation
73 GTO 00  and quit
74 ...     start f(x) here
nn GTO 09  end f(x) with this line

The program iterates until the last two approximation agree to display precision. So set FIX 6 if a result with six digit accuracy is desired. Please note that the 12C handles the display different from other HP calculators, for instance in cases where a number is so close to zero that other calculators would switch to scientific notation. So if the program does not stop, press [R/S] and start over with a different display setting.

Usage:

Enter program
Enter code for f(x) as usual, starting at line 74
End f(x) code with GTO 09

Example function: e^x – 3x = 0

Code:

g [GTO] 73
f [P/R]       73-43,33 00

74 e^x
75 LstX
76 3
77 x
78 -
79 GTO 09

f [P/R]

Set the desired accuracy via FIX mode.
For instance, for six decimals set FIX 6.
Enter 0 and 1 as initial guesses.

Code:

f 6
0 [ENTER] 1  [R/S]
=> 0,619061

Now also find the root between 1 and 2.

Code:

1 [ENTER] 2  [R/S]
=> 1,512135

Note that f(a) and f(b) should have opposite signs. If this condition is not met the program generates an "Error 0" message. If you want to continue anyway, clear the error display and press 0 [R/S] (any other non-negative value will do as well).

The last two approximations are returned in X and Y so that these can be viewed with [X<>Y].

Code:

f 9
[X<>Y]   1,512134546
[X<>Y]   1,512134552

So the last approximation is correct in 9 digits, probably even in all 10.
Another [R/S] restarts the iteration with these two values as initial guesses, and in fact this yields the same result 1,512134552 in both X and Y.

Addendum: the mentioned Wikipedia article includes a different way of calculating the new approximation c which is said to have some numerical advantages. If you want to try this in the above program simply replace line 30...40 with the code below. Since the line count remains the same the is no need to adjust any GTOs.

Code:

.. ...
30 RCL 1
31 RCL 5
32 x
33 RCL 2
34 RCL 4
35 x
36 -
37 RCL 5
38 RCL 4
39 -
40 /
.. ...

Dieter