Table of contents

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

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

where a\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.

{F=Fi=m  r¨TranslationM=Mi+Me=Σri×Fi+Me=Iˉθ¨Rotation \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 Iˉ\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

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

where Rˉ\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 [Fx,Fy,Fz]t=mr¨{\left\lbrack F_x ,F_y ,F_z \right\rbrack }^t =m\ddot{\bm r} and [Mx,My,Mz]T=Iˉθ¨{\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 R0\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.,

F=ddt(m(t)  r˙(t))    F=m˙r˙+mr¨\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.

Step 1 - Free body diagram

Step 2 - Equations of motion

From the figure we get

Fext=P[cosθ,sinθ]TFN=[0,N]TFW=[0,mg]TFf=[μkN,0]T \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

a=[ax,0]T \bm a = [a_x, 0]^T

Using Newtons second law of motion we can formulate

F=maP[cosθsinθ]+[0N]+[0mg]+[μkN0]=[max0] \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
(Pcos(θ)Nμk=axmNgmPsin(θ)=0) \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 axa_x, of course, and NN 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 =
Pμksin(θ)Pcos(θ)+gmμkm -\dfrac{P\,\mu_k \,\sin \left(\theta \right)-P\,\cos \left(\theta \right)+g\,m\,\mu_k }{m}
N =
gm+Psin(θ) g\,m+P\,\sin \left(\theta \right)

From the expression we can see that the vertical component of PP 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 PP during a short time. Let P={4000t0.30t>0.3P=\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 =
{4000  if    t<3100  if    310<t \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[0,0.3]t\in[0,0.3]s and another differential equation with initial values for the range t(0.3,tend]t\in(0.3,t_\text{end}].

The strategy:

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

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

⚠ Note
After tendt_\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

x¨=d2dt2x(t)=ax \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)
ax=P501473500 a_x = \dfrac{P}{50}-\dfrac{1473}{500}

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

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

syms x(t)
DE1 = diff(x,t,2) == a_x
502t2  x(t)=P501473500 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)
t2(10P1473)1000 \dfrac{t^2 \,{\left(10\,P-1473\right)} }{1000}
x1(t) = subs(x1, P, 4000)
38527t21000 \dfrac{38527\,t^2 }{1000}
v1(t) = diff(x1,t)
38527t500 \dfrac{38527\,t}{500}

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

DE2 = vpa(diff(x,t,2) == subs(a_x,P,0) ,2)
2t2  x(t)=2.9 \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)
(x(310)=0.069((t  x(t))t=310)=0.46) \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.5t2+24.0t3.6 -1.5\,t^2 +24.0\,t-3.6
v2(t) = diff(x2,t);
v2 = vpa(v2,2)
24.02.9t 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')]
u1(tt2)t1t2u2(tt1)t1t2 \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)
{80000t3  if    t<320800080000t3  if    t<3100  if    310<t \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))
x1(t)=dt2(800000t31473)1000 x_1(t) = d\frac{t^2 \,{\left(\dfrac{800000\,t}{3}-1473\right)}}{1000}
v1(t) = diff(x1,t)
x1(t)=t(800000t31473)500+800t23 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))
x2(t)=t(80t6)6tt2(800t3785271000) 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)
v2(t)=160t2t(800t3785271000)800t2312 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)
x3(t)=12t1473t21000 x_3(t) = 12\,t-\dfrac{1473\,t^2 }{1000}
v3(t) = diff(x3,t)
v3(t)=121473t500 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 PP. 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.

F=maPN1=m1a1N1N2=m2a2N2=m3a3 \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 =
(PN1=m12t2  x1(t)N1N2=m22t2  x2(t)N2=m32t2  x3(t)) \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 =
t2(N1P)2m1 -\dfrac{t^2 \,{\left(N_1 -P\right)} }{2\,m_1 }
x2 = s.x2
t2(N1N2)2m2 \dfrac{t^2 \,{\left(N_1 -N_2 \right)} }{2\,m_2 }
x3 = s.x3
N2t22m3 \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.

x1=x2t2(N1P)2m1=t2(N1N2)2m2x2=x3t2(N1N2)2m2=N2t22m3 \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 =
P(m2+m3)m1+m2+m3 \dfrac{P\,{\left(m_2 +m_3 \right)} }{m_1 +m_2 +m_3 }
N2 =
Pm3m1+m2+m3 \dfrac{P\,m_3 }{m_1 +m_2 +m_3 }

