Table of contents

Transient and state based solutions
Example 1 - Linear friction
Example 2 - Contact force
Example 3 - Constrained motion
Example 4 - Throwing particles through the air
1. Finding the optimal launch angle
2. Wind resistance
Example 5
Example 6 - The Damper
Example 7 - Spring - Damper
Example 7 - The pendulum
Example 8 - The parachute
Example 9 - The bungee cord
Example 10 - The robot arm
Example 11 - The banked track
Example 12 - The roller coaster

Particle kinetics

The method for solving problems in kinetics is an extension from statics and the methods used in the corresponding sections. Kinetics is the study of relations between forces and accelerations and a starting point is Newton's equation of motion

\sum {\bm F}_i =m\;\bm a

where $\bm a$ is the center of mass (COM) of a body (or particle). For completeness we shall here state the complete generalized equations of motion extended by Euler et al.

\left\lbrace \begin{array}{ll} \bm F=\sum {\bm F}_i =m\;\ddot{\bm r} & \textrm{Translation}\\ \bm M=\sum {\bm M}_i +{\bm M}_e =\Sigma {\bm r}_i \times {\bm F}_i +{\bm M}_e =\bar{I} \ddot{\theta} & \textrm{Rotation} \end{array}\right.

where $\bar{I}$ is the angular moment of inertia for the body at the center of mass, otherwise known as the rotational inertia of a rigid body. We shall give a proper introduction of this quantity in the section of rigid body kinetics. Note that the moment of inertia is given by the geometric relation

\bar{I} =\int_V \rho |\bar{\bm R} |^2 \;dV

where $\bar{\bm R}$ is the vector from the center of rotation (which in this case is the COM) and all infinitesimal points on the body.

Notice that in this general formulation we have ${\left\lbrack F_x ,F_y ,F_z \right\rbrack }^t =m\ddot{\bm r}$ and ${\left\lbrack M_x ,M_y ,M_z \right\rbrack }^T =\bar{I} \ddot{\theta}$ , thus six equations (ODEs) per component, where for each component we derive these equations by the use of a free body diagram, exactly as in statics. In order to get a specific solution to the system of equations we need corresponding initial conditions and/or boundary conditions.

Suffice to say for now, if the body has near zero size, then $\bm R\to 0$ and we can assume particle dynamics and disregard the rotational part of the equation of motion, resulting in only three degrees of freedom instead of six. Similarly if the rotational acceleration $\ddot{\theta}$ is sufficiently low, we can assume particle dynamics. We will revisit the generalized equations of motion in rigid body kinetics.

Note that Newtons second law of motion does not explicitly relate force to acceleration, but rather momentum, paraphrased: "The time rate of change of a bodies momentum equals the force acting upon it" i.e.,

\bm F=\dfrac{d}{dt}\left(m\left(t\right)\;\dot{\bm r} \left(t\right)\right)\iff \bm F=\dot{m} \dot{\bm r} +m\ddot{\bm r}

Thus, we can create models where the body has a change of mass, e.g., rockets. When a body has constant mass the expression reduces to what we are familiar with.

The solution of the equations of motion is typically done by solving a differential equation using either a symbolic or numerical solver. In traditional courses on classical mechanics, countless chapters are spent on deriving various analytical methods and formulas for simplifying problems or solving problems specifically designed to be solved using simple hand calculations. We shall here instead focus on solving problems completely, and be explicit with the visualization of the whole motion.

Transient and state based solutions

We can apply two main methods for solving and analyzing dynamics problems. Transient, meaning we formulate an ODE and solve it either analytically of numerically. This gives us the position, velocity and acceleration vectors as a function of time over the complete time frame of interest. For situations where we are interested in solutions only at a particular time (or state) we can formulate the problem as an energy equation instead using Lagrangian or Hamiltonian formulation of mechanics. This can simplify the calculations drastically, since we do not need to iterate to the state from the initial state (using a numerical method) nor determine the accelerations or include forces that do no work.

Example 1 - Linear friction

A box is being pushed from rest.

How much force $P$ is required to accelerate the box to $a_0 \;m/s^2$ ?
What is the velocity of the box after $t_1$ seconds?
How far does the box glide? At which time does it come to a stop?

Step 1 - Free body diagram

Step 2 - Equations of motion

From the figure we get

\begin{align*} & \bm F_{\text{ext}}=P[\cos\theta,-\sin\theta]^{T}\\ & \bm F_{N}=[0,N]^{T}\\ & \bm F_{W}=[0,-mg]^{T}\\ & \bm F_{f}=[-\mu_{k}N,0]^{T} \end{align*}

and the the acceleration vector is

\bm a = [a_x, 0]^T

Using Newtons second law of motion we can formulate

\begin{gather*} \sum\bm{F}=m\bm{a} \Rightarrow\\ P\begin{bmatrix}\cos\theta\\ -\sin\theta \end{bmatrix}+\begin{bmatrix}0\\ N \end{bmatrix}+\begin{bmatrix}0\\ -mg \end{bmatrix}+\begin{bmatrix}-\mu_{k}N\\ 0 \end{bmatrix}=\begin{bmatrix}ma_{x}\\ 0 \end{bmatrix} \end{gather*}

clear

syms mu_k m g v_0 P theta N a_x
F_ext = P*[cos(theta); -sin(theta)];
F_N = [0; N];
F_W = [0; -m*g];
F_f = [-mu_k*N; 0];
a = [a_x; 0];

ekv = F_ext + F_N + F_W + F_f == m*a

\left(\begin{array}{c} P\,\cos \left(\theta \right)-N\,\mu_k =a_x \,m\\ N-g\,m-P\,\sin \left(\theta \right)=0 \end{array}\right)

Step 3 - Solve for acceleration

From this equation system we need to solve for the same amount of unknowns as we have number of equations, we choose $a_x$ , of course, and $N$ for the additional unknown. This is typical in mechanics to get additional unknowns, usually contact forces, or reaction forces. The choice is based on careful analysis of what is considered input parameter and what is an unknown quantity.

[a_x, N] = solve(ekv, [a_x, N])

a_x =

-\dfrac{P\,\mu_k \,\sin \left(\theta \right)-P\,\cos \left(\theta \right)+g\,m\,\mu_k }{m}

N =

g\,m+P\,\sin \left(\theta \right)

From the expression we can see that the vertical component of $P$ causes additional friction in the negative x-direction whereas the horizontal component gives rise to an acceleration in the positive x-direction. Analyzing symbolic expressions before substitution is a good idea in the modeling stage, since we can verify the model much easier.

Step 4 - Formulating a set of differential equations (kinematics)

Let us assume that the box is given a push from $P$ during a short time. Let $P=\left\lbrace \begin{array}{ll} 4000 & t\le 0\ldotp 3\\ 0 & t>0\ldotp 3 \end{array}\right.$

syms t
P_const = piecewise(t<0.3, 4000, ...
                    t>0.3 , 0)

P_const =

