Applied Mechanics

by Mirza Cenanovic & Carl-Johan Jonsson

Table of contents

  1. Transient and state based solutions
  2. Example 1 - Linear friction
    1. Free body diagram
  3. Example 2 - Contact force
  4. Example 3 - Constrained motion
  5. Example 4 - Throwing particles through the air
    1. Finding the optimal launch angle
    2. Wind resistance
  6. Example 5
  7. Example 6 - The Damper
  8. Example 7 - Spring - Damper
  9. Example 7 - The pendulum
  10. Example 8 - The parachute
  11. Example 9 - The bungee cord
  12. Example 10 - The robot arm
  13. Example 11 - The banked track
  14. 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

∑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 R→0\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.

Free body diagram

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
(P cos⁡(θ)−N μk=ax mN−g m−P sin⁡(θ)=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) 
[a_x, N] = solve(ekv, [a_x, N])
a_x =
−P μk sin⁡(θ)−P cos⁡(θ)+g m μkm -\dfrac{P\,\mu_k \,\sin \left(\theta \right)-P\,\cos \left(\theta \right)+g\,m\,\mu_k }{m}
N =
g m+P sin⁡(θ) 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.

Let us assume that the box is given a push from PP during a short time. Let P={4000t≤0.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.

Now, we shall use some numerical values to solve the problem.

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

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).

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

a_x = subs(a_x)
P50−1473500 \dfrac{P}{50}-\dfrac{1473}{500} 
syms x(t)
DE1 = m*diff(x,t,2) == a_x
50 ∂2∂t2  x(t)=P50−1473500 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 (10 P−1473)50000 \dfrac{t^2 \,{\left(10\,P-1473\right)} }{50000} 
x1(t) = subs(x1, P, 4000)
38527 t250000 \dfrac{38527\,t^2 }{50000} 
v1(t) = diff(x1,t)
38527 t25000 \dfrac{38527\,t}{25000} 
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) 
DE1 = vpa(m*diff(x,t,2) == subs(a_x,P,0) ,2)
50.0 ∂2∂t2  x(t)=−2.9 50.0\,\dfrac{\partial^2 }{\partial t^2 }\;x\left(t\right)=-2.9 
x2(t) = dsolve(DE1,BV2);
x2 = vpa(x2,2)
−0.029 t2+0.48 t−0.072 -0.029\,t^2 +0.48\,t-0.072 
v2(t) = diff(x2,t);
v2 = vpa(v2,2)
0.48−0.059 t 0.48-0.059\,t

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 = 1.8832

Now we can plot!

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')

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 (t−t2)t1−t2−u2 (t−t1)t1−t2 \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)
{80000 t3  if    t<3208000−80000 t3  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 = m*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)=t2 (800000 t3−1473)50000 x_1(t) = \dfrac{t^2 \,{\left(\dfrac{800000\,t}{3}-1473\right)} }{50000} 
v1(t) = diff(x1,t)
t (800000 t3−1473)25000+16 t23 \dfrac{t\,{\left(\dfrac{800000\,t}{3}-1473\right)} }{25000}+\dfrac{16\,t^2 }{3} 
DE2 = m*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 (8 t5−325)−3 t25−t2 (16 t3−7852750000) x_2(t) = t\,{\left(\dfrac{8\,t}{5}-\dfrac{3}{25}\right)}-\dfrac{3\,t}{25}-t^2 \,{\left(\dfrac{16\,t}{3}-\dfrac{78527}{50000}\right)} 
v2(t) = diff(x2,t)
16 t5−2 t (16 t3−7852750000)−16 t23−625 \dfrac{16\,t}{5}-2\,t\,{\left(\dfrac{16\,t}{3}-\dfrac{78527}{50000}\right)}-\dfrac{16\,t^2 }{3}-\dfrac{6}{25} 
DE3 = m*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)=6 t25−1473 t250000 x_3(t) = \dfrac{6\,t}{25}-\dfrac{1473\,t^2 }{50000} 
v3(t) = diff(x3,t)
625−1473 t25000 \dfrac{6}{25}-\dfrac{1473\,t}{25000} 
t3 = double(solve(v3==0, t))
t3 = 4.0733
x_max = double(x3(t3))
x_max = 0.4888
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.

