Relay-Version: version B 2.10 5/3/83; site utzoo.UUCP
Posting-Version: version B 2.10.2 9/18/84; site lanl.ARPA
Path: utzoo!watmath!clyde!burl!ulysses!unc!mcnc!philabs!cmcl2!lanl!jp
From: jp@lanl.ARPA
Newsgroups: net.math
Subject: Re: matrix diagonalization
Message-ID: <25628@lanl.ARPA>
Date: Thu, 9-May-85 14:04:24 EDT
Article-I.D.: lanl.25628
Posted: Thu May  9 14:04:24 1985
Date-Received: Sat, 11-May-85 03:49:16 EDT
References: <1995@burdvax.UUCP>
Sender: newsreader@lanl.ARPA
Organization: Los Alamos National Laboratory
Lines: 40

> 
>   I am looking for a utility program that diagonalizes an NxN matrix
> (complex), ( also giving the eigenvalues and vectors).
> 
> Ed Dougherty
> allegra ! psuvax ! burdvax ! edwardd


Somewhere I saw a message that implied that your matrix was tri-diagonal.
If this is the case there is no need to get out the cannon to kill an ant.
Given a set of equations of the form Y = M X where M is tri diagonal, solve
the first equation for x2 in terms of x1, solve the second equation for x3
in terms of x2 and x1, then substitute the value of x2 in terms of x1 from
the previous equation to get x3 in terms of x1.  Continue in this fashion
until you get to the last two equations.  Both of these relate xn to x1.
Solve them for xn and x1, then substitute the value for x1 in all the expressions
for xi  1<i<n and you have it.   I have used this simple technique for linear
chains of oscillators (coupled microwave cavity accelerators).  The technique
can be simply extended to the case of next-nearest neighbor coupling (5-diag-
onal matrices) by carrying thru the calculation in terms of x1 and x2 and
solving a 4X4 matrix at the end.  The problems I calculated with this method
were high Q resonators Q on the order of 10000 to 30000 with coupling factors
(spectra relative bandwidth) on the order of .02 to .1 depending on the
problem.  

To get the eigenvectors it is easiest to eliminate the "losses" and solve
for the squares of the natural frequencies.  This can be done by Newton's
method.  It is relatively easy to calculate the determinant of a tridiadgonal
matrix (this yields the characteristic polynomial which must be zero at the
eigenvalue)  and I once had an algorithm for similarly calculating the derivative
of the determinant.  The only problem I had with this technique was when the
matrix rank was about 8000.  I ran into some overflow problems because of the
magnitude of some of the intermediate coefficients.  This was solved by a trick
which I remember involved renormalizing the results occasionally.

If you are interested in a further discussion of this topic please contact me
by mail or telephone.  I have quite a bit of experience on this particular
problem and variations of it.

Jim Potter  jp@lanl.arpa   (505) 667-9615