\left\lbrace \begin{array}{cl} 4000 & \;\textrm{if}\;\;t<\dfrac{3}{10}\\ 0 & \;\textrm{if}\;\;\dfrac{3}{10}< t \end{array}\right.

In order to solve this problem we need to tackle this piecewise, we need to set up a differential equation with initial values for the range $t\in[0,0.3]$ s and another differential equation with initial values for the range $t\in(0.3,t_\text{end}]$ .

The strategy:

We can then formulate the differential equation and solve it to get the position $x\left(t\right)$ which we can differentiate to get $v\left(t\right)$ .

Then we can find out at what time the velocity is zero, $t_\text{end}$ , it is at this time that the box has come to stop.

⚠ Note

After $t_\text{end}$ our friction model is no longer physical! That would mean that the box starts sliding backwards, which is absurd.

We formulate the differential equation in general by

\ddot{x} = \dfrac{d^2}{dt^2}x(t)=a_{x}

Now, we shall use some numerical values for brevity.

m=50; mu_k=0.3; g = 9.82; theta = 0*pi/180;

a_x = subs(a_x)

a_x = \dfrac{P}{50}-\dfrac{1473}{500}

We split the solution into before and after $t=0\ldotp 3\mathrm{s}$ .

For $t\in[0,0.3]$ we have

syms x(t)
DE1 = diff(x,t,2) == a_x

50\,\dfrac{\partial^2 }{\partial t^2 }\;x\left(t\right)=\dfrac{P}{50}-\dfrac{1473}{500}

x_dot = diff(x,t);
BV1 = [x(0) == 0, x_dot(0) == 0];
x1(t) = dsolve(DE1,BV1)

\dfrac{t^2 \,{\left(10\,P-1473\right)} }{1000}

x1(t) = subs(x1, P, 4000)

\dfrac{38527\,t^2 }{1000}

v1(t) = diff(x1,t)

\dfrac{38527\,t}{500}

For $t>0.3$ we have that $P=0$

DE2 = vpa(diff(x,t,2) == subs(a_x,P,0) ,2)

\dfrac{\partial^2 }{\partial t^2 }\;x\left(t\right)=-2.9

t1 = 0.3;
BV2 = vpa([x(t1) == x1(t1), x_dot(t1) == v1(t1)],2)

\left(\begin{array}{cc} x\left(\dfrac{3}{10}\right)=0.069 & {\left({ {\left(\dfrac{\partial }{\partial t}\;x\left(t\right)\right)}\left|\right.}_{t=\dfrac{3}{10} } \right)}=0.46 \end{array}\right)

x2(t) = dsolve(DE1,BV2);
x2 = vpa(x2,2)

-1.5\,t^2 +24.0\,t-3.6

v2(t) = diff(x2,t);
v2 = vpa(v2,2)

24.0-2.9t

To get the time where the box stops gliding

t2 = double(solve(v2==0, t))

t2 = 8.1466

x_max = double(x2(t2))

x_max = 94.1597

Now we can plot the situation!

figure; hold on
fplot(x1,[0,t1],'Linewidth',2)
fplot(x2,[t1,t2],'Linewidth',2)
xlabel('$t$','interpreter','latex')
ylabel('$x(t)$','interpreter','latex')
title('Position of the box')

figure; hold on
fplot(v1,[0,t1],'Linewidth',2)
fplot(v2,[t1,t2],'Linewidth',2)
xlabel('$t$','interpreter','latex')
ylabel('$v(t)$','interpreter','latex')
title('Velocity of the box')

Modifying the force to be piecewise linear instead

We can model this push as a linear force instead. Let us create a linear function for the load.

syms t1 t2 t
A = [1 t1
     1 t2];
C = A\eye(2);
phi = simplify([1 t]*C);
P1Interp = phi*[sym('u1');sym('u2')]

\dfrac{u_1 \,{\left(t-t_2 \right)} }{t_1 -t_2 }-\dfrac{u_2 \,{\left(t-t_1 \right)} }{t_1 -t_2 }

P1Interp = matlabFunction(P1Interp)

P1Interp = 
    @(t,t1,t2,u1,u2)(u1.*(t-t2))./(t1-t2)-(u2.*(t-t1))./(t1-t2)

P_lin = piecewise(t<0.15, P1Interp(t,0,0.15,0,4000), ...
                  t<0.3 , P1Interp(t,0.15,0.3,4000,0), ...
                  t>0.3 , 0)

\def\arraystretch{2.2} \left\lbrace \begin{array}{cl} \dfrac{80000\,t}{3} & \;\textrm{if}\;\;t< \dfrac{3}{20}\\ 8000-\dfrac{80000\,t}{3} & \;\textrm{if}\;\;t< \dfrac{3}{10}\\ 0 & \;\textrm{if}\;\;\dfrac{3}{10}< t \end{array}\right.

figure
fplot(P_lin,[0,2],'b')
xlabel('$t$ [s]','interpreter','latex')
ylabel('$P(t)$ [N]','interpreter','latex')
title('Piecewise linear load')

syms P
a_x = subs(a_x);

t1 = 0.15; t2 = 0.3;

DE1 = diff(x,t,2) == a_x;
x_dot = diff(x,t);
BV1 = [x(0) == 0, x_dot(0) == 0];
x1(t) = dsolve(DE1,BV1);
x1(t) = subs(x1, P, P1Interp(t,0,t1,0,4000))

x_1(t) = d\frac{t^2 \,{\left(\dfrac{800000\,t}{3}-1473\right)}}{1000}

v1(t) = diff(x1,t)

x_1(t) = \dfrac{t\,{\left(\dfrac{800000\,t}{3}-1473\right)}}{500}+\dfrac{800\,t^2 }{3}

DE2 = diff(x,t,2) == a_x;
BV2 = [x(t1) == x1(t1), x_dot(t1) == v1(t1)];
x2(t) = dsolve(DE2,BV2);
x2(t) = subs(x2, P, P1Interp(t,t1,t2,4000,0))

x_2(t) = t\,{\left(80\,t-6\right)}-6\,t-t^2 \,{\left(\dfrac{800\,t}{3}-\dfrac{78527}{1000}\right)}

v2(t) = diff(x2,t)

v_2(t) = 160\,t-2\,t\,{\left(\dfrac{800\,t}{3}-\dfrac{78527}{1000}\right)}-\dfrac{800\,t^2 }{3}-12

DE3 = diff(x,t,2) == subs(a_x,P,0);
BV3 = [x(t2) == x2(t2), x_dot(t2) == v2(t2)];
x3(t) = dsolve(DE3,BV3)

x_3(t) = 12\,t-\dfrac{1473\,t^2 }{1000}

v3(t) = diff(x3,t)

v_3(t) = 12-\dfrac{1473\,t}{500}

t3 = double(solve(v3==0, t))

t3 = 4.0733

x_max = double(x3(t3))

x_max = 24.4399

