HP 11C real root finder [Newton Method]

[Carlos CM (Mexico)]

LBL 0
RCL 1
GSB 1
RCL 1
GSB 2
/
CHS
RCL 1
+
STO 3
RCL 1
-
ABS
RCL 2
X>Y?
RTN
RCL 3
STO 1
GTO 0

LBL 1
f(x)=0 code

LBL 2
f'(x)=0 code

R1: old value (or initial value)
R2: Tolerance of error (0.001, 0.0001, 0.000000001 etc)

THE RESULT:
R3: new value (root)

EXAMPLE

Find the root of f(x)=x^3 - 3x^2-6x+8

rewritten as f(x)= [(x-3)x-6]x +8

f(x) code:

LBL1
ENTER
ENTER
ENTER
3
-
*
6
-
*
8
+
RTN


f'(x)= 3x(x-6)-6
f'(x) code

LBL 2
ENTER
ENTER
6
-
*
3
*
6
-
RTN

The roots are [-2, 1, 4]

give a initial value, for example: -3

-3 STO 1

give a tolerance of error, for example: 0.0001

0.0001 STO 2

BEGIN THE PROGRAM...

GSB 0

runnning...

when the error < TOL the program stop and display:

0.0001

you can find the root in the register 3

RCL 3

DISPLAY: -2.0001

try with initial value = 2

















[/size][/font][/size][/size][/font][/font]

[Thomas Klemm]

Using storage arithmetic allows us to get rid of a few steps:

LBL 0
RCL 1
GSB 1
RCL 1
GSB 2
/
STO- 1
ABS
RCL 2
X<=Y?
GTO 0
RCL 1
RTN


You made a mistake calculating the 1st derivative.
This is the corrected program:

f'(x)= 3x^2 - 6x - 6 = (3x - 6)x - 6
f'(x) code

LBL 2
ENTER
ENTER
3
*
6
-
*
6
-
RTN

Kind regards
Thomas

[Carlos CM (Mexico)]

Thanks to reply Thomas and for your help in the code.

Best regards

[Namir]

Here is a version that requires coding f(x) only since it approximate f'(x) as:

f'(x) = (f(x+h) - f(x))/h

Where h = 0.001*(ABS(X)+1)

The new version uses registers R1 through R4.

Code:

LBL 0
RCL 1
ABS
1
+
EEX
3
CHS
*
STO 3            # calculate and store increment h
RCL 1
GSB 1
STO 4            # calculate and store f(x)
RCL 1
RCL 3
+
GSB 1            # calculate f(x+h)
RCL 4
-
1/X
RCL 4
*
RCL 3
*                    # calculate diff = h *f(x)/(f(x+h) - f(x))
STO- 1
ABS
RCL 2
X<=Y?
GTO 0
RCL 1
RTN

[Dieter]

(01-12-2014 08:31 AM)Namir Wrote:  ...h = 0.001*(ABS(X)+1)

First of all, instead of multiplying with \(10^{-3}\), dividing by \(10^{3}\) is one step shorter. ;-)

This method for determining h will work in most cases, but not for very small arguments. Consider \(x = 10^{-4}\) or even \(x = 10^{-40}\). That's why I prefer \(h = x/10^4\). On the 34s, the result can be easily rounded to 1 or 2 significant digits (RSD 1) to prevent slight roundoff errors. As usual, \(x=0\) is handled as \(x=1\).

Dieter

[Thomas Klemm]

If you only want to solve polynomials with real coefficients Bairstow's Method can be used.
In a nutshell: instead of a root a quadratic factor is found in each iteration. To find the roots of this factor the classic formula is used.

Example:

\(x^3-3x^2-6x+8=0\)

First enter the coefficients of the polynomial. Always start with register 9:

1 STO 9
-3 STO .0
-6 STO .1
8 STO .2


Then specify the order of the polynomial with a loop-control value. It defines the registers you used for the coefficients:

9.012 STO 6


And now give an initial guess. Probably {1, 1} will do in all cases:

1 STO 7
STO 8


Start the program B:

