# What's the difference of using different state-space forms for solving eigenvalue problem?

I have a set of equations of motion describing a planetary gear train of 18 DoF (sun, 3 planets, carrier and ring), they have the general form of:

$$\mathbf{M{\ddot q}+{\Omega_c}G{\dot q}+{Kq}=F(t)}$$

• $K$: stiffness symmetrical matrix
• $M$: Mass symmetrical matrix
• $G$: Gyroscopic skew symmetrical matrix
• $\Omega_c$: angular speed
• $F$: excitation
• $q$: generalized coordinates

While trying to solve the eigenvalue problem with Matlab, I noticed there´s a whole bunch of ways to represent the equations in state space form and I don´t understand what the main difference is, or which one I should use, since when calculating I get different results in the eigenvectors.

But for the free vibrations I have the equation without damping:

$$\mathbf{M{\ddot q}+{Kq}=0}$$ 1. and to solve the standard eigenvalue problem of:

$$\mathbf {\omega^2M\phi_i-{K}\phi_i=0}$$ I use the Matlab code:

[v,D]=eig(K,M)
root=sqrt(D)
natural_frq=root/(2*pi)


which gives me all real natural frequencies and eigenvector per frequency.

2. I also tried calculating them in the form of: $$\mathbf{{\dot x}=Ax+Bu}$$

$$\begin{bmatrix} \mathbf{\dot q} \\ \mathbf{\ddot q} \end{bmatrix} = \begin{bmatrix} \mathbf{0} & \mathbf{I} \\ \mathbf{-M^{-1}{K}} & \mathbf{-M^{-1}{G}} \end{bmatrix} \begin{bmatrix} \mathbf{ q} \\ \mathbf{\dot q} \end{bmatrix}+ \begin{bmatrix} \mathbf{0} \\ \mathbf{-M^{-1}} \end{bmatrix}u$$

and by using:

[v,D]=eig(A)


I get similar results but in complex numbers (there´s no need for the square root), so with some extra code I choose only the imaginary parts and select the positive values of the eigenvalues vector, which instead of 18 elements is now 36. The frequencies I get with this method against the first ones are paractically the same, but the eigenvectors are not, they are also complex and the results don´t look correct.

3. Then there's another state space form:

$$\begin{bmatrix} \mathbf{I} & \mathbf{0} \\ \mathbf{0} & \mathbf{M} \end{bmatrix} \begin{bmatrix} \mathbf{\dot q} \\ \mathbf{\ddot q} \end{bmatrix} + \begin{bmatrix} \mathbf{0} & \mathbf{-I} \\ \mathbf{K} & \mathbf{G} \end{bmatrix}\begin{bmatrix} \mathbf{ q} \\ \mathbf{\dot q} \end{bmatrix} = 0$$

If I use the code: [v,D]=eig(-B,A) I get similar results as the last method, I'm not sure which form of state space this form is, I would guess it's $$A{\dot x} = -Bu$$.

The main question is that I don't know which is the best solution to my problem and why I get complex frequencies in the second method, if my angular speed $$\Omega_c$$ is zero, all damping terms should also go to zero and thus not affect the system results , right?

EDIT: why do I get very different values when calculating the eigenvectors through the standard eigenproblem form $$\mathbf {\omega^2M\phi_i-{K}\phi_i=0}$$ and with the state space?

Regarding the different state space formulations. First, M might not be invertible. In many cases it will be invertible, but sometimes it will not be. In your specific case, this would correspond to a gear with zero mass. Obviously that never happens in the real world as all gears would have some finite mass. But there are situations, where the modeling choices and assumptions lead to a DOF with zero mass. In that case, you cannot use the first form of the state space equations. You must use the second one.

Now, if M is invertible, then analytically, both formulations will give the same answer. But numerically, they will invariably give slightly different answers due to floating point roundoff error. Depending on the solver that you use, one form or the other might end up being more accurate. Speed differences are also possible. This depends on the exact structure of your matrices and the details of your solver. If you want more details, pick up a book on numerical methods that describes the inner workings of the various eigenvalue solvers.

