Solve a 2x2 Eigensystem

Question

For those with a little linear algebra background, the challenge is as simple as this: determine the eigenvalues and eigenvectors of a given complex 2x2 matrix. You may skip ahead to The Challenge for I/O details, etc. For those who need a little refresher on eigensystems, read on.

Background

The characteristic equation of a matrix A is defined by

det| A - λI | = 0

where λ is a complex (scalar) parameter, I is the identity matrix and det|...| is the determinant. The left-hand side evaluates to a polynomial in λ, the characteristic polynomial, which is quadratic in the case of 2x2 matrices. The solutions of this characteristic equation are the eigenvalues of A, which we will denote as λ₁ and λ₂.

Now the eigenvectors v_i of A satisfy

A vi = λi vi

For each λ_i, this gives you a system of two equations in two unknowns (the components of v_i), which can be solved quite easily. You will notice that the system is actually underspecified, and the magnitude of the eigenvectors are not determined by the equations. We will usually want the eigenvectors to be normalised, that is √(|x|² + |y|²) = 1, where x and y are the vector components, |x|² is x multiplied by its complex conjugate.

Note that the eigenvalues may be degenerate, i.e. λ₁ = λ₂. In this case, you may or may not be able to satisfy the single system of equations with two linearly independent eigenvectors.

The Challenge

Given a 2x2 matrix with complex elements, determine its two (possibly identical) eigenvalues and a normalised eigenvector for each eigenvalue. The resulting numbers must be accurate to at least 3 (decimal) significant digits. You may assume that the real and imaginary parts of any matrix element is in the range [-1,1].

You may write a function or a program, taking input via STDIN, command-line argument, prompt or function argument. You may output the result to STDOUT, a dialog box or as the function return value.

You may use any convenient (but unambiguous) string or list format for input and output. You can also choose between pairs of floats or complex types to represent the individual numbers.

You must not use built-in functions for solving eigensystems (like Mathematica's Eigenvectors or Eigensystem) or equation solvers.

This is code golf, so the shortest answer (in bytes) wins.

Examples

Each example is three lines: the input, the eigenvalues and the corresponding eigenvectors in the same order. Note that the eigenvectors are only determined up to their phase, and that in the case of degenerate eigenvalues, the eigenvectors may actually be arbitrary (as in the first example).

[[1.0, 0.0], [0.0, 1.0]]
[1.0, 1.0]
[[1.0, 0.0], [0.0, 1.0]]

[[0.0, 0.4], [-0.1, -0.4]]
[-0.2, -0.2]
[[0.894427, -0.447214], [0.894427, -0.447214]]

[[0.3, 0.1], [0.4, -0.9]]
[-0.932456, 0.332456]
[[-0.0808731, 0.996724], [0.951158, 0.308703]]

[[0.5, -1.0], [0.8, -0.5]]
[0.74162i, - 0.74162i]
[[0.745356, 0.372678 - 0.552771i], [0.745356, 0.372678 + 0.552771i]]

[[-0.0539222 + 0.654836i, -0.016102 + 0.221334i], [0.739514 - 0.17735i, -0.0849216 + 0.77977i]]
[0.238781 + 0.984333i, -0.377625 + 0.450273i]
[[0.313668 + 0.322289i, 0.893164], [-0.236405 - 0.442194i, 0.865204]]

[[-0.703107 - 0.331792i, 0.286719 - 0.587305i], [-0.418476 + 0.396347i, -0.885934 + 0.50534i]]
[-1.13654 - 0.32678i, -0.4525 + 0.500329i]
[[0.833367, -0.248208 - 0.493855i], [-0.441133 - 0.408236i, 0.799215]]

[[-0.156312 + 0.788441i, 0.045056 - 0.579167i], [0.130741 - 0.97017i, 0.049183 - 0.590768i]]
[-0.181759 + 1.11738i, 0.0746298 - 0.919707i]
[[0.86955, -0.493846 + 0.000213145i], [0.318856 - 0.0181135i, 0.94763]]

Answer 1

MATLAB, 91

A standard technique to obtain a normalized vector and remove the useless degree of freedom is representing the elements of the vector as the cosine and sine of some angle.

I originally tried to code in Python, but its math handling proved to be too brain-damaged. Its math functions refused to accept complex values, and it does not understand that floating-point division by zero is OK.

function[]=f(a,b,c,d)
L=(a+d+[1,-1]*((a-d)^2+4*b*c)^.5)/2
t=atan((L-a)/b);v=[cos(t);sin(t)]

First the two eigenvalues are printed under the heading L =. Then two column vectors are printed under the corresponding values of L, under v =. The code may fail to give linearly independent vectors in cases where it is possible to do so (such a program normally would be considered broken), but Martin said it is not required.

Answer 2

Python 2, 198 bytes

a,b,c,d=input()
H=(a+d)/2
D=(H*H-a*d+b*c)**.5
X,Y=H+D,H-D
p,q,r,s=[[1,0,0,1],[b,X-a,b,Y-a],[X-d,c,Y-d,c]][2*(c!=0)or(b!=0)]
A=abs
V=A(A(p)+A(q)*1j)
W=A(A(r)+A(s)*1j)
print[X,Y],[[p/V,q/V],[r/W,s/W]]

Input is a flat list of 4 complex numbers via STDIN, e.g.

[0.0+0j, 0.4+0j, -0.1+0j, -0.4+0j]

Note that Python uses j instead of i for complex numbers.

Output is two lists, the first containing the eigenvalues and the second containing the eigenvectors, e.g.

[(-0.2+0j), (-0.2+0j)]
[[(0.8944271909999159+0j), (-0.4472135954999579+0j)], [(0.8944271909999159+0j), (-0.4472135954999579+0j)]]

(newline inserted for clarity)

Answer 3

Lua, 462 455 431 427 bytes

There is no built-in complex math in Lua. No vector operations either. All had to be rolled by hand.

a,b,c,d,e,f,g,h=...x=math.sqrt z=print i=a-g j=b-h
k=(i^2-j^2)/2+2*(c*e-d*f)m=x(k^2+(i*j+2*(c*f+d*e))^2)n=x(m+k)o=x(m-k)i=(a+g+n)/2
j=(b+h+o)/2 k=(a+g-n)/2 l=(b+h-o)/2 z(i,j,k,l)q=c^2+d^2 r=e^2+f^2 s=q+r if s==0
then z(1,0,0,0,0,0,1,0)else if r==0 then m,n,o,p=c,d,c,d c,d=i-a,j-b e,f=k-a,l-b
u=x(q+c^2+d^2)v=x(q+e^2+f^2)else m,n=i-g,j-h o,p=k-g,l-h c,d=e,f
u=x(r+m^2+n^2)v=x(r+o^2+p^2)end z(m/u,n/u,o/v,p/v,c/u,d/u,e/v,f/v)end

Run from the command-line with the following arguments:

lua eigen.lua Re(a) Im(a) Re(b) Im(b) Re(c) Im(c) Re(d) Im(d)

Produces the following output:

Re(lambda1) Im(lambda1) Re(lambda2) Im(lambda2)
Re(v11) Im(v11) Re(v12) Im(v12) Re(v21) Im(v21) Re(v22) Im(v22)

...for a,b,c,d the 4 components of the input matrix, lambda1 and lambda2 the two eigenvalues, v11,v21 the first unit eigenvector, and v12,v22 the second unit eigenvector. For example,

lua eigen.lua 1 0  1 0  1 0  0 0

...produces...

1.6180339887499 0   -0.61803398874989   0
0.85065080835204    0   -0.52573111211913   0   0.52573111211913    0   0.85065080835204    0