figure;hold on
fplot(x1,[0,t1],'Linewidth',2)
fplot(x2,[t1,t2],'Linewidth',2)
fplot(x3,[t2,t3],'Linewidth',2)
xlabel('$t$','interpreter','latex')
ylabel('$x(t)$','interpreter','latex')

figure; hold on
fplot(v1,[0,t1],'Linewidth',2)
fplot(v2,[t1,t2],'Linewidth',2)
fplot(v3,[t2,t3],'Linewidth',2)
xlabel('$t$','interpreter','latex')
ylabel('$v(t)$','interpreter','latex')
title('Velocity of the box')

Comparing the two models, we note the significant difference in both velocity and sliding distance and for the case where the load is piece-wise linear, the box is starting to de-accelerate 0.15 seconds earlier.

Example 2 - Contact force

The three boxes are accelerated on the smooth surfaces by the force $P$ . Determine the motion and contact forces between the boxes.

1 - Free body diagram

Do determine the contact forces between the boxes we need to split the system up in three subsystems where we can apply Newton II och each individual box. Using Newton III on each subsystem we get the following situation with forces.

2 - Equations of motion

From these three subsystems, we formulate three kinetic equations as differential equations with the initial conditions stating that all three boxes have zero initial velocity and position.

\begin{align*} \sum F & = ma \\ & \Rightarrow \\ P-N_{1} & = m_{1}a_{1}\\ N_{1}-N_{2} & = m_{2}a_{2}\\ N_{2} & = m_{3}a_{3} \end{align*}

3 - Adding differential equations (kinematics)

clear
syms x1(t) x2(t) x3(t)
syms N1 N2 P m1 m2 m3
v1 = diff(x1,t); v2 = diff(x2,t); v3 = diff(x3,t);
a1 = diff(x1,t,2); a2 = diff(x2,t,2); a3 = diff(x3,t,2);

DE = [P - N1  == m1*a1
      N1 - N2 == m2*a2 
      N2      == m3*a3];

BV =

\def\arraystretch{2.4} \left(\begin{array}{c} P-N_1 =m_1 \,\dfrac{\partial^2 }{\partial t^2 }\;x_1 \left(t\right)\\ N_1 -N_2 =m_2 \,\dfrac{\partial^2 }{\partial t^2 }\;x_2 \left(t\right)\\ N_2 =m_3 \,\dfrac{\partial^2 }{\partial t^2 }\;x_3 \left(t\right) \end{array}\right)

BV = [v1(0)==0; x1(0)==0;
      v2(0)==0; x2(0)==0;
      v3(0)==0; x3(0)==0;];
s = dsolve(DE,BV);
x1 = s.x1

x1 =

-\dfrac{t^2 \,{\left(N_1 -P\right)} }{2\,m_1 }

x2 = s.x2

\dfrac{t^2 \,{\left(N_1 -N_2 \right)} }{2\,m_2 }

x3 = s.x3

\dfrac{N_2 \,t^2 }{2\,m_3 }

4 - Geometric coupling

To get the forces we need to couple the positions, creating two equations, from which we can solve for the forces.

\def\arraystretch{2.5} \begin{align*} x_{1} & =x_{2}\Rightarrow-\dfrac{t^{2}\left(N_{1}-P\right)}{2m_{1}}=\dfrac{t^{2}\left(N_{1}-N_{2}\right)}{2m_{2}}\\ x_{2} & =x_{3}\Rightarrow\dfrac{t^{2}\left(N_{1}-N_{2}\right)}{2m_{2}}=\dfrac{N_{2}t^{2}}{2m_{3}} \end{align*}

[N1,N2] = solve([x1==x2, x2==x3], [N1,N2])

N1 =

\dfrac{P\,{\left(m_2 +m_3 \right)} }{m_1 +m_2 +m_3 }

N2 =

\dfrac{P\,m_3 }{m_1 +m_2 +m_3 }

Example 3 - Constrained motion

A log with a mass $m_1 =200$ kg is to be pulled up a ${30}^{\circ }$ ramp with a friction coefficient of $\mu_k =0\ldotp 5$ using a pulley system according to the figure and a counter-mass of $m_2 =125$ kg. The counter-mass travels 6 m, determine the velocity as it reaches the ground, starting from a resting position.

1 - Free body diagram

2 - Equations of Motion (Kinetics)

For the log we have:

\underset{\boldsymbol{F}_{W1}}{\underbrace{m_{1}g\begin{bmatrix}-\sin\theta\\ -\cos\theta \end{bmatrix}}}+\underset{\boldsymbol{F}_{N}}{\underbrace{\begin{bmatrix}0\\ N \end{bmatrix}}}+\underset{\boldsymbol{F}_{T}}{\underbrace{\begin{bmatrix}2T\\ 0 \end{bmatrix}}}+\underset{\boldsymbol{F}_{f}}{\underbrace{\begin{bmatrix}-\mu_{k}N\\ 0 \end{bmatrix}}}=\underset{m\boldsymbol{a}_{c}}{\underbrace{\begin{bmatrix}m_{1}a_{c}\\ 0 \end{bmatrix}}}

For the counterweight we have:

F_{W2}-T=m_{2}a_{B}

clear
theta = 30;
syms m1 m2 g N mu_k T a_C a_A
F_W1 = m1*g*[-sind(theta); -cosd(theta)];
F_f = [-mu_k*N; 0];
F_N = [0; N];
F_T = [2*T; 0];
eq = F_W1 + F_f + F_T + F_N == m1*[a_C; 0]
eq = [eq;
      m2*g - T == m2*a_A]

eq =

\left(\begin{array}{c} 2\,T-N\,\mu_k -\dfrac{g\,m_1 }{2}=a_C \,m_1 \\ N-\dfrac{\sqrt{3}\,g\,m_1 }{2}=0 \\ g\,m_2 -T=a_A \,m_2 \end{array}\right)

3 - Geometric constraints - Kinematics

We begin by modeling the relation between the point $C$ and $A$ .

Note that the rope length is constant, let $L = 2s_C + s_A$

syms s_C(t) s_A(t) L
eq = L == 2*s_C + s_A

eq(t) =

L=s_A \left(t\right)+2\,s_C \left(t\right)

We need to get ${\ddot{s} }_C$ and ${\ddot{s} }_A$ to use in the kinetics formulations. Let's differentiate L twice

diff(eq,t,2)

0=\dfrac{\partial^2 }{\partial t^2 }\;s_A \left(t\right)+2\,\dfrac{\partial^2 }{\partial t^2 }\;s_C \left(t\right)

We note that as $a_A$ increases, $a_C$ needs to decrease according to the geometric relation in the figure, thus we get:

a_A = 2 a_C

Solving for the unknowns - formulating ODE and solving

The complete system of equations we need to deal with:

eq = [eq;
      a_A == -2*a_C]

\left(\begin{array}{c} 2\,T-N\,\mu_k -\dfrac{g\,m_1 }{2}=a_C \,m_1 \\ N-\dfrac{\sqrt{3}\,g\,m_1 }{2}=0\\ g\,m_2 -T=a_A \,m_2 \\ a_A = 2 a_C \end{array}\right)

