| Title: | Mathematics at DEC |
| Moderator: | RUSURE::EDP |
| Created: | Mon Feb 03 1986 |
| Last Modified: | Fri Jun 06 1997 |
| Last Successful Update: | Fri Jun 06 1997 |
| Number of topics: | 2083 |
| Total number of notes: | 14613 |
This is a really interesting conference. (Un)fortunately it's like
looking up a word in the dictionary in that you start reading the
meanings of all the other words and forget the one you wanted in the
first place.
I am looking to be apprised of the state of the art and the standard
of the art in finding the zeroes of polynomials in a single complex
variable where the polynomial has real coefficients.
I did a DIR/TITLE="ROOTS" and a DIR/TITLE="POLYNOMIALS" but did not
get the desired information. I have information on a Graeffe's root
squaring method and a zillion variants of Newton's and other
gradient based methods. It seems that closed form solutions are
done at degree 4 and MAPLE and and a couple of other methods are
done at degree 5 (besides MAPLE is not reliable, repeatable, etc.).
There's a method in the HP15C "Advanced Functions Handbook" which I
haven't read yet and I just obtained a copy of "Intro to Applied
Numerical Analysis" by Hamming (McGraw-Hill Computer Science Series
c 1971) which has some stuff on this subject in it.
So far it looks like most algorithms have certain "singularities"
which are dependent on the polynomial and/or the last estimate of
the root. I desire an algorithm for root finding in electrical
system analysis. The algorithm needs to have the following
properties:
1. Finds ALL roots.
2. Conjugate pairs can be indentified so that they can be
kept exactly numerically conjugate.
3. User specified trial roots are not necessary (maybe I could
get around this but I want to avoid convergence properties
dependent on the trial root choice so that the trial root
choice can be automated and made transparent at some level)
4. 30-40th degree polynomials are possible.
OK people, you can stop laughing. With those specifications I'm
sure I don't have to go to great lengths to emphasize my
mathematical naivete'. If this information exists elsewhere in this
conference, I apologize. Simply redirect me and I will delete this
note if possible.
Thanks,
Coop
PS. Ultimately I would like to be able to use roots generated from
such an algorithm as input into a program that would give me partial
fraction expansion. But I don't have that algorithm yet either and
there are some sticky problems to resolve there.
| T.R | Title | User | Personal Name | Date | Lines |
|---|---|---|---|---|---|
| 866.1 | Not solvable...? | CHOVAX::YOUNG | Dumb, Expensive, Dumb ... (Pick Two) | Sat Apr 23 1988 23:19 | 7 |
> I am looking to be apprised of the state of the art and the standard
> of the art in finding the zeroes of polynomials in a single complex
> variable where the polynomial has real coefficients.
I thought that the general case for this problem was unsolvable
for polynomials of order 5 and above?
| |||||
| 866.2 | Good rootfinders are available | CTCADM::ROTH | If you plant ice you'll harvest wind | Mon Apr 25 1988 10:48 | 67 |
One of the standard black box routines is the Jenkens Traub rootfinder.
It's fairly complicated, but quite reliable and pretty fast. It's
described in Jenkens' Stanford PHD Thesis, along with working code;
also it is described in one of the numerical analysis journals, either
a SIAM journal, or possibly Mathematics of Computation or Numeriche
Mathematic. I'd have to check this for you. ALGOL code is available
in the published references, and I think FORTRAN code is in his thesis.
The ACM collected algorithms has an old FORTRAN version of Jenksns-Traub
which uses complex coefficients - more general than you need.
However, another competetive method is to populate the companion matrix
(a matrix whose characteristic equation happens to be the polynomial
in question) with the coefficients and solve that with an eigenvalue
finder. This turns out to work very well in practice, and I've had
no trouble with 100'th order (!) polynomials using this method.
There exists a very efficient (cubically convergent) eigenvalue finder
for matrices in so-called "upper Hessenberg" form - where the matrix
looks like this:
x x x x x
x x x x x
x x x x
x x x
x x
No more than one diagonal below the major diagonal is nonzero. This
luckily happens to be the form the companion matrix is in, so no trickey
preprocessing has to be done.
The eigenvalue finder you want is in EISPACK (the public domain package
of eigenroutines), and is called HQR - (Hessenberg QR). It returns an
array of real and imaginary eigenvalues, with conjugate pairs adjacent
in the list. The QR algorithm is the state of the art in eigenvalue
finders for upper Hessenberg matrices.
An HQR routine is in the book "Numerical Recipes". That book also has
a Laguerre root finder, but be warned that it probably won't work for
extremely high order polynomials - you have to play tricks like keeping
your own exponents when doing trial evaluations of high degree polynomials
because the dynamic range of numbers encountered can be very high.
The reason the eigen-approach is good for electrical network analysis
is that you can often cast the problem directly in such a form
rather than the round about way of making up a characteristic
polynomial first. And there are generalized eigenvalue routines just for
such situations.
Good references on this are:
"Numerical Recipes"
Press, et al - Cambridge
"EISPACK Guide"
Dongarra, et al - Springer Verlag
"Matrix Computations"
Golub and Van Loan - Johns Hopkins
I'll post a routine you can use below.
Hope this helps...
- Jim
| |||||
| 866.3 | HQR and demo program | CTCADM::ROTH | If you plant ice you'll harvest wind | Mon Apr 25 1988 11:21 | 335 |
PROGRAM DEMO_DHQR
C
C Show how to call the HQR eigenvalue finder to get complex roots of a
C high degreee polynomial
C
IMPLICIT REAL*8 (A-H, O-Z)
C
C Array of polynomial coefficients in decreasing powers of Z
C
REAL*8 C(51)
C
C Work array for N square matrix passed by reference
C
REAL*8 H(2500)
C
C Resulting real and imaginary parts of eigenvalues
C
REAL*8 WR(50), WI(50)
REAL*8 SCALE
REAL*4 MTH$RANDOM
INTEGER SEED
SEED = 12345678
DO 500 NDEG=6,51
C
C Setup a random polynomial for quick timing test
C
C F(Z) = C[1]*Z^DEG + C[2]*Z^(DEG-1) + ... + C[DEG+1]
C
DO 100 I=1,NDEG+1
C(I) = MTH$RANDOM(SEED)+1.0E-6*MTH$RANDOM(SEED)
100 CONTINUE
C
C Zero the N by N matrix and eigenvalue work arrays
C
DO 110 I=1,NDEG*NDEG
H(I) = 0.0
110 CONTINUE
DO 120 I=1,NDEG
WR(I) = 0.0
WI(I) = 0.0
120 CONTINUE
TYPE *,'NDEG = ',NDEG
C
C Pass an upper Hessenberg matrix whose characteristic polynomial is
C the polynomial we wish to factor to DHQR and recover the roots as the
C eigenvalues.
C
SCALE = -1.0/C(1)
J = 1
DO 130 I=1,NDEG
H(J) = C(I+1)*SCALE
IF (I.NE.NDEG) H(I+J) = 1.0
J = J+NDEG
130 CONTINUE
TYPE *,'Calling DHQR'
CALL LIB$INIT_TIMER
CALL DHQR(NDEG, NDEG, 1, NDEG, H, WR, WI, IERR)
CALL LIB$SHOW_TIMER
TYPE *,'DHQR RETURNED IERR = ', IERR
DO 200 I=1,NDEG
TYPE *, I, WR(I), WI(I)
200 CONTINUE
500 CONTINUE
END
SUBROUTINE DHQR(NM, N, LOW, IGH, H, WR, WI, IERR)
C DATE WRITTEN 760101 (YYMMDD)
C REVISION DATE 830518 (YYMMDD)
C CATEGORY NO. D4C2B
C KEYWORDS EIGENVALUES,EIGENVECTORS,EISPACK
C AUTHOR SMITH, B. T., ET AL.
C PURPOSE Computes eigenvalues of a real upper Hessenberg matrix
C using the QR method.
C
C This subroutine is a translation of the ALGOL procedure HQR,
C NUM. MATH. 14, 219-231(1970) by Martin, Peters, and Wilkinson.
C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 359-371(1971).
C
C This subroutine finds the eigenvalues of a REAL
C UPPER Hessenberg matrix by the QR method.
C
C On INPUT
C
C NM must be set to the row dimension of two-dimensional
C array parameters as declared in the calling program
C dimension statement.
C
C N is the order of the matrix.
C
C LOW and IGH are integers determined by the balancing
C subroutine BALANC. If BALANC has not been used,
C set LOW=1, IGH=N.
C
C H contains the upper Hessenberg matrix. Information about
C the transformations used in the reduction to Hessenberg
C form by ELMHES or ORTHES, if performed, is stored
C in the remaining triangle under the Hessenberg matrix.
C
C On OUTPUT
C
C H has been destroyed. Therefore, it must be saved
C before calling HQR if subsequent calculation and
C back transformation of eigenvectors is to be performed.
C
C WR and WI contain the real and imaginary parts,
C respectively, of the eigenvalues. The eigenvalues
C are unordered except that complex conjugate pairs
C of values appear consecutively with the eigenvalue
C having the positive imaginary part first. If an
C error exit is made, the eigenvalues should be correct
C for indices IERR+1,...,N.
C
C IERR is set to
C Zero for normal return,
C J if the J-th eigenvalue has not been
C determined after a total of 30*N iterations.
C
C Questions and comments should be directed to B. S. Garbow,
C APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
C ------------------------------------------------------------------
C REFERENCES B. T. SMITH, J. M. BOYLE, J. J. DONGARRA, B. S. GARBOW,
C Y. IKEBE, V. C. KLEMA, C. B. MOLER, *MATRIX EIGEN-
C SYSTEM ROUTINES - EISPACK GUIDE*, SPRINGER-VERLAG,
C 1976.
C
IMPLICIT REAL*8 (A-H, O-Z)
INTEGER I, J, K, L, M, N
INTEGER EN, LL, MM, NA, NM, IGH, ITN, ITS, LOW, MP2, ENM2, IERR
REAL*8 H(NM,N), WR(N), WI(N)
REAL*8 P, Q, R, S, T, W, X, Y, Z, NORM, S1, S2
LOGICAL NOTLAS
IERR = 0
NORM = 0.0
K = 1
C
C Store roots isolated by BALANC and compute matrix norm
C
DO 50 I = 1,N
DO 40 J = K, N
40 NORM = NORM + ABS(H(I,J))
K = I
IF (I .GE. LOW .AND. I .LE. IGH) GO TO 50
WR(I) = H(I,I)
WI(I) = 0.0
50 CONTINUE
EN = IGH
T = 0.0
ITN = 30*N
C
C Search for next eigenvalues
C
60 IF (EN .LT. LOW) GO TO 1001
ITS = 0
NA = EN - 1
ENM2 = NA - 1
C
C Look for single small sub-diagonal element
C FOR L=EN STEP -1 UNTIL LOW DO
C
70 DO 80 LL = LOW, EN
L = EN + LOW - LL
IF (L .EQ. LOW) GO TO 100
S = ABS(H(L-1,L-1)) + ABS(H(L,L))
IF (S .EQ. 0.0) S = NORM
S2 = S + ABS(H(L,L-1))
IF (S2 .EQ. S) GO TO 100
80 CONTINUE
C
C Form shift
C
100 X = H(EN,EN)
IF (L .EQ. EN) GO TO 270
Y = H(NA,NA)
W = H(EN,NA) * H(NA,EN)
IF (L .EQ. NA) GO TO 280
IF (ITN .EQ. 0) GO TO 1000
IF (ITS .NE. 10 .AND. ITS .NE. 20) GO TO 130
C
C Form exceptional shift
C
T = T + X
DO 120 I = LOW, EN
120 H(I,I) = H(I,I) - X
S = ABS(H(EN,NA)) + ABS(H(NA,ENM2))
X = 0.75 * S
Y = X
W = -0.4375 * S * S
130 ITS = ITS + 1
ITN = ITN - 1
C
C Look for two consecutive small sub-diagonal elements
C FOR M=EN-2 STEP -1 UNTIL L DO
C
DO 140 MM = L, ENM2
M = ENM2 + L - MM
Z = H(M,M)
R = X - Z
S = Y - Z
P = (R * S - W) / H(M+1,M) + H(M,M+1)
Q = H(M+1,M+1) - Z - R - S
R = H(M+2,M+1)
S = ABS(P) + ABS(Q) + ABS(R)
P = P / S
Q = Q / S
R = R / S
IF (M .EQ. L) GO TO 150
S1 = ABS(P) * (ABS(H(M-1,M-1)) + ABS(Z) + ABS(H(M+1,M+1)))
S2 = S1 + ABS(H(M,M-1)) * (ABS(Q) + ABS(R))
IF (S2 .EQ. S1) GO TO 150
140 CONTINUE
150 MP2 = M + 2
DO 160 I = MP2, EN
H(I,I-2) = 0.0
IF (I .EQ. MP2) GO TO 160
H(I,I-3) = 0.0
160 CONTINUE
C
C Double qr step involving rows l to EN and columns M to En
C
DO 260 K = M, NA
NOTLAS = K .NE. NA
IF (K .EQ. M) GO TO 170
P = H(K,K-1)
Q = H(K+1,K-1)
R = 0.0
IF (NOTLAS) R = H(K+2,K-1)
X = ABS(P) + ABS(Q) + ABS(R)
IF (X .EQ. 0.0) GO TO 260
P = P/X
Q = Q/X
R = R/X
170 S = SIGN(SQRT(P*P+Q*Q+R*R),P)
IF (K .EQ. M) GO TO 180
H(K,K-1) = -S * X
GO TO 190
180 IF (L .NE. M) H(K,K-1) = -H(K,K-1)
190 P = P+S
X = P/S
Y = Q/S
Z = R/S
Q = Q/P
R = R/P
C
C Row modification
C
DO 210 J = K, EN
P = H(K,J) + Q * H(K+1,J)
IF (.NOT. NOTLAS) GO TO 200
P = P + R * H(K+2,J)
H(K+2,J) = H(K+2,J) - P * Z
200 H(K+1,J) = H(K+1,J) - P * Y
H(K,J) = H(K,J) - P * X
210 CONTINUE
J = MIN0(EN,K+3)
C
C Column modification
C
DO 230 I = L, J
P = X * H(I,K) + Y * H(I,K+1)
IF (.NOT. NOTLAS) GO TO 220
P = P + Z * H(I,K+2)
H(I,K+2) = H(I,K+2) - P * R
220 H(I,K+1) = H(I,K+1) - P * Q
H(I,K) = H(I,K) - P
230 CONTINUE
260 CONTINUE
GO TO 70
C
C One root found
C
270 WR(EN) = X + T
WI(EN) = 0.0
EN = NA
GO TO 60
C
C Two roots found
C
280 P = (Y - X) / 2.0
Q = P * P + W
Z = SQRT(ABS(Q))
X = X + T
IF (Q .LT. 0.0) GO TO 320
C
C Real pair
C
Z = P + SIGN(Z,P)
WR(NA) = X + Z
WR(EN) = WR(NA)
IF (Z .NE. 0.0) WR(EN) = X - W / Z
WI(NA) = 0.0
WI(EN) = 0.0
GO TO 330
C
C Complex pair
C
320 WR(NA) = X + P
WR(EN) = X + P
WI(NA) = Z
WI(EN) = -Z
330 EN = ENM2
GO TO 60
C
C Set error -- no convergence to an eigenvalue after 30*n iterations
C
1000 IERR = EN
1001 RETURN
END
| |||||