(12C) Newton's Method

[Gamo]

This program is an adaptation of the HP 25 Solver program from our Member name Dwight Sturrock under the topic (12C) Numeric Solver
I make update to this program from Dieter advice as well.

I have been looking for the Solver Program for HP-12C for a long time, I didn't notice before that Dwight post this since 2013
This is a very good all around Solver for HP-12C and this one work very well.

Newton's Method Solver work extreamly fast on HP-12C+ and I'm very satisfy with this program.

The program uses registers 0, 1, 2, 3 and 4.

>Put f(x) program to line 41
>End with GTO 19

Procedure:
1. Store initial guess to R1 (STO 1)
2. Store tolerance to R2 (STO 2)

Program:

Code:


CLx
STO 0
RCL 1
GTO 41
Rv
STO 4
1
STO 0
RCL 1
RCL 1
X=0
e^x
EXX
5
/
STO 3
+
GTO 41
X=0
GTO 39
RCL 0
X=0
GTO 05
Rv
RCL 4
/
1
-
1/x
RCL 3
x
STO-1
ENTER
x
SQRT
RCL 2
X<=Y
GTO 01
RCL 1
R/S
Line 41 (Start Program)
.
.
.
GTO 19 (End)

Example:
Solve for X^X=Y (X to power of X for the given Y)

X^X=1000

Formula use to solve this equation: LN(X)*X – LN(Y) = 0

Press [g] [GTO] 40 and Switch to Program Mode: Program start at Line 41
LN
LSTx
x
RCL 5
LN
-
GTO 19 (End)

Setup equation:
1. Initial Guess for 4 > STO 1
2. Tolerance > EEX CHS 6 > STO 2
3 1000 > STO 5 (This number can be change to any number you like to solve)
4. R/S to solve. Result 4.555535705

RCL 1 is the store answer of the problem.
RCL 4 to check for accuracy. [0.000000000] is the best result with ± estimation.

Gamo

[Dieter]

(03-24-2018 01:10 PM)Gamo Wrote:  This is a very good all around Solver for HP-12C and this one work very well.

...if you remove line 33. ;-)

The sequence 2 ENTER x SQRT does not make sense, this always yields 2. Remove the line with the "2" and the program should work as intended (compare the absolute value of the last correction term with the tolerance in R2). This also means that all references to line 42 should be changed to line 41.

There is an error in the instructions for the example: the tolerance has to be stored in R2. So it's EEX 6 CHS STO 2.
I think you should also mention that x^x=y here is transformed to x · ln x – ln y = 0, which is how the function is coded.

The program may have an error in that it stops the iteration not only if f(x)=0 but also if f(x+h)=0. OK. this can be fixed easily.
And maybe the final step should better be GTO 00 instead of R/S so that the program can be restarted with a simple R/S.
Here is my attempt at a corrected version:

Code:

01 STO 1
02 1
03 STO 0
04 RCL 1
05 GTO 42
06 RCL 0
07 x=0?
08 GTO 25
09 Rv
10 STO 4
11 x=0?
12 GTO 40
13 0
14 STO 0
15 RCL 1
16 x=0?
17 e^x
18 EEX
19 4
20 /
21 STO 3
22 RCL 1
23 +
24 GTO 42
25 Rv
26 RCL 4
27 -
28 RCL 3
29 /
30 RCL 4
31 x<>y
32 /
33 STO-1
34 ENTER
35 x
36 SQRT
37 RCL 2
38 x<=y?
39 GTO 02
40 RCL 1
41 GTO 00
42 .       start f(x) here
   .
   .
   .
nn GTO 06  end with this line

- store the tolerance in R2
- enter a guess, start with R/S

By the way, I have also tried a Regula Falsi / Illinois version for the 12C. This is slightly more effort as the 12C does not support subroutines and there are three different cases (as opposed to two for Newton) that have to be handled by the code that accepts the return value of the function call. But it is possible – in 70 steps.

Dieter

[Gamo]

I took off the 2 and rearrange GTO lines and now the program is updated.

Retry the same X^X=Y equation.
Now work flawlessly for this particular problem.

I like to see the Regula Falsi Method if it possible to make it to 12C

HP-12C should include SOLVE like the higher model of the HP Financial Calculator.

Gamo

