35s - find roots of 3rd and higher order equations « Next Oldest | Next Newest »

 ▼ Chris C Junior Member Posts: 21 Threads: 6 Joined: Feb 2012 02-19-2012, 08:42 PM 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 Posting Freak Posts: 1,755 Threads: 112 Joined: Jan 2005 02-20-2012, 04:10 AM 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. ``` ``` ▼ Bart (UK) Posting Freak Posts: 850 Threads: 10 Joined: Mar 2009 02-20-2012, 07:01 PM Nice one, thanks Valentin :-) Peter Murphy (Livermore) Member Posts: 167 Threads: 33 Joined: Jul 2011 02-20-2012, 07:45 PM 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 Posting Freak Posts: 3,229 Threads: 42 Joined: Jul 2006 02-20-2012, 09:00 PM How about a 34s port for inclusion as an internal complex solve command? :-) - Pauli ▼ Paul Dale Posting Freak Posts: 3,229 Threads: 42 Joined: Jul 2006 02-21-2012, 06:16 AM 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 Senior Member Posts: 735 Threads: 34 Joined: May 2007 02-21-2012, 08:06 PM 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 ``` ▼ bill platt Posting Freak Posts: 2,448 Threads: 90 Joined: Jul 2005 02-21-2012, 08:53 PM This is so cool, it deserves its own permanent article. Peter Murphy (Livermore) Member Posts: 167 Threads: 33 Joined: Jul 2011 02-21-2012, 09:31 PM 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 Senior Member Posts: 735 Threads: 34 Joined: May 2007 02-21-2012, 10:14 PM 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 Senior Member Posts: 260 Threads: 0 Joined: Oct 2008 02-22-2012, 06:31 PM 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. ▼ Thomas Klemm Senior Member Posts: 735 Threads: 34 Joined: May 2007 02-25-2012, 02:55 PM 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 2,384 11-21-2013, 08:25 PM Last Post: Chris Pem10 HP prime: logs and roots Alberto Candel 5 1,585 11-20-2013, 12:31 PM Last Post: Manolo Sobrino My 2nd and 3rd Prime Programs Jeff O. 7 2,032 10-31-2013, 08:18 AM Last Post: Jeff O. 3rd order ODE diff eq Richard Berler 0 814 10-23-2013, 09:53 PM Last Post: Richard Berler Prime cant find Memory available Dougggg 5 1,699 10-07-2013, 07:24 PM Last Post: Han HP Prime can't find how to use partfrac function ! dg1969 3 1,286 10-04-2013, 09:25 PM Last Post: Joe Horn Rearanging the "Fixed" sort order of Prime's Apps icons Joe Horn 3 1,375 10-02-2013, 10:24 AM Last Post: Han HP Prime Solving Nonlinear System of Equations for Complex Results Helge Gabert 11 3,471 09-30-2013, 03:44 AM Last Post: From Hong Kong OT - A lucky find - Busicom Handy-LE LE-120A Cristian Arezzini 2 1,485 09-26-2013, 04:43 AM Last Post: Cristian Arezzini New article on a new type of neo linear equations Namir 0 1,108 08-11-2013, 10:27 AM Last Post: Namir

Forum Jump: 