Matrices form a fundamental part of electrical networks. Components such as capacitors and inductors introduce complex number arithmetic. Other fields where complex eigenvalues are common place are control engineering and mechanical engineering. Even when the matrix only has real numbers, there might be complex eigenvalues. One example is a clock, where the imaginary part of the eigenvalue is the frequency of the pendulum.
History
This program was inspired discussions with Michiel Niemijer [1] and two publications. First John Ferman’s article “Eigenvalues & Eigenvectors by Brute Force using the MathLex ROM” [2], that was written for the HP71 and determines the characteristic polynomial of the matrix, and then uses a SOLVE. The other publication was Ruckdeschel’s book “Basic Scientific Subroutines” [3] where he describes a program (ALLROOT) to calculates complex roots from a polynomial with real coefficients.
This article describes an implementation on the HP41 using functionality from the Advantage ROM. Given the limitation of the Advantage ROM compared to the MathLex ROM, and the limited amount of memory on the HP41, this program only determines eigenvalues.
Method
The program first determines the characteristic polynomial. An eigenvalue is determined and removed from the equation. Remaining eigenvalues are determined until all are found.
Because the algorithm also finds complex roots, it takes a bit longer compared to routines that only find real roots. Matrices with double eigenvalues pose some additional challenges, but the algorithm will always finds the roots.
The program consists of two parts. The first part (KV) determines the characteristic polynomial. The second part (ALLROOT) determines the roots of the polynomial (the eigenvalues).
Determining the Characteristic Polynomial
Given a square matrix A, we want to find a polynomial whose roots are precisely the eigenvalues of A. A scalar λ is an eigenvalue of A if and only if there is an eigenvector v≠0 such that
Av=λv
This can be rearranged to
(AλI)v=0
Here I denotes the identity matrix.
If there exists an inverse
(AλI)^{1}
then both sides can be leftmultiplied by it, to obtain v=0. Therefor, if λ is such that AλI is invertible, λ cannot be an eigenvalue. It can be shown that the converse holds too. If AλI is not invertible, λ is an eigenvalue.
Linear algebra teaches us that a matrix (here AλI) is noninvertible if and only if its determinant is zero, thus leading to the characteristic equation.
The lefthand size of this equation can be seen to be a polynomial function in λ, whose coefficients depend on the entries of A. This polynomial is called the characteristic equation. Its degree is n. It also can be shown that any n×n matrix has at most n eigenvalues. [7]
—
Since v is nonzero, this means that the matrix λIA is singular (noninvertible). This in turn means that its determinant is 0. Thus the roots of the function det(λIA) are the eigenvalues of A, and it is clear that this determinant is a polynomial in λ. [5]
Trace tr(A) is by definition the sum of the diagonal entries of A (a_{1,1}+a_{2,2}+..+a_{n,n}) and also equals the sum of the eigenvalues. Thus, for complex matrix A
det(exp(A))=exp(tr(A))
Here exp() denotes the matrix exponential of A, because every eigenvalue λ corresponds to the eigenvalue exp(λ) of exp(A) [6].
The Leibniz formula for determinants expresses the determinant of a square matrix
xxx see http://en.wikipedia.org/wiki/Leibniz_formula_for_determinants
—
http://en.wikipedia.org/wiki/Companion_matrix
The function “KV” determines the characteristic polynomial of the square matrix “MAT”.
Line 135..142 will determine MAT .. MAT^{n} where n is the number of columns (=rows).
Line 088..103, determine the sum of all diagonal elements of MAT^{I}.
Line 143..155, stores the sums in trace matrix T. These are called the traces T_{1} .. T_{n}.
B is removed and the coefficients are calculated and stored in matrix C (C is redimensioned to N×2).
Matrix W to store the roots is created (dimension N×2).
The characteristic polynomial of the matrix follows as
This polynomial has N solutions for x (the eigenvalues). a_{n1} .. a_{0} can now be calculated as
with a_{n}=1.
The coefficients a_{n1} .. a_{0} are stored in matrix C
Line 166..215 calculate a_{n1} t/m a_{0} from T_{1} t/m T_{n}.
Upon completion
 R00: “MAT”
 R01: 0.000
 R02: scratch
 R03: scratch
 Alpha: “POL”
 EM: de matrices MAT, C, T, W
 Matrix MAT is unchanged and can be discarded
 Matrix C contains the coefficients a_{0} t/m a_{n1} (dim: N×2)
 Matrix T contains T_{1} t/m T_{n} (dim: N×2) and can be discharged
 Matrix W is still empty (dim: N×2), will be filled by ALLROOT
Determine eigenvalue and remove it from the equation
ALLROOT is an iterative algorithm that attempts to calculate the roots of a given series or function by repeatedly using the Muller’s Method [4] in the complex plane, and removing the roots already found by division. It will iterate until one of the following conditions is satisfied:
 The result of two successive iterations matches within 9 digits. Line 765..778 use SCI 8 and RND.
 There were over 100 iterations. This may occur with large characteristic polynomials or when the eigenvalues are close such as triple eigenvalues. The program will stop with a BEEP and display “TE MOEILIJK”. The X and Y register will have the approximation of the root. If you are satisfied with this, you can press R/S and the program with continue as if that is the exact root. Otherwise, you can find the results so far. In that R00 has the total number of roots found. The roots themselves can be found viewed by using XEQ “MEDIT” on matrix W.
