The Museum of HP Calculators

HP Forum Archive 20

 35s - find roots of 3rd and higher order equationsMessage #1 Posted by Chris C on 19 Feb 2012, 8:42 p.m. 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.

 Re: 35s - find roots of 3rd and higher order equationsMessage #2 Posted by Valentin Albillo on 20 Feb 2012, 4:10 a.m.,in response to message #1 by Chris C 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: Best regards from V. ``` ```

 Re: 35s - find roots of 3rd and higher order equationsMessage #3 Posted by Bart (UK) on 20 Feb 2012, 7:01 p.m.,in response to message #2 by Valentin Albillo Nice one, thanks Valentin :-)

 Re: 35s - find roots of 3rd and higher order equationsMessage #4 Posted by Peter Murphy (Livermore) on 20 Feb 2012, 7:45 p.m.,in response to message #2 by Valentin Albillo 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!

 Re: 35s - find roots of 3rd and higher order equationsMessage #5 Posted by Paul Dale on 20 Feb 2012, 9:00 p.m.,in response to message #4 by Peter Murphy (Livermore) How about a 34s port for inclusion as an internal complex solve command? :-) - Pauli

 Re: 35s - find roots of 3rd and higher order equationsMessage #6 Posted by Paul Dale on 21 Feb 2012, 6:16 a.m.,in response to message #5 by 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

 Re: 35s - find roots of 3rd and higher order equationsMessage #7 Posted by Thomas Klemm on 21 Feb 2012, 8:06 p.m.,in response to message #4 by Peter Murphy (Livermore) 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: $\frac{\sqrt{1-\frac{2f(x)f''(x)}{f'(x)^2}}-1}{\frac{f''(x)}{f'(x)}}$ Compare this with the 2nd formula in Halley-Verfahren but don't ask me why it's missing in the English version: $\frac{-\frac{2f(x)}{f'(x)}}{1+\sqrt{1-2\frac{f(x)f''(x)}{f'(x)^2}}}$ To see that both formulas are equivalent use the following abbreviations: $u=\frac{f(x)}{f{(x)}'}$ $v=\frac{f{(x)}''}{f{(x)}'}$ The trick is to use the binomial theorem (a + b)(a - b) = a2 - b2: $\frac{-2u}{1+\sqrt{1-2uv}}\cdot\frac{\sqrt{1-2uv}-1}{\sqrt{1-2uv}-1}=\frac{-2u(\sqrt{1-2uv}-1)}{(1-2uv)-1}=\frac{\sqrt{1-2uv}-1}{v}$ 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 ```

 Re: 35s - find roots of 3rd and higher order equationsMessage #8 Posted by bill platt on 21 Feb 2012, 8:53 p.m.,in response to message #7 by Thomas Klemm This is so cool, it deserves its own permanent article.

 Re: 35s - find roots of 3rd and higher order equationsMessage #9 Posted by Peter Murphy (Livermore) on 21 Feb 2012, 9:31 p.m.,in response to message #7 by Thomas Klemm 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!

 Re: 35s - find roots of 3rd and higher order equationsMessage #10 Posted by Thomas Klemm on 21 Feb 2012, 10:14 p.m.,in response to message #9 by Peter Murphy (Livermore) 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

 Re: 35s - find roots of 3rd and higher order equationsMessage #11 Posted by C.Ret on 22 Feb 2012, 6:31 p.m.,in response to message #10 by Thomas Klemm 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.

 Re: 35s - find roots of 3rd and higher order equationsMessage #12 Posted by Thomas Klemm on 25 Feb 2012, 2:55 p.m.,in response to message #11 by C.Ret 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

Go back to the main exhibit hall