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



#2

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.


#3

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.


#4

Nice one, thanks Valentin :-)

#5

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!


#6

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


- Pauli


#7

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

#8

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) = a2 - b2:

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


#9

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

#10

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!


#11

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


#12

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 xn+1 = xn – Hal(xn) converge to one of the f function roots depending of the initial vaue x0.

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 xx=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 x3-3.x2+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 X3 – (1+10i).X2 - (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.


#13


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


Possibly Related Threads...
Thread Author Replies Views Last Post
  [HP Prime] Tips for Solving Differential Equations More Efficiently Chris Pem10 8 380 11-21-2013, 08:25 PM
Last Post: Chris Pem10
  HP prime: logs and roots Alberto Candel 5 256 11-20-2013, 12:31 PM
Last Post: Manolo Sobrino
  My 2nd and 3rd Prime Programs Jeff O. 7 382 10-31-2013, 08:18 AM
Last Post: Jeff O.
  3rd order ODE diff eq Richard Berler 0 122 10-23-2013, 09:53 PM
Last Post: Richard Berler
  Prime cant find Memory available Dougggg 5 276 10-07-2013, 07:24 PM
Last Post: Han
  HP Prime can't find how to use partfrac function ! dg1969 3 245 10-04-2013, 09:25 PM
Last Post: Joe Horn
  Rearanging the "Fixed" sort order of Prime's Apps icons Joe Horn 3 225 10-02-2013, 10:24 AM
Last Post: Han
  HP Prime Solving Nonlinear System of Equations for Complex Results Helge Gabert 11 563 09-30-2013, 03:44 AM
Last Post: From Hong Kong
  OT - A lucky find - Busicom Handy-LE LE-120A Cristian Arezzini 2 196 09-26-2013, 04:43 AM
Last Post: Cristian Arezzini
  New article on a new type of neo linear equations Namir 0 195 08-11-2013, 10:27 AM
Last Post: Namir

Forum Jump: