Quantcast

Documentation Center

  • Trial Software
  • Product Updates

Eigenvalue Trajectories

This example applies several numeric, symbolic, and graphic techniques to study the behavior of matrix eigenvalues as a parameter in the matrix is varied. This particular setting involves numerical analysis and perturbation theory, but the techniques illustrated are more widely applicable.

In this example, you consider a 3-by-3 matrix A whose eigenvalues are 1, 2, 3. First, you perturb A by another matrix E and parameter . As t increases from 0 to 10-6, the eigenvalues , , change to , , .

This, in turn, means that for some value of , the perturbed matrix A(t) = A + tE has a double eigenvalue . The example shows how to find the value of t, called τ, where this happens.

The starting point is a MATLAB® test example, known as gallery(3).

A = gallery(3)
A =
  -149      -50     -154
   537      180      546
   -27       -9      -25

This is an example of a matrix whose eigenvalues are sensitive to the effects of round-off errors introduced during their computation. The actual computed eigenvalues may vary from one machine to another, but on a typical workstation, the statements

format long
e = eig(A)

produce

e =
   1.000000000010722
   1.999999999991790
   2.999999999997399

Of course, the example was created so that its eigenvalues are actually 1, 2, and 3. Note that three or four digits have been lost to round-off. This can be easily verified with the toolbox. The statements

syms x;
B = sym(A);
e = eig(B)'
p = charpoly(B, x)
f = factor(p)

produce

e =
[1,  2,  3]

p =
x^3 - 6*x^2 + 11*x - 6

f =
(x - 3)*(x - 1)*(x - 2)

Are the eigenvalues sensitive to the perturbations caused by round-off error because they are "close together"? Ordinarily, the values 1, 2, and 3 would be regarded as "well separated." But, in this case, the separation should be viewed on the scale of the original matrix. If A were replaced by A/1000, the eigenvalues, which would be .001, .002, .003, would "seem" to be closer together.

But eigenvalue sensitivity is more subtle than just "closeness." With a carefully chosen perturbation of the matrix, it is possible to make two of its eigenvalues coalesce into an actual double root that is extremely sensitive to round-off and other errors.

One good perturbation direction can be obtained from the outer product of the left and right eigenvectors associated with the most sensitive eigenvalue. The following statement creates the perturbation matrix:

E = [130,-390,0;43,-129,0;133,-399,0]
E =
   130  -390     0
    43  -129     0
   133  -399     0

The perturbation can now be expressed in terms of a single, scalar parameter t. The statements

syms x t
A = A + t*E

replace A with the symbolic representation of its perturbation:

A = 
[130*t - 149, - 390*t - 50, -154]
[ 43*t + 537,  180 - 129*t,  546]
[ 133*t - 27,  - 399*t - 9,  -25]

Computing the characteristic polynomial of this new A

p = charpoly(A, x)

gives

p =
x^3 + (- t - 6)*x^2 + (492512*t + 11)*x - 1221271*t - 6

p is a cubic in x whose coefficients vary linearly with t.

It turns out that when t is varied over a very small interval, from 0 to 1.0e–6, the desired double root appears. This can best be seen graphically. The first figure shows plots of p, considered as a function of x, for three different values of t: t = 0, t = 0.5e–6, and t = 1.0e–6. For each value, the eigenvalues are computed numerically and also plotted:

x = .8:.01:3.2;
for k = 0:2
  c = sym2poly(subs(p,t,k*0.5e-6));
  y = polyval(c,x);
  lambda = eig(double(subs(A,t,k*0.5e-6)));
  subplot(3,1,3-k)
  plot(x,y,'-',x,0*x,':',lambda,0*lambda,'o')
  axis([.8 3.2 -.5 .5])
  text(2.25,.35,['t = ' num2str( k*0.5e-6 )]);
end
Warning: Imaginary parts of complex X and/or Y arguments ignored 

The bottom subplot shows the unperturbed polynomial, with its three roots at 1, 2, and 3. The middle subplot shows the first two roots approaching each other. In the top subplot, these two roots have become complex and only one real root remains.

The next statements compute and display the actual eigenvalues

e = eig(A);
ee = subexpr(e)
sigma = 
(1221271*t)/2 + (t + 6)^3/27 - ((492512*t + 11)*(t + 6))/6 + (((492512*t)/3 - (t + 6)^2/9 + 11/3)^3...
                    + ((1221271*t)/2 + (t + 6)^3/27 - ((492512*t + 11)*(t + 6))/6 + 3)^2)^(1/2) + 3
 
