35s - find roots of 3rd and higher order equations

[Chris C]

Are there any programs (or built-in methods I'm overlooking) for the 35s to find the roots of equations higher than second order?

Thank you.

[Valentin Albillo]

Quote:

---

Are there any programs (or built-in methods I'm overlooking) for the 35s to find the roots of equations higher than second order?

Thank you.

---

Sure. Have a look at this article of mine, it features a short program for the HP35s (plus assorted examples) which will allow you to find any roots, real or complex, of most any equation, polynomial or not:

HP35s - Going back to the roots

Best regards from V.

[Bart (UK)]

Nice one, thanks Valentin :-)

[Peter Murphy (Livermore)]

Valentín,

Is there a reference for the algorithm on which you based this superb program? I'd like to start there, rather than reverse-engineer your program. Thanks!

[Paul Dale]

How about a 34s port for inclusion as an internal complex solve command? :-)

- Pauli

[Paul Dale]

It is in the library now :-)

Not sure if I ported it over completely properly or not but it handles the first example from the PDF article okay.

- Pauli

[Thomas Klemm]

In order to analyze Valentin's program I've translated it to Python while trying to keep it as close to the original as possible. You will find the whole program at the end. The examples are added as well.

Let me go through it step by step.

First we need the definition of some constants:

Z = 1 + 0j

S = 1E-4

T = S**2   # i.e. 1E-8

Y = 0.5

Now we use the function:

    U = F(X) / Y

This leaves us with: U = 2F(X)

For the next step some approximations are used:

• F"(X) ~ (F(X+S) - 2F(X) + F(X-S))/S^2
  
• F'(X) ~ (F(X+S) - F(X-S))/(2S)
  

    V, W = (W + V - U) / T, Y * (V - W) / S

So we end up with:

• V ~ F"(X)
  
• W ~ F'(X)
  

Next step is to divide F"(X) by F'(X):

    V /= W

Now we have: V ~ F"(X)/F'(X)

Since I don't have a stack I use a variable D which is not used in Valentin's program:

    D = ( (Z - V * U / W) ** Y - Z ) / V

This gives the following expression:

Compare this with the 2nd formula in Halley-Verfahren but don't ask me why it's missing in the English version:

To see that both formulas are equivalent use the following abbreviations:

The trick is to use the binomial theorem (a + b)(a - b) = a² - b²:

Kind regards

Thomas

---

from math import pi, copysign

from cmath import sin, cos

from cmath import phase as arg

# 1. Find a root of : (a) x^x = pi , (b) x^x = i

def F(X):

    return X**X - pi

X = 2

def F(X):

    return X**X - 1j

X = 1j

# 2. Find all roots of: ( 2 + 3i ) x^3 - (1 + 2i ) x^2 - ( 3 + 4i ) x - ( 6 + 8i ) = 0

def F(X):

    return (2+3j)*X**3 - (1+2j)*X**2 - (3+4j)*X - (6+8j)

X = 1

X = -1

X = 1+1j

# 3. Attempt to find a complex root of: x^3 - 6x - 2 = 0

def F(X):

    return X**3 - 6*X - 2

X = 2+3j

# 4. Solve Leonardo di Pisa's equation: x^3 + 2 x^2 + 10 x - 20 = 0

def F(X):

    return X**3 + 2*X**2 + 10*X - 20

X = 1

X = -6

# 5. Find several complex roots of : Sin ( 2 x - 4 i ) + 3 x^2 - ( 1 + 5 i ) = 0

def F(X):

    return sin(2*X - 4j) + 3*X**2 - (1+5j)

X = 0

X = -1

X = pi

# here starts the program

Z = 1 + 0j

X *= Z

S = 1E-4

T = S**2

Y = 0.5

D = 1

while (T < abs(D)):

    U = F(X) / Y

    X += S

    V = F(X)

    X -= S

    X -= S

    W = F(X)

    V, W = (W + V - U) / T, Y * (V - W) / S

    X += S

    V /= W

    D = ( (Z - V * U / W) ** Y - Z ) / V

    X += D

    D /= X

W = arg(X)

if abs(sin(W)) < T:

    X = copysign(abs(X), cos(W))

print "X = %r" % X

[bill platt]

This is so cool, it deserves its own permanent article.

[Peter Murphy (Livermore)]

Thomas,