Solving for the unknown variables

[T,N,a_C,a_A] = solve(eq, [T,N,a_C,a_A] )

T =

-\frac{\sqrt{3}\,g\,m_1 \,m_2 \,\mu_k }{m_1 -4\,m_2 }

N =

\dfrac{\sqrt{3}\,g\,m_1 }{2}

a_C =

-\dfrac{g\,m_1 -4\,g\,m_2 +\sqrt{3}\,g\,m_1 \,\mu_k }{2\,{\left(m_1 -4\,m_2 \right)} }

a_A =

\dfrac{g\,m_1 -4\,g\,m_2 +\sqrt{3}\,g\,m_1 \,\mu_k }{m_1 -4\,m_2 }

m1 = 200; m2 = 125; mu_k = 0.5; g=9.81;

T = double(subs(T))

T = 707.9758

N = double(subs(N))

N = 1.6991e+03

a_C = double(subs(a_C))

a_C = -2.0731

a_A = double(subs(a_A))

a_A = 4.1462

Now we can formulate the ODE

syms s_C(t)
v_A= diff(s_A,t);
s_A(t) = vpa(dsolve(diff(s_A,t,2) == a_A, v_A(0)==0, s_A(0)==0), 2)

s_A(t) =

2.0730969295836985111236572265625\,t^2

v_A(t) = diff(s_A,t)

v_A(t) =

4.146193859167397022247314453125\,t

t1 = double(max(solve(s_A==6)))

t1 = 1.7012

v_A_end = vpa(v_A(t1),3)

7.05

figure; hold on
fplot(s_A,[0,t1],'Linewidth',2)
plot(t1, s_A(t1),'o','MarkerFaceColor','k')
xlabel('$t$ [s]','interpreter','latex')
ylabel('$s_C(t)$ [m]','interpreter','latex')
title('Distance of the log')
grid on;

Example 4 - Throwing particles through the air

Plot the path of the projectile until it reaches the ground.
How long does the projectile hurl through the air?
How far away does the projectile impact ground?
What is the optimal angle $\theta$ to maximize the distance?
Add wind resistance according to $F_W = c{\bm v}_p^2 + c{\bm v}_w^2$ , where $c=0\ldotp 03$ is some geometric constant, ${\bm v}_p$ is the velocity of the particle and ${\bm v}_w$ is the velocity of the wind.

clear
v0 = 24;
theta = 35;
syms theta v0
V0 = v0*[cosd(theta); sind(theta)]

V0 =

\left(\begin{array}{c} v_0 \,\cos \left(\dfrac{\pi \,\theta }{180}\right)\\ v_0 \,\sin \left(\dfrac{\pi \,\theta }{180}\right) \end{array}\right)

syms x(t) y(t)
r(t) = [x(t); y(t)];
r_dot(t) = diff(r,t);
r_ddot = diff(r,t,2);

syms m g
DE = r_ddot*m == [0; -m*g]

DE(t) =

\def\arraystretch{2.3} \left(\begin{array}{c} m\,\dfrac{\partial^2 }{\partial t^2 }\;x\left(t\right)=0\\ m\,\dfrac{\partial^2 }{\partial t^2 }\;y\left(t\right)=-g\,m \end{array}\right)

BV = [r_dot(0) == V0;
      r(0) == [0;2]]

\left(\begin{array}{c} {\left({ {\left(\dfrac{\partial }{\partial t}\;x\left(t\right)\right)}\left|\right.}_{t=0} \right)}=v_0 \,\cos \left(\dfrac{\pi \,\theta }{180}\right)\\ {\left({ {\left(\dfrac{\partial }{\partial t}\;y\left(t\right)\right)}\left|\right.}_{t=0} \right)}=v_0 \,\sin \left(\dfrac{\pi \,\theta }{180}\right)\\ x\left(0\right)=0\\ y\left(0\right)=2 \end{array}\right)

[x(t),y(t)] = dsolve(DE,BV)

x(t) =

t\,v_0 \,\cos \left(\dfrac{\pi \,\theta }{180}\right)

y(t) =

-\dfrac{g\,t^2 }{2}+v_0 \,\sin \left(\dfrac{\pi \,\theta }{180}\right)\,t+2

theta = 35; g = 9.82; v0 = 24; m = 4;
xn = subs(x); yn = subs(y);

t1 = double(max(solve(yn == 6)))

t1 = 2.4744

x1 = double(subs(xn,t,t1))

x1 = 48.6457

figure; hold on
fplot(xn,yn,[0,t1],'Displayname', '$\theta=35^\circ$');
hp = plot(NaN, NaN,'--r','Linewidth',2,'Displayname', ['$\theta=',num2str(0),'^\circ$']);
plot([0,30,40,250],[0,0,6,6],'-k','Handlevisibility','off')
plot([40,40],[6,8.4],'-k','Handlevisibility','off')
axis equal
axis([0,70,0,25])
xlabel('$x$ [m]','Interpreter','latex');
ylabel('$y$ [m]','Interpreter','latex');

Finding the optimal launch angle

% syms g v0 m theta
% xn = subs(x,[g,v0,m],[9.82, 24, 4])

24\,t\,\cos \left(\dfrac{\pi \,\theta }{180}\right)

% yn = subs(y,[g,v0,m],[9.82, 24, 4])

-\dfrac{491\,t^2 }{100}+24\,\sin \left(\dfrac{\pi \,\theta }{180}\right)\,t+2

% xn = matlabFunction(xn);
% yn = matlabFunction(yn);
% tn = linspace(0,10,100);
% legend('show','Interpreter','latex')
% for thetan = [0:1:90, 89:-1:0]
%     X = xn(tn, thetan);
%     Y = yn(tn, thetan);
%     set(hp,'XData',X,'YData',Y,'Displayname', ['$\theta=',num2str(thetan),'^\circ$'])
%     drawnow
% end

Implementing a visualization we can quickly get a very good overview over the situation and graphically read the maximum height, distance, valid path etc.

To get the optimal angle to reach the maximum distance analytically we can get an expression for the time $t\left(\theta \right)$ when the height is 6. Note that for some angles the particle never reaches this height, mathematically this is equivalent to not finding any roots, we need to be aware. We then substitute the expression for time into the expression for $x\left(t\left(\theta \right)\right)$ and plot the graph.

t2 = solve(yn==6, t)

\left(\begin{array}{c} \dfrac{1200\,\sin \left(\dfrac{\pi \,\theta }{180}\right)}{491}-\dfrac{50\,\sqrt{576\,{\sin \left(\dfrac{\pi \,\theta }{180}\right)}^2 -\dfrac{1964}{25} } }{491}\\ \dfrac{1200\,\sin \left(\dfrac{\pi \,\theta }{180}\right)}{491}+\dfrac{50\,\sqrt{576\,{\sin \left(\dfrac{\pi \,\theta }{180}\right)}^2 -\dfrac{1964}{25} } }{491} \end{array}\right)

