Closed Form Solution

%
%   Mass-Spring-Damper system 
%   Closed form solution 
%

%
%   constant inputs
%

m = 1;                      %  mass of the system in meter
k = 50;                     %  stiffness of the spring (N/m)
c = 20;                     %  damping coefficient
wn = sqrt(k/m);             %  natural frequency (rad/sec)
zeta = c/(2*wn*m)           %  viscous damping factor

tf = 2;                     %  calculation time
n = 500;                    %  number of data points

s1 = (-zeta+sqrt(zeta^2-1))*wn;
s2 = (-zeta-sqrt(zeta^2-1))*wn;

%
%   initial conditions
%

x_0 = 0;                    %  initial displacement
x_dot_0 = 1;                %  initial velocity
x0 = [x_0,x_dot_0];

%
%   calculation of amplitudes & x(t)
%

t = 0:t_final/n:t_final;

if zeta ~= 1                       % underdamped and overdamped cases

  A = [[1,1];[s1,s2]];
  X = A^(-1)*x0';

  X_1 = X(1);
  X_2 = X(2);

  x = real(X_1*exp(s1*t)+X_2*(exp(s2*t)))+...
      imag(X_1*exp(s1*t)+X_2*(exp(s2*t)));               %  displacement x(t)
  x_dot = real(s1*X_1*exp(s1*t)+s2*X_2*(exp(s2*t)))+...
          imag(s1*X_1*exp(s1*t)+s2*X_2*(exp(s2*t)));     %  velocity x_dot(t)

elseif zeta == 1                    % critically damped case

  X_1 = x_0;
  X_2 = x_dot_0+X_1*wn;
  x = (X_1+X_2.*t).*exp(-wn*t);                          %  displacement x(t)
  x_dot = (-wn*(X_1+X_2.*t)+X_2).*exp(-wn*t);            %  velocity x_dot(t)

end
%
%   plotting the results 
%

figure(1); clf; orient tall;
subplot(2,1,1),plot(t,x);
xlabel('Time (sec.)');
ylabel('Displacement (m)');
title('Mass-Spring-Damper system : Closed form result (underdamped case)')

subplot(2,1,2),plot(t,x_dot);
xlabel('Time (sec.)');
ylabel('Velocity (m/sec.)');