[Dieter]

(03-25-2018 03:05 AM)Gamo Wrote:  I like to see the Regula Falsi Method if it possible to make it to 12C

[x] Done. (12C) Solve f(x)=0 with modified Regula Falsi method.

(03-25-2018 03:05 AM)Gamo Wrote:  HP-12C should include SOLVE like the higher model of the HP Financial Calculator.

I don't think it has to. Look, the 12C is not a universal tool, it is a specialized calculator for financial applications. That's why it has a very limited set of scientific functions: square root, power, natural log and e^x – and that's it. There is not even an x² key, let alone sin, cos, tan, common log and antilog and others. With the limited ressources of 1981 you could either have a financial calculator with a quite complete function set for such applications, or you could get a more universal scientific device like the 11C. Looking at 37 years of success this seems to have been a quite good decision. ;-)

Dieter

[Dieter]

(03-24-2018 01:31 PM)Dieter Wrote:  Here is my attempt at a corrected version:
...
- store the tolerance in R2
- enter a guess, start with R/S

Gamo, in the Regula Falsi thread you said that using the Newton solver is a bit cumbersome as – in the original version you posted – both the tolerance and an initial guess have to be prestored. The latter can be fixed easily (cf. the version I posted above), but also the exit condition can be simplified so that no user interaction is required.

The above versions use a fixed tolerance that is prestored by the user. On the other hand my Regula Falsi program quits as soon at the last two approximation agree to display precision. Here is another option: the following Newton solver iterates until the relative difference between the last two approximations drops below 5 E–9. This should provide 8 correct digits, but since the Newton method converges quadratically, in most cases the result should be close to exact.

Here is the code:

Code:

01 STO 1
02 1
03 STO 0
04 RCL 1
05 GTO 47
06 RCL 0
07 x=0?
08 GTO 25
09 Rv
10 STO 4
11 x=0?
12 GTO 45
13 0
14 STO 0
15 RCL 1
16 x=0?
17 e^x
18 EEX
19 4
20 /
21 STO 3
22 RCL 1
23 +
24 GTO 47
25 Rv
26 RCL 4
27 -
28 RCL 3
29 /
30 RCL 4
31 x<>y
32 /
33 STO-1
34 RCL 3
35 /
36 2
37 EEX
38 5
39 +
40 LstX
41 -
42 x=0?
43 GTO 45
44 GTO 02
45 RCL 1
46 GTO 00
47 .       start f(x) here
   .
   .
nn GTO 06  end f(x) with this line

Try it with your example function 3^x – x^3 = 0.

Code:

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

47 3
48 x<>y
49 y^x
50 LstX
51 3
52 y^x
53 -
54 GTO 06

f [P/R]

Now try this with an initial guess of 2:

Code:

 2  [R/S]  =>  2,478052683

You can always add one more iteration by pressing [R/S] again.
In the above case the result remains the same, with 10 digit precision f(2,478052683) yields exactly 0.

Now start with a guess of 2,5 which is even closer to the solution:

Code:

2,5 [R/S]  =>  2,478052686

The result differs in the last digit.
Now press [R/S] again and see what you get:

Code:

    [R/S]  =>  2,478052680
    [R/S]  =>  2,478052686
    [R/S]  =>  2,478052680
    [R/S]  =>  2,478052686
     ...

See? The iteration oscillates between 2,478052680 and ...86. That's why it is not a good idea to continue until the last two approximations match exactly (which would have caused an infinite loop here). The method in the above program stops if they agree in at least 8 digits, which usually means that it will not get more accurate anyway. Nevertheless you can always press [R/S] again and see if the result improves.

This way an infinite loops is avoided (well, at least in most cases). Try the function with a guess of 6 and you see that it would otherwise oscillate endlessly between 3,000000002 and 2,999999998.

In detail: line 34-35 divides the last correction term by h, and since h = x / 1E+4 this yields the current relative error times 1E+4. Dividing by h instead of x automatically handles the x=0 case for which h is set to 1E–4. Then 2E+5 is added and subtracted again, which shifts a sufficiently small error beyond the 10th digit so that the result is zero if the error is small enough. Selecting a constant with more or less digits (2E+6 or 2E+4 etc.) changes the error threshold, but I think 2E+5 is just right for a 10-digit calculator while h = x / 1E+4.