t2 = t2(2)

\dfrac{1200\,\sin \left(\dfrac{\pi \,\theta }{180}\right)}{491}+\dfrac{50\,\sqrt{576\,{\sin \left(\dfrac{\pi \,\theta }{180}\right)}^2 -\dfrac{1964}{25} } }{491}

xn = subs(x,[g,v0,m],[9.82, 24, 4])

24\,t\,\cos \left(\dfrac{\pi \,\theta }{180}\right)

xn = formula(subs(xn, t, t2))

24\,\cos \left(\dfrac{\pi \,\theta }{180}\right)\,{\left(\dfrac{1200\,\sin \left(\dfrac{\pi \,\theta }{180}\right)}{491}+\dfrac{50\,\sqrt{576\,{\sin \left(\dfrac{\pi \,\theta }{180}\right)}^2 -\dfrac{1964}{25} } }{491}\right)}

figure
fplot(xn,[0,90])
xlabel('$\theta$ [deg]','Interpreter','latex');
ylabel('$x$ [m]','Interpreter','latex');

theta_max = double(vpasolve(diff(xn)==0,45))

theta_max = 47.0985

x_max = double(subs(xn,theta, theta_max))

x_max = 54.5092

The maximum distance is at 54 m and is reached by setting the launch angle to ${47}^{\circ }$ .

Wind resistance

Adding wind resistance according to $F_W =c{\bm v}_p^2 +c{\bm v}_w^2$ , where $c=0\ldotp 03$ is some constant, ${\bm v}_p$ is the velocity of the particle and ${\bm v}_w$ is the velocity of the wind. We formulate the ODE

\begin{array}{l} m\ddot{x} =-c\;{ {\dot{x} }_p }^2 +c\;{v_W }^2 \\ m\ddot{y} =-m\;g-c{ {\dot{y} }_p }^2 +c\;{v_W }^2 \end{array}

This becomes a non-linear ODE since we have to deal with the squared terms ${\dot{x} }_p$ and ${\dot{y} }_p$ which are the particles velocity. The sign is negative since the wind acts against the particle velocity. In addition to this, the surrounding air might actually have a velocity on its own, modeled by the vector ${\bm v}_W$ . Here $c$ denotes a constant which actually is a combination of the drag coefficient, air density and the projectile cross-sectional area.

Let's set the wind velocity to ${\bm v}_W =\left\lbrack -2,0\right\rbrack$

clear

syms theta v0
c = 0.03;
V0 = v0*[cosd(theta); sind(theta)]

V0 =

\left(\begin{array}{c} v_0 \,\cos \left(\dfrac{\pi \,\theta }{180}\right)\\ v_0 \,\sin \left(\dfrac{\pi \,\theta }{180}\right) \end{array}\right)

syms x(t) y(t)
r(t) = [x(t); y(t)];
r_dot(t) = diff(r,t);
r_ddot = diff(r,t,2);

v_W = [-2;0];
syms m g
DE = r_ddot*m == [0; -m*g] + c* -r_dot.^2 + v_W

\def\arraystretch{2.3} \left(\begin{array}{c} m\,\dfrac{\partial^2 }{\partial t^2 }\;x\left(t\right)=-\dfrac{3\,{ {\left(\dfrac{\partial }{\partial t}\;x\left(t\right)\right)} }^2 }{100}-2\\ m\,\dfrac{\partial^2 }{\partial t^2 }\;y\left(t\right)=-\dfrac{3\,{ {\left(\dfrac{\partial }{\partial t}\;y\left(t\right)\right)} }^2 }{100}-g\,m \end{array}\right)

BV = [r_dot(0) == V0;
      r(0) == [0;2]]

\def\arraystretch{2.3} \left(\begin{array}{c} {\left({ {\left(\dfrac{\partial }{\partial t}\;x\left(t\right)\right)}\left|\right.}_{t=0} \right)}=v_0 \,\cos \left(\dfrac{\pi \,\theta }{180}\right)\\ {\left({ {\left(\dfrac{\partial }{\partial t}\;y\left(t\right)\right)}\left|\right.}_{t=0} \right)}=v_0 \,\sin \left(\dfrac{\pi \,\theta }{180}\right)\\ x\left(0\right)=0\\ y\left(0\right)=2 \end{array}\right)

[x(t),y(t)] = dsolve(DE,BV)

Warning: Unable to find symbolic solution.
 
x(t) =
 
[ empty sym ]
 
 
y(t) =
 
[ empty sym ]

As we can see, dsolve is not able to find a symbolic solution, we need to tackle this numerically.

Applying Euler forward approximation we get

{\ddot{r} }_x =a_x =\dfrac{d}{dt}\left(v_x \right)\approx \dfrac{\Delta v_x }{\Delta t}=\dfrac{v_{x,i+1} -v_{x,i} }{\Delta t}\Rightarrow v_{x,i+1} =v_{x,i} +a_x \Delta t

and similarly

{\dot{r} }_x =v_x =\dfrac{d}{dt}\left(x\right)\approx \dfrac{\Delta x}{\Delta t}=\dfrac{x_{i+1} -x_i }{\Delta t}\Rightarrow x_{i+1} =x_i +v_x \Delta t

From

\def\arraystretch{1.3} \begin{array}{l} m\ddot{x} =-c\;{ {\dot{x} }_p }^2 +c\;{v_W }^2 \\ m\ddot{y} =-m\;g-c{ {\dot{y} }_p }^2 +c\;{v_W }^2 \end{array}

we have

a_x =\dfrac{-c\;{ {\dot{x} }_p }^2 +c\;{v_W }^2 }{m}

and

a_y =\dfrac{-m\;g-c{ {\dot{y} }_p }^2 +c\;{v_W }^2 }{m}

clear
v0 = 24; c = 0.03; theta = 47; m = 4; g = 9.82;
V0 = v0*[cosd(theta); sind(theta)];
v_W = [-20;22];

figure; hold on
hpath = plot(NaN, NaN,'-r','Linewidth',2,'Displayname', ['$\theta=',num2str(0),'^\circ$']);
hp = plot(NaN, NaN, 'ok','MarkerFaceColor','k');
plot([0,30,40,250],[0,0,6,6],'-k','Handlevisibility','off')
plot([40,40],[6,8.4],'-k','Handlevisibility','off')
axis equal; axis([0,70,0,25])
xlabel('$x$ [m]','Interpreter','latex');
ylabel('$y$ [m]','Interpreter','latex');
title(['$\theta=',num2str(theta),'^\circ, v_0=',num2str(v0), ...
    ' \mathrm{ m/s}, v_W=[',num2str(v_W(1)),',',num2str(v_W(2)),'] \mathrm{ m/s}$'], ...
    'interpreter','latex')

