Small Solver Program

[Gamo]

Recently I try to make a small but versatile solver program specifically
for HP Voyager Series calculator.

This solver solve for f(x) = 0 and use some type of counter for iteration.
Due to the iteration process this program will run slow on Original Calculator.
With modern 12C, 15C, and all other emulator this will run extremely fast.

This program is intended for HP-11C but this example will use
HP-15C emulator so that this solver program result can be compare
with the 15C build-in SOLVE function.

---------------------------------------------------------
To run program.

Enter 2 guesses closest to the root.

X1 [ENTER] X2 [A] display answer // Use this Solver

X1 [ENTER] X2 f [SOLVE] [E] display answer // Use 15C build-in SOLVE
---------------------------------------------------------
Example for X^X = Y // X to the power of X equal to Y

X on R1 // Register Storage 1 // For this Sover
Y on R5 // Register Storage 5

X on R2 // Register Storage 2 // For 15C SOLVE
Y on R5 // Register Storage 5
--------------------------------------------------------
Program:

Code:


LBL A   // This Solver
STO 0
Rv
STO 1
90   
STO I
CLx
LBL 1
RCL 0
/
CHS
RCL 1
+
STO 1
DSE I
GTO B
RCL 1
RTN
------------------------------
LBL B  // f(x) program for X^X = Y
RCL 1
LN
RCL 1
x
RCL 5
LN
-
GTO 1
-------------------------------
LBL E  // 15C SOLVE
STO 2 // f(x) program for X^X = Y
LN
RCL 2
x
RCL 5
LN
-
RTN

-------------------------------------------
Example: X^X = 1000

STORE 1000 to R5 // [STO] 5

This Solver

4 [ENTER] 6 f [A] display 4.555535704
-------------------------------------------------------
15C SOLVE

4 [ENTER] 6 f [SOLVE] [E] display 4.555535705
-------------------------------------------------------
Remark:
For this solver answer is in R1
For 15C SOLVE answer is in R2


Gamo

[Thomas Klemm]

(02-14-2019 05:25 AM)Gamo Wrote:  Example: X^X = 1000

You can rewrite this equation such that \(x\) is a fixed point of \(f(x)\), that is \(x=f(x)\):

\(x=\frac{3}{LOG(x)}\)

Then you can iterate the following program with a starting value, say \(x=4\):

Code:

001-  43  13 :    LOG
002-       3 :    3
003-      34 :    x<>Y
004-      10 :    ÷

Example:

4
R/S
4.982892143
R/S
4.301189432
R/S
4.734933901
R/S
4.442378438
R/S
4.632377942
R/S
4.505830646
R/S
4.588735607
R/S
4.533824357
R/S
4.569933524
R/S
4.546075273
R/S
4.561789745
R/S
4.551417837
R/S
4.558254212
R/S
4.553744141
R/S
4.556717747
R/S
4.554756405
R/S
4.556049741
R/S
4.555196752
R/S
4.555759258
R/S
4.555388285
R/S
4.555632930
R/S
4.555471588
R/S
4.555577990
R/S
4.555507820
R/S
4.555554095
R/S
4.555523578
R/S
4.555543703
R/S
4.555530431
R/S
4.555539183
R/S
4.555533412
R/S
4.555537217
R/S
4.555534708
…
4.555535704
R/S
4.555535706
R/S
4.555535704
R/S
4.555535706
…

From then on the result flips between these last two values.

Cheers
Thomas

[Thomas Klemm]

If you don't want to press R/S repeatedly:

Code:

001-  43  13 :    LOG
002-       3 :    3
003-      34 :    x<>Y
004-      10 :    ÷
005-  42  31 :    PSE
006-  43  40 :    x=0

[Albert Chan]

(02-14-2019 07:06 AM)Thomas Klemm Wrote:  \(x=\frac{3}{LOG(x)}\)

4
R/S
4.982892143
R/S
4.301189432...

It would be nice if we can temper the oscillation, or slow convergence.
Let x0 = 4, x1, x2 = the first two iterated values.