ee =
                                                                                       t/3...
                       + sigma^(1/3) - ((492512*t)/3 - (t + 6)^2/9 + 11/3)/sigma^(1/3) + 2
 t/3 - sigma^(1/3)/2 - (3^(1/2)*(sigma^(1/3) + ((492512*t)/3 - (t + 6)^2/9 + 11/3)/sigma^(1/3))*i)/2...
                                           + ((492512*t)/3 - (t + 6)^2/9 + 11/3)/(2*sigma^(1/3)) + 2
 t/3 - sigma^(1/3)/2 + (3^(1/2)*(sigma^(1/3) + ((492512*t)/3 - (t + 6)^2/9 + 11/3)/sigma^(1/3))*i)/2...
                                           + ((492512*t)/3 - (t + 6)^2/9 + 11/3)/(2*sigma^(1/3)) + 2
pretty(ee)

showing that e(2) and e(3) form a complex conjugate pair:

/    t        1/3            \
|    - + sigma    - #3 + 2   |
|    3                       |
|                            |
|          1/3               |
| t   sigma                  |
| - - -------- - #2 + #1 + 2 |
| 3       2                  |
|                            |
|          1/3               |
| t   sigma                  |
| - - -------- + #2 + #1 + 2 |
\ 3       2                  /

where

                           2
         492512 t   (t + 6)    11
         -------- - -------- + --
             3          9       3
   #1 == ------------------------
                       1/3
                2 sigma

                       1/3
         sqrt(3) (sigma    + #3) i
   #2 == -------------------------
                     2

                           2
         492512 t   (t + 6)    11
         -------- - -------- + --
             3          9       3
   #3 == ------------------------
                      1/3
                 sigma

Next, the symbolic representations of the three eigenvalues are evaluated at many values of t to produce a plot of their trajectories.

tvals = (2:-.02:0)' * 1.e-6;
r = size(tvals,1);
c = size(e,1);
lambda = zeros(r,c);
for k = 1:c
   lambda(:,k) = double(subs(e(k),t,tvals));
end
plot(lambda,tvals)
xlabel('\lambda'); ylabel('t');
title('Eigenvalue Transition')
Warning: Imaginary parts of complex X and/or Y arguments ignored 

Above t = 0.8e-6, the graphs of two of the eigenvalues intersect, while below t = 0.8e–6, two real roots become a complex conjugate pair. What is the precise value of t that marks this transition? Let τ denote this value of t.

One way to find the exact value of τ involves polynomial discriminants. The discriminant of a quadratic polynomial is the familiar quantity under the square root sign in the quadratic formula. When it is negative, the two roots are complex.

There is no discrim function in the toolbox, but there is the polylib::discrim function in the MuPAD® language.

Use these commands

syms a b c x
evalin(symengine,'polylib::discrim(a*x^2+b*x+c, x)')

to show the generic quadratic's discriminant, b2 - 4ac:

ans =
b^2 - 4*a*c

The discriminant for the perturbed cubic characteristic polynomial is obtained, using

discrim = feval(symengine,'polylib::discrim',p,x)

which produces

discrim =
242563185060*t^4 - 477857003880091920*t^3 +...
 1403772863224*t^2 - 5910096*t + 4

The quantity τ is one of the four roots of this quartic. You can find a numeric value for τ with the following code.

s = solve(discrim);
tau = vpa(s)
tau =
 1970031.04061804553618913725474883634597991201389
 0.000000783792490596794010485879469854518820556090553664
 0.00000107692481604921513807537160160597784208236311263 -...
 0.00000308544636502289065492747*i
 0.00000308544636502289065492746538275636180217710757295*i +...
 0.00000107692481604921513807537160160597784249167873707

Of the four solutions, you know that

tau = tau(2)

is the transition point

tau =
0.00000078379249059679401048084

because it is closest to the previous estimate.

A more generally applicable method for finding τ is based on the fact that, at a double root, both the function and its derivative must vanish. This results in two polynomial equations to be solved for two unknowns. The statement

sol = solve(p,diff(p,'x'))

solves the pair of algebraic equations p = 0 and dp/dx = 0 and produces

sol = 
    t: [4x1 sym]
    x: [4x1 sym]

Find τ now by

format short
tau = double(sol.t(2))

which reveals that the second element of sol.t is the desired value of τ:

tau =
  7.8379e-007

Therefore, the second element of sol.x

sigma = double(sol.x(2))

is the double eigenvalue

sigma =
    1.5476

To verify that this value of τ does indeed produce a double eigenvalue at , substitute τ for t in the perturbed matrix A(t) = A + tE and find the eigenvalues of A(t). That is,

e = eig(double(subs(A, t, tau)))
e =
    1.5476
    1.5476
    2.9048

confirms that is a double eigenvalue of A(t) for t = 7.8379e–07.

Was this topic helpful?