Undamped forced vibration for a multi dof system

I'm trying to determine the amplitude vs the excitation frequency of a multi dof dynamic system. I'm very confused from reading different methods, although I understand there's a very general one, which is not giving me the results that I think should be when plotting the response. Considering a system with no damping $\mathbf{M{\ddot x}+{Kx}=F(t)}$, and by assuming the harmonic excitation of the force is a cosine the system becomes $\mathbf{-MX{cos(\omega t)}+{KX}cos(\omega t)=fcos(\omega t)}$. From here the cosines get cancelled out so only $\mathbf{[K-M{\omega^2}]X=f}$ remains. By solving for the $X$ vector in the end I get $X = \frac{F/M}{|\omega_n^2-\omega^2|}$.

So with Matlab I wrote the following code, with parameters I found on an example from brown.edu, in which the explanation is not complete, my system is actually a 24 DoF system, so I need the code a bit automated:

    k1 = 2;
k2 = 1;
k3 = 1;
m1 = 1;
m2 = 1;

M=[m1 0; ...
0 m2];
K=[k1+k2 -k2; ...
-k2 k3+k3];
F=[1;...
1];

[V,D] = eig(K,M);
D     = diag(sqrtm(D));
frequencies   = (D);
[rows,~] = size(V);

f=linspace(0,0.7,100);
omega=2*pi*f;
Mag = M\F;
Xt = zeros(rows, length(omega));
for i = 1:rows
for j = 1:length(omega)
Xt(i,j) =  abs((frequencies(i)^2 - omega(j)^2)\Mag(i)) ;
end
end
figure
for m = 1:length(frequencies)
plot(omega, Xt(m,:));
set(gca,'YScale','log')
hold on
axis tight
end


Where the main part is the Xt(i,j) = (abs (frequencies(i)^2 - omega(j)^2)\Mag(i) ); , which follows the determined formula for the displacements of each degree of freedom.

When plotting what I get is:

which shows the masses vibrations being at resonance, but from what I understand there should be more peaks and some antiresonance.

What am I missing in the general math?

EDIT_1: I changed the bit of the code to calculate the absolute values of the equation like: $$X(\omega) = \left|\sum_n \frac{F/M}{\omega_n^2 - \omega^2}\right|$$ but I get the same when plotting. I found another solution but I don't understand how its determined or why choose it. From ${M{\ddot x_i}+{Kx_i}=f(t)sin\omega t}$ the steady state solution is $x_{iss} = x_{is} sin \omega t + x_{ic} cos \omega t$ , where $i = 1,2$ , then $$\begin{bmatrix} \mathbf{-(\omega)^2M+{K}} & \mathbf{0} \\ \mathbf{0} & \mathbf{-(\omega)^2M+{K}} \end{bmatrix} \begin{bmatrix} \mathbf{x_{ic}} \\ \mathbf{x_{is}} \end{bmatrix} = \begin{bmatrix} \mathbf{0} \\ \mathbf{f} \end{bmatrix}$$ which is then solved to get ${X_{i}} = (x_{ic} + x_{is} )^{1/2}$ First, I don´t understand very well if just $x_{iss}$ is substituted in ${M{\ddot x_i}+{Kx_i}=f(t)sin\omega t}$ for $x_i$. And is this just assuming a sum of $cos$ and $sin$ instead of my first approach of only $cos$? Then, do I need to consider adding the transient response also, which I´m actually not sure how it would be.

EDIT_2: I think I figured out what I was doing wrong. I changed the equation from $$X(\omega) = \left|\sum_n \frac{F/M}{\omega_n^2 - \omega^2}\right|$$ to $$X = [K-M\omega^2]^{-1}F$$ which mathematically I though they were the same??? The code part that changed is as follows:

for j = 1:length(omega)
Xt(:,j)  =  abs((K-M*omega(j)^2)\F);
end
figure
for m = 1:length(frequencies)
plot(omega,Xt(m,:))
hold on
axis tight
end


And the result plot is:

For my extended MDoF system I got:

In which I excited all DoF, even when I only excite a single one, it is still very hard to notice what's going on, is there a better way of analysing the response?