Dieter

[Gamo]

The new updated user friendly version of the Newton's Method now work even better than the previous version just start your guess and go R/S !!!

I try the same X^X=Y

Store Y in Register 5 (STO 5)

X^X=1000

1000 [STO 5]
2 R/S > 4,555535706
3 R/S > 4,555535705
4 R/S > 4,555535705

150.25 [STO 5]
2 R/S > 3.773964520
3 R/S > 3.773964519
4 R/S > 3.773964520
5 R/S > 3.773964519
6 R/S > 3.773964519

123456789 [STO 5]
4 R/S > 8.640026440
5 R/S > 8.640026441
6 R/S > 8.640026437
7 R/S > 8.640026440
8 R/S > 8.640026440
9 R/S > 8.640026439
10 R/S > 8.640026440

-----------------------------------------------------------------------------------------------
Solve: LN(X) + 3X = 10.8074

Change to: LN(X) + 3X - 10.8074 = 0

Program Code:
LN
LSTx
3
x
+
RCL 5
-

Solve Equation by first store 10.8074 to Register 5 [STO 5]

4 > R/S > 3.213360869 check tolerant gave exact 0.000000000
3 > R/S > 3.213360870 check tolerant gave -0.000000020


So far this version is very accurate and no need to manually change tolerant.

Gamo

[Carsen]

I like it too. Gamo said it all. Key in your guess and press [R/S]. Well done

[Dieter]

(03-29-2018 05:11 AM)Carsen Wrote:  I like it too. Gamo said it all. Key in your guess and press [R/S]. Well done

Thank you. But there is always a way to squeeze out a few bytes. Here is an optimized version with some changes. These include:

• Instead of R0, R1, R3 and R4 the program now uses four consecutive registers, R0...R3.
• The f(x)=0 test was removed. It is not required since in this case the correction term is zero as well so that the iteration will exit anyway. This way one more call of f(x+h) may be required, but hey: who cares if you can save two steps? ;-)
• The exit condition now is set to a relative error < 1E–8. This is slightly higher than before but it makes sure that even in the worst case the program finishes if the last change in x affected only the last one ore two digits. You can always add one more iteration with a simple [R/S]. Also the test of the last correction term (is it below the threshold or not?) has been changed to a slightly more effective method.
• Here and there the code was optimized to save the one or other step.
  

All this lead to the version below which now is down to 40 steps (plus the f(x) code).
Here is a commented listing:

Code:

01 STO 1   // store initial guess
02 1
03 STO 0   // flag = 1
04 RCL 1   // call f(x)
05 GTO 41
06 RCL 0   // function calls return here
07 x=0?    // is this a return from f(x+h)?
08 GTO 22  // then continue with Newton correction
09 STO-0   // flag = 0
10 Rv
11 STO 2   // save f(x)
12 RCL 1   // calculate h
13 x=0?
14 e^x     // set h = x/10000, or 1/10000 if x=0
15 EEX
16 4
17 /
18 STO 3   // store h
19 RCL 1
20 +
21 GTO 41  // call f(x+h)
22 Rv
23 RCL 2   // calculate correction term
24 -       // according to Newton's method
25 RCL 3
26 /       // f'(x)
27 STO/2
28 RCL 2   // correction term = f(x)/f'(x)
29 STO-1   // adjust x
30 RCL 3
31 /       // = relative correction * 10000
32 EEX
33 4
34 x       // = relative correction * 1E+8
35 INTG
36 x=0?    // was |relative correction| < 1E-8 ?
37 GTO 39  // then exit
38 GTO 02  // else do another loop
39 RCL 1   // return x
40 GTO 00  // and quit
41 .       // start f(x) here
   .
   .
nn GTO 06  // end f(x) with this line

As usual, comments and suggestions are welcome.

Dieter

[Gene]

Great job as always.

I should turn you loose on my lifetime HP 41 Yahtzee game code.

First wrote it back in 1985 and always looking to make it shorter, faster, better.

;-)

[Gamo]

The latest update to only 40 line Newton's Method program work very good.

Find a root of a polynomial of f(x) = X^4 - 4X^3 + 8X^2 + 20X - 65

