(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

[Gamo]

I test this to find a root of a polynomial of this equation.

f(x)=X^4 - 4X^3 + 8X^2 + 20X - 65

Program Method:

ENTER
ENTER
ENTER
4
-
x
8
+
x
20
+
x
65
-

FIX 9
Initial Guess used: 0 ENTER 3 R/S > instant answer > 2.236067977

My favorite X^X=Y also work flawlessly.

For this equation somehow not working.

Solve X^3 = 3^X
Rewrite equation to X^3 - 3^X = 0

Program Method:

ENTER
ENTER
3
Y^X
X<>Y
3
X<>Y
Y^X
-

The answer is 2.478052679

I try initial guess 1 ENTER 3 R/S > return 3 press R/S again > Error 0

Did I do something wrong?

I try on the Newton's Method and work OK

Gamo

[Dieter]

(03-26-2018 06:26 AM)Gamo Wrote:  Solve X^3 = 3^X
Rewrite equation to X^3 - 3^X = 0

...

The answer is 2.478052679

One answer is 2,47805. But there are two.

(03-26-2018 06:26 AM)Gamo Wrote:  I try initial guess 1 ENTER 3 R/S > return 3

Sure. X=3 is a solution: 3^3 = 3^3.

(03-26-2018 06:26 AM)Gamo Wrote:  press R/S again > Error 0

Why do you press R/S again? What do you want to do this way?

If you press R/S again you restart the solver with two identical guesses 3 and 3. So the calculated secant is a horizontal line that will never cross the x-axis and a division by zero will occur. That's why the two guesses must be different.

However, the program can be changed so that it does accept two identical guesses if these are already a root of f(x). See below.

(03-26-2018 06:26 AM)Gamo Wrote:  Did I do something wrong?

You simply didn't notice that 3 is one of the two solutions. ;-)

- If you enter 1 and 3 as initial guesses the program correctly returns X=3.
- Of course the program can also calculate the other root: simply start with different guesses, for instance 1 and 2,5 => 2,478052679.

BTW, the exact root is 2,4780526802883..., but with 10 digit precision everything between 2,478052678 and ...683 evaluates to f(x)=0. If you start with 1 and 2,8 or 1 and 2,9 you will get these results.

Edit: here is a modified version that works even with two identical guesses if these are a root of f(x).
f(x) now starts at line 76 and ends with GTO 10.

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

Additional feature: in the f(x) code, x now can be addressed with "RCL 3".
So your example can be coded like this:

3
Y^X
3
RCL 3
Y^X
–

Dieter

[Gamo]

The modify version also used Register 3 in the program, is that OK to use in the
f(x) program with RCL 3 ?

I try both Newton's and Regula Falsi methods on Modern 12C+ and on Official 12C App on Android OS. Both devices run both methods equally fast of about one second.

For Newton's Method I find that the program is shorter is only about 41 lines of code versus 75 lines of the Regula Falsi Method program.

Newton's Method is more cumbersome to input all the necessary data for regular user such as Store Initial Guess, Store tolerance

Regula Falsi Method is very easy to input data, user only put two initial guesses then R/S that's all.

I like to know how long it take to compute with the older version of the HP-12C For comparison Newton's Method run slow on HP-11C but with Regula Falsi this method run faster and very reliable it can take almost any equation with ease.

Gamo

[Dieter]

(03-26-2018 01:04 PM)Gamo Wrote:  The modify version also used Register 3 in the program, is that OK to use in the
f(x) program with RCL 3 ?

Yes, x is always stored in R3 before the function is called. So you can use "RCL 3" for every occurence of x. And of course x it is also passed in the X-register. ;-)

(03-26-2018 01:04 PM)Gamo Wrote:  Newton's Method is more cumbersome to input all the necessary data for regular user such as Store Initial Guess, Store tolerance

Take a look at the Newton version I posted in the respective thread. It's the same 41 lins but the user only has to store the tolerance. Of course this program can also be modified so that it uses the same method as the Regula Falsi program: compare the last two rounded approximations. This way there is no need to enter a tolerance. In other words: all this can also be done with the Newton program.

(03-26-2018 01:04 PM)Gamo Wrote:  I like to know how long it take to compute with the older version of the HP-12C For comparison Newton's Method run slow on HP-11C but with Regula Falsi this method run faster and very reliable it can take almost any equation with ease.

I can't say anything about execution times on a classic hardware 12C but beware: for every method there are cases where the algorithm will be slow or even fail. There is no universal method that will always return an exact result reliably and fast.

Dieter

[Carsen]

This is part of the reason why I love the HP-12C so much. Programming. This program really made my day!

[Dieter]

(03-26-2018 04:37 PM)Dieter Wrote:  

(03-26-2018 01:04 PM)Gamo Wrote:  Newton's Method is more cumbersome to input all the necessary data for regular user such as Store Initial Guess, Store tolerance

Take a look at the Newton version I posted in the respective thread.

Gamo and Carsen: There is a new version of the Newton solver with a modified exit condition I just posted in the respective thread. You do not have to prestore anything and the danger of an infinite loop (which may happen sooner as one might expect) is (mostly) avoided – take a look at the examples. I like this approach better than the previous ones.

Dieter