• I don't know what your last solution is trying to do. It's obvious that $x_{ic} = 0$, at least when there is no damping...??? Aug 11 '17 at 2:46
• "mathematically I thought they were the same???" One is in physical coordinates, the other is supposed to be in modal coordinates - but you are probably still doing something wrong in modal coordinates. For a 24 DOF model, it doesn't make much difference which you choose. For a 1,000,000 DOF model, reducing it to say 1,000 modes and ignoring the rest makes a very big practical difference (i.e. it saves hours of computing time). Aug 11 '17 at 6:53
• "In which I excited all DoF, even when I only excite a single one" If you apply a force at only one position, the all the other DOFs should be responding! Note, if you convert to modal coordinates correctly, your single DOF physical force will be applied to all the modes in general (unless you applied the force at a nodal position for some modes). Aug 11 '17 at 6:59
• I think you really need to review multi degree of freedom analysis in modal coordinates - you just seem to be blundering from one misunderstanding to another. See any good dynamics textbook. My answer is a summary of most of the mathematics in theory, but there's a limit to how much of a whole textbook chapter that I'm going to write on SE for free (and you have already reached that limit!) Aug 11 '17 at 7:02
• @alephzero I followed your suggestion and reviewed the mdof chapters, I understand better now and see where you were going, although from the various textbooks I read , none of them show any practical examples, which is why I have trouble understanding all these concepts, I'll just keep trying my best to understand, thank you for your great help, I really do appreciate it Aug 13 '17 at 22:54

The code in your OP seems to be for a 2 DOF system not a 24 DOF system. You could check that by printing the "frequencies" variable from your eigensolution.

The reason you are not getting anti-resonances is because you are taking the absolute values in the wrong place. You should be calculating $$X(\omega) = \left|\sum_n \frac{F/M}{\omega_n^2 - \omega^2}\right|$$ not $$X(\omega) = \sum_n \frac{F/M}{\left|\omega_n^2 - \omega^2\right|}$$

The contribution to $X$ from each mode changes its phase by 180 degrees, and therefore changes sign, when the frequency passes through the corresponding $\omega_n$. The antiresonances occur when the contributions from different modes cancel out.

Note, you will probably see the antiresonances better on a plot if you use a log scale for $X$.

EDIT following comments by the OP:

You also need to decompose $F$ into its modal components and apply them to the corresponding terms in the sum. If $$M\ddot x + K x = F.$$ Pre-multiply by the eigenvectors $\Phi$: $$\Phi^T M\ddot x + \Phi^TK x = \Phi^T F$$ But you know $\ddot x = -\omega^2 x$, so $$(-\omega^2 \Phi^T M + \Phi^TK) x = \Phi^T F$$ Also $\Phi^T M \Phi = I$ for mass-normalized eigenvectors, and $\Phi^T K \Phi = \text{diag}\,(\omega_n^2)$, so inserting $\Phi \Phi^{-1}$ into the left hand side gives $$(-\omega^2 + \text{diag}\,(\omega_n^2)) \Phi^{-1}X = \Phi^T F.$$ Decomposing $X$ into its modal components gives $$X = \sum_n \Phi_n\xi_n,$$ so you end up with the decoupled equations $$(-\omega^2 + \omega_n^2) \xi_n = \Phi_n^T F$$.

In other words, if $\mathbf{f}$ is the vector applied forces, the $F$ in each term of your sum should be $\Phi_n^T \mathbf{f}$ where $\Phi_n$ is the $n$th column of the eigenvector matrix.

Finally, you need to recombine the $\xi_n$ to get $X$, and then take its absolute value.

• Thanks! yes, I actually meant to say that this example is for 2 DoF but my system is 24, which I´ll just replace the matrices after getting this one right. I changed my code to : $$X(\omega) = \left|\sum_n \frac{F/M}{\omega_n^2 - \omega^2}\right|$$ but I still don't get the right result, I´ll edit my main question. Aug 11 '17 at 0:11
• See the edit to my answer. Aug 11 '17 at 1:24