Getting natural frequencies, damping ratios and modes of vibration from the state-space format of equations

My question is fairly simple. I have a highly complex nonlinear model dynamic model, and I want to linearize it around a working point so I get the matrices A,B,C and D for the state-space format of ODEs.

From that (linearized system), I would like to extract the natural frequencies, the damping ratios, and the modes of vibration for each degree of freedom. Is it the eigenvalues and eigenvectors for the ss(A,B,C,D) that give me information about it? If I do:

sys = ss(A,B,C,D);
[s,v]=eig(sys);

s would be my eigenvalues and v my eigenvectors. From this matrices s and v, I get the natural frequencies and the modes of vibration, respectively?

Another question is, my model has 7DoF, so I have 14 states to represent its dynamics. I though I would have only 7 eigenvalues of the system, but if I procceed in this way, I'll get an eigenvalue for all the displacements and the velocities (so 14 eigenvalues, thus 14 natural frequencies) Does this make physical sense? Is this correct?

[Pedro Calorio]

Also, what would be the different between the following:

%I have a given M, C and K matrix for n DoF
[s,v] =eig(inv(M)*K)
wn = sqrt(v) %square root of the eigenvalue
modes = v

And

sys = ss(A,B,C,D) %state space format of my dynamical system
[s,v]=eig(sys)

In the first method I get n natural frequencies, while in the last one I'll obtain 2*n natural frequencies (all second order ODEs...)

What is right what is wrong? Does existis a different natural frequency and damping ratio for displacement and velocity?

[David Goodmanson]

Hi Pedro, the short answer is, there are two possible signs for the square root of the eigenvalue and both of them count, so things work out all right.

[Pedro Calorio]

David, could you explain with a little bit more details?

The order I get my eigenvalues from eig is the order of the states vector?

For example:

X = [x1 x2 x3 x4]' state vector

I get from eig(sys) a column vector of:

eigvalues = [a b c d e f g h]'

a,b is refering to x1 pair

c,d is refering to x2 pair

and so on?

Answer 1

[David Goodmanson]

Hi Pedro,

here is an example, two masses and two springs, with dash pots in parallel with the springs so there is a force equal to -c*v = -c*x' as well as -k*x from the spring. The springs have unstretched length zero, and the masses are allowed to pass through each other and through the attachment point on the left.

 |   k1,c1      k2,c2
 |----------m1----------m2
 |
x=0

There are two displacements and two velocities, and the state space has four dimensions. The equations are

m1*x1'' = -k1*x1 -c1*x1' + k2(x2-x1) + c2*(x2'-x1')
m2*x1'' =                  k2(x1-x2) + c2*(x1'-x2')

For convenience the state vector is in the order [x1; x2; x1'; x2']. Four dimensions mean there are four eigenvalues alpha. Each solution is of the form exp(alpha*t) * eigenvector. As you say the first eigenvalue goes with the first column of v (first eigenvector) and so forth.

alpha =
  -0.2094 + 1.6475i
  -0.2094 - 1.6475i
  -0.0239 + 0.4910i
  -0.0239 - 0.4910i

The displacements of the four independent solutions are shown in the plots (no velocities are plotted). You can take linear combinations of these four to satisfy four boundary conditions, usually positions and velocities at t=0.

The first two solutions are complex conjugates of each other. You can take the sum and difference of these to get two independent real solutions, or you can take the real and imaginary parts of the first solution as is done below.

Same idea for the third and fourth solutions.

In he first two solutions m1 and m2 move opposite each other, and in the third and fourth solutions the two masses move in the same direction. The k2 spring is more compressed in the first two solutions, leading to a much higher natural frequency than in the other case.

k1 = 1;
k2 = 2;
m1 = 3;
m2 = 1;
c1 = .2;
c2 = .3;

K = [       0     0        1   0  
            0     0        0   1
      -(k1+k2)   k2  -(c1+c2) c2
           k2   -k2       c2 -c2];
M = [1 0 0 0
     0 1 0 0 
     0 0 m1 0 
     0 0 0 m2];

[v alpha] = eig(K,M,'vector');
% same as [v alpha] = eig(inv(M)*K,'vector')
alpha

t = 0:.01:40;
xs = v(:,1)*exp(alpha(1)*t);  % higher frequency, column 1 of v
figure(1)
plot(t,real(xs(1:2,:))); grid on
figure(2)
plot(t,imag(xs(1:2,:))); grid on

xs = v(:,3)*exp(alpha(1)*t);  % lower frequency, column 3 of v
figure(3)
plot(t,real(xs(1:2,:))); grid on
figure(4)
plot(t,imag(xs(1:2,:))); grid on