## Documentation Center |

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, *b*^{2} -
4*ac*:

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?