In my course of dynamics of structures, I am struggling with some basic principles. As a practice, I am plotting (or attempting to plot) the curves that are shown in the course book, with the help of Matlab. For example, I want to plot the phase angle theta as a function of beta, the dimensionless frequency, for a damped SDOF. The code I have is the following:

%Stiffness of the spring in [N/m]
k = 1;
%Mass in [kg]
m = 1;
%Damping constant in [Ns/m]
c = 0.1;
%Initial displacement in [m]
u_0 = 0;
%Initial velocity in [m/s]
v_0 = 0;
%Resonance frequency in [rad/s]
omega_res = sqrt(k/m);
%Resonance frequency in [Hz]
f_res = (omega_res)/(2*pi);
%Eigenperiod in [s]
T_res = 1/f_res;
%Critical damping in [Ns/m]
c_c = 2*sqrt(k*m);
%Damping ratio in [/]
xi = c/c_c;
%Damped eigenfrequency in [rad/s]
omega_d = omega_res * sqrt(1-xi^2);
%Damped eigenfrequency for strongly damped system in [rad/s]
omega_hat = omega_res * sqrt(xi^2-1);
t = 0:0.1:100;
beta_v = 0:0.001:4;
%D as a function of beta
D = 1./sqrt((1-beta_v.^2).^2+(2*beta_v.*xi).^2);

axis([0 4 0 10])
xlabel('Dimensionless frequency [-]')
ylabel('Dynamic amplification factor [-]')

%theta as a function of beta
theta = atan((2*beta_v*xi)./(1-beta_v.^2));

axis([0 4 0 pi])
grid on
xlabel('Dimensionless frequency [-]')
ylabel('Phase angle [rad]')

But plotting gives me the curve of an undamped SDOF. The curve I get:

enter image description here

The plot I am supposed to get:

enter image description here

Does anyone have a clue what is wrong in my approach?

A second struggle occurs when you want to plot the the response for a particular beta. The curve I achieved was not what I expected. My code is the following:

beta_a = 0.2;
D_a = 1./sqrt((1-beta_a.^2).^2+(2*beta_a.*xi).^2);
theta_a = atan((2*beta_a*xi)./(1-beta_a.^2));
p_hat = 1;
omega = beta_a.*omega_res;
u_p = p_hat/k.*D_a.*exp(i.*(omega.*t-theta_a));
B = u_0 - p_hat./k.*D_a.*exp(-i.*theta_a);
A = (v_0 + xi*omega_res*B - p_hat/k*D_a.*i*omega*exp(-i.*theta_a))/omega_d;
u_h = exp(-xi.*omega_res.*t).*(A.*sin(omega_d.*t) + B.*cos(omega_d.*t));
u = u_p + u_h;
grid on
xlabel('Time [s]')
ylabel('Displacement [m]')

Can anyone help me with this line of thought? What's wrong?

  • 1
    $\begingroup$ It would help if you include pictures of what you expect the responses to be and what you're actually getting from your code. $\endgroup$
    – Paul
    Feb 21, 2016 at 14:38
  • $\begingroup$ I am not allowed to post more than two pictures, so I included the two of the phase angle. $\endgroup$
    – Beunhaas
    Feb 21, 2016 at 21:56
  • 1
    $\begingroup$ Your code for plotting the phase angle seems to be wrong. Could you add that to your question? Otherwise the rest of the code is OK. See my answer below. $\endgroup$ Feb 22, 2016 at 22:24
  • $\begingroup$ I added the missing code for plotting the phase angle. $\endgroup$
    – Beunhaas
    Feb 23, 2016 at 20:19

1 Answer 1


With your code, and with $c = 0.3$, I get the following results: enter image description here enter image description here enter image description here

The amplitude and phase look OK. But the displacement does not show any damping. To see why read on below.

I'm not sure about your notation. So the equations below may differ from those in your textbook.


The ODE you are trying to solve is $$ \ddot{u} + 2\xi\omega_{\text{res}}\dot{u} + \omega_{\text{res}}^2 u = \frac{\hat{p}}{m}\text{Re}\left[\exp(i\omega t)\right] $$ where $$ \xi = \frac{c}{c_c} ,~~ c_c = 2 m \omega_{\text{res}}, ~~ \omega_{\text{res}} = \sqrt{\frac{k}{m}} $$ and $c$ is the damping factor, $k$ is the stiffness, $m$ is the mass, $\omega$ is the forcing frequency.

Homogeneous solution

The homogeneous solution is $$ u_h = A\,\exp(-\xi\omega_{\text{res}} t)\sin(\omega_d t + \phi_d) $$ where $\phi_d$ is a phase and $$ \omega_d = \omega_{\text{res}}\sqrt{1 - \xi^2} $$

Notice that the exponential decay term occurs only in the homogeneous solution.

Particular solution

The particular solution is $$ u_p = \frac{\hat{p}}{k} \frac{\text{Re}\left[\exp[i(\omega t - \phi)]\right]}{\sqrt{\left[1 - \beta_a)^2\right]^2 + \left[2\xi\beta_a\right]^2}} $$ where $$ \beta_a = \frac{\omega}{ \omega_{\text{res}}}, ~~\tan\phi = \frac{2\xi\beta_a}{1 - \beta_a^2} \,. $$

Your solution

In your solution the homogeneous term is zero and therefore there is no decay term in your solution. Otherwise your implementation looks correct.

  • $\begingroup$ I have some questions about your explanation: - In the critical damping coeffcient, you have an extra term m for the mass which is unknown for me. - Is there a reason you are increasing the damping to 0.3? - Due to the initial condition of u_0 = 0, the curve has to start at a displacement that equals 0? $\endgroup$
    – Beunhaas
    Feb 23, 2016 at 20:29
  • $\begingroup$ And where is the term phi_d coming from? $\endgroup$
    – Beunhaas
    Feb 23, 2016 at 20:43
  • 1
    $\begingroup$ Your critical damping coeff $c_c = 2\sqrt(k m) = 2 m\sqrt{k/m} = 2 m \omega_{\text{res}}$. The damping coeff is higher just to make the damping more obvious in the figures. When using the initial condition to find $A$ and $B$ make sure that you use $u = u_h + u_p$. The $\phi_d$ term comes if you write your equation in the compact form in the updated answer. $\endgroup$ Feb 23, 2016 at 20:52
  • $\begingroup$ It is obvious to search A and B, using the total solution instead of only the homogeneous one. I changed it in my code. Unfortunately, I am still not getting the results I expect. $\endgroup$
    – Beunhaas
    Feb 24, 2016 at 18:19

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.