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 by a discussions with Michiel Niemijer [1], John Ferman’s article [2] and the routine “ALLROOT” from Ruckdeschel’s book [3].
We implement the algorithm on the HP41 to calculate eigenvalues. To do so, we use functionality from the Advantage ROM. We will not calculate the eigenvectors, due to limitations of this ROM and the limited amount of memory on the HP41.
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 \(\lambda\) is an eigenvalue of A if and only if there is an eigenvector v≠0 such that
??
$$Av=\lambda v Av=\lambda v$$
This can be rearranged to
$$(A\lambda I)v=0$$
here \(I\) denotes the identity matrix.
If there exists an inverse
$$(A\lambda I)^{1}$$
then both sides can be leftmultiplied by it, to obtain v=0. Therefor, if \(\lambda\) is such that \(A\lambda I\) is invertible, \(\lambda\) cannot be an eigenvalue. It can be shown that the converse holds too. If \(A\lambda I\) is not invertible, \(\lambda\) is an eigenvalue.
Linear algebra teaches us that a matrix (here \(A\lambda 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 \(\lambda\), 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 \(\lambda IA\) is singular (noninvertible). This in turn means that its determinant is 0. Thus the roots of the function \(det(\lambda IA)\) are the eigenvalues of \(A\), and it is clear that this determinant is a polynomial in \(\lambda\). [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 \(\lambda\) corresponds to the eigenvalue \(exp(\lambda\) of \(exp(A)\) [6].
The Leibniz formula for determinants expresses the determinant of a square matrix
xxx see https://en.wikipedia.org/wiki/Leibniz_formula_for_determinants
—
https://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
$$x^n+a_{n1}x^{n1}+\ldots a_1x+a_0=0$$
This polynomial has N solutions for x (the eigenvalues). a_{n1} .. a_{0} can now be calculated as
$$a_i=\frac{1}{1n}\sum{a_k\ T_{ki}}$$
with \(a_n=1\)
The coefficients \(a_{n1}\ldots a_0\) are stored in matrix \(C\).
$$\begin{aligned}
C&=\left(\begin{array}{cc}
\Re({a_0}) & \Im({a_0})\\
\vdots & \vdots\\
\Re({a_{n1}}) & \Im({a_{n1}})\\
\end{array}\right) \Rightarrow\\
\\
T&=\left(\begin{array}{cc}
\Re({T_1}) & \Im({T_1})\\
\vdots & \vdots\\
\Re({T_n}) & \Im({T_n})\\
\end{array}\right)
\end{aligned}$$
Line 166..215 calculate \(a_{n1} \ldots a_0\) from \(T_1\ldots 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 – 14x + 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 \(\left(\frac{dp}{dz}\right)^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 \(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 \(2N\).
ALLROOT
will look for roots until \(W\) is filled. If the dimension of \(W\) is less then \(2N\), 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
Example 1
Input matrix MAT using XEQ “MEDIT”
$$\left(\begin{array}{cc}
8 & 4\\
3 & 6
\end{array}\right)$$
Key strokes  Display 

“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
$$\left(\begin{array}{cc}
3.394 & 2.04\times10^{20}\\
10.606 & 4.717\times 10^{18}
\end{array}\right)$$
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.
Example2
Input matrix MAT in
$$\left(\begin{array}{rrr}
8 & 4 & 3\\
4 & 8 & 1 \\
5 & 5 & 7
\end{array}\right)$$
Key strokes  Display 

“MAT”  
CF 02  
XEQ “KV”  0,000 
R/S  0,0000000 00 
XEQ “MEDIT”  4,0010000 00 
FIX 3  4,001 
$$\left(\begin{array}{rrr}
3.446 & 2.505\times10^{21}\\
9.782&3.015\\
9.7825&3.015
\end{array}\right)$$
That implies that the eigenvalues are
$$\left\{\begin{array}{c}
3.436\\
9.782 – 3.015j\\
9.782 + 3.015j
\end{array}\right.$$
Source code
And now for the bad news .. I have the program on magnetic tape, replaced the gummy wheels in the tape reader, but still it will not read the tapes. I will look for the hard copy, and post it later.
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 ALLROOT calculates complex roots from a polynomial with real coefficients 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 https://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.