Rate = (x2-x1)/(x1-x0) ≈ (4.30 - 4.98) / (4.98 - 4) = -0.694

If the same trend continued, we expect final % = 1/(1-r) ~ 60%

x_∞ ≈ x0/(1-r) = (x1 - (x1-x0))/(1-r) = (x1 - r x0) / (1-r)

Use weighted fixed-point equation x = 0.6 * 3/log10(x) + 0.4 x

4.6
4.555924149
4.555537395
4.555535712
4.555535705 (converged)

Edit: compare with Newton's method, x = (ln(1000) + x) / (ln(x) + 1)
4
4.571001573
4.555546101
4.555535705 (converged)

[Thomas Klemm]

(02-15-2019 12:07 AM)Albert Chan Wrote:  It would be nice if we can temper the oscillation, or slow convergence.

We can also use Aitken's delta-squared process to accelerate the speed of convergence:

\(a_{n}=x_{n+2}-\frac {(x_{n+2}-x_{n+1})^{2}}{(x_{n+2}-x_{n+1})-(x_{n+1}-x_{n})}\)

This leads to the following program for the HP-11C:

Code:

001-42,21, 1 :    ▸LBL 01
002-      36 :     ENTER
003-    32 0 :     GSB 00
004-      30 :     -
005-  43  36 :     LSTx
006-      36 :     ENTER
007-    32 0 :     GSB 00
008-      30 :     -
009-  43  36 :     LSTx
010-      34 :     x<>y
011-  43  33 :     R↑
012-      34 :     x<>y
013-      30 :     -
014-  43  36 :     LSTx
015-  43  11 :     x²
016-      34 :     x<>y
017-      10 :     ÷
018-      30 :     -
019-  43  32 :     RTN
020-42,21, 0 :    ▸LBL 00
021-  43  12 :     LN
022-    45 0 :     RCL 00
023-      34 :     x<>y
024-      10 :     ÷
025-  43  32 :     RTN

Example:

Solve: \(x^x=1000\)

1000 LN
STO 0

4.6
GSB 1
4.555665779
R/S
4.555535706
R/S
4.555535705
R/S
Error 0

Feel free to add a check for 0 after line 13 to prevent the error and loop back to label 1 instead of returning in line 19.

Cheers
Thomas

[Albert Chan]

(02-15-2019 06:26 PM)Thomas Klemm Wrote:  

(02-15-2019 12:07 AM)Albert Chan Wrote:  It would be nice if we can temper the oscillation, or slow convergence.

We can also use Aitken's delta-squared process to accelerate the speed of convergence:

\(a_{n}=x_{n+2}-\frac {(x_{n+2}-x_{n+1})^{2}}{(x_{n+2}-x_{n+1})-(x_{n+1}-x_{n})}\)

Amazingly, my rate formula is same as Aitken extrapolation formula !
Assuming we have 3 values, x0,x1,x2 and tried to guess where it should end up.

My rate analysis: r = (x2-x1)/(x1-x0) = Δx1 / Δx0

We need this for the proof:
(Δx0)²
= ((Δx0 - Δx1) + Δx1)²
= (Δx0 - Δx1)² + 2 * Δx1 (Δx0 - Δx1) + (Δx1)²

If same trend continue, where it will ends up
= x0 + Δx0 * 1/(1-r)
= x0 + (Δx0)² / (Δx0 - Δx1)
= x0 + (Δx0 - Δx1) + 2 * Δx1 + (Δx1)² / (Δx0 - Δx1)
= x0 + (x1-x0) + (x1-x2) + 2*(x2-x1) − (Δx1)² / (Δx1 - Δx0)
= x2 − (Δx1)² / (Δx1 - Δx0)

[Thomas Klemm]

(02-15-2019 09:16 PM)Albert Chan Wrote:  Amazingly, my rate formula is same as Aitken extrapolation formula !

That's not really a surprise, is it?
It's essentially the same idea: geometric series

\(1\,+\,r\,+\,r^{2}\,+\,r^{3}\,+\,\cdots \;=\;{\frac {1}{1-r}}\)

[Csaba Tizedes]

