|
Loading LOADSUB ALL FROM "PROOTS.HTS"
or LOADSUB FROM "MATHLIB.HTS"
or LOADSUB Proots FROM "MATHLIB.HTS"
Usage INTEGER M,N
REAL P(*)
COMPLEX Z(*)
CALL Proots(N,P(*),M,Z(*))
Description
Proots attempts to finds all values of z that satisfy the relation p(z) = 0, where p is a polynomial of degree n, which must be 10 or less. The roots will usually be returned in Z in order of increasing modulus. Since the coefficients of p are real, complex roots returned in Z will occur in conjugate pairs.
If a root is not found after m iterations of the algorithm, that root is set to CMPLX(MAXREAL,0.0) to indicate an error and Proot returns without attempting to find more roots.
The first element in the array P represents the constant term in the polynomial; the second element, the linear term; the third the quadratic term, etc. The interpretation of the elements in P is without regard to the OPTION BASE in effect or any lower bound specified when P was declared.
The array Z must contain enough elements to hold all the solutions to the expression p(z) = 0, that is, Z should contain at least n elements.
Proots uses LaGuerre's method to find a real root or a conjugate pair of roots to the equation p(z) = 0. It then reduces the equation by dividing p(z) by the term z - zn for a real root zn or z – 2ℜ(zn) + |zn|2 for a complex root and repeats the procedure for the reduced polynomial. This stops when all the roots are found or when the algorithm fails to find a new root after m iterations.
Typical values for m might be 50, 100, or 200.
Errors
Proots causes an HTBasic error if n < 2, n > 10, if P contains fewer than n + 1 elements, if Z contains fewer than n elements, or if an evaluation of the polynomial being used results in a value greater in magnitude than MAXREAL.
Example
The following BASIC program calculates the roots of the function
x3 - 3x2 + 3x - 2 = 0.
10 LOADSUB ALL FROM "PROOTS.HTS"
20 INTEGER I
30 REAL C(0:3)
40 COMPLEX Z(1:3)
50 READ C(*)
60 DATA -2.0,3.0,-3.0,1.0
70 CALL Proots(3,C(*),100,Z(*))
80 FOR I=1 TO 3
90 PRINT USING """("",MZ.6D,"","",MZ.6D,"")""";Z(I)
100 NEXT I
110 END
When run, it prints
( 0.500000,-0.866025)
( 0.500000, 0.866025)
( 2.000000, 0.000000).
The roots of the equation can be found by hand. They are
See Also
Froot, Paderiv, Pderiv, Pinteg, Poly
| 
