Jacobi's rotation method

Directly to the input form

• In order to find all eigenvalues (and eigenvectors) of a real symmetric matrix A one may use Jacobi's rotation method, which is one of the oldest methods in use nowadays.
• The idea of the method is quite simple: by using a sequence of unitary similarity transforms one reduces the sum of squares of off diagonal elements of the matrix, such that in the limit a unitary similarity transformation to diagonal form is obtained.
• These transformations are two dimensional rotations in a (p,q)- plane of Rⁿ which annihilate the element (p,q) and, of course, (q,p). Previously generated zero elements are overwritten by this, necessarily, but the sum of squares of off diagonal elements decreases in the strict sense.
• By construction each accumulation point of the sequence of products of these rotations is unitary and a complete eigensystem of the given matrix. A.
• Convergence is of second order in general, more precisely
  t_k+1 <= C t_k²
  for some constant C depending on the eigenvalue distribution, where
  t= ( ( ∑_{over all p < q} (a_p,q² ) / (n(n-1)) ) ^1/2.
  if one counts one full sweep over all off diagonal elements as one step.
• Here we use the so called threshhold value version where over one sweep only those elements are annihilated, which are in absolute value larger than a fraction of t_k. Afterwards, the treshhold value is reduced by some factor < 1.
• Due to the application of unitary transformations only the error in the computed eigenvalues will be limited rouhgly by ''number of steps times matrix norm times machine precision ''.
• This is code translated from Rutishausers code published in the Handbook Linear Algebra (Wilkinson , Reinsch)

Your input

• You may take one from four cases:
  1. A 8x8 case with a cluster of large eigenvalues (exactly known).
  2. A quindiagonal matrix which occurs in discretizing the differential equation of the bending of a bar. Here the dimension is variable. The eigenvalues are known exactly.
  3. You may choose a set of eigenvalues. With the help of an internally defined unitary matrix a full symmetric matrix is computed then and used as input. So you can check the result and the effort.
  4. You may type in the lower triangle of a matrix (thought as symmetric). Input is then rowwise, the entries separated by blank or comma. The code then computes the eigenvalues and eigenvectors for you. Since the exact eigenvalues are not known, you get no error bound then . In principle the value (n(n-1))^1/2 × t_k provides such an error bound, but this neglects the roundoff error in performing the transformations, which might be larger than this.
• The dimension n of the matrix in the cases 2,3,4 above. For case 2 there must hold 3 <= n <= 30, otherwise 2 <= n <= 30 !
• You can have a printed output of A .
• You can have a printed output of the eigensystem.

Output

• If desired the input A and/or the eigensystem.
• The number of rotations applied.
• The eigenvalues found and their errors, as far as the exact eigenvalues are known
• A plot which shows the decay of the sum of squares of off diagonal elements over the step number.

Questions

• How influences the eigenvalue distribution the speed of convergence ?
• If you vary some matrix element (symmetrically), how change the eigenvalues and how the eigenvectors?

The input form