GSB B


Now the coefficients of the quadratic factor can be found in registers 7 and 8:

RCL 7
1.000000

RCL 8
-2.000000


The coefficient of the first term is always 1. Thus the factor is \(x^2+x-2\).
Now we can solve this quadratic equation using program A:

1
RCL 7
RCL 8
GSB A

1.000000
X<>Y
-2.000000


Thus: \(x^2+x-2=(x-1)(x+2)\).

What is left? Let's have a look at it:

RCL 6
9.010

RCL 9
1.000000

RCL .0
-4.000000


This is a linear factor. Thus we end up with the following factorization: \(x^3-3x^2-6x+8=(x-1)(x+2)(x-4)\)
Therefore the solutions are: {1, -2, 4}


Example from Bunuel66's solution to the crossed ladders problem

\(C^4-30C^3+700C^2-21,000C+157,500=0\)

Enter the coefficients of the polynomial:

1 STO 9
-30 STO .0
700 STO .1
21,000 STO .2
157,500 STO .3


Specify the loop-control value:

9.013 STO 6


Use the initial guess {1, 1}:

1 STO 7
STO 8


Start the program:

GSB B


The coefficients of the quadratic factor are in registers 7 and 8:

RCL 7
-35.178025

RCL 8
248.596895


Solve this quadratic equation using program A:

1
RCL 7
RCL 8
GSB A

25.384938
X<>Y
9.793087


What about the rest?

RCL 6
9.011

RCL 9
1.000000

RCL .0
5.178025

RCL .1
633.555781


If we try to solve this quadratic equation we get an Error 0.
That's because we're trying to calculate the square root of a negative value.
However we can still get the desired result:

RCL 9
RCL .0
RCL .1
GSB B
Error 0

<-
CHS
626.852796

\(\sqrt{x}\)
25.037029

X<>Y
-2.589012

Thus the final list of solutions is:

• 25.384938
• 9.793087
• -2.589012\(\pm\)25.037029i
  

Cheers
Thomas

[Namir]

(01-12-2014 01:26 PM)Dieter Wrote:  

(01-12-2014 08:31 AM)Namir Wrote:  ...h = 0.001*(ABS(X)+1)

First of all, instead of multiplying with \(10^{-3}\), dividing by \(10^{3}\) is one step shorter. ;-)

This method for determining h will work in most cases, but not for very small arguments. Consider \(x = 10^{-4}\) or even \(x = 10^{-40}\). That's why I prefer \(h = x/10^4\). On the 34s, the result can be easily rounded to 1 or 2 significant digits (RSD 1) to prevent slight roundoff errors. As usual, \(x=0\) is handled as \(x=1\).

Dieter

Sure, dividing by 1000 is a step shorter. Originially, I had learned to calculate h using an If statement:

Code:

If |x| >= 1 Then
  h= x/100
else
  h = 1/100
end if

Until I realized one day that .01*(|x|+1) does the job and eliminates the need for labels and GOTOs. I recently started using .001 instead of 0.01. Using the expression for h ensures that if x=0, h is not zero.

Namir

[Dieter]

Namir, sorry that I was not able to express myself clearly. I wanted to point out that I do not think it's a good idea to use h = 0,01 or 0,001 for any value below 1. This may lead to significant errors since h may be much, much larger than x. For instance, if x = 1E-10, h = 1E-3 is not recommended. Here, h = 1E-13 should be better.

That's why I prefer to set h = 0,001 x or similar. With the only exception x=0 where x may be 0,001.

(01-15-2014 05:53 AM)Namir Wrote:  Until I realized one day that .01*(|x|+1) does the job and eliminates the need for labels and GOTOs. I recently started using .001 instead of 0.01. Using the expression for h ensures that if x=0, h is not zero.

This can be coded this way, for instance:

Code:

 X=0?
 e^x
 EEX
  3
  /

Look, no labels or gotos required. ;-)

Or even more elegant on the 34s:

Code:

 X=0?
 INC X
 SDR 3
 RSD 1   ' optional

Dieter