Why the results are different while using eig() to solve syms and double Jacobian matrix?

Dearl all,

When I use eig to find the eigenvalues for a Jacobian matrix, the outputs are different if I define the Jacobian matrix as syms and double types. The order of each values are different in particular. I appreciate any clues for this difference:))

Many thanks,

JL

Here is an example:

% Create a syms J matrix and a dounble J matrix
JacobianSyms = vpa([1, 2, 3; 2, 1 8; 3, 8 9]);
JacobianDoub = double(JacobianSyms);

% Calculate their eigenvalues
[evs, eigenvaluesfromJacobianSyms] = eig(JacobianSyms)

evs = 

outPutfromJacobianSyms = 

[evd, eigenvaluesfromJacobianDoub] = eig(JacobianDoub)

evd = 3×3

-0.0294 0.9675 0.2513 0.8534 -0.1066 0.5103 -0.5205 -0.2294 0.8225

outPutfromJacobianDoub = 3×3

-3.9479 0 0 0 0.0681 0 0 0 14.8798

Answer 1

[Bruno Luong]

If I ask a set of 3 integer numbers > 1 so that 30 is the product, what do you tell me?

(5,3,2) or
(5,2,3) or
(3,5,2) or
(3,2,5) or
(2,5,3) or
(2,3,5)

?

It's the same thing with eigen value and vectors, the order is just arbitrary set. So when you ask two different implementations for eigen values/vectors, you'll get the order that is arbitrary set by this routine.

[Juntong Lai]

Hi Bruno,

Thank you very much for your answer!

I understand what you mean but I am still curious about what routine results in the arbitrary order, and I am also wondering if it is possible to make the orders consistent.

ATB,

JL

[Bruno Luong]

" it is possible to make the orders consistent."

There is no such thing as orders consistent, since the consistent is not universally defined (yet, or never).

Is (5,3,2) more consistent than (3,2,5)?

Up to you, you can sort eigen value in increasing order for example (after getting them from both routines) and shufle accordingly eigen vectors with the same permutation. If you are unlucky there is multiple-order eigen values (same value but appear multiple times) in this case eigen vectors are arbitrary (even up to scaling).

[Juntong Lai]

Hi Bruno,

In your case, it indeed makes sense. However, I was doing a linearisation of nonlinear equations so it is important that the order is consistent with the the oder of Jacobian matrix. Each eigenvalue represents the convergency of each variable in its dimension. Back to the way to calculate the eigenvalues, if I use the double type for the jacobian matrix then the results make more sense as they are consistent. In turn, the arbitrary order of the eigenvalues makes the solutions confused. Thus, I was managing to see if I can find the mechanism that leads to the arbitrary order when using syms matrix.

I appreciate your clarification but I am wondering if you know what is the exact routine shufling eigenvalues when using syms matrix compared to using double matrix.

Many thanks,

JL

[Torsten]

However, I was doing a linearisation of nonlinear equations so it is important that the order is consistent with the the oder of Jacobian matrix. Each eigenvalue represents the convergency of each variable in its dimension.

You want to say that the eigenvalues are somehow connected with the variables of the nonlinear equations ? Can you elaborate on this ? Because I don't understand what you mean.

% Create a syms J matrix and a dounble J matrix
JacobianSyms = vpa([1, 2, 3; 2, 1 8; 3, 8 9]);
JacobianDoub = double(JacobianSyms);
% Calculate their eigenvalues
[evs, eigenvaluesfromJacobianSyms] = eig(JacobianSyms);
evs*eigenvaluesfromJacobianSyms*evs.'

ans = 

[evd, eigenvaluesfromJacobianDoub] = eig(JacobianDoub);
evd*eigenvaluesfromJacobianDoub*evd.'

ans = 3×3

1.0000 2.0000 3.0000 2.0000 1.0000 8.0000 3.0000 8.0000 9.0000

Thus order of eigenvalues and eigenvectors does not matter to reproduce the original matrix.

[Bruno Luong]

Usualy the symbolic solves the roots of the characteristic polynomial, in case of third order (3 x 3) matrix the "arbitrary" order is the order the formal cardano roots formulae is written down like for second order polynomial we select for 2x2 matrix reduced to 2nd order polynomial

-b-sqrt(delta)/(2*a) first then
-b+sqrt(delta)/(2*a) second