It'll take me some time to work through your analysis of Valentin's program, but I must thank you now; I'm very grateful for the effort you've gone to, and I hope I can do it justice by understanding it well enough to write an RPL program (for 48G or 50g) that will do all that the 35S program does.

Many thanks!

[Thomas Klemm]

Meanwhile I found another reference:
Halley's Irrational Formula (cf. formula #4)

Let me know if something is not clear in my analysis and please post your RPL program.

Cheers

Thomas

[C.Ret]

The numerical analysis Halley’s method is a root-finding algorithm used for functions of one variable with a continuous second derivative.

Where f is the function were are investigating for roots.

Posing Hal(x)= 2.F(x)*F’(x) / ( 2 * [F’(x)]^2 – F(x).F(“(x)) )
The series x_n+1 = x_n – Hal(x_n) converge to one of the f function roots depending of the initial vaue x₀.

In RPL, the code to follow such an iterative method will be as simple as the following:
Imagine the initial value for the root-finding is in first stack level.

«  (1,0) *                            @  Converting x0 into complex z0 

   0                                  @     

   DO

      -                               @  x(n+1) <-- x(n) – Hal(x(n))

      DUP 1 DISP                      @  Display actual estimate

      DUP Hal                         @  Hal(x)

   UNTIL DUP ABS eps < END            @  if Ha(x) small enough, stop

   IF – DUP ARG SIN ABS eps <         @  if arg(x) small enough

        THEN RE END                   @  convert it to real

   CLMF

»

‘R.HALLEY’ STO

In RPL, invoking R.Halley may return the closest root from the inputted initial. It is quite like the LBL A program on HP-35C.

Nut this will be true only if the user function Hal is the appropriate one.

From the above formula, one may calculate the Hal function from any given secondly derivative continuous function and input it in the RPL calculator.

In RPL, a general formulae of the Hal user’function may be :

«  ->  x

   '2*F(x)*dF(x)/(2*SQ(dF(x))-F(x)*d2F(x))’

»

‘Hal’ STO

Where F(x), dF(x) and d2F(x) stand respectively for function, first derivate and second derivate of the function under investigation.

One major advantage of RPL is that they all are able to formally derivate any given function.
So the user have not to calculate, write down the first and second derivates, nor the ‘Hal’ function, its RPL system can easily made it :

One way is to use one initialization program of the form :

«  DUP ->STRU                @  copy ‘x’ as unquote string “x”

   -> x Sx           

      « 

        "'F("  Sx + STR->    @  put ‘F(x)’ on stack

        OVER = DEF           @  define user function ‘F(x)=...’

        x d/dx COLCT         @  derivate expression

        "'dF("  Sx + STR->   @  put ‘dF(x)’ on stack

        OVER = DEF           @ define user function  F’ =dF()

        x d/dx COLCT         @  derivate expression

        "'d2F(" Sx + STR->   @  put ‘d2F(x)’ on stack

        SWAP = DEF           @  define user function F” = d2F()

        { F dF d2F Hal R.Halley INIT } MENU

                             @ create dedicated custom menu

        35 CF                @ on HP-28S numeric evaluation of constant

»

‘INIT’ STO

The user have to enter formal definition of the function/expression to investigate in level 2: and the name of the variable to scan in level 1:. Then pressing INIT will built up first and second derivate and create a soft menu with dedicated key.
Next steps are to seek for root by trying different initial guess.
As for the HP-35C, true real root are put in stak as real number and complex roots as complex number.

Examples:

To solve x^x=PI :

>> x x ^ PI – x INIT (takes a few secondes on an HP-28S)

3:                                            3:                                         

2:                          'x^x-ð'           2:                                         

1:                              'x'           1:                                         

 [INIT ]                                      [  F  ][ dF  ][ d2F ][ Hal ][R.HALL][INIT ]

>> 2 [R.HALL]

(2,0)                                         3:                                         

1:                                            2:                                         

2                                             1:                            1.85410596792

[  F  ][ dF  ][ d2F ][ Hal ][R.HALL][INIT ]   [  F  ][ dF  ][ d2F ][ Hal ][R.HALL][INIT ]

The soft key can easily be use to check results or to show formulae :

>> x [ dF  ] returns ‘x^x+LN(x)*x^x’

And

>> 'dF' RCL return « -> x  'x^x+LN(x)*x^x' »

