.

Complex Eigenvalues for HP-41cv/cx

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 HP-71 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 HP-41 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 HP-41, 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 left-multiplied 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 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 λ, 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 λI-A is singular (non-invertible).  This in turn means that its determinant is 0.  Thus the roots of the function det(λI-A) 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 (a1,1+a2,2+..+an,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 .. MATn where n is the number of columns (=rows).

Line 088..103, determine the sum of all diagonal elements of MATI.

Line 143..155, stores the sums in trace matrix T. These are called the traces T1 .. Tn.

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). an-1 .. a0 can now be calculated as
with an=1.

The coefficients an-1 .. a0 are stored in matrix C

Line 166..215 calculate an-1 t/m a0 from T1 t/m Tn.

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 a0 t/m an-1 (dim: N×2)
  • Matrix T contains T1 t/m Tn (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: x2 – 14*x + 36 = 0 and the eigenvalue found is 3,394 then it is removed using:

Because x2 – 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 root-seeking algorithm is restricted to real coefficients.  That implies that its complex roots come in conjugate pairs, such that if z1 is a root, its complex conjugate conj(z1) 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  (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 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)=zn.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 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 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.
      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

  1. 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.
    :

  2. 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
PPC-Journal V13N1, page 22-23
[3] Basic Scientific Subroutines, Vol. II, chapter 7
F.R. Ruckdeschel, 1981
McGraw-Hill Publications, New York
[4] Muller’s Method, a generalized secant method of root finding by using quadratic 3-point 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

Coert Vonk

Independent Firmware Engineer at Los Altos, CA
Welcome to the things that I couldn’t find.This blog shares some of the notes that I took while deep diving into various fields.Many such endeavors were triggered by curious inquiries from students. Even though the notes often cover a broader area, the key goal is to help the them adopt, flourish and inspire them to invent new technology.
Coert Vonk

Latest posts by Coert Vonk (see all)

Leave a Reply

You can use these HTML tags

<a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <s> <strike> <strong>

  

  

  

Protected with IP Blacklist CloudIP Blacklist Cloud