t = 0; dt = 0.1; i = 1;
vx = V0(1); vy = V0(2);
x = 0; y = 2;
X = NaN(100,1); Y = X;
while 1
    ax = (-c*vx^2 + c*sign(v_W(1))*v_W(1)^2)/m;
    ay = (-m*g -c*vy^2 + c*sign(v_W(2))*v_W(2)^2)/m;
    
    vx = vx + ax*dt;
    vy = vy + ay*dt;
    
    x = x + vx*dt;
    y = y + vy*dt;
    
    X(i) = x; Y(i) = y;
    set(hp,'XData',x,'YData',y)
    set(hpath,'XData',X,'YData',Y)
    % pause(0.1)
    drawnow

    t = t + dt; i = i + 1;
    if (x > 10 && y <= 0)
        break
    end
end

We note that even with no wind at all we can no longer make it over the obstacle, no matter what launch angle we choose, with the original initial velocity. We either need to make design changes to increase the initial velocity or hope for favorable winds.

Note that the parameter $c$ significantly affects the path, play around with the parameters!

Note how in general the path does not form a parabolic path which is why the swedish word "kastparabel" is very wrong to use.

Example 5

Small packages are released from rest at $A$ and slide down the smooth circular surface of radius $R$ to a conveyor at $B$ . Determine the expression for the the normal contact force $N$ between the guide and each package in terms of $\theta$ and specify the angular velocity $\omega$ of the conveyor wheel of radius $r$ to prevent any sliding on the belt as the object transitions from the circular path to the conveyor.

Solution: We start with the kinematics and describe the position of the package. With the origin according to the figure, we express the angle as a function of time and express the positin in Cartesion coordinates.

e_R ={\left\lbrack -\cos \theta \left(t\right),-\sin \theta \left(t\right)\right\rbrack }^T

We state the normal and tangential directions next, $e_n =-e_R$ and the tangential we can easily get by $e_t =e_n \times e_z$ according to the right hand rule.

Then we can formulate the objects velocity and acceleration by differentiating the position vector.

clear
syms theta(t)
syms R r 

e_R = [-cos(theta); -sin(theta); 0];
e_n = -e_R;
e_t = cross(e_n,[0;0;1])

e_t(t) =

\left(\begin{array}{c} \sin \left(\theta \left(t\right)\right)\\ -\cos \left(\theta \left(t\right)\right)\\ 0 \end{array}\right)

We get the position by

p = R*e_R

p(t) =

\left(\begin{array}{c} -R\,\cos \left(\theta \left(t\right)\right)\\ -R\,\sin \left(\theta \left(t\right)\right)\\ 0 \end{array}\right)

Then the velocity

p_dot = subs(diff(p,t), diff(theta), sym('theta_dot'))

p_dot(t) =

\left(\begin{array}{c} R\,\dot{\theta} \,\sin \left(\theta \left(t\right)\right)\\ -R\,\dot{\theta} \,\cos \left(\theta \left(t\right)\right)\\ 0 \end{array}\right)

And the acceleration

diff(p,t,2)

ans(t) =

\def\arraystretch{2.5} \left(\begin{array}{c} R\,\sin \left(\theta \left(t\right)\right)\,\dfrac{\partial^2 }{\partial t^2 }\;\theta \left(t\right)+R\,\cos \left(\theta \left(t\right)\right)\,{ {\left(\dfrac{\partial }{\partial t}\;\theta \left(t\right)\right)} }^2 \\ R\,\sin \left(\theta \left(t\right)\right)\,{ {\left(\dfrac{\partial }{\partial t}\;\theta \left(t\right)\right)} }^2 -R\,\cos \left(\theta \left(t\right)\right)\,\dfrac{\partial^2 }{\partial t^2 }\;\theta \left(t\right)\\ 0 \end{array}\right)

p_ddot = subs(diff(p,t,2), [diff(theta,2), diff(theta)], ...
                           [sym('theta_ddot'), sym('theta_dot')])

p_ddot(t) =

\left(\begin{array}{c} R\,\cos \left(\theta \left(t\right)\right)\,{\dot{\theta} }^2 +R\,\ddot{\theta} \,\sin \left(\theta \left(t\right)\right)\\ R\,{\dot{\theta} }^2 \,\sin \left(\theta \left(t\right)\right)-R\,\ddot{\theta} \,\cos \left(\theta \left(t\right)\right)\\ 0 \end{array}\right)

Now we can just project these onto the natural coordinates

v_n = p_dot.' * e_n

0