Example 3 - Constrained motion

A log with a mass m1=200m_1 =200 kg is to be pulled up a 30{30}^{\circ } ramp with a friction coefficient of μk=0.5\mu_k =0\ldotp 5 using a pulley system according to the figure and a counter-mass of m2=125m_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:

m1g[sinθcosθ]FW1+[0N]FN+[2T0]FT+[μkN0]Ff=[m1ac0]mac \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:

FW2T=m2aB 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 =
(2TNμkgm12=aCm1N3gm12=0gm2T=aAm2) \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 CC and AA.

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

syms s_C(t) s_A(t) L
eq = L == 2*s_C + s_A
eq(t) =
L=sA(t)+2sC(t) L=s_A \left(t\right)+2\,s_C \left(t\right)

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

diff(eq,t,2)
0=2t2  sA(t)+22t2  sC(t) 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 aAa_A increases, aCa_C needs to decrease according to the geometric relation in the figure, thus we get:

aA=2aC 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]
(2TNμkgm12=aCm1N3gm12=0gm2T=aAm2aA=2aC) \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 =
3gm1m2μkm14m2 -\frac{\sqrt{3}\,g\,m_1 \,m_2 \,\mu_k }{m_1 -4\,m_2 }
N =
3gm12 \dfrac{\sqrt{3}\,g\,m_1 }{2}
a_C =
gm14gm2+3gm1μk2(m14m2) -\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 =
gm14gm2+3gm1μkm14m2 \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.0730969295836985111236572265625t2 2.0730969295836985111236572265625\,t^2
v_A(t) = diff(s_A,t)
v_A(t) =
4.146193859167397022247314453125t 4.146193859167397022247314453125\,t
t1 = double(max(solve(s_A==6)))
t1 = 1.7012
v_A_end = vpa(v_A(t1),3)
7.05 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

clear
v0 = 24;
theta = 35;
syms theta v0
V0 = v0*[cosd(theta); sind(theta)]
V0 =
(v0cos(πθ180)v0sin(πθ180)) \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) =
(m2t2  x(t)=0m2t2  y(t)=gm) \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]]
(((t  x(t))t=0)=v0cos(πθ180)((t  y(t))t=0)=v0sin(πθ180)x(0)=0y(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) =
tv0cos(πθ180) t\,v_0 \,\cos \left(\dfrac{\pi \,\theta }{180}\right)
y(t) =
gt22+v0sin(πθ180)t+2 -\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])
24tcos(πθ180) 24\,t\,\cos \left(\dfrac{\pi \,\theta }{180}\right)
% yn = subs(y,[g,v0,m],[9.82, 24, 4])
491t2100+24sin(πθ180)t+2 -\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(θ)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(t(θ))x\left(t\left(\theta \right)\right) and plot the graph.

t2 = solve(yn==6, t)
(1200sin(πθ180)49150576sin(πθ180)21964254911200sin(πθ180)491+50576sin(πθ180)2196425491) \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)
1200sin(πθ180)491+50576sin(πθ180)2196425491 \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])
24tcos(πθ180) 24\,t\,\cos \left(\dfrac{\pi \,\theta }{180}\right)
xn = formula(subs(xn, t, t2))
24cos(πθ180)(1200sin(πθ180)491+50576sin(πθ180)2196425491) 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{47}^{\circ }.

Wind resistance

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

