Mass-Spring system

Here is the mass-spring system example. The purpose of this example is to solve the free vibration, which means no excitation exists on the system. The procedure for deriving equation of motion for 1 degree of freedom mass-spring system will be described. The purpose of this example is to help you become familiar with setting up the equation of motion and solving it using Matlab.
Problem Statement


The below is the drawing of the mass-spring system. We will solve for the homogeneous solution of the system, which means the external excitation force is zero. As you all know, this homogeneous solution corresponds to the transient response of the system. We will define the natural frequency of the system. After that, we will set up equation of motion (FBD, force balance), get the closed form of the solution and try to solve it using Matlab.

In this example, we will use both closed form solution and ODE45 solver to compare the results. (If you have no idea what ODE45 solver is, then click here.)

Free Body Diagram and Equation of Motion


From the free body digram we can write the force balance equations as follow. Consider the right-going direction as the positive x-direction to get
\begin{displaymath}-kx+F(t)=m \ddot{x}\end{displaymath} (1)

where m is the mass of the system, k is the stiffness of the spring, and F(t) is the external force acting on the mass.

However, in the current problem we are focusing on the homogeneous solution, thus the external force is zero, i.e. F(t)=0. Then the above equation becomes
\begin{displaymath}m \ddot{x} +kx = 0\end{displaymath} (2)

and we have an equation of motion for the mass-spring system for the homogeneous solution. We can define the natural frequency of the system,$\omega_n$, as
\begin{displaymath}\omega_n^2 = \frac{k}{m}\end{displaymath} (3)

Then the equation of motion becomes
\begin{displaymath}\ddot{x} + \omega_n^2 x = 0\end{displaymath} (4)

Closed Form Solution


Now we need to think about the closed form solution of the equation above. The solution form can be written as a sinusoidal function. Therefore,
\begin{displaymath}x(t) = Acos(\omega_n t - \phi)\end{displaymath} (1)
\begin{displaymath}\dot{x}(t) = -\omega_n Asin(\omega_n t - \phi)\end{displaymath} (2)

where A and $\phi$ are amplitude and phase angle respectively. (Cosine was an arbitrary choice and sine would work equally well.)

A and $\phi$ are unknown and must be determined using the initial conditions x(0)=x0 and $\dot{x}(0)=\dot{x}_0$. In terms of these initial conditions, A and $\phi$ become
\begin{displaymath}A = \sqrt{x_0^2+\left(\frac{\dot{x}_0}{\omega_n}\right)^2}\end{displaymath} (3)
\begin{displaymath}\phi = cos^{-1} \left(\frac{x_0}{A}\right)\end{displaymath} (4)

Therefore, with the initial conditions given, we have a complete closed form solution for the response of the system.
 
 

We can code the closed form solution into MATLAB. The following is the MATLAB code for closed form solution of mass-spring system.

At the beginning of the program, we need to write some constant inputs. n is the number of the data points. 

m = 1;                      %  mass of the system in meter
k = 50;                     %  stiffness of the spring (N/m)
wn = sqrt(k/m);             %  natural frequency (rad/sec)

t_final = 2;                %  calculation time
n = 500;                    %  number of data points
Then, initial conditions are followed. 
x_0 = 1;                    %  initial displacement
x_dot_0 = 0;                %  initial velocity
Calculate the amplitude and the phase angle with given initial conditions. 
X_0 = sqrt(x_0^2+(x_dot_0/wn)^2);
phi = acos(x_0/X_0);
And we can determine the displacement x(t) and velocity x_dot(t). 
t = 0:t_final/n:t_final;
x = X_0*cos(wn*t+phi);           %  displacement x(t)
x_dot = -wn*x_0*sin(wn*t+phi);   %  velocity x_dot(t)
Click here for the entire MATLAB source code for closed form solution and plotting.

Integrating EOM with Matlab


Here we will use the ODE45 solver in Matlab to integrate the equations of motion. The following is the MATLAB code for ODE45 solver solution of mass-spring system. If you need to see the basic operation of ODE45 solver, see Matlab Assistance.

Some constant inputs need to be fixed before we proceed. 

m = 1;                      %  mass of the system in meter
k = 50;                     %  stiffness of the spring (N/m)
wn = sqrt(k/m);             %  natural frequency (rad/sec)
t_final = 2;                %  calculation time
And the initial conditions were shown. tol is the tolerance for ODE45 solver. Higher it is, more accurate the answer will be, more time the calculation takes. 
x_0 = 1;                    %  initial displacement
x_dot_0 = 0;                %  initial velocity
X_0 = [x_0,x_dot_0];        %  form a vector(array) of initial conditions
tol = 1e-6;                 %  tolerance of error for ODE45 solver
Using ODE45 solver, we can solve for the displacement x(t) and velocity x_dot(t). y is used as a dummy variable. 
[t,y] = ode45('pms_sol',0,t_final,X_0,tol);

x = y(:,1);                 %  displacement x(t)
x_dot  = y(:,2);            %  velocity x_dot(t)
Here's the subprogram pms_sol.m. 
function ydot = pms_sol(t,y)

m = 1;                      %  mass of the system in meter
k = 50;                     %  stiffness of the spring (N/m)

yd1 = y(2);
yd2 = -k*y(1)/m;

ydot = [yd1;yd2];
Here are the complete list of the main program (PC,Mac), main program (Unix) and subprogram. The main programs also include the plotting commands.

Plotting Results


Results are given in 2 parts. The first gives the plotting commands and results using the closed form solution, the second gives the plotting commands and results using the ODE solver. Just as you may have guessed, both methods should give the same answers. Check them out.

Return to Example Table