Edit based on your latest edit, I'm not that my original answer really answers your question. These two methods really should give you similar answers. The fact that they are not implies that you are probably doing something wrong, but I don't know what. Suggestion: don't start with an 18x18 system. Start with a 2x2 system, something small enough that you can work out the exact answer by hand on one sheet of paper, or even better look up the answer in a textbook. That should give you insight into where your matlab code is going wrong. If not, a 2x2 system is small enough that you can post the entire thing here (all the inputs, all the outputs, all the code), and we can probably figure it out. As is, without seeing your input matrices, or the whole code, and only portions of some of the eigenvectors, I'm not sure if I can be of more help.

• in my case the M matrix is always invertible, so I should be getting "slightly" same results between solving by the standard solution eig(K,M) and the state space form eig(A) right? I don't understand why I get such different results. Jul 25 '17 at 13:51
• Unless the M matrix has a simple form (for example it is diagonal) finding its inverse and then multiplying by K and G is just a waste of computer time in method (2), and probably less well conditioned numerically than the eigenproblem in method (3). If M is large (say 10,000 degrees of freedom or more) and sparse, you never want to invert it numericallly to get a large, non-sparse, and poorly conditioned matrix! Finding a few thousand eigenpairs of a million-degree-of-freedom system is perfectly possible, but only if you do it the right way, not an inefficient wrong way! Jul 25 '17 at 16:19

For any sinusoidal motion, the $\mathbf{q}$ really represents the motion $\Re(\mathbf{q}\,e^{i \omega t})$ in the time domain. $\mathbf{q}$ is complex because in general, the different elements of the eigenvector have different phase angles.

This will be the case for your problem since your first equation seems to contain a gyroscopic term - if the rotation is about the $x$ axis, the Eigenvector components if the $y$ and $z$ directions will be 90 degrees out of phase to represent the "whirling" motion of the structure.

Also, $\mathbf{\dot q} = i \omega\, \mathbf{q}$, so the vector $\begin{bmatrix} \mathbf{q} & \mathbf{\dot q} \end{bmatrix}^T$ can not have all its terms real, except for the special case where $\omega = 0$.

Remember that the eigenvectors are to some extent arbitrary even for a real eigenproblem. Even if you insist they are mass normalized, you can still multiply them by $-1$. For a complex eigenproblem, you can multiply them by any complex value on the unit circle, of the form $\cos \theta + i \sin \theta$ for an arbitrary $\theta$.

If you want to compare the solutions for two complex eigenproblems, you have to fix these arbitrary differences - for example, by "rotating" the vectors around the unit circle so that one element of the eigenvector is real and positive, and "scaling" the vectors so that the real positive value is $+1.0$.

These arbitrary factors don't matter if you just want to use the eigensolution to compute something which is real for physical reasons. When you take a linear combination of the eigenvectors, the scaling factors multiplying each vector are also complex, and they "clean up" the arbitrariness in the eigensolution, in exactly the same way that it doesn't matter if you calculate $\mathbf{+x}$ or $\mathbf{-x}$ for a real eigenvector.

• that's great, thanks. I don't understand what you mean about multiplying the vector by -1 or by any complex number, what effect would this have on the eigenvectors? and if matlab's command eig gives me the results normalized, is there still need to normalize them? Jul 25 '17 at 14:01
• If you calculate a real eigenvector, there is no difference between, [+1.0, -0.5] and [-1.0, +0.5]. For a complex Eigen vector, there is no difference between [1,-0.5i] or [i, 0.5] or [-1, 0.5i], or [-i, -0.5], or even something like [0.6+0.4i, -0.2-0.3i]. Draw all these examples on an Argand diagram, and you will see they are "the same pattern" but rotated through different angles. Jul 25 '17 at 16:16
• There are several mathematical ways to "normalize" eigenvectors, and if your mathematical eigenvector contains both $\mathbf{q}$ and $\mathbf{\dot q}$ the normalization doesn't necessarily mean anything physically. You might be better off just normalizing based on the $\mathbf{q}$ part of the vector and doing the same transformation to the $\mathbf{\dot q}$ part, but Matlab probably won't do that automatically. Jul 25 '17 at 17:26

You basically want to solve for the homogeneous solutions of the system