(02-14-2019 05:25 AM)Gamo Wrote:  a small but versatile solver program specifically for HP Voyager Series calculator.
Example for X^X = Y // X to the power of X equal to Y

Yes, the classic problem: which number you can power to itself to get the largest number on your calculator: x^x=10^100. Do LOG() of both sides: x×LOG(x)=100, then do the iteration: x_0:=50, then x_i+1:=100/log(x_i). The simple fixed point iteration and converges.

By the way, if you interested, here is my SOLVE(i) for HP-12C, which can solve any equation for ANY variable like on 15C, without rewriting the equation. Looks like indirect addressing without indirect addressing. 28 steps only.

Yes, my English is terrible, but more convenient than read the full material in Hungarian.
At the end of pdf you can find lots of applications which solved with this little solver (like ODE solve with EULER+SOLVE(i) on 12C (!!!), pressure vessel pneumatic conveying pipe sizing program, humidity and dew point solving, outer temperature calculation of a concrete silo wall or comparison of present values of educated and not educated people salaries with the little SOLVE combined with 12C's cash flow calculation capabilities - I guess this is shows clearly the power of this little routine possibilities and why I say everywhere, the SOLVER is MUST to implement ALL current calculators, I hope somebody read and understand what I want to say).

The small secant method in this paper for CASIO fx-50F meantime reduced to 13 steps, I guess the smallest solver for these blind machines. You can find attached.

Csaba

[Thomas Klemm]

(02-16-2019 12:24 PM)Csaba Tizedes Wrote:  By the way, if you interested, here is my SOLVE(i) for HP-12C

Steps of secant method from 03 to 18:

Code:

03-   45 2 :    RCL 2
04-   45 0 :    RCL 0
05-     30 :    -
06-  43 36 :    LSTx
07-   44 2 :    STO 2
08-     35 :    CLx
09-   45 3 :    RCL 3
10-   45 1 :    RCL 1
11-     30 :    -
12-  43 36 :    LSTx
13-   44 3 :    STO 3
14-     33 :    R↓
15-     10 :    ÷
16-   45 1 :    RCL 1
17-     20 :    ×
18-44 30 0 :    STO+ 0

This can be shortened to:

Code:

03-   45 2 :    RCL 2
04-   45 0 :    RCL 0
05-   44 2 :    STO 2
06-     30 :    -
07-   45 3 :    RCL 3
08-   45 1 :    RCL 1
09-   44 3 :    STO 3
10-     30 :    -
11-     10 :    ÷
12-   45 1 :    RCL 1
13-     20 :    ×
14-44 30 0 :    STO+ 0

Don't hesitate to publish your programs in this forum. You can use Google to translate to English.

Cheers
Thomas

[Csaba Tizedes]

(02-16-2019 01:42 PM)Thomas Klemm Wrote:  Steps of secant method from 03 to 18 can be shortened to.

Thanks, this was eons before, I am sure this is only the first iteration and required to fine tuning.

(02-16-2019 01:42 PM)Thomas Klemm Wrote:  Don't hesitate to publish your programs in this forum. You can use Google to translate to English.

OK, sometimes I feel myself as Don Quixote... Frankly speaking, who want to discuss about calculators and calculator programming now? 15 years before I wrote something useful, today I am only (an angry) user...

Csaba

[Gamo]

I just fine tune this solver program for HP-11C

Added the iterations limit to 28 loops count
Added the condition test to end program when root is found.

If iteration is over 28 loops count this indicate that user must use a closer guess.

-------------------------------------------

Procedure:

Use R1 for unknown X in your equation. End your formula with [RTN]

X1 ENTER X2 [A] display Answer

Answer is in R1

------------------------------------------------------------------------
Program:

Code:


LBL A
STO 0
Rv
STO 1
28
STO I  // Store Loop Count Limit
CLx
---------------------
LBL 0
GSB B
GSB 1
STO 1
GSB B
GSB 1
RCL 1
X=Y  // End if Root is found
GTO 2
DSE  // Loop Limit Counter
GTO 0
CLx
FIX 9  // 0.000000000 indicate that Maximum Loops is use up.
PSE
PSE
FIX 4 
GTO 2
---------------------
LBL 1  // Solver Routine
RCL 0
÷
CHS
RCL 1
+
RTN
--------------------
LBL 2
RCL 1  // Answer 
RTN
--------------------
LBL B
.  // f(x) equation start here
.
.  // X = Register 1 [R1]
.
.
RTN

Gamo

[Thomas Klemm]

(02-17-2019 02:57 AM)Gamo Wrote:  I just fine tune this solver program for HP-11C

Your program doesn't work for this simple example: \(x^2+x=x(x+1)=0\)
I've entered the following program:

Code:

LBL B
1
RCL 1
+
LSTx
×
RTN

Example:

-2
ENTER
-0.5
A

9.999999 99

The reason is that the function value is divided by -0.5, making it larger and larger:

Code:

LBL 1  // Solver Routine
RCL 0
÷
CHS
RCL 1
+
RTN

---

Another thing is this in this loop:

Code:

LBL 0
GSB B
GSB 1
STO 1
GSB B
GSB 1
RCL 1
X=Y  // End if Root is found
GTO 2
DSE  // Loop Limit Counter
GTO 0

When you return to label 0, the same value in register 1 will be used as before.
This function evaluation can be avoided if you save the result of the previous evaluation.

---

You could move this block:

Code:

LBL 2
RCL 1  // Answer 
RTN

Over here and avoid the jump GTO 2:

Code:

CLx
FIX 9  // 0.000000000 indicate that Maximum Loops is use up.
PSE
PSE
FIX 4 
LBL 2
RCL 1  // Answer 
RTN

---

Take a look at Csaba's solution for the HP-12C using the secant method.

Cheers
Thomas

[Gamo]

Thanks Thomas Klemm

I try single guess and it work.

CHS 2 [A] display -1

If using two negative guesses and result is
Overflow then use Single guess instead.

Formula Routine is

RCL 1
1
+
RCL 1
X

Gamo

[Thomas Klemm]

(02-17-2019 02:33 PM)Gamo Wrote:  If using two negative guesses and result is Overflow then use Single guess instead.

Then try this function \(f(x)=4x(x+1)\):

Code:

LBL B
RCL 1
1
+
RCL 1
×
4
×
RTN

Examples:

-0.5
ENTER
0.5
A

9.999999 99

---

-0.5
ENTER
A

9.999999 99

---

0.5
ENTER
A

9.999999 99


Even initial guesses very close to the solutions lead to the same result:

-1.00001
ENTER
-0.99999
A

9.999999 99

---

-0.00001
ENTER
0.00001
A

9.999999 99

Your algorithm is fundamentally broken.

Cheers
Thomas

[Gamo]

Your're right that only work for certain situations.

Gamo

[Thomas Klemm]

If you compare your algorithm with Newton's method:

\(x_{n+1}=x_{n}-\frac {f(x_{n})}{f'(x_{n})}\)

you may notice that the number in register 0 should in fact be \(f'(x_{n})\) or at least a good approximation thereof.

And with this in mind indeed the solution can now be found:

\(f'(x) = \frac{d}{dx}4 x (x + 1) = 8 x + 4\)

\(f'(-1) = -4\)

\(f'(0) = 4\)

Examples:

-1.5
ENTER
-4
A

-1.00000

---

0.5
ENTER
4
A

0.00000

---

Thus your 2nd parameter is an estimate of the slope at the root.

---

And then you don't really have to call the function twice with the same value:

Code:

LBL A
STO 0
Rv
STO 1
28
STO I  // Store Loop Count Limit
LBL 0
GSB B
RCL 0
÷
STO - 1
X=0  // End if Root is found
GTO 1
DSE  // Loop Limit Counter
GTO 0
CLx
FIX 9  // 0.000000000 indicate that Maximum Loops is use up.
PSE
PSE
FIX 4
LBL 1
RCL 1  // Answer 
RTN
LBL B
.  // f(x) equation start here
.
.  // X = Register 1 [R1]
.
.
RTN

So the question remains how to approximate the derivation at the roots.
If you keep record of the previous function evaluation you can get an estimate of the slope using the secant method.

Cheers
Thomas

[Dieter]

(02-18-2019 05:20 AM)Thomas Klemm Wrote:  If you compare your algorithm with Newton's method:

\(x_{n+1}=x_{n}-\frac {f(x_{n})}{f'(x_{n})}\)

you may notice that the number in register 0 should in fact be \(f'(x_{n})\) or at least a good approximation thereof.

Yes. Gamo already posted this method some time ago in the General Software Library. I have mentioned several times that the first input it NOT a guess for the root but an estimate for the function's DERIVATIVE.

For more details take a look at this post: http://www.hpmuseum.org/forum/thread-123...#pid111680

Gamo, you really should correct and update the instructions. Here it still says:

Quote:Enter 2 guesses closest to the root.

But again: that's NOT true. The first input is NOT a guess for the root. You now have seen a variety of cases that prove it.

Dieter

[Thomas Klemm]

(02-18-2019 07:46 PM)Dieter Wrote:  Yes. Gamo already posted this method some time ago in the General Software Library. I have mentioned several times that the first input it NOT a guess for the root but an estimate for the function's DERIVATIVE.

You could think him stubborn.
But it made Csaba post his program which made me translate some of his papers from Hungarian to English and stumble totally unrelated upon one of his previous posts:
Quadratic fit without linear algebra (32SII) - NOT only for statistics lovers - LONG
Which I still do not understand in detail.

So all in all I'm happy with the outcome.

Cheers
Thomas

[Albert Chan]

Quadratic fit without linear algebra (32SII) - NOT only for statistics lovers - LONG

Linear regression on slope had another benefit, which my Casio cannot do.
Doing Quadratic Regression on Casio, you don't really know if it is any good.
There is no "correlation coefficient" to lookup. (Why?)

Linear regression on slope can show whether a Quadratic curve is warranted.

I think the code add slope data, like this:

Data points (0,90) and (30, 90.2), add to linear regression: (30+0)/2, (90.2-90)/(30-0)
Next point (60, 86.4), add to linear regression: (60+30)/2, (86.4-90.2)/(60-30)
...
N points (sorted) thus generate N-1 slopes to regress.

[Csaba Tizedes]

(02-18-2019 10:22 PM)Thomas Klemm Wrote:  

(02-18-2019 07:46 PM)Dieter Wrote:  Yes. Gamo already posted this method some time ago in the General Software Library. I have mentioned several times that the first input it NOT a guess for the root but an estimate for the function's DERIVATIVE.

You could think him stubborn.
But it made Csaba post his program which made me translate some of his papers from Hungarian to English and stumble totally unrelated upon one of his previous posts:
Quadratic fit without linear algebra (32SII) - NOT only for statistics lovers - LONG
Which I still do not understand in detail.

So all in all I'm happy with the outcome.

Cheers
Thomas

I feel myself like a star (at least for 15 minutes...)

Yes, this is my "linear regression on slope", as Chan wrote. When I was in my "Summer Practice" at a company I have modelled lots of water network flows between pump houses (sources), consumers and reservoirs, and sometimes the characteristics curves given only in paper form. Thats why I wrote this little routine.

As you can see, if the measured points has significant error, the slope will be very inaccurate then the regression also inaccurate and the coefficients can be meaningless.

But as a programming question this was an interesting "game" for me.

To the original question from Gamo: I will check carefully the lists here, but I have another idea: what if, if we just simple walking along on the function's curve until the sign is changing, then we turning back, decreasing the stepsize and walking again. If the sign of the function changed, we are turning back again, decreasing the size of the steps and so on... If the stepsize is less than a small value, we found the root.

This "marching method" not so elegant, not so fast, but good to play with it, just for a try - maybe it can be shorter than the other methods. (I have an exactly same method for finding a minimum point of a function for CASIO fx-3650P/Sencor SEC103; perfectly works, like a ball at the bottom of a hole, slowly find the deepest point, when swings around the extremum).

Csaba