How to determine the dynamic response of a spring-mass system?

Question

The Problem

I have the following system for which I would like to determine its dynamic response, taken from these lecture notes pg.20:

where the forcing functions are $u_1 = 1-U(t - t_0)$, $u_2=0$, and $u_3=0$, $U(\cdot)$ is the unit step function and $t_0 = 2\pi/\omega_1$ where $\omega_1$ is the first undamped natural frequency.

---

My Attempt

I know the time domain representation of the system to be of the typical format as shown in Eq.(1), and therefore the Fourier transform is as shown in Eq.(2):

\begin{gather} \mathbf{M}\ddot{\mathbf{x}} + \mathbf{C}\dot{\mathbf{x}} + \mathbf{Kx} = \mathbf{u}(t) \tag{1} \\ % \Rightarrow [-\omega^2\mathbf{M} + i\omega\mathbf{C} + \mathbf{K}]\mathbf{X}(i\omega) = \mathbf{U}(i\omega) \tag{2} \end{gather}

where since we know that $\mathcal{L}\{1 - U(t - t_0)\} = 1 - \frac{\exp{(-st_0})}{s}$, then the Fourier transform of the right hand side of Eq.(2) is:

$$\mathbf{U}(i\omega) = \left[\left(1 - \frac{e^{-i\omega t_0}}{i\omega}\right),\, 0,\, 0,\right]^T$$

As a result, the dynamic response would be:

$$\mathbf{X}(i\omega) = [-\omega^2\mathbf{M} + i\omega\mathbf{C} + \mathbf{K}]^{-1}\mathbf{U}(i\omega) \tag{3}$$

Since the response by which I am comparing my answer is given in generalised modal coordinates $\mathbf{q}$, then the modal transformation of Eq.(3) follows where $\mathbf{V}$ is the eigenvector matrix thus we know that the modal mass, stiffness, damping matrices should be $\bar{\mathbf{M}} = \mathbf{I}$ (i.e. the identity matrix), $\bar{\mathbf{K}} = \mathbf{\Omega}^2$ (i.e. a diagonal marix containing the square of the natural frequencies), and $\bar{\mathbf{C}} = \mathbf{V}^T\mathbf{C}\mathbf{V}$ respectively. So Eq.(3) becomes:

$$\mathbf{q}(i\omega) = \left[-\omega^2\mathbf{I} + i\omega\bar{\mathbf{C}} + \mathbf{\Omega}^2 \right]^{-1}\left[\mathbf{V}^T\mathbf{U}(i\omega)\right] \tag{4}$$

Thefore, the dynamic response would be the magnitude of the imaginary numbers obtained from Eq.(4) above, for all positive $\omega$.

However, when I follow the above steps and solve the problem in MATLAB, the answer is far from the one provided in the excercise. Admitedly, the code is very simple as well (see below) so I struggle to see where I could have gone wrong:

clear ; clc

% Input
n = 3 ;         % Number of DOFs
m = 1 ;         % Mass [kg]
k = 1 ;         % Stiffnes [N/m]
c = 0.2 ;       % Viscous damping [Ns/m]
nv = [1 0 0]' ; % Load application vector (load applied to 1st DOF)

% Matrices
M = diag(m*ones(n, 1)) ;
K = diag(2*k*ones(n,1)) - diag(k*ones(n-1,1),-1) - diag(k*ones(n-1,1),1) ;
C = diag(c*ones(n,1)) ;

% Undamped Natural Frequencies & Modes
[V, Omsq] = eig(K, M) ;
Om = diag(sqrt(abs(Omsq))) ; % Natural Frequencies

% Frequency Aanalysis
t0 = 2*pi/Om(1) ;           % Time where unit step is applied
omega = linspace(0, 3, 100) ; % Frequency vector

for i = 1:numel(omega)
    
    s = 1i*omega(i) ;
    
    Us = V'*((1 - exp(-s*t0)/s)*nv) ;
    FRF(:,i) = (s^2*eye(n) + s*(V'*C*V) + Omsq)\Us ;
end

% Plotting
figure ; plot(omega, abs(FRF)) ; grid on ; set(gca, 'YScale', 'log') ; ylim([0 10]) ; 

Answer 1

As TimWescott indicated, this problem is more suitable for control theory for a fuller understanding than for digital signal processing, yet there's no control.se as far as I know. Hence you can give it a try here too since such damped-spring systems is not uncommon subject for audio / acoustic engineers etc.

First of all you have a programming error: Your for-loop variable is set to i , a potential error because there's already the complex unit $j$ in your body code in the line $$s = j \omega \tag{1}$$which you have implemented as :

s = 1i*omega(i)

which is wrong: you are using the same letter i , both for indexing integer into the angular frequency vector omega, and also for the imaginary unit $j$. I've replaced the for-loop variable with k instead.

Second, you have a signal processing mistake. Your nonzero input function is $$ u_1(t) = 1 - u(t-t_0) \tag{2}$$whose (bilateral - double sided) Laplace transform is

$$ - \frac{ e^{ -s t_0} }{s} \tag{3} $$

You can see this by noting that the input is actually: $$u_1(t) = u(t_0 - t) \tag{4}$$

And the third issue is about the diagonalized computation. I'have replaced the diagonalised frequency response calculation, with the undiagonalised one with the matrcices M,K, C. The result is improved,as shown on the plot below, but still not exactly what you have been given as the answer plot. I guess you may correct it further.

I put the corrected code and the corresponding MATLAB output below :

clc; clear variables ; close all

% Input
n = 3 ;         % Number of DOFs
m = 1 ;         % Mass [kg]
k = 1 ;         % Stiffnes [N/m]
c = 0.2 ;       % Viscous damping [Ns/m]
nv = [1 0 0]' ; % Load application vector (load applied to 1st DOF)

% Matrices
M = diag(m*ones(n, 1)) ;
K = diag(2*k*ones(n,1)) - diag(k*ones(n-1,1),-1) - diag(k*ones(n-1,1),1) ;
C = diag(c*ones(n,1)) ;

% Undamped Natural Frequencies & Modes
[V, Omsq] = eig(K,M); 
Om = diag(sqrt(abs(Omsq))) ; % Natural Frequencies

% Frequency Aanalysis
t0 = 2*pi/Om(1) ;           % Time where unit step is applied
omega = linspace(0, 2.5, 100) ; % Frequency vector

for k = 1:numel(omega)
    
    s = 1i*omega(k) ;
           
    Us = -(exp(-s*t0)/s)*nv; 
    FRF(:,k) = (s^2*M + s*C + K)\Us ;

    %Us = V'*((1 - exp(-s*t0)/s)*nv) ;
    %FRF(:,k) = (s^2*eye(n) + s*(V'*C*V) + Omsq)\Us ;

end

% Plotting
figure ; plot(omega, abs(FRF)) ; grid on ; set(gca, 'YScale', 'log') ; ylim([0.01 10]) ;