Re: 35s - find roots of 3rd and higher order equations Message #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.
|