wscott@EN.ECN.PURDUE.EDU (Wayne H Scott) (03/14/90)
Posted for a friend:
From ryoder Mon Mar 12 21:22:34 1990
Date: Mon, 12 Mar 90 17:33:11 -0500
From: ryoder (Robert W Yoder)
To: wscott
Here are new, improved versions of the root finder routines I sent before.
They all use the same input format and are no longer constrained to use
a specific variable.
I also included a cubic root finder that uses ROOT, although it
is subject to the error of a poor first guess.
Last is the new improved version of Knapik's CUBIC.
NOTE: Programs containing the variable 's1' assume flag 34 is set,
i.e. 'General' solutions mode.
These programs are good for NUMERICAL roots only.
If ROOT2 returns a root of 0, it may cause a crash in CUBE2.
CUBE2 calls ROOT2 which in turn calls NWTN. Each can also be used alone.
example:
input: level 3: cubic equation S^3+8*S^2+6*S+48.
level 2: variable S
level 1: first guess for root 2
output: level 3: root 1 -8
level 2: root 2 (0,2.44948974278)
level 1: root 3 (0,-2.44948974278)
CUBE2 [7A18]
<< ROOT2 -> f x r
<< r f x r - / COLCT
x 3 TAYLR x QUAD DUP
1 's1' STO EVAL SWAP
-1 's1' STO EVAL
's1' PURGE
>>
>>
ROOT2 calls NWTN repeatedly until it gets the same answer twice.
example:
input: level 3: equation (any degree) S^3+8*S^2+6*S+48.
level 2: variable S
level 1: first guess for root 2
output: level 3: original equation S^3+8*S^2+6*S+48.
level 3: original variable S
level 2: one root -8
ROOT2 [244B]
<<
DO DUP 'a' STO
NWTN DUP
UNTIL a SAME
END 'a' PURGE
>>
NWTN is set up so you can use it repeatedly simply by pressing 'NWTN'
again and again. Each iteration takes you closer to the nearest root,
and usually 4 to 5 iterations will find it exactly.
example:
input: level 3: equation (any degree) S^3+8*S^2+6*S+48.
level 2: variable S
level 1: first guess for root 2
output: level 3: original equation S^3+8*S^2+6*S+48.
level 3: original variable S
level 2: root approximation 0
for this example, the next iteration gives the exact value of -8.
NWTN [800F] (Newton's Method) '&' represents 'derivative'
<< -> f x g
<< f x g f DUP x &
/ - g x STO EVAL x
PURGE
>>
>>
The next program is a stand-alone cubic root finder using the built-in
ROOT function instead of NWTN. The nice part is that it is short and
takes advantage of built-in functions. The bad part is that it is susceptible
to all the quirks inherent in ROOT. For example, if you give it the same
input as in my example for CUBE2 above, it gives the wrong answer.
example:
input: level 3: cubic equation S^3+8*S^2+6*S+48.
level 2: variable S
level 1: first guess for root -7
output: level 3: root 1 -8
level 2: root 2 (0,2.44948974278)
level 1: root 3 (0,-2.44948974278)
CUBE [B65]
<< -> f x g
<< f x g ROOT x
PURGE DUP f SWAP x
SWAP - / x 3 TAYLR x
QUAD DUP 1 's1' STO
EVAL SWAP -1 's1'
STO EVAL 's1' PURGE
>>
>>
I decided to work on David Knapik's CUBIC to see if I could find the problem.
I located the same book he used and found he had made no mistakes. When I
checked another book, I found the same algorithm, (Cardin's Method), but this
one just happened to mention that when taking the cube roots, the REAL ones
are what must be used!
That is, if the value of which you are taking the cube root is negative,
it will have 3 roots: 1 real and 2 complex. And for reasons known only to HP,
the '28 happily gives you one of the complex roots.
I solved the problem by saving the SIGN of the value; taking the ABS;
taking the cube root; then restoring the sign.
Here is an improved version using 5% less memory and running 5 times FASTER!
I eliminated all global variables; shortened variable names to one character;
and eliminated redundant computations by using 3 more intermediate variables.
Now for the bad news: after trying this program successfully on all the
equations for which the original CUBIC had given me garbage answers, I
found one for which it doesn't work; x^3+23*x^2+167*x+385. If anyone can
tell me why, let me know. (the roots for the above are -5, -7, -11)
example: to solve: x^3+A*x^2+B*x+C (first coefficient must be 1!)
input: level 3: A
level 2: B
level 1: C
output: level 3: root 1
level 2: root 2
level 1: root 3
CUBIC [2C3]
<< -> a b c
<< a 3 / -> d
<< b 3 / d SQ - 3 ^ -> f
<< b a * 3 c * - 6 / d 3 ^ - -> g
<< f g SQ + SQRT -> h
<< 0 1 FOR p g h -1 p ^ + DUP SIGN -> s
<< ABS 3 INV ^ s * >> NEXT -> j k
<< j k + -> m
<< m d - ;first root
m -2 / d -
-3 SQRT 2 / j k - * -> n o
<< n o + ;second root
n o - ;third root
>> >> >> >> >> >> >> >> >>
Written by: Robert W Yoder
ryoder@en.ecn.purdue.edu
_______________________________________________________________________________
Wayne Scott | INTERNET: wscott@en.ecn.purdue.edu
Electrical Engineering | BITNET: wscott%ea.ecn.purdue.edu@purccvm
Purdue University | UUCP: {purdue, pur-ee}!en.ecn.purdue.edu!wscott