Complex eigenvalues on HP-41

Scientific calculator software to compute complex eigenvalues on HP-41. Commonly used in control, electrical and mechanical engineering.

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 HP-41 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 HP-41.

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

This can be rearranged to

here \(I\) denotes the identity matrix.

If there exists an inverse

then both sides can be left-multiplied 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 non-invertible if and only if its determinant is zero, thus leading to the characteristic equation.

The left-hand 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 non-zero, this means that the matrix \(\lambda I-A\) is singular (non-invertible). This in turn means that its determinant is 0. Thus the roots of the function \(det(\lambda I-A)\) 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\)

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. 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ⁿ 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₁ .. T_n.

\(B\) is removed and the coefficients are calculated and stored in matrix C (\(C\) is re-dimensioned 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_n-1 .. a₀ can now be calculated as

with \(a_n=1\)

The coefficients \(a_{n-1}\ldots a_0\) are stored in matrix \(C\).

Line 166..215 calculate \(a_{n-1} \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₀ t/m a_n-1 (dim: N×2)
• Matrix T contains T₁ 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:

1. The result of two successive iterations matches within 9 digits. Line 765..778 use SCI 8 and RND.
2. 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² – 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 root-seeking algorithm is restricted to real coefficients. That implies that its complex roots come in conjugate pairs, such that if z₁ 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 Cauchy-Riemann equations.

There are two important properties to finding the roots of a polynomial. [2]

1. 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 real-valued 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 Cauchy-Riemann differential equations. These reduce the derivative analysis to a single function, say μ(x,y). It also leads to the orthogonality property.
2. Mapping z→1/z prevents overflow errors by mapping the polynomial function to its inverse:
   q(z)=zⁿ.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 HP-41cv/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

1. Setup the matrix \(C\) and \(W\)
   • By hand
     1. 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).
     2. 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.
     3. Put the text POL in the alpha register
   • By calling KV
     1. Enter the matrix \(MAT\)
     2. XEQ “KV”
2. XEQ “ALLROOT” (when using “KV”, you can just press R/S instead)
3. 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”

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.

Example 2

Input matrix MAT

That implies that the eigenvalues are

Source code

And now for the bad news .. I have the program on magnetic cards, replaced the gummy wheels in the card reader, but still it will not read the tapes.

References

Around the same time, the routine allroot.bas applying Mueller’s parabolic fitting method in the complex plane was published [2]. It was one of those rare moments that I drew inspiration from a piece of Basic code.