Program for f(x)
[ENTER] [ENTER] [ENTER] 4 [-] [x] 8 [+] [x] 20 [+] [x] 65 -

5 > R/S > 2.236067978
4 > R/S > 2.236067978
3 > R/S > 2.236067977
2 > R/S > 2.236067977
1 > R/S > 2.236067978

Answer is 2.236067977

Gamo

Can I check for accuracy on Register 2

[Dieter]

(03-30-2018 01:15 PM)Gamo Wrote:  5 > R/S > 2.236067978
4 > R/S > 2.236067978
3 > R/S > 2.236067977
2 > R/S > 2.236067977
1 > R/S > 2.236067978

Answer is 2.236067977

The true result is right in the middle. Even with 13-digit precision the answer is 2,236067977500. See below.

(03-30-2018 01:15 PM)Gamo Wrote:  Can I check for accuracy on Register 2

In a way, yes. But let's be precise for a moment now. What is accuracy? How do you define it? The term may relate to two different things.

On the one hand you may refer to the difference between the returned result and the actual true answer. For instance, if the true result is 2 and the program returns 1,999999999 this is accurate to 9 places, or it has an error of 1 ULP.

On the other hand you may take a look at f(x) for the returned value. Does it return exactly zero, or 0,000000001 or something else?

Both definitions may come to different conclusions. For instance, there are cases where a whole range of results yields f(x) = exactly 0. Remember the 3^x–x^3 case? Here anything between 2,478052679 and 2,478052683 will return exactly zero on a correctly working 10-digit calculator. The true result is 2,4780526802883...

On the other hand even a function result that is not zero does not mean that the result has a slight error. It may be exact to all 10 digits. Consider the function x^2–2=0. The exact 10-digit result here is 1,414213562, but f(1,414213562) yields –0,000000001. The next higher answer 1,414213563 returns +0,000000002. So there simply is no 10-digit result that returns exactly 0.

You can also check by how much the final value of x has been corrected. This can be a measure for its accuracy. But again, even if the last two approximations agree and the correction term is zero, this does not neccessarily mean that the result is dead on.

Back to the programs now. You said that in the original version of the program (cf. your initial post) RCL 4 could be used to check the accuracy of the result. In fact this is not possible. The original version calculates f(x) and stores this in R4. But then a "better x" is calculated and finally displayed. So R4 holds the function value of the second-to-last x, and not that of the finally returned x (!).

In the last version of my program you may recall R2 to see the last correction term. If R2 is zero, f(x) returns zero as well. Nonzero values mean that f(x) of the previous (!) approximatino was not zero, but x should be at least as accurate as this. So if x is for example 3,5 and RCL 2 returns 3E–9 you know that the final correction was 3 units in the last digit. The program terminates if this update only affected the last one or two digits. In most cases this means that the new result is as good as it gets with 10-digit precision. Due to the quadratic convergence of Newton's method, in a perfect world even 15 digits should be fine.

There is an easy way to get another impression of the accuracy: press [R/S] again for an update of x. In some cases this may show how x oscillates between two values. Remember the example where the last two digits changed between ...80 and ...86? Here the calculator simply cannot do better, for instance because of limited accuracy in the function itself, and the last digit remains uncertain.

He is an example for the polynomial in your last post.

5 [R/S] => 2,236067978
RCL 2 => 4,938E–10

So f(x) is not exactly zero, but the last correction was less than 5 units in the 11th digit, thus not affecting the 10-digit result. In fact the true result is right in the middle between 2,236067977 and ...78. The program tried to correct it to 2,2360679775 but due to 10-digit precision it could not do so.

Dieter

[Gamo]

New equation: X = (e^-X) solve for X

Change equation to (e^-X) - X = 0

f(x) Program:
STO 5
CHS
e^x
RCL 5
-

3 R/S > 0,567143290
2 R/S > 0,567143290
1 R/S > 0,567143290

X = 0,567143290

e^-0,567143290 = 0,567143290

Gamo

[Dieter]

(03-31-2018 01:15 PM)Gamo Wrote:  f(x) Program:
STO 5
CHS
e^x
RCL 5
-

It is one of the advantages of RPN that you rarely need data registers for such cases.

ENTER
CHS
e^x
X<>Y
–

Or, even simpler:

CHS
e^x
LstX
+

Dieter