When a root is found it will sound a TONE 90. When ALLROOT is finished it will sound a BEEP.
 line 220..224 are for the iteration.
 line 225..261 find an appropriate complex starting value.
 line 262..279 are again for the iteration.
 LBL 99 searches for a root and uses POL.
 POL evaluates the quotient of the characteristic polynomial and the polynomial of the roots found. Because prior roots are not removed from the original polynomial, it greatly reduces error accumulation.
Upon completion
 R00: number of roots found
 R01..25: scratch
 R26: “POL” (value from alpha at the time ALLROOT was called).
 F00: set, when not all roots were found by satisfying criteria 1.
 Alpha: “W”, so you can view the roots by calling XEQ “MEDIT”.
 The display mode is SCI] 8
 Matrix matrix C is unchanged.
Assume the characteristic polynomial is: x^{2} – 14*x + 36 = 0 and the eigenvalue found is 3,394 then it is removed using:
Because x^{2} – 14*x + 36 is also zero for the second eigenvalue, the error in the prior eigenvalue barely matters. Refer to the examples for more details.
The rootseeking algorithm is restricted to real coefficients. That implies that its complex roots come in conjugate pairs, such that if z_{1} is a root, its complex conjugate conj(z_{1}) is also a root. Polynomial with complex coefficients can be rewritten as two independent polynomials with real number coefficients using CauchyRiemann equations.
There are two important properties to finding the roots of a polynomial. [2]
 Polynomials are always analytic in the entire complex plane, because all its derivatives (dp/dz)^{n} exist. Therefore, it can be represented by a Taylor series from which the real and imaginary parts can be separated
p(z)=μ(x,y)+j.ν(x,y), where μ(x,y) and ν(x,y) are realvalued functions.
Now instead of finding a value of complex value z so that p(z) is zero, we can find real values of x and y that simultaneously make μ(x,y) and ν(x,y) zero. Applying the partial derivatives, we get the CauchyRiemann differential equations. These reduce the derivative analysis to a single function, say μ(x,y). It also leads to the orthogonality property.  Mapping z→1/z prevents overflow errors by mapping the polynomial function to its inverse:
q(z)=z^{n}.p(1/z)
It can be shown that q(z) is also a polynomial of degree n, and that its coefficients the same as those of p(z), but with the order reversed.
Instructions
You may want to experiment with real eigenvalues using “Eigenvalues for HP41cv/cx“.
The input matrix should be square. Matrices with real coefficients should be entered using MEDIT, after which F02 should be cleared. For complex coefficients, CMEDIT should be used and F02 be set. The name of the matrix should be placed in the alpha register. The names can not be “b”, “c”, “t” or “w” because those are used internally by the program. The matrix may be stored in regular register above R26.
The text “POL” in the alpha register is also needed when calling ALLROOT. When MAT has only real coefficients, its size should be be N*N. For complex coefficients, it should be 2N*2N.
To find the roots of a polynomial, you can
 Setup the matrix C and W
 By hand
 Create a matrix C size N*2, and fill in the coefficients. Even if the coefficients are only real number, you should enter it as if it were complex numbers (imaginary part 0).
 Create a matrix W with a size between 1*2 and N*2. ALLROOT will look for roots until W is filled. If the dimension of W is less then N*2, not all roots will be determined.
 Put the text “POL” in the alpha register
 By calling “KV”
 Enter the matrix “MAT”
 XEQ “KV”
 XEQ “ALLROOT” (when using “KV”, you can just press R/S instead)
 View the roots using XEQ “MEDIT” (with “W” in alpha).
To view the matrices C, T or W you should always use XEQ “MEDIT” (not XEQ “CMEDIT”).
Examples
 Input matrix MAT using XEQ “MEDIT”
“MAT” CF 02 XEQ “KV” 0.000 (in alpha “POL”) R/S 0,000000000 (in alpha “W”) The MAT, C, T, en W are now in place. By calling XEQ “MEDIT”, we can inspect the roots
The two eigenvalues are 3,394 en 10,606 when the imaginary parts are neglected given that they are a magnitudes smaller as the real parts.
:  Input matrix MAT in
“MAT” CF 02 XEQ “KV” 0,000 R/S 0,0000000 00 XEQ “MEDIT” 4,0010000 00 FIX 3 4,001
That implies that the eigenvalues are 3,436
 9,782 – 3,015j
 9,782 + 3,015j
Listing
And now for the bad news .. I printed the program at the time, but lost the printout. The magnetic cards are a 10 hour flight away, but once I get them, I will include the code here.
References
[1]  Correspondence met Michiel Niemijer  
[2]  Eigenvalues & Eigenvectors by brute force using the Math Mex ROM John Ferman PPCJournal V13N1, page 2223 

[3]  Basic Scientific Subroutines, Vol. II, chapter 7 F.R. Ruckdeschel, 1981 McGrawHill Publications, New York 

[4]  Muller’s Method, a generalized secant method of root finding by using quadratic 3point interpolation David Muller, 1958 Wikipedia 

[5]  Characteristic polynomial Wikipedia 

[6]  Determinant, relation to eigenvalues and trace Wikipedia 

[7]  Eigenvalues and eigenvectors, matrices Wikipedia 
see http://en.wikipedia.org/wiki/Characteristic_polynomial
Around the same time, the routine allroot.bas applying Meuller’s parabolic fitting method in the complex plane was published [2]. Perhaps it was one of those rare moments that a piece of basic code inspired me.
Coert Vonk
Latest posts by Coert Vonk (see all)
 Arduino Stuff  20161001
 ESP8266 reads Google Calendar  20160920
 Soldering with Kids  20160908
Leave a Reply