$$\textbf{M}\, \ddot{\textbf{q}} + \textbf{D}\, \dot{\textbf{q}} + \textbf{K}\, \textbf{q} = \textbf{0}.$$

It can be shown that each solution will be of the form

$$\textbf{q}_i(t) = \textbf{u}_i\,e^{\mu_i\,t},$$

where $\textbf{u}_i$ is a vector of the same dimension as $\textbf{q}$ (usually normalized to a length of one) and $\lambda_i$ is a scalar. In general any homogeneous solution can be constructed by taking a linear combination these solution (I will ignore for now that you can also have Jordan blocks of bigger then one by one, which will increase the possible solutions).

Substituting the solution into the homogeneous equation yields

$$\textbf{M}\, \textbf{u}_i\,\mu_i^2\,\,e^{\mu_i\,t} + \textbf{D}\, \textbf{u}_i\,\mu_i\,\,e^{\mu_i\,t} + \textbf{K}\, \textbf{u}_i\,e^{\mu_i\,t} = \left(\textbf{M}\, \mu_i^2 + \textbf{D}\, \mu_i + \textbf{K}\right) \textbf{u}_i\,e^{\mu_i\,t} = \textbf{0}.$$

Since $e^{\mu_i\,t}$ is always non-zero (for finite time) it can be ignored since the remaining terms on the left hand side of the equation still has equal to zero in order to satisfy the right hand side. Doing so already starts to make the equation look more like a generalized eigenvalue problem

$$\textbf{A}\, \textbf{v} = \lambda\, \textbf{B}\, \textbf{v},$$

however $\mu_i$ appears both linearly and quadratically. In the special case when $\textbf{B} = 0$ then the generalized eigenvalue problem can be formulated using $\textbf{A} = \textbf{K}$, $\textbf{B} = \textbf{M}$, $\textbf{v} = \textbf{u}_i$ and $\lambda = -\mu_i^2$. Another special case is when $\textbf{M}$, $\textbf{D}$ and $\textbf{K}$ are all simultaneously diagonalizable by pre and post multiplying them by the same matrices. Namely in that case the differential equation can be written as a system of decoupled second order differential equations, like a single mass-spring-damper system. One way such transformation could be constructed when $\textbf{M}$, $\textbf{D}$ and $\textbf{K}$ are symmetric is

$$\mathbf{U} = \begin{bmatrix}\mathbf{u}_1 & \mathbf{u}_2 & \cdots & \mathbf{u}_n\end{bmatrix},$$

where $\mathbf{u}_i$ are the eigenvectors of the eigenvalue problem with $\textbf{B} = 0$. The matrices with which you pre and post multiply would then be $\mathbf{U}^\top$ and $\mathbf{U}$ respectively. Often if the system is not simultaneously diagonalizable and the damping is small, then assuming it is can be a good approximation (so setting the non-diagonal elements of $\mathbf{U}^\top \mathbf{B}\, \mathbf{U}$ to zero). A specific form for which this is the case is called Rayleigh damping, when $\mathbf{B} = \alpha\,\mathbf{M} + \beta\,\mathbf{K}$.

The advantage of all the previous mentioned methods is that the generalized eigenvalue problem you have to solve are of the same dimension as $\mathbf{q}$. However if none of cases mentioned previously apply to your system you would have to solve the generalized eigenvalue problem of twice the dimension of $\mathbf{q}$. There are multiple ways of formulating this. A few of them are

$$\begin{bmatrix} \mathbf{I} & \mathbf{0} \\ \mathbf{0} & \mathbf{M} \end{bmatrix} \begin{bmatrix} \dot{\mathbf{q}} \\ \ddot{\mathbf{q}} \end{bmatrix} = \begin{bmatrix} \mathbf{0} & \mathbf{I} \\ -\mathbf{K} & -\mathbf{D} \end{bmatrix} \begin{bmatrix} \mathbf{q} \\ \dot{\mathbf{q}} \end{bmatrix},$$

$$\begin{bmatrix} \mathbf{I} & \mathbf{0} \\ \mathbf{0} & \mathbf{I} \end{bmatrix} \begin{bmatrix} \dot{\mathbf{q}} \\ \ddot{\mathbf{q}} \end{bmatrix} = \begin{bmatrix} \mathbf{0} & \mathbf{I} \\ -\mathbf{M}^{-1}\mathbf{K} & -\mathbf{M}^{-1}\mathbf{D} \end{bmatrix} \begin{bmatrix} \mathbf{q} \\ \dot{\mathbf{q}} \end{bmatrix},$$