mx¨=c  x˙p2+c  vW2my¨=m  gcy˙p2+c  vW2\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 x˙p{\dot{x} }_p and y˙p{\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 vW{\bm v}_W. Here cc 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 vW=[2,0]{\bm v}_W =\left\lbrack -2,0\right\rbrack

clear

syms theta v0
c = 0.03;
V0 = v0*[cosd(theta); sind(theta)]
V0 =
(v0cos(πθ180)v0sin(πθ180)) \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
(m2t2  x(t)=3(t  x(t))21002m2t2  y(t)=3(t  y(t))2100gm) \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]]
(((t  x(t))t=0)=v0cos(πθ180)((t  y(t))t=0)=v0sin(πθ180)x(0)=0y(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

r¨x=ax=ddt(vx)ΔvxΔt=vx,i+1vx,iΔtvx,i+1=vx,i+axΔt{\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

r˙x=vx=ddt(x)ΔxΔt=xi+1xiΔtxi+1=xi+vxΔt{\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

mx¨=c  x˙p2+c  vW2my¨=m  gcy˙p2+c  vW2 \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

ax=c  x˙p2+c  vW2ma_x =\dfrac{-c\;{ {\dot{x} }_p }^2 +c\;{v_W }^2 }{m}

and

ay=m  gcy˙p2+c  vW2ma_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 cc 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 AA and slide down the smooth circular surface of radius RR to a conveyor at BB. Determine the expression for the the normal contact force NN between the guide and each package in terms of θ\theta and specify the angular velocity ω\omega of the conveyor wheel of radius rr 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.

eR=[cosθ(t),sinθ(t)]Te_R ={\left\lbrack -\cos \theta \left(t\right),-\sin \theta \left(t\right)\right\rbrack }^T

We state the normal and tangential directions next, en=eRe_n =-e_R and the tangential we can easily get by et=en×eze_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) =
(sin(θ(t))cos(θ(t))0) \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) =
(Rcos(θ(t))Rsin(θ(t))0) \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) =
(Rθ˙sin(θ(t))Rθ˙cos(θ(t))0) \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) =
(Rsin(θ(t))2t2  θ(t)+Rcos(θ(t))(t  θ(t))2Rsin(θ(t))(t  θ(t))2Rcos(θ(t))2t2  θ(t)0) \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) =
(Rcos(θ(t))θ˙2+Rθ¨sin(θ(t))Rθ˙2sin(θ(t))Rθ¨cos(θ(t))0) \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 0
v_t = simplify(p_dot.' * e_t)
Rθ˙ R\,\dot{\theta}
a_n = simplify(p_ddot.' * e_n)
Rθ˙2 R\,{\dot{\theta} }^2
a_t = simplify(p_ddot.' * e_t)
Rθ¨ R\,\ddot{\theta}

Kinetics - From the free body diagram we get the forces

syms m g N
F_W = m*g*[0;-1;0]
(0gm0) \left(\begin{array}{c} 0\\ -g\,m\\ 0 \end{array}\right)
F_N = N*e_n
(Ncos(θ(t))Nsin(θ(t))0) \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]
(gmsin(θ(t))gmcos(θ(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]
(N0) \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'))
(Ngmsin(θ)=Rmθ˙2gmcos(θ)=Rmθ¨) \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 =
Rmθ˙2+gmsin(θ) R\,m\,{\dot{\theta} }^2 +g\,m\,\sin \left(\theta \right)
theta_ddot =
gcos(θ)R \dfrac{g\,\cos \left(\theta \right)}{R}

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

0vvdv=0θatds\int_0^v v\,dv=\int_0^{\theta } a_t \,ds

turns into:

0vvdv=0sRθ¨ds\int_0^v {vdv}=\int_0^s R\ddot{\theta} ds

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

0vvdv=0θRθ¨Rdθ\int_0^v {vdv}=\int_0^{\theta } R\ddot{\theta} {Rd}\theta

where v=Rθ˙v=R\dot{\theta}

syms v theta theta_dot
eq = int(v,v,0,v) == int(R*theta_ddot*R, theta, 0, theta)
v22=Rgsin(θ) \dfrac{v^2 }{2}=R\,g\,\sin \left(\theta \right)
eq = subs(eq, v, R*theta_dot)
R2θ˙22=Rgsin(θ) \dfrac{R^2 \,{\dot{\theta} }^2 }{2}=R\,g\,\sin \left(\theta \right)
theta_dot = solve(eq, theta_dot);
theta_dot = theta_dot(1)
2gsin(θ)R \dfrac{\sqrt{2}\,\sqrt{g}\,\sqrt{\sin \left(\theta \right)} }{\sqrt{R} }
N = subs(N, 'theta_dot', theta_dot)
3gmsin(θ) 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 ωr=θ˙θ=π/2R\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)
2gR \dfrac{\sqrt{2}\,\sqrt{g} }{\sqrt{R} }
omega = solve(omega*r == theta_dot_2*R, omega )
2Rgr \dfrac{\sqrt{2}\,\sqrt{R}\,\sqrt{g} }{r}

Example 6 - The Damper

A box with a mass of mm is gliding without friction with a speed of v0v_0 when it hits a damper with a damping coefficient of cc. 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: Fd=cx˙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 mx¨=Fm\ddot{x} =F and the IVP becomes

{mx¨=cx˙x˙(0)=v0x(0)=0\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 =
(x(0)=0((t  x(t))t=0)=v0) \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)
mv0cmv0ectmc \dfrac{m\,v_0 }{c}-\dfrac{m\,v_0 \,{e }^{-\dfrac{c\,t}{m} } }{c}
x_dot = diff(x,t)
v0ectm 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 x˙(x)\dot{x} \left(x\right) this can be done using the phase-plane. We can plot the parametric plot [x˙(t),x(t)]\left\lbrack \dot{x} \left(t\right),x\left(t\right)\right\rbrack for a range, e.g., t[0,5]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 tt in x(t)x\left(t\right) to get t(x)t\left(x\right) and plug this into x˙(t(x))=x˙(x)\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)
mlog(cxmv0mv0)c -\dfrac{m\,\log \left(-\dfrac{c\,x-m\,v_0 }{m\,v_0 }\right)}{c}
v_x = simplify(subs(x_dot))
(cxmv0mv0)m/c { {\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)
1x 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

F=kxcx˙+Fext=mx¨\sum F=-\textrm{kx}-c\dot{x} +F_{\textrm{ext} } =m\ddot{x}

where FextF_{\textrm{ext}} is an externally added force, kk is the spring coefficient and cc is the damping coefficient. The equation of motion is usually written on the form

mx¨+cx˙+kx=Fextm\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

{mx¨+cx˙+kx=0x(0)=δx˙(0)=0\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 =
δet(cσ1)2m(c+σ1)2σ1δet(c+σ1)2m(cσ1)2σ1where    σ1=c24km \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

{mx¨+cx˙+kx=0x(0)=0x˙(0)=v0\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 =
mv0et(cc24km)2mc24kmmv0et(c+c24km)2mc24km \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

x¨+2ζωx˙+ω2x=u\ddot{x} +2\zeta \omega \dot{x} +\omega^2 x=u

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

{x(0)=δx˙(0)=0\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 =
δeωt(ζζ21)(ζ+ζ21)2ζ21δeωt(ζ+ζ21)(ζζ21)2ζ21 \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.1cos(t) 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 ω=1\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)
et(0.250.968i)(0.05+0.0129i)+et(0.25+0.968i)(0.050.0129i) {\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.0e1.0t158.0e1.0t 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.105e0.209t0.00455e4.79t 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 =
eωt(ζ+ζ21)(ζζ21)2ζ21eωt(ζζ21)(ζ+ζ21)2ζ21 \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

Example 7 - The pendulum

t:F=mg  sinθ=matt:\sum F=-mg\;\sin \theta ={ma}_t

Thus

at=g  sinθa_t =-g\;\sin \theta

The arc length gives us

s=Lθ,  v=s˙=Lθ˙,at=s¨=Lθ¨s=L\theta ,\;v=\dot{s} =L\dot{\theta} ,a_t =\ddot{s} =L\ddot{\theta}

Thus we get the ODE for the simple pendulum

d2θdt2=gLsinθ\dfrac{d^2 \theta }{ {dt}^2 }=-\dfrac{g}{L}\sin \theta

which is a non-linear ODE with no closed form solution. This needs to be solved using a numerical method for a given initial condition, e.g.,

{θ(0)=θ0θ˙(0)=0\left\lbrace \begin{array}{ll} \theta \left(0\right)=\theta_0 & \\ \dot{\theta} \left(0\right)=0 & \end{array}\right.

For very small angles, θ<<1\theta << 1, we can utilize the small angle approximation, sinθθ\sin \theta \approx \theta and get a modified linear ODE

d2θdt2=gLθfor  θ<<1\dfrac{d^2 \theta }{ {dt}^2 }=-\dfrac{g}{L}\theta \quad\textrm{for}\;\theta <<1

which for the above initial condition has the closed form solution

θ(t)=θ0cos(tgL)for  θ0<<1\theta \left(t\right)=\theta_0 \cos \left(t\sqrt{\dfrac{g}{L} }\right) \quad\textrm{for}\;\theta_0 <<1
Your browser does not support html 5. Cannot display graphics.

Example 8 - The parachute

Example 9 - The bungee cord

Example 10 - The robot arm

Example 11 - The banked track

Example 12 - The roller coaster