and traditionaly mathematicians decide by convention that sqrt(delta) or delta^(1/n) returns the solution with smallest complex argument.

I don't know how the express the order numerical method in simple term. Most of numerical method uses QR method of Hessenberg reduction, usualy it converge to eigen values in decreasing order of their abs values, but some of the accelerating procedure can raise out different order https://article.sciencepublishinggroup.com/html/10.11648.j.pamj.20160504.15.html

So you see that the eigen order is not completely arbitrary, there is some dependency of method used by the specific algorithm, but it is quite complexe and fragile that user should NOT assume any consistent order between methods used.

We are not here to assist you to interpret the eigen values if they have some physical meaing related to your Jacobian ODE, so the order lust reflect this physical interpretation.

We are here to explain to you why you should not expect two eig functions give the same consistent order. There is no such thing mathematically.

Regarless for a given matrix A all the algorithms give D, V such that

A*V = V*D

or

A = V*D*inv(V)

For any permutation vector p (of 1:n)

Vp = V(:,p)

Dp = D(p,p)

is still a valid decomposition of A

A = Vp*Dp*inv(Vp)

In your case if |V,D] is resp [evs, eigenvaluesfromJacobianSyms], p = [3 2 1], Vp,Dp match [evd, eigenvaluesfromJacobianDoub]. They match up too the permutation.

[Juntong Lai]

Hi Bruno,

Thank you very much for your kind explaination!

I think I am going to stick with using the jacobian matrix of the double type other than syms, but your further clarification clears my confusion.

In terms of Each eigenvalue represents the convergency of each variable in its dimension, you are right, Torsten. I am sorry for a bit of confusion above, please let me make the statement clearer.

In a linearisation of nonlinear ordinary differential equations, the eigenvalues can indicate the convergency of each variable in the corresponding directions (eigenvectors) in a phase plane, i.e. Ideally, the value of each variable will be convergent to zero if the eigenvalue is negative, and it will be away from zero if the eigenvalue is positive. So far, it does not matter how the order of eigenvalues is. However, it becomes important for me to ensure the order consistent if there is an additional dimension that is independently convergent. Namely, Torsten you are completely right in the case I came up with, but there will be a noise eigenvalue (always negative) to interfere with indexing the other eigenvalues to see whether they are negative/positive/complex values. I can traverse all of them but it will be much easier to analise the others if I keep the noise eigenvalue at the end of the matrix. That is the 'order' I meant essentially.

For instance:

% Create a syms J matrix and a dounble J matrix
JacobianSyms = vpa([1, 2, 3, 9, 9; ...
    2, 1, 8, 9, 9; ...
    3, 8, 9, 9, 9; ...
    9, 8, 3, 9, 9; ...
    0, 0, 0, 0, -1,]);
JacobianDoub = double(JacobianSyms);

% Calculate their eigenvalues
[evs, eigenvaluesfromJacobianSyms] = eig(JacobianSyms)

evs = 

eigenvaluesfromJacobianSyms = 

[evd, eigenvaluesfromJacobianDoub] = eig(JacobianDoub)

evd = 5×5

-0.3297 -0.3199 0.5400 0.7083 0.1188 -0.4520 0.2928 0.5967 -0.4953 -0.0198 -0.6208 0.7674 -0.0680 0.3967 0.0000 -0.5492 -0.4723 -0.5896 -0.3093 -0.7129 0 0 0 0 0.6909

eigenvaluesfromJacobianDoub = 5×5

24.3799 0 0 0 0 0 5.2621 0 0 0 0 0 -6.9940 0 0 0 0 0 -2.6480 0 0 0 0 0 -1.0000

In this example, -1.0 would be the noise eigenvalue, and I do not want to analyse it as I know it will be always a negative real value. Thus, the variable that may not be convergent at a point will be possibly regarded as convergent to the point if -1.0 appears in the middle of the eigenvalue matrix, though the above example can not mislead the judgement as there is no complex value here and the stability is quite obvious.

Again, thank you both for answering my questions and it is now well solved.

ATB,

JL

[Bruno Luong]

BTW if you want to take the control and impose your order, for example sorted on real part of eigen value just do this

[V,D] =eig(A)
[~,p] = sort(real(diag(D)));
V = V(:,p)
D = D(p,p)