Free body diagram

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

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 =
(P−N1=m1 ∂2∂t2  x1(t)N1−N2=m2 ∂2∂t2  x2(t)N2=m3 ∂2∂t2  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 (N1−P)2 m1 -\dfrac{t^2 \,{\left(N_1 -P\right)} }{2\,m_1 } 
x2 = s.x2
t2 (N1−N2)2 m2 \dfrac{t^2 \,{\left(N_1 -N_2 \right)} }{2\,m_2 } 
x3 = s.x3
N2 t22 m3 \dfrac{N_2 \,t^2 }{2\,m_3 }

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

[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 =
P m3m1+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.

Free body diagram

Kinematics - we being by modeling the relation between the point CC and BB.

Note that the rope length is constant, let L=2sC+sAL={2s}_C +s_A

clear
syms s_C(t) s_A(t) L 
theta = 30;
eq = L == 2*s_C + s_A
eq(t) =
L=sA(t)+2 sC(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=∂2∂t2  sA(t)+2 ∂2∂t2  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 aCa_C increases, aAa_A needs to decrease according to the second figure, we thus have have 0=−aA+2aC0={-a}_A +2a_C

Kinetics

syms m1 m2 g N mu_k T a_C a_A
F_W1 = m1*g*[-sind(theta); -cosd(theta)]
(−g m12−3 g m12) \def\arraystretch{2.3} \left(\begin{array}{c} -\dfrac{g\,m_1 }{2}\\ -\dfrac{\sqrt{3}\,g\,m_1 }{2} \end{array}\right) 
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 =
(2 T−N μk−g m12=aC m1N−3 g m12=0) \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 \end{array}\right)

For the component AA we have

m2g−T=m2aBm_2 g-T=m_2 a_B 
eq = [eq;
      0 == -2*a_C + a_A
      m2*g - T == m2*a_A]
(2 T−N μk−g m12=aC m1N−3 g m12=00=aA−2 aCg m2−T=aA m2) \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\\ 0=a_A -2\,a_C \\ g\,m_2 -T=a_A \,m_2 \end{array}\right)

Solving for the unknown variables

[T,N,a_C,a_A] = solve(eq, [T,N,a_C,a_A] )
T =
m1 (2 g m2+3 g m2 μk)m1+4 m2 \dfrac{m_1 \,{\left(2\,g\,m_2 +\sqrt{3}\,g\,m_2 \,\mu_k \right)} }{m_1 +4\,m_2 }
N =
3 g m12 \dfrac{\sqrt{3}\,g\,m_1 }{2}
a_C =
−g m1−4 g m2+3 g m1 μk2 (m1+4 m2) -\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 =
−g m1−4 g m2+3 g m1 μkm1+4 m2 -\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 = 1.0041e+03
N = double(subs(N))
N = 1.6991e+03
a_C = double(subs(a_C))
a_C = 0.8885
a_A = double(subs(a_A))
a_A = 1.7769

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) =
0.88847011269535869359970092773438 t2 0.88847011269535869359970092773438\,t^2 
v_A(t) = diff(s_A,t)
v_A(t) =
1.7769402253907173871994018554688 t 1.7769402253907173871994018554688\,t 
t1 = double(max(solve(s_A==6)))
t1 = 2.5987
v_A_end = vpa(v_A(t1),3)
4.62 4.62 
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 =
(v0 cos⁡(π θ180)v0 sin⁡(π θ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) =
(m ∂2∂t2  x(t)=0m ∂2∂t2  y(t)=−g m) \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)=v0 cos⁡(π θ180)((∂∂t  y(t))∣t=0)=v0 sin⁡(π θ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) =
t v0 cos⁡(π θ180) t\,v_0 \,\cos \left(\dfrac{\pi \,\theta }{180}\right)
y(t) =
−g t22+v0 sin⁡(π θ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])
24 t cos⁡(π θ180) 24\,t\,\cos \left(\dfrac{\pi \,\theta }{180}\right) 
% yn = subs(y,[g,v0,m],[9.82, 24, 4])
−491 t2100+24 sin⁡(π θ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)
(1200 sin⁡(π θ180)491−50 576 sin⁡(π θ180)2−1964254911200 sin⁡(π θ180)491+50 576 sin⁡(π θ180)2−196425491) \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)
1200 sin⁡(π θ180)491+50 576 sin⁡(π θ180)2−196425491 \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⁡(π θ180) 24\,t\,\cos \left(\dfrac{\pi \,\theta }{180}\right) 
xn = formula(subs(xn, t, t2))
24 cos⁡(π θ180) (1200 sin⁡(π θ180)491+50 576 sin⁡(π θ180)2−196425491) 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  g−cy˙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 =
(v0 cos⁡(π θ180)v0 sin⁡(π θ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
(m ∂2∂t2  x(t)=−3 (∂∂t  x(t))2100−2m ∂2∂t2  y(t)=−3 (∂∂t  y(t))2100−g m) \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)=v0 cos⁡(π θ180)((∂∂t  y(t))∣t=0)=v0 sin⁡(π θ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+1−vx,iΔt⇒vx,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+1−xiΔt⇒xi+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  g−cy˙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  g−cy˙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) =
(−R cos⁡(θ(t))−R sin⁡(θ(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) =
(R sin⁡(θ(t)) ∂2∂t2  θ(t)+R cos⁡(θ(t)) (∂∂t  θ(t))2R sin⁡(θ(t)) (∂∂t  θ(t))2−R cos⁡(θ(t)) ∂2∂t2  θ(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) =
(R cos⁡(θ(t)) θ˙2+R θ¨ sin⁡(θ(t))R θ˙2 sin⁡(θ(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]
(0−g m0) \left(\begin{array}{c} 0\\ -g\,m\\ 0 \end{array}\right) 
F_N = N*e_n
(N cos⁡(θ(t))N sin⁡(θ(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]
(−g m sin⁡(θ(t))g m cos⁡(θ(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'))
(N−g m sin⁡(θ)=R m θ˙2g m cos⁡(θ)=R m θ¨) \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 θ˙2+g m sin⁡(θ) R\,m\,{\dot{\theta} }^2 +g\,m\,\sin \left(\theta \right)
theta_ddot =
g cos⁡(θ)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

∫0vv dv=∫0θat ds\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=R g sin⁡(θ) \dfrac{v^2 }{2}=R\,g\,\sin \left(\theta \right) 
eq = subs(eq, v, R*theta_dot)
R2 θ˙22=R g sin⁡(θ) \dfrac{R^2 \,{\dot{\theta} }^2 }{2}=R\,g\,\sin \left(\theta \right) 
theta_dot = solve(eq, theta_dot);
theta_dot = theta_dot(1)
2 g sin⁡(θ)R \dfrac{\sqrt{2}\,\sqrt{g}\,\sqrt{\sin \left(\theta \right)} }{\sqrt{R} } 
N = subs(N, 'theta_dot', theta_dot)
3 g m sin⁡(θ) 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)
2 gR \dfrac{\sqrt{2}\,\sqrt{g} }{\sqrt{R} } 
omega = solve(omega*r == theta_dot_2*R, omega )
2 R gr \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)
m v0c−m v0 e−c tmc \dfrac{m\,v_0 }{c}-\dfrac{m\,v_0 \,{e }^{-\dfrac{c\,t}{m} } }{c} 
x_dot = diff(x,t)
v0 e−c tm 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)
−m log⁡(−c x−m v0m v0)c -\dfrac{m\,\log \left(-\dfrac{c\,x-m\,v_0 }{m\,v_0 }\right)}{c} 
v_x = simplify(subs(x_dot))
(−c x−m v0m v0)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)
1−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

∑F=−kx−cx˙+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 =
δ e−t (c−σ1)2 m (c+σ1)2 σ1−δ e−t (c+σ1)2 m (c−σ1)2 σ1where    σ1=c2−4 k m \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 =
m v0 e−t (c−c2−4 k m)2 mc2−4 k m−m v0 e−t (c+c2−4 k m)2 mc2−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} }-\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 (ζ−ζ2−1) (ζ+ζ2−1)2 ζ2−1−δ e−ω t (ζ+ζ2−1) (ζ−ζ2−1)2 ζ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} }-\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⁡(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.25−0.968 i) (0.05+0.0129 i)+et (−0.25+0.968 i) (0.05−0.0129 i) {\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 e−1.0 t−158.0 e−1.0 t 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 e−0.209 t−0.00455 e−4.79 t 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 (ζ+ζ2−1) (ζ−ζ2−1)2 ζ2−1−e−ω t (ζ−ζ2−1) (ζ+ζ2−1)2 ζ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} }-\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

© Mirza Cenanovic. Last modified: April 05, 2024. Website built with Franklin.jl and the Julia programming language.