v_t = simplify(p_dot.' * e_t)

R\,\dot{\theta}

a_n = simplify(p_ddot.' * e_n)

R\,{\dot{\theta} }^2

a_t = simplify(p_ddot.' * e_t)

R\,\ddot{\theta}

Kinetics - From the free body diagram we get the forces

syms m g N
F_W = m*g*[0;-1;0]

\left(\begin{array}{c} 0\\ -g\,m\\ 0 \end{array}\right)

F_N = N*e_n

\left(\begin{array}{c} N\,\cos \left(\theta \left(t\right)\right)\\ N\,\sin \left(\theta \left(t\right)\right)\\ 0 \end{array}\right)

We want to convert these into natural coordinates as to state the equations of motion in natural coordinates

F_W_n = [F_W.' * e_n
         F_W.' * e_t]

\left(\begin{array}{c} -g\,m\,\sin \left(\theta \left(t\right)\right)\\ g\,m\,\cos \left(\theta \left(t\right)\right) \end{array}\right)

F_N_n = [simplify(F_N.' * e_n)
         F_N.' * e_t]

\left(\begin{array}{c} N\\ 0 \end{array}\right)

We can now formulate the equations of motion

eq = F_W_n + F_N_n == m*[a_n
                         a_t];

And solve for the normal force as well as the angular acceleration

eq = formula(subs(eq, theta, 'theta'))

\left(\begin{array}{c} N-g\,m\,\sin \left(\theta \right)=R\,m\,{\dot{\theta} }^2 \\ g\,m\,\cos \left(\theta \right)=R\,m\,\ddot{\theta} \end{array}\right)

[N,theta_ddot] = solve(eq, [N, 'theta_ddot'])

N =

R\,m\,{\dot{\theta} }^2 +g\,m\,\sin \left(\theta \right)

theta_ddot =

\dfrac{g\,\cos \left(\theta \right)}{R}

Now to eliminate $\dot{\theta}$ , we need to find out the angular velocity at $B$ , using $v\;dv=a\;ds$ we get

\int_0^v v\,dv=\int_0^{\theta } a_t \,ds

turns into:

\int_0^v {vdv}=\int_0^s R\ddot{\theta} ds

and with $ds={Rd}\theta$ and we have

\int_0^v {vdv}=\int_0^{\theta } R\ddot{\theta} {Rd}\theta

where $v=R\dot{\theta}$

syms v theta theta_dot
eq = int(v,v,0,v) == int(R*theta_ddot*R, theta, 0, theta)

\dfrac{v^2 }{2}=R\,g\,\sin \left(\theta \right)

eq = subs(eq, v, R*theta_dot)

\dfrac{R^2 \,{\dot{\theta} }^2 }{2}=R\,g\,\sin \left(\theta \right)

theta_dot = solve(eq, theta_dot);
theta_dot = theta_dot(1)

\dfrac{\sqrt{2}\,\sqrt{g}\,\sqrt{\sin \left(\theta \right)} }{\sqrt{R} }

N = subs(N, 'theta_dot', theta_dot)

3\,g\,m\,\sin \left(\theta \right)

Now, the geometric constraint means that the box needs to have the same speed as the conveyor belt, thus $\omega r=\dot{\theta} |_{\theta =\pi /2} R$ must hold! From this we can solve for $\omega$

syms omega
theta_dot_2 = subs(theta_dot,theta,pi/2)

\dfrac{\sqrt{2}\,\sqrt{g} }{\sqrt{R} }

omega = solve(omega*r == theta_dot_2*R, omega )

\dfrac{\sqrt{2}\,\sqrt{R}\,\sqrt{g} }{r}

Example 6 - The Damper

A box with a mass of $m$ is gliding without friction with a speed of $v_0$ when it hits a damper with a damping coefficient of $c$ . Determine the position and speed as a function of time.

Solution: We set the origin at the moment of impact. The only force acting of the box in the horizontal direction is the drag force for low velocity flow: $F_d =-c\dot{x}$ . For high velocity the drag is proportional to the velocity squared. The distinction between low and high speed is given by the Reynolds number. Thus with Newton we have $m\ddot{x} =F$ and the IVP becomes

\left\lbrace \begin{array}{ll} m\ddot{x} =-c\dot{x} & \\ \dot{x} \left(0\right)=v_0 & \\ x\left(0\right)=0 & \end{array}\right.

We use dsolve to get the position as a function of time

clear
syms x(t) v0 m c

x_dot = diff(x,t);
DE = diff(x,2,t)*m == -c*x_dot;
BV = [x(0)==0; x_dot(0)==v0]

BV =

\left(\begin{array}{c} x\left(0\right)=0\\ {\left({ {\left(\dfrac{\partial }{\partial t}\;x\left(t\right)\right)}\left|\right.}_{t=0} \right)}=v_0 \end{array}\right)

x = dsolve(DE,BV)

\dfrac{m\,v_0 }{c}-\dfrac{m\,v_0 \,{e }^{-\dfrac{c\,t}{m} } }{c}

x_dot = diff(x,t)

v_0 \,{e}^{-\dfrac{c\,t}{m} }

m = 1; c = 1; v0 = 1;

figure; hold on;
fplot(subs(x),[0,5],'LineWidth',2,'DisplayName','$x(t)$')
fplot(subs(x_dot),[0,5],'LineWidth',2,'DisplayName',"$x'(t)$")
legend('show','interpreter','latex')
xlabel('$t$','interpreter','latex');

It is many times useful to express the velocity as a function of the position $\dot{x} \left(x\right)$ this can be done using the phase-plane. We can plot the parametric plot $\left\lbrack \dot{x} \left(t\right),x\left(t\right)\right\rbrack$ for a range, e.g., $t\in \left\lbrack 0,5\right\rbrack$

figure
fplot(subs(x_dot),subs(x),[0,5],'LineWidth',2,'DisplayName','$\dot{x}(x)$')
legend('show','interpreter','latex')
xlabel('$x(t)$','interpreter','latex'); 
ylabel('$\dot x(t)$','interpreter','latex')
title('Phase plane')

The alternative would be to first solve for $t$ in $x\left(t\right)$ to get $t\left(x\right)$ and plug this into $\dot{x} \left(t\left(x\right)\right)=\dot{x} \left(x\right)$ which of course is more cumbersome, but we get an expression which can be analyzed.

t = solve(x==sym('x'), t)

-\dfrac{m\,\log \left(-\dfrac{c\,x-m\,v_0 }{m\,v_0 }\right)}{c}

v_x = simplify(subs(x_dot))

{ {\left(-\dfrac{c\,x-m\,v_0 }{m\,v_0 }\right)} }^{m/c}

substitution gives

m = 1; c = 1; v0 = 1; syms x
subs(v_x)

1-x

Then we plot

figure
fplot(subs(v_x),[0,1],'LineWidth',2,'DisplayName','$\dot{x}(x)$')
legend('show','interpreter','latex')
xlabel('$x(t)$','interpreter','latex'); 
ylabel('$\dot x(x)$','interpreter','latex')

Example 7 - Spring - Damper

A spring and damper system is an extremely common model in dynamics.

Schematically the system can be drawn like this

Free body diagram gives

\sum F=-\textrm{kx}-c\dot{x} +F_{\textrm{ext} } =m\ddot{x}

where $F_{\textrm{ext}}$ is an externally added force, $k$ is the spring coefficient and $c$ is the damping coefficient. The equation of motion is usually written on the form

m\ddot{x} +c\dot{x} +\textrm{kx}=F_{\textrm{ext} }

For the following we shall assume zero external force for brevity. We have two options for initial values, either a nudge is given to the system in terms of an initial non-zero displacement

\left\lbrace \begin{array}{ll} m\ddot{x} +c\dot{x} +\textrm{kx}=0 & \\ x\left(0\right)=\delta & \\ \dot{x} \left(0\right)=0 & \end{array}\right.

clear
syms x(t) c k m delta
v = diff(x,t);

DE = m*diff(x,2,t) + c*v + k*x == 0;
BV = [x(0)==delta; v(0)==0];
x = dsolve(DE,BV)

x =

\begin{array}{l} \dfrac{\delta \,{e }^{-\dfrac{t\,{\left(c-\sigma_1 \right)} }{2\,m} } \,{\left(c+\sigma_1 \right)} }{2\,\sigma_1 }-\dfrac{\delta \,{e }^{-\dfrac{t\,{\left(c+\sigma_1 \right)} }{2\,m} } \,{\left(c-\sigma_1 \right)} }{2\,\sigma_1 }\\ \mathrm{}\\ \textrm{where}\\ \mathrm{}\\ \;\;\sigma_1 =\sqrt{c^2 -4\,k\,m} \end{array}

Another option is to use an initial velocity

\left\lbrace \begin{array}{ll} m\ddot{x} +c\dot{x} +\textrm{kx}=0 & \\ x\left(0\right)=0 & \\ \dot{x} \left(0\right)=v_0 & \end{array}\right.

syms x(t) v0 

BV = [x(0)==0; v(0)==v0];
x = dsolve(DE,BV)

x =

\dfrac{m\,v_0 \,{\mathrm{e} }^{-\dfrac{t\,{\left(c-\sqrt{c^2 -4\,k\,m}\right)} }{2\,m} } }{\sqrt{c^2 -4\,k\,m} }-\dfrac{m\,v_0 \,{\mathrm{e} }^{-\dfrac{t\,{\left(c+\sqrt{c^2 -4\,k\,m}\right)} }{2\,m} } }{\sqrt{c^2 -4\,k\,m} }

As we see we get a closed form solution for the position of the mass. This model is common enough to have special names and meaning attached to parts of the equation.

The standard form of the model is

\ddot{x} +2\zeta \omega \dot{x} +\omega^2 x=u

where $u=\frac{F_{\textrm{ext}} }{m}$ , $\omega =\sqrt{\frac{k}{m}}$ is known as the natural frequency, $\zeta =\frac{c}{2m\omega }$ is known as the damping ratio (zeta) and has special names associated with certain values, $\zeta =0$ is undamped, $\zeta <1$ is underdamped, $\zeta =1$ is critically damped and $\zeta >1$ is overdamped.

\left\lbrace \begin{array}{ll} x\left(0\right)=\delta & \\ \dot{x} \left(0\right)=0 & \end{array}\right.

syms omega zeta u x(t) delta
v = diff(x,t);

DE = diff(x,2,t) + 2*zeta*omega*v + omega^2*x == 0;
BV = [x(0)==delta; v(0)==0];
x = simplify(dsolve(DE,BV))

x =

\dfrac{\delta \,{\mathrm{e} }^{-\omega \,t\,{\left(\zeta -\sqrt{\zeta^2 -1}\right)} } \,{\left(\zeta +\sqrt{\zeta^2 -1}\right)} }{2\,\sqrt{\zeta^2 -1} }-\dfrac{\delta \,{\mathrm{e} }^{-\omega \,t\,{\left(\zeta +\sqrt{\zeta^2 -1}\right)} } \,{\left(\zeta -\sqrt{\zeta^2 -1}\right)} }{2\,\sqrt{\zeta^2 -1} }

Setting some parameters to numerical values we can plot the oscillation

First undamped:

m = 1; k = 1; c = 0.0; delta = 0.1;
omega = sqrt(k/m)

omega = 1

zeta = c/(2*m*omega)

zeta = 0

syms t
vpa(simplify(subs(x)),3)

0.1\,\cos \left(t\right)

Here we can see that the solution is a cosine function with an constant amplitude of $\delta$ and the frequency of $\omega =1$ .

figure; hold on
fplot(subs(x),[0,10],'LineWidth',2,'Displayname',['$\zeta=',sprintf('%0.2f',zeta),'$'])
xlabel('$t$','interpreter','latex')
ylabel('$x(t)$','interpreter','latex')
legend('show','interpreter','latex')

Then we have underdamped

m = 1; k = 1; c = 0.5; delta = 0.1;
omega = sqrt(k/m)

omega = 1

zeta = c/(2*m*omega)

zeta = 0.2500

vpa(simplify(subs(x)),3)

{\mathrm{e} }^{t\,{\left(-0.25-0.968\,\mathrm{i}\right)} } \,{\left(0.05+0.0129\,\mathrm{i}\right)}+{\mathrm{e} }^{t\,{\left(-0.25+0.968\,\mathrm{i}\right)} } \,{\left(0.05-0.0129\,\mathrm{i}\right)}

fplot(subs(x),[0,10],'LineWidth',2,'Displayname',['$\zeta=',sprintf('%0.2f',zeta),'$'])

Next we have a critically damped system

m = 1; k = 1; 
omega = sqrt(k/m)

omega = 1

c = 2.0000001; % Due to numerical instability
delta = 0.1;
zeta = c/(2*m*omega)

zeta = 1.0000

vpa(simplify(subs(x)),3)

158.0\,{\mathrm{e} }^{-1.0\,t} -158.0\,{\mathrm{e} }^{-1.0\,t}

fplot(subs(x),[0,10],'LineWidth',2,'Displayname',['$\zeta=',sprintf('%0.2f',zeta),'$'])

Finally an overdamped system

m = 1; k = 1; 
omega = sqrt(k/m)

omega = 1

c = 5; 
delta = 0.1;
zeta = c/(2*m*omega)

zeta = 2.5000

vpa(simplify(subs(x)),3)

0.105\,{\mathrm{e} }^{-0.209\,t} -0.00455\,{\mathrm{e} }^{-4.79\,t}

fplot(subs(x),[0,10],'LineWidth',2,'Displayname',['$\zeta=',sprintf('%0.2f',zeta),'$'])

% This will plot the spring
% L = 5; b = 1; N = 20; delta = -1;
% y = linspace(L,0-delta,N+2);
% x = y*0+0;
% x(3:2:end-3) = x(3:2:end-3)-b/2;
% x(4:2:end-3) = x(4:2:end-3)+b/2;
% 
% figure; hold on
% hSpring = plot(x,y,'k-');
% hMass = plot(x(end),y(end),'ko','MarkerFaceColor','k','MarkerSize',10);
% axis equal
% axis([-1,1,-L/2,L])
% 
% 
% 
% syms omega zeta u x(t)
% v = diff(x,t);
% DE = diff(x,2,t) + 2*zeta*omega*v + omega^2*x == 0;
% BV = [x(0)==delta; v(0)==0];
% x = simplify(dsolve(DE,BV))

x =

\dfrac{ {\mathrm{e} }^{-\omega \,t\,{\left(\zeta +\sqrt{\zeta^2 -1}\right)} } \,{\left(\zeta -\sqrt{\zeta^2 -1}\right)} }{2\,\sqrt{\zeta^2 -1} }-\dfrac{ {\mathrm{e} }^{-\omega \,t\,{\left(\zeta -\sqrt{\zeta^2 -1}\right)} } \,{\left(\zeta +\sqrt{\zeta^2 -1}\right)} }{2\,\sqrt{\zeta^2 -1} }

% 
% m = 1; k = 1; c = 5;
% omega = sqrt(k/m)

omega = 1

% zeta = c/(2*m*omega)

zeta = 2.5000

% syms t
% f = matlabFunction(simplify(subs(x)));
% 
% title(['$\omega=1,\; \zeta=',sprintf('%0.2f',zeta),',\;\delta=1$'],'Interpreter','latex')
% xlabel('$t$','Interpreter','latex')
% 
% t1 = 10; dt = 0.1;
% axis([-1,t1+1,-L/2,L])
% hPath = plot(NaN,NaN,'b-','LineWidth',2);
% X = NaN(length(t1/dt)); Y = X;
% i = 1;
% for t = 0:dt:t1
% 
%     y1 = f(t);
%     y = linspace(L,0-y1,N+2);
%     x = y*0+t;
%     x(3:2:end-3) = x(3:2:end-3)-b/2;
%     x(4:2:end-3) = x(4:2:end-3)+b/2;
%     set(hSpring,'YData',y,'XData',x)
%     set(hMass,'Ydata',y(end),'XData',x(end))
%     X(i) = x(end); Y(i) = y(end);
%     set(hPath,'Ydata',Y,'XData',X)
%     drawnow
% 
%     append = true;
%     if i == 1
%         append = false;
%     end
%     exportgraphics(gca,"overdamped.gif","Append",append)
%     i = i + 1;
% 
% end