$$\begin{bmatrix} \mathbf{D} & \mathbf{M} \\ \mathbf{M} & \mathbf{0} \end{bmatrix} \begin{bmatrix} \dot{\mathbf{q}} \\ \ddot{\mathbf{q}} \end{bmatrix} = \begin{bmatrix} -\mathbf{K} & \mathbf{0} \\ \mathbf{0} & \mathbf{M} \end{bmatrix} \begin{bmatrix} \mathbf{q} \\ \dot{\mathbf{q}} \end{bmatrix}.$$

For this extended state space it can easily be seen that the solutions will be of the form

$$\begin{bmatrix} \mathbf{q}_i(t) \\ \dot{\mathbf{q}}_i(t) \end{bmatrix} = \begin{bmatrix} \mathbf{u}_i \\ \mu_i\,\mathbf{u}_i \end{bmatrix} e^{\mu_i\,t}.$$

Substituting this into the three equations above and factoring out $e^{\mu_i\,t}$ yields

$$\begin{bmatrix} \mathbf{I} & \mathbf{0} \\ \mathbf{0} & \mathbf{M} \end{bmatrix} \begin{bmatrix} \mathbf{u}_i \\ \mu_i\,\mathbf{u}_i \end{bmatrix} \mu_i = \begin{bmatrix} \mathbf{0} & \mathbf{I} \\ -\mathbf{K} & -\mathbf{D} \end{bmatrix} \begin{bmatrix} \mathbf{u}_i \\ \mu_i\,\mathbf{u}_i \end{bmatrix},$$

$$\begin{bmatrix} \mathbf{I} & \mathbf{0} \\ \mathbf{0} & \mathbf{I} \end{bmatrix} \begin{bmatrix} \mathbf{u}_i \\ \mu_i\,\mathbf{u}_i \end{bmatrix} \mu_i = \begin{bmatrix} \mathbf{0} & \mathbf{I} \\ -\mathbf{M}^{-1}\mathbf{K} & -\mathbf{M}^{-1}\mathbf{D} \end{bmatrix} \begin{bmatrix} \mathbf{u}_i \\ \mu_i\,\mathbf{u}_i \end{bmatrix},$$

$$\begin{bmatrix} \mathbf{D} & \mathbf{M} \\ \mathbf{M} & \mathbf{0} \end{bmatrix} \begin{bmatrix} \mathbf{u}_i \\ \mu_i\,\mathbf{u}_i \end{bmatrix} \mu_i = \begin{bmatrix} -\mathbf{K} & \mathbf{0} \\ \mathbf{0} & \mathbf{M} \end{bmatrix} \begin{bmatrix} \mathbf{u}_i \\ \mu_i\,\mathbf{u}_i \end{bmatrix}.$$

These are all just generalized eigenvalue problems. I am not an expert on numerically solving these problems, however I have few things that can be noted. Namely the second problem formulation requires a matrix inverse, which might not exist, but does reduce it to a normal eigenvalue problem (not sure if this requires any less computation time in general). The third problem formulation will yield symmetric matrices if $\textbf{M}$, $\textbf{D}$ and $\textbf{K}$ are symmetric as well, which might have some attractive properties when doing numerical calculations. So depending on your system and solver a different formulation might give more accurate results.

• thanks for your detailed explanation. Is there a specific reason you are using the $\mu_i$ instead of $\omega_i$ ? Also, I got the results I wanted for now through the generalised eigenvalue problem, so would there be a reason to try and solve it through the state space form? from what I understand on my limited knowledge of the topic, this would only be useful for control purposes? thanks Jul 31 '17 at 16:17
• @quo I just wanted to used a different symbol then $\lambda$, so there is no meaning behind the use of $\mu_i$. Also any generalized eigenvalue problem is probably derived from a state space formulation of the homogeneous differential equation. That is also what I tried to give more insights into in my answer. Jul 31 '17 at 16:36