>> ‘d2F’ RCL return « -> x  '(x^x+LN(x)*x^x)*LN(x)+x^x+LN(x)*x^x+x^(x-1)' »

To solve x³-3.x²+4.x-2 = 0 :

>> 'x*x*x-3*x*x+4*x-2' COLCT x [INIT] 0 [R.HALL] (2,3) [R.HALL] (2,3)[CHS] [R.HALL] 5 FIX

3:                                            3:                                  1.00000

2:               '-2-3*x^2+x^3+4*x'           2:                        (1.00000,1.00000)

1:                              'x'           1:                       (1.00000,-1.00000)

[  F  ][ dF  ][ d2F ][ Hal ][R.HALL][INIT ]   [  F  ][ dF  ][ d2F ][ Hal ][R.HALL][INIT ]

To solve x To solve X³ – (1+10i).X² - (33-6i).X + (5-40i) = 0 :

>> 'X^3-(1,10)*X*X-(33,-6)*X+(5,-40)' COLCT X [INIT] 1 [R.HALL] (1,1) [R.HALL] (0,1)[CHS] [R.HALL]

3:                                            3:                                   (-1,2)

2:     '(5,-40)+(-1,-10)*X^2+X^3+(-33,6)*X'   2:                                    (2,3)

1:                                      'X'   1:                                    (0,5)

[  F  ][ dF  ][ d2F ][ Hal ][R.HALL][INIT ]   [  F  ][ dF  ][ d2F ][ Hal ][R.HALL][INIT ]

Edited: 22 Feb 2012, 6:48 p.m.

[Thomas Klemm]

As I'm not familiar with the HP-28 I made a few modifications to you programs which allows me to use them with an HP-48:

%%HP: T(3)A(D)F(.);

DIR

  eps .00000001

  Hal

    \<< \-> x '2*F(x)*dF(x)/(2*SQ(dF(x))-F(x)*d2F(x))'

    \>>

  R.Halley

    \<< (1,0) * 0

      DO

        -

        DUP 1 DISP

        DUP Hal

      UNTIL DUP ABS eps <

      END

      IF - DUP ARG SIN ABS eps <

      THEN RE

      END

    \>>

  INIT

    \<< \-> x

      \<<

        "'F(" x + STR\->

        OVER = DEFINE

        x \.d COLCT

        "'dF(" x + STR\->

        OVER = DEFINE

        x \.d COLCT

        "'d2F(" x + STR\->

        SWAP = DEFINE

        { F dF d2F Hal R.Halley INIT } MENU

        -2 SF

      \>>

    \>>

END

However I encountered a problem when I tried to solve Valentin's 4th problem:

Quote:

---

# 4. Solve Leonardo di Pisa's equation: x^3 + 2 x^2 + 10 x - 20 = 0

 1, XEQ A -> X =  1.36880810782

-6, XEQ A -> X = -1.68440405391 i -3.4313313502

---

When I use the same intital values I always get the same result for 1 as well as for -6:

>> 'x^3+2*x^2+10*x-20' x [INIT ]

3:

2:                     'x^3+2*x^2+10*x-20'

1:                                     'x'

[  F  ][ DF  ][ D2F ][ HAL ][R.HAL][INIT ]

>> 1 [R.HAL] -6 [R.HAL]

4:

3:

2:                           1.36880810782

1:                           1.36880810782

[  F  ][ DF  ][ D2F ][ HAL ][R.HAL][INIT ]

It is obvious that from real initial values you will never get complex roots. That's the benefit of using Halley's Irrational Formula because a square root is used. But thanks to your clever design all I had to do was to replace the program Hal.

In your version of Hal both functions F(x) and dF(x) are called twice with the same value. That's why I save the result of the function calls in local variables u and v (well not exactly that, but their quotients):

  Hal

    \<< \-> x

      \<<

        x dF

        x F

        OVER /

        x d2F

        ROT /

        \-> u v

        '2*u/(1+\v/(1-2*u*v))'

      \>>

    \>>

With this program I get the same results as with the HP-35S:

>> 1 [R.HAL] -6 [R.HAL]

3:

2:                           1.36880810782

1:                     (-1.68440405391,

                       -3.4313313502)

[  F  ][ DF  ][ D2F ][ HAL ][R.HAL][INIT ]

Thanks for posting your program: I've learned a lot from the ideas of your INIT routine. Very nice!

Kind regards

Thomas