6.3  Particle Kinematics

What is kinematics?

Kinematics is the study of motion without regard to what causes it. We describe how a body moves, its position, velocity and acceleration, without asking why it moves. The forces and moments responsible for the motion belong to kinetics, which we treat in later chapters. This separation is useful because the geometry of motion obeys its own rules independently of the forces involved: a particle tracing a circular arc has a centripetal acceleration whether it is a planet in orbit or a ball on a string.

We begin with the simplest case: particle kinematics.

A body whose rotation is insignificant compared to the path its center of gravity travels is called a particle. More precisely, if either the rotational acceleration or the physical size of the body is small enough, we can neglect changes in angular momentum and treat the body as a point mass. The equation of motion then reduces to Newton’s second law,

\[ \sum F = m a \]

and the rotational equation of motion \(\sum M = I \alpha\) can be ignored. Later we will generalize to rigid bodies where rotation matters, but for now these two examples illustrate the distinction:

Position vector

The position of a body’s center of gravity at time \(t\) is described by the position vector \(\boldsymbol{r}(t) = [x(t), y(t), z(t)]^\mathsf{T}\), where each component is a function of the same independent variable \(t\). These components are called coordinate functions, and the curve traced out by \(\boldsymbol{r}(t)\) as \(t\) varies is the path of the particle.

Velocity

The mean velocity between two positions \(\boldsymbol{r}(t)\) and \(\boldsymbol{r}(t+\Delta t)\) can be described as \(\bar{\boldsymbol{v}} = \frac{\Delta \boldsymbol{r}}{\Delta t}\). In the limit where \(\Delta t \to 0\) we get the time derivate of the position vector, known as the velocity vector

\[ \lim_{\Delta t \to 0} \frac{\Delta \boldsymbol{r}}{\Delta t} = \frac{d\boldsymbol{r}}{dt} =: \dot{\boldsymbol{r}} = \boldsymbol{v} \tag{6.3.1}\]

As can be seen in the animation, the velocity vector \(\dot{\boldsymbol{r}}\) at \(t\) is always tangential to the path.

The dot notation in \(\dot{r}\) or in \(\dot{\ddot{\ddot{\ddot{\ddot{r}}}}} = \frac{d^9}{dt^9}r\) (a.k.a. Newton’s abomination) comes from Newton who studied dynamic problems and needed to develop the mathematics to deal with these problems.
The mathematics came to be known as calculus which was expanded and is really made up by a great many contributors (Leibnitz, Lagrange). see e.g., this video for an historical overview.

The magnitude of the velocity vector, \(|\dot{\boldsymbol{r}}|\), is known as the speed. The difference between velocity and speed being that velocity also contains information about the direction of the object.

Acceleration

Similarly, the acceleration is the time derivative of the velocity,

\[ \lim_{\Delta t \to 0} \frac{\Delta \dot{\boldsymbol{r}}}{\Delta t} = \frac{d\dot{\boldsymbol{r}}}{dt} = \frac{d^2\boldsymbol{r}}{dt^2} =: \ddot{\boldsymbol{r}} = \boldsymbol{a} \tag{6.3.2}\]

From the animation we can see that the acceleration points inward, toward the center of curvature, and that it generally differs in direction from both the position and velocity vectors. The acceleration vector can be decomposed into a tangential and a normal component, \(\boldsymbol{a} = a_t \boldsymbol{e}_t + a_n \boldsymbol{e}_n\). The normal component \(a_n\) is responsible for changes in direction, while the tangential component \(a_t\) is responsible for changes in speed. If the particle speeds up through a curve, so that \(|\dot{\boldsymbol{r}}_1| > |\dot{\boldsymbol{r}}_0|\), a tangential component is present and the resulting acceleration is no longer perpendicular to the velocity. Only when the speed is constant (\(a_t = 0\)) is the acceleration purely normal, that is, perpendicular to the path.

Rectilinear motion

In the special case of rectilinear motion (motion without curvature), such as described in Figure 6.3.3

the position vector reduces to a single component, for example \(x(t)\), and the velocity and acceleration reduce to the time derivatives of this single component.

The proposed working method is described in what follows.

Differential equations

Classical mechanics textbooks often present ready-made formulas for position, velocity and acceleration under specific assumptions (typically constant acceleration). It is usually more tedious to figure out when these formulas are valid than to set up and solve the differential equation directly. We shall therefore work with \(\boldsymbol{r}, \dot{\boldsymbol{r}}\) and \(\ddot{\boldsymbol{r}}\) throughout and treat kinematics problems as ordinary differential equations (ODEs).

Differentiation and integration act componentwise on vectors,

\[ \dot{\boldsymbol{r}}(t) = [\dot x(t), \dot y(t), \dot z(t)]^\mathsf T \quad \text{and} \quad \int \boldsymbol{r}(t) \, dt = \left[\int x(t) \, dt, \int y(t) \, dt, \int z(t) \, dt\right]^\mathsf T \]

so we can formulate our problems in vector form and solve them using either a symbolic ODE solver such as sympy.dsolve() or a numerical integrator such as the Euler method.

Coordinate systems

Different coordinate systems align their basis vectors with different aspects of the motion. This alignment often simplifies the resulting expressions, but more importantly, it reveals the physical structure of the problem: polar coordinates naturally separate radial and angular motion, natural coordinates separate speed from curvature, and so on. Whichever system we choose, we should commit to it for the duration of the analysis.

We shall here derive the most common coordinate systems and their basis vectors such that this projection can be done. The three systems we consider are:

• Cartesian, our typical xyz-system with fixed basis vectors \(\boldsymbol{e}_x = [1, 0, 0]^\mathsf{T}\), \(\boldsymbol{e}_y = [0, 1, 0]^\mathsf{T}\) and \(\boldsymbol{e}_z = [0, 0, 1]^\mathsf{T}\).
• Natural (tangent, normal, bi-normal), \(\boldsymbol{e}_t, \boldsymbol{e}_n, \boldsymbol{e}_b\), attached to the particle and aligned with its path.
• Polar (2D) or cylindrical (3D), \(\boldsymbol{e}_r, \boldsymbol{e}_\theta, \boldsymbol{e}_z\), aligned with the radial direction from a fixed origin.

The Cartesian basis is constant and needs no further discussion. The natural and polar bases, however, rotate as the particle moves, which means their basis vectors are functions of time. This has important consequences when we differentiate position, velocity and acceleration, as we shall see below.

The time derivative of a unit vector

In both the natural and polar coordinate systems we encountered the time derivative of a unit vector, \(\dot{\boldsymbol{e}}\). We now show that this derivative is always perpendicular to the unit vector itself.

A unit vector satisfies \(\boldsymbol{e} \cdot \boldsymbol{e} = 1\) by definition. Differentiating both sides with respect to time and applying the product rule gives

\[ \frac{d}{dt}(\boldsymbol{e} \cdot \boldsymbol{e}) = \dot{\boldsymbol{e}} \cdot \boldsymbol{e} + \boldsymbol{e} \cdot \dot{\boldsymbol{e}} = 2\,\boldsymbol{e} \cdot \dot{\boldsymbol{e}} = 0 \]

and since \(\boldsymbol{e} \cdot \dot{\boldsymbol{e}} = 0\) we conclude that \(\boldsymbol{e} \perp \dot{\boldsymbol{e}}\). Note that \(|\dot{\boldsymbol{e}}| \neq 1\) in general.

The figure confirms this geometrically: as \(\Delta t \to 0\), the change \(\Delta\boldsymbol{e}\) becomes perpendicular to \(\boldsymbol{e}\).

import sympy as sp
from mechanicskit import la, ltx
t = sp.symbols('t', real=True, positive=True)
theta = sp.Function('theta', real=True)(t)

rr = sp.Matrix([sp.cos(theta), sp.sin(theta)])

ltx("\\bm r =", rr)

ee_r = sp.trigsimp(rr.normalized())

ltx("\\bm e_r =", ee_r)

ee_r_dot = ee_r.diff(t)
ltx(r"\dot{\bm e}_r =", ee_r_dot, arraystretch=2.5)

theta_dot = sp.symbols(r'\dot{\theta}', real=True)
ee_theta = ee_r_dot.normalized().subs(sp.diff(theta, t), theta_dot)
ee_theta = sp.trigsimp(ee_theta)

ltx(r"{\bm e}_\theta =", ee_theta, arraystretch=2.5)

Polar coordinates (\(r-\theta\))

We now derive the kinematic equations in polar coordinates. In practice we can formulate \(\dot{\boldsymbol r}\) and \(\ddot{\boldsymbol r}\) in Cartesian coordinates and project onto any basis, so there is no need to memorize the expressions below. The purpose here is to explore the connections between these quantities and their physical meaning.

We begin by recalling the basis vectors established earlier,

\[ \boldsymbol{e}_r = \frac{\boldsymbol{r}}{|\boldsymbol{r}|}, \qquad \boldsymbol{e}_\theta = \frac{\dot{\boldsymbol{e}}_r}{|\dot{\boldsymbol{e}}_r|} \]

and writing the position vector as

\[ \boldsymbol{r} = r\boldsymbol{e}_r \]

where \(r = |\boldsymbol{r}|\) is the radial distance from the origin.

Natural coordinates (\(t-n\))

We now derive the kinematic equations in natural coordinates. The acceleration is particularly revealing in this basis since it splits cleanly into a component along the path (change of speed) and a component perpendicular to it (change of direction). As before, we can always obtain these by projecting the total acceleration onto any basis, so the expressions below need not be memorized.

We recall the basis vectors established earlier,

\[ \boldsymbol{e}_t = \frac{\dot{\boldsymbol{r}}}{|\dot{\boldsymbol{r}}|}, \qquad \boldsymbol{e}_n = \frac{\dot{\boldsymbol{e}}_t}{|\dot{\boldsymbol{e}}_t|} \]

A tiny view on differential geometry

The theory of curves is an extensive topic within the field of differential geometry. Here, we will only introduce some basic concepts, such as the arc length. For a general curve, the arc length may not have a closed-form solution and is found by integrating the speed \(|\dot{\boldsymbol{r}}(t)|\):

\[ s(t) = \int_0^t |\dot{\boldsymbol{r}}(t)| dt \]

Using the arc length \(s\) as the parameter instead of time \(t\) can sometimes lead to simpler expressions for derivatives. However, since \(s\) is itself a function of \(t\), this change of variable does not necessarily simplify computer-based calculations.

Another important concept is curvature, \(\kappa(t) = \frac{1}{R(t)}\). As expected, the curvature approaches zero as the radius of curvature \(R(t)\) tends to infinity, which corresponds to a straight path. The center of curvature, \(\boldsymbol{r}_C(t)\), is given by:

\[ \boldsymbol{r}_C(t) = \boldsymbol{r}(t) + R(t)\boldsymbol{e}_n(t) \]

This point is the center of the osculating circle (a term coined by Leibniz as Circulus Osculans), which is the circle that best approximates the curve at a given point. Natural coordinates are particularly useful for this type of analysis.

The torsion of a curve, \(\tau(t)\), measures how sharply it twists out of its plane of curvature. A simple analogy is the pitch of a screw. The derivations for both curvature and torsion can be complex, see e.g., the Frenet-Serret formulas, so the formulas are stated here without proof:

\[ \kappa(t) = \frac{|\dot{\boldsymbol{r}}(t) \times \ddot{\boldsymbol{r}}(t)|}{|\dot{\boldsymbol{r}}(t)|^3} \]

\[ \tau(t) = - \frac{(\dot{\boldsymbol{r}}(t) \times \ddot{\boldsymbol{r}}(t)) \cdot \dot{\ddot{\boldsymbol{r}}}(t)}{|\dot{\boldsymbol{r}}(t) \times \ddot{\boldsymbol{r}}(t)|^2} \]

where \(\dot{\ddot{\boldsymbol{r}}}\) denotes the third-order time derivative of the position vector, also known as jerk or jolt. For the names of higher-order derivatives, see this list.

import sympy as sp
from mechanicskit import la, ltx

t = sp.symbols('t', real=True, positive=True)
R = sp.symbols('R', real=True, positive=True)

r = R * sp.Matrix([sp.cos(t), sp.sin(t), 0])

ltx("\\bm r =", r)

r_dot = r.diff(t)
r_ddot = r.diff(t, 2)
r_dddot = r.diff(t, 3)

# kappa = |r_dot x r_ddot| / |r_dot|^3
kappa = sp.simplify((r_dot.cross(r_ddot)).norm() / r_dot.norm()**3)

ltx("\\kappa =", kappa)

# tau = - (r_dot x r_ddot) . r_dddot / |r_dot x r_ddot|^2
# The cross product r_dot.cross(r_ddot) is not zero, so we can proceed.
tau = sp.simplify(- (r_dot.cross(r_ddot)).dot(r_dddot) / (r_dot.cross(r_ddot)).norm()**2)

ltx("\\tau =", tau)

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation, PillowWriter

t = sp.symbols('t', real=True, positive=True)
rr = sp.Matrix([t, sp.sqrt(t) * sp.sin(2*t)])

rr_dot = rr.diff(t)
rr_ddot = rr.diff(t, 2)

kappa = (rr_dot[0]*rr_ddot[1] - rr_dot[1]*rr_ddot[0]) / rr_dot.norm()**3
R_curv = 1 / kappa

ee_t = rr_dot / rr_dot.norm()
ee_n = sp.Matrix([-ee_t[1], ee_t[0]])

rr_C = rr + R_curv * ee_n




r_func = sp.lambdify(t, rr, 'numpy')
r_C_func = sp.lambdify(t, rr_C, 'numpy')
R_func = sp.lambdify(t, R_curv, 'numpy')
e_t_func = sp.lambdify(t, ee_t, 'numpy')

fig, ax = plt.subplots(figsize=(8, 6))
fig.patch.set_facecolor('#212121')
ax.set_facecolor('#212121')
ax.set_aspect('equal')
ax.tick_params(axis='x', colors='white')
ax.tick_params(axis='y', colors='white')
for spine in ax.spines.values():
    spine.set_edgecolor('white')

t_vals = np.linspace(0.01, 10, 400)
path = r_func(t_vals)
ax.plot(path[0].flatten(), path[1].flatten(), 'c-')

particle, = ax.plot([], [], 'ro')
circle_line, = ax.plot([], [], 'g-')
tangent_arrow = ax.quiver(0, 0, 0, 0, color='b', scale=1, angles='xy', scale_units='xy')
radius_arrow = ax.quiver(0, 0, 0, 0, color='r', scale=1, angles='xy', scale_units='xy')
ax.set_xlim(0, 10)
ax.set_ylim(-3, 3)
ax.grid(False)

def update(frame):
    pos = r_func(frame)
    center = r_C_func(frame)
    radius = R_func(frame)
    et = e_t_func(frame)
    particle.set_data(pos[0], pos[1])
    tangent_arrow.set_offsets(pos.flatten())
    tangent_arrow.set_UVC(et[0], et[1])
    rv = center - pos
    rv = rv / np.linalg.norm(rv)
    radius_arrow.set_offsets(pos.flatten())
    radius_arrow.set_UVC(rv[0], rv[1])
    theta = np.linspace(0, 2*np.pi, 100)
    circle_line.set_data(center[0] + radius*np.cos(theta), center[1] + radius*np.sin(theta))
    return particle, circle_line, tangent_arrow, radius_arrow

ani = FuncAnimation(fig, update, frames=t_vals, blit=False, interval=50)
ani.save('graphics/kinematics_osculating_circle.gif', writer=PillowWriter(fps=20))
plt.close(fig)

Circular motion

A special case of curvilinear motion is planar circular motion. We study this using polar coordinates and the rotating arm shown below.

We express the particle’s position as \(\boldsymbol{r} = r(t)[\cos\theta(t), \sin\theta(t)]^\mathsf{T}\).

theta = sp.Function('theta', real=True, positive=True)(t)
r = sp.Function('r', real=True, positive=True)(t)

rr = r * sp.Matrix([sp.cos(theta), sp.sin(theta)])
ltx(r"\bm r =", rr)

rr_dot = rr.diff(t)
ltx(r"\dot{\bm r} =", rr_dot, arraystretch=2.5)

We clean up by substituting Newton notation for the derivatives.

dot_subs = {
    sp.diff(theta, t, 2): sp.symbols(r'\ddot{\theta}'),
    sp.diff(theta, t):    sp.symbols(r'\dot{\theta}'),
    sp.diff(r, t, 2):     sp.symbols(r'\ddot{r}'),
    sp.diff(r, t):        sp.symbols(r'\dot{r}'),
}

rr_dot = rr_dot.subs(dot_subs)
ltx(r"\dot{\bm r} =", rr_dot, arraystretch=1.5)

The acceleration is

rr_ddot = rr.diff(t, 2).subs(dot_subs)
ltx(r"\ddot{\bm r} =", rr_ddot, arraystretch=1.5)

SymPy handles the tedious differentiation for us. We now define the polar basis vectors and project onto them.

ee_r = sp.Matrix([sp.cos(theta), sp.sin(theta)])
ee_theta = sp.Matrix([-sp.sin(theta), sp.cos(theta)])

where \(\boldsymbol{e}_\theta\) is perpendicular to \(\boldsymbol{e}_r\), rotated 90 degrees counter-clockwise. Projecting onto the polar basis using the dot product (valid since the basis is orthonormal),

rr_dot_polar = sp.simplify(sp.Matrix([ee_r.dot(rr_dot), ee_theta.dot(rr_dot)]))

ltx(r"\dot{\bm r}_{r\theta} =", rr_dot_polar, arraystretch=1.5)

rr_ddot_polar = sp.simplify(sp.Matrix([ee_r.dot(rr_ddot), ee_theta.dot(rr_ddot)]))

ltx(r"\ddot{\bm r}_{r\theta} =", rr_ddot_polar, arraystretch=1.5)

We recover the same expressions derived by hand in Section 6.3.9.

circ_subs = {sp.symbols(r'\dot{r}'): 0, sp.symbols(r'\ddot{r}'): 0}

ltx(r"\dot{\bm r}_{r\theta} =", rr_dot_polar.subs(circ_subs), arraystretch=1.5)

ltx(r"\ddot{\bm r}_{r\theta}  =", rr_ddot_polar.subs(circ_subs), arraystretch=1.5)

Relative motion

It is often easier to express position, velocity and acceleration relative to some other point rather than the origin. This simplifies the modeling and reflects how engineers typically work: relative measurements are often more practical and more accurate than absolute ones. This approach, called relative motion analysis, will be used extensively in rigid body dynamics.

If we know the position of particle \(B\) and the position of particle \(A\) relative to \(B\), written \(\boldsymbol{r}_{A/B}\) (read: “\(A\) seen from \(B\)”), the absolute position of \(A\) is

\[ \boldsymbol{r}_A = \boldsymbol{r}_B + \boldsymbol{r}_{A/B} \]

It may be more convenient to express \(\boldsymbol{r}_{A/B}\) and \(\boldsymbol{r}_B\) separately than to formulate \(\boldsymbol{r}_A\) directly. Differentiating with respect to time gives

\[ \dot{\boldsymbol{r}}_A = \dot{\boldsymbol{r}}_B + \dot{\boldsymbol{r}}_{A/B} \quad \text{or} \quad \boldsymbol{v}_A = \boldsymbol{v}_B + \boldsymbol{v}_{A/B} \]

and

\[ \ddot{\boldsymbol{r}}_A = \ddot{\boldsymbol{r}}_B + \ddot{\boldsymbol{r}}_{A/B} \quad \text{or} \quad \boldsymbol{a}_A = \boldsymbol{a}_B + \boldsymbol{a}_{A/B} \]

Example 1: Rectilinear motion

The position of a particle moving along a straight line is given by \(s = 2t^3 - 24t + 6\), where \(s\) is in meters and \(t\) in seconds. Determine

1. the time for the particle to reach a velocity of 72 m/s from its initial condition at \(t = 0\),

2. the acceleration when \(v = 30\) m/s, and

3. the net displacement from \(t = 1\) to \(t = 4\).

import sympy as sp
from mechanicskit import la, ltx
t = sp.symbols('t', real=True, positive=True)
s = 2*t**3 - 24*t + 6
v = sp.diff(s, t)
a = sp.diff(v, t)

ltx(r"s(t) =", s, r"\\v(t) =", v, r"\\ a(t) =", a)

(a) Find \(t\) when \(v = 72\) m/s.

t1 = sp.solve(sp.Eq(v, 72), t)
ltx(r"v(t) = 72 \Rightarrow t =", t1)

(b) Find the acceleration when \(v = 30\) m/s.

t2 = sp.solve(sp.Eq(v, 30), t)
a30 = a.subs(t, t2[0])
ltx(r"v(t) = 30 \Rightarrow t =", t2[0],
    r"\\ a(t_2) =", a30, r"\text{ m/s}^2")

(c) Net displacement from \(t = 1\) to \(t = 4\).

s4 = s.subs(t, 4)
s1 = s.subs(t, 1)
ds = s4 - s1

ltx(r"s(4) - s(1) =", s4, r"-", s1, r"=", ds, r"\text{ m}")

Example 2: Relative motion

Passengers in an A380 flying east at 800 km/h observe a JAS 39 Gripen passing underneath with a true heading of \(45^\circ\) (north-east). To the passengers, the Gripen appears to move away at \(60^\circ\) as shown. Determine the true speed of the Gripen.

Solution

We have two unknowns: the absolute speed of the Gripen, \(v_B\), and the relative speed observed from the A380, \(v_{B/A}\). The relative motion equation is

\[ \boldsymbol{v}_B = \boldsymbol{v}_A + \boldsymbol{v}_{B/A} \]

Using the known directions,

\[ v_{B}\begin{bmatrix}\cos 45^\circ\\ \sin 45^\circ\end{bmatrix} = \begin{bmatrix}800\\ 0\end{bmatrix} + v_{B/A}\begin{bmatrix}-\cos 60^\circ\\ \sin 60^\circ\end{bmatrix} \]

import sympy as sp
from mechanicskit import la, ltx
v_B, v_BA = sp.symbols('v_B v_{B/A}', real=True, positive=True)

vv_B = v_B * sp.Matrix([sp.cos(sp.rad(45)), sp.sin(sp.rad(45))])
vv_A = sp.Matrix([800, 0])
vv_BA = v_BA * sp.Matrix([-sp.cos(sp.rad(60)), sp.sin(sp.rad(60))])

eq = sp.Eq(vv_B, vv_A + vv_BA)
eq | la(arraystretch=2.5)

sol = sp.solve(eq, (v_B, v_BA))
sol[v_B] | la.evalf(3)

sol[v_BA].evalf(3) | la

The Gripen’s true speed is 717 km/h, and it moves away from the passengers at 586 km/h.

Example 3: Data fitting

The acceleration performance of a new solar car is being tested at the Jönköping airfield. The JTH students are measuring time against velocity and the data is presented as two vectors. Create a suitable continuous function of the speed and determine how far the car has traveled at the last measurement point. Plot the distance, time and acceleration curves.

Given data:

T=[1.1486, 2.0569, 3.1485, 4.0443, 5.0165, 6.2552, 7.1682, 8.2789, 9.209, 10.175]

V=[9.1023, 16.245, 22.732, 26.446, 33.052, 36.889, 41.304, 43.496, 45.861, 48.319]

Solution

We fit a quadratic polynomial to the data. In ascending powers of \(t\),

\[ v_i = c_0 + c_1 t_i + c_2 t_i^2 \quad \forall\, i \in [1, m] \]

where \(m = 10\) in our case. On matrix form this becomes

\[ \underbrace{\left[\begin{array}{ccc} 1 & t_1 & t_1^2 \\ 1 & t_2 & t_2^2 \\ \vdots & \vdots & \vdots \\ 1 & t_m & t_m^2 \end{array}\right]}_{\bm{X}} \underbrace{\left[\begin{array}{c} c_0 \\ c_1 \\ c_2 \end{array}\right]}_{\bm{c}} = \underbrace{\left[\begin{array}{c} v_1 \\ v_2 \\ \vdots \\ v_m \end{array}\right]}_{\bm{v}} \]

The matrix \(\bm{X}\) is called a Vandermonde matrix. Note that the columns are ordered by ascending powers of \(t\): the first column is \(t^0 = 1\), the second is \(t^1\), the third is \(t^2\). This means \(c_0\) is the constant term, \(c_1\) the linear coefficient, and \(c_2\) the quadratic coefficient, that is, c[i] multiplies \(t^i\).

Since we have more equations (\(m = 10\)) than unknowns (\(n = 3\)), the system is overdetermined and has no exact solution. Instead we seek the \(\bm{c}\) that minimizes the squared error,

\[ E(\bm{c}) = ||\bm{X}\bm{c} - \bm{v}||^2 = (\bm{X}\bm{c} - \bm{v})^\mathsf{T}(\bm{X}\bm{c} - \bm{v}) \]

Expanding and differentiating with respect to \(\bm{c}\) gives

\[ \frac{\partial E}{\partial \bm{c}} = 2\bm{X}^\mathsf{T}\bm{X}\bm{c} - 2\bm{X}^\mathsf{T}\bm{v} = \bm{0} \]

which yields the normal equations,

\[ \boxed{\bm{c} = (\bm{X}^\mathsf{T} \bm{X})^{-1} \bm{X}^\mathsf{T} \bm{v}} \]

import numpy as np
from mechanicskit import la, ltx


# Given data:
T = np.array([1.1486, 2.0569, 3.1485, 4.0443, 5.0165, 6.2552, 7.1682, 8.2789, 9.209, 10.175])
V = np.array([9.1023, 16.245, 22.732, 26.446, 33.052, 36.889, 41.304, 43.496, 45.861, 48.319])

# np.vander builds the Vandermonde matrix
# increasing=True gives columns [1, t, t^2], matching c_0 + c_1*t + c_2*t^2
X = np.vander(T, 3, increasing=True)

ltx(r"\bm X =", X, show_shape=True)

The following python code implements \(\mathbf{c}=\left(\mathbf{X}^T \mathbf{X}\right)^{-1} \mathbf{X}^T \mathbf{v}\) literally

c = np.linalg.solve(X.T @ X, X.T @ V)
ltx(r"\bm c =", c, arraystretch=1.0, show_shape=True)

But doing this manually is tedious and error prone, instead we should use the least squares solver in numpy, numpy.linalg.lstsq

c = np.linalg.lstsq(X, V)[0]

c | la(show_shape=True)

import sympy as sp
t = sp.Symbol('t')

# Since X has columns [1, t, t^2], c[i] multiplies t^i
v = c[0] + c[1]*t + c[2]*t**2

ltx(r"v(t) =", v, precision=3)

We now have a continuous function \(v(t)\) that we can differentiate and integrate symbolically.

import matplotlib.pyplot as plt

v_func = sp.lambdify(t, v, 'numpy')
tn = np.linspace(0, max(T), 100)
vn = v_func(tn)

fig, ax = plt.subplots()
ax.plot(T, V, 'ro', label='Data')
ax.plot(tn, vn, 'b-', label='Fit')
ax.set(xlabel='Time, t [s]', ylabel='Velocity, v [m/s]', title='Velocity vs. Time')
ax.grid(True)
ax.legend()
ax.fill_between(tn, vn, color='b', alpha=0.2)
plt.show()

# Acceleration is the derivative of velocity
a = sp.diff(v, t)
ltx(r"a(t) =", a, precision=3)

# Plot acceleration
an = sp.lambdify(t, a, 'numpy')(tn)
fig, ax = plt.subplots()
ax.plot(tn, an, 'm-')
ax.set(xlabel='Time, t [s]', ylabel='Acceleration, a [m/s^2]', title='Acceleration vs. Time')
ax.grid(True)
plt.show()

# Total distance is the integral of velocity
s_tot = sp.integrate(v, (t, 0, max(T)))

ltx(r"s_{tot} = \int_0^{10.175} v(t)\, dt =", s_tot, precision=4)

s = sp.integrate(v)

ltx(r"s(t) =", s, precision=3)

sn = sp.lambdify(t, s, 'numpy')(tn)
fig, ax = plt.subplots()
ax.plot(tn, sn, 'k-')
ax.set(xlabel='Time, t [s]', ylabel='Distance, s [m]', title='Distance vs. Time')
ax.grid(True)
plt.show()

Example 4: Oscillating slider

A spring-mounted slider moves along the \(x\)-axis with initial velocity \(v_0 = 2\) m/s as it crosses the mid-position where \(s = 0\) at \(t = 0\). The two springs exert a retarding force proportional to the displacement, giving an acceleration \(a = -k^2 s\) with \(k = 4\). Determine the displacement, velocity and acceleration of the slider.

k, v0 = sp.symbols('k v0', positive=True)
x = sp.Function('x')(t)

DE = sp.Eq(x.diff(t, 2), -k**2 * x)
ICs = {x.subs(t, 0): 0, x.diff(t).subs(t, 0): v0}

DE | la

x_sol = sp.simplify(sp.dsolve(DE, ics=ICs).rhs)
ltx(r"x(t) =", x_sol)

v_sol = x_sol.diff(t)
a_sol = v_sol.diff(t)

ltx(r"v(t) =", v_sol)

ltx(r"a(t) =", a_sol)

data = {k: 4, v0: 2}
x_n = x_sol.subs(data)
v_n = v_sol.subs(data)
a_n = a_sol.subs(data)

ltx(r"x(t) =", x_n, r",\quad v(t) =", v_n, r",\quad a(t) =", a_n)

t1 = 2 * np.pi
t_vals = np.linspace(0, t1, 400)

x_func = sp.lambdify(t, x_n, 'numpy')
v_func = sp.lambdify(t, v_n, 'numpy')
a_func = sp.lambdify(t, a_n, 'numpy')

fig, ax = plt.subplots(figsize=(10, 6))
ax.plot(t_vals, v_func(t_vals), linewidth=2, label='Velocity [m/s]')
ax.fill_between(t_vals, v_func(t_vals), alpha=0.2, label='Position [m]')
ax.spines['bottom'].set_position('zero')
ax.spines['right'].set_color('none')
ax.spines['top'].set_color('none')
ax.set_xlabel(r'Time $t$ [s]', loc='right')
ax.legend()
ax.grid(True)
plt.show()

The total area under the velocity curve is the distance the slider has traveled, which for \(t = 2\pi\) is 2 m. But the signed area is zero, since the slider returns to its starting position after \(2\pi\) seconds, \(x(2\pi) = 0\) m.

signed_area = sp.integrate(v_n, (t, 0, 2*sp.pi))
total_distance = sp.integrate(sp.Abs(v_n), (t, 0, 2*sp.pi))

ltx(r"\int_0^{2\pi} v\, dt =", signed_area,
    r",\qquad \int_0^{2\pi} |v|\, dt =", total_distance)

Finally, plotting \(v(t)\) against \(x(t)\) gives the phase portrait: an ellipse, characteristic of simple harmonic motion.

fig, ax = plt.subplots()
ax.plot(x_func(t_vals), v_func(t_vals), linewidth=2)
ax.spines['bottom'].set_position('zero')
ax.spines['left'].set_position('zero')
ax.spines['right'].set_color('none')
ax.spines['top'].set_color('none')
ax.set_xlabel(r'$x(t)$ [m]', loc='right')
ax.set_ylabel(r'$v(t)$ [m/s]', loc='top')
# ax.set_aspect('equal')
ax.set_xlim(-1, 1)
ax.set_ylim(-2.1, 2.1)
ax.grid(True)
plt.show()

Example 5: Piecewise functions

We have noticed the position of a particular nice motorcycle,

\[ x(t) = \begin{cases} t^2 & 0 \le t < 10 \\ 20t - 100 & t \ge 10 \end{cases} \]

Plot its velocity and acceleration curves for \(t \in [0, 30]\) s.

Solution

We define the piecewise position function and differentiate to obtain velocity and acceleration. Note that at \(t = 10\) the two pieces meet: \(10^2 = 100 = 20 \cdot 10 - 100\), so the position is continuous. The velocity transitions from \(2t = 20\) m/s to a constant 20 m/s, which is also continuous. The acceleration, however, jumps from 2 m/s\(^2\) to 0.

x = sp.Piecewise((t**2, t < 10), (20*t - 100, True))
ltx(r"x(t) =", x)

x_dot = sp.diff(x, t)
x_ddot = sp.diff(x, t, 2)

ltx(r"\dot x(t) =", x_dot, r",\quad\ddot x(t) =", x_ddot)

from mechanicskit import fplot

fig, axes = plt.subplots(3, 1, figsize=(8, 8), sharex=True)

fplot(x, (t, 0, 30), ax=axes[0], label=r'$x(t)$')
axes[0].set_ylabel('Position [m]')
axes[0].legend()
axes[0].grid(True)

fplot(x_dot, (t, 0, 30), ax=axes[1], label=r'$\dot x(t)$')
axes[1].set_ylabel('Velocity [m/s]')
axes[1].legend()
axes[1].grid(True)

fplot(x_ddot, (t, 0, 30), ax=axes[2], label=r'$\ddot x(t)$')
axes[2].set_ylabel(r'Acceleration [m/s$^2$]')
axes[2].set_xlabel('Time $t$ [s]')
axes[2].legend()
axes[2].grid(True)

fig.tight_layout()
plt.show()

Example 6: Solving an ODE using Euler’s method

A sports car manufacturer is limiting the power for one of the models by limiting the fuel consumption when the throttle is fully engaged such that the acceleration at any moment is proportional with the constant \(k = 0.1\) s\(^{-1}\) with respect to the difference in wanted top speed \(v_1 = 80\) m/s and the current speed \(v\). Assume the car starts from rest and the driver puts the pedal to the metal.

• Draw the distance and velocity curves
• How long does it take to reach top speed and how far has the car traveled?
• What is the 0-100 km/h time and distance?

Solution

From the problem description we formulate the ODE with initial conditions,

\[ \ddot{x}(t) = k(v_1 - \dot{x}(t)), \quad \dot{x}(0) = 0, \quad x(0) = 0 \]

x = sp.Function('x')(t)
v1, k = sp.Rational(80), sp.Rational(1, 10)

DE = sp.Eq(x.diff(t, 2), k * (v1 - x.diff(t)))
ICs = {x.subs(t, 0): 0, x.diff(t).subs(t, 0): 0}

x_sol = sp.dsolve(DE, ics=ICs).rhs
v_sol = x_sol.diff(t)

ltx(r"x(t) =", x_sol, r",\quad v(t) =", v_sol)

from mechanicskit import fplot

fig, ax = plt.subplots()
fplot(v_sol * sp.Rational(36, 10), (t, 0, 70), ax=ax, label=r'$v(t)$')
ax.axhline(100, color='b', linestyle='--', label='100 km/h')
ax.set(xlabel='Time [s]', ylabel='Velocity [km/h]')
ax.set_xlim(0, 70)
ax.set_ylim(0, 300)
ax.legend()
ax.grid(True)
plt.show()

t_99 = sp.solve(sp.Eq(v_sol, sp.Rational(99, 100) * v1), t)[0]
x_99 = x_sol.subs(t, t_99)

ltx(r"t_{99\%} =", t_99, r"\text{ s} \approx", t_99.evalf(4), r"\text{ s}")
ltx(r"x(t_{99\%}) =", x_99.evalf(4), r"\text{ m}")

# 100 km/h = 100/3.6 m/s
v_100kmh = sp.Rational(1000, 36)
t_100 = sp.solve(sp.Eq(v_sol, v_100kmh), t)[0]
x_100 = x_sol.subs(t, t_100)

ltx(r"t_{0-100} =", t_100.evalf(4), r"\text{ s}")
ltx(r"x_{0-100} =", x_100.evalf(4), r"\text{ m}")

dt = 0.1
t_end = 70
n = int(t_end / dt) + 1

t_euler = np.zeros(n)
x_euler = np.zeros(n)
v_euler = np.zeros(n)

k_val, v1_val = 0.1, 80

for i in range(n - 1):
    v_euler[i+1] = v_euler[i] + k_val * (v1_val - v_euler[i]) * dt
    x_euler[i+1] = x_euler[i] + v_euler[i+1] * dt
    t_euler[i+1] = t_euler[i] + dt

x_exact = sp.lambdify(t, x_sol, 'numpy')
v_exact = sp.lambdify(t, v_sol, 'numpy')

fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 8), sharex=True)

ax1.plot(t_euler, v_euler, 'r--', lw=3, label='Euler')
ax1.plot(t_euler, v_exact(t_euler), 'b-', label='Exact')
ax1.set_ylabel('Velocity [m/s]')
ax1.legend()
ax1.grid(True)

ax2.plot(t_euler, x_euler, 'r--', lw=3, label='Euler')
ax2.plot(t_euler, x_exact(t_euler), 'b-', label='Exact')
ax2.set_ylabel('Position [m]')
ax2.set_xlabel('Time [s]')
ax2.legend()
ax2.grid(True)

fig.tight_layout()
plt.show()

Example 7: Curvilinear motion in polar coordinates

A particle \(B\) is confined to the radial groove of a rotating arm. The angular position and radial distance are given by

\[ \theta(t) = 0.2t + 0.02t^3, \qquad r(t) = 0.2 + 0.04t^2 \]

Analyse the motion of the particle in Cartesian and polar coordinates for \(t \in [0, 5]\) s.

Solution

We express the position in Cartesian form, differentiate, then project onto the polar basis to obtain the radial and angular components.

theta = sp.Rational(1, 5)*t + sp.Rational(1, 50)*t**3
r = sp.Rational(1, 5) + sp.Rational(1, 25)*t**2

rr = r * sp.Matrix([sp.cos(theta), sp.sin(theta)])
ltx(r"\theta(t) =", theta, r",\qquad r(t) =", r)

rr_dot = rr.diff(t)
rr_ddot = rr.diff(t, 2)

We define the polar basis and project the Cartesian velocity and acceleration onto it.

ee_r = sp.Matrix([sp.cos(theta), sp.sin(theta)])
ee_theta = sp.Matrix([-sp.sin(theta), sp.cos(theta)])

v_polar = sp.trigsimp(sp.Matrix([ee_r.dot(rr_dot), ee_theta.dot(rr_dot)]))
a_polar = sp.trigsimp(sp.Matrix([ee_r.dot(rr_ddot), ee_theta.dot(rr_ddot)]))

ltx(r"\dot{\bm r}_{r\theta} =", v_polar, r",\quad \ddot{\bm r}_{r\theta} =", a_polar, arraystretch=2.0)

We can verify: \(\dot{r} = 0.08t\) and \(r\dot{\theta} = (0.2 + 0.04t^2)(0.2 + 0.06t^2)\), consistent with the general polar velocity formula \(\dot{\boldsymbol{r}} = \dot{r}\boldsymbol{e}_r + r\dot{\theta}\boldsymbol{e}_\theta\).

rr_func = sp.lambdify(t, rr, 'numpy')
rr_dot_func = sp.lambdify(t, rr_dot, 'numpy')
rr_ddot_func = sp.lambdify(t, rr_ddot, 'numpy')

t_vals = np.linspace(0, 5, 300)
path = np.array([rr_func(ti).flatten() for ti in t_vals])

fig, ax = plt.subplots(figsize=(8, 6))
ax.plot(path[:, 0], path[:, 1], 'b-', linewidth=2)
ax.set_aspect('equal')
ax.set_xlabel(r'$x(t)$ [m]')
ax.set_ylabel(r'$y(t)$ [m]')
ax.set_title('Path of the particle')
ax.grid(True)

# Mark integer time points
for ti in [1, 2, 3,4,5]:
    pos = rr_func(ti).flatten()
    ax.plot(pos[0], pos[1], 'ko', markersize=5)
    ax.annotate(f'$t={ti}$', pos, textcoords='offset points',
               xytext=(8, 8), fontsize=11,
               bbox=dict(boxstyle='round,pad=0.2', fc='white', alpha=0.8))

plt.show()

from matplotlib.animation import FuncAnimation, FFMpegWriter

# Precompute component curves
t_plot = np.linspace(0.1, 5, 300)
vx = np.array([rr_dot_func(ti).flatten()[0] for ti in t_plot])
vy = np.array([rr_dot_func(ti).flatten()[1] for ti in t_plot])
ax_vals = np.array([rr_ddot_func(ti).flatten()[0] for ti in t_plot])
ay = np.array([rr_ddot_func(ti).flatten()[1] for ti in t_plot])

fig = plt.figure(figsize=(14, 6))
ax_path = fig.add_subplot(121)
ax_comp = fig.add_subplot(122)

# Left: path + arrows
ax_path.plot(path[:, 0], path[:, 1], 'b-', linewidth=1, alpha=0.4)
ax_path.set_aspect('equal')
ax_path.set_xlabel(r'$x$ [m]')
ax_path.set_ylabel(r'$y$ [m]')
ax_path.grid(True)
ax_path.set_xlim(-1.5, 0.7)
ax_path.set_ylim(-0.5, 1.3)
ax_path.set_title('Path')

dot, = ax_path.plot([], [], 'ko', markersize=6)
v_arrow = ax_path.annotate('', xy=(0, 0), xytext=(0, 0),
    arrowprops=dict(arrowstyle='->', color='blue', lw=2))
a_arrow = ax_path.annotate('', xy=(0, 0), xytext=(0, 0),
    arrowprops=dict(arrowstyle='->', color='red', lw=2))

# Right: component curves (faint background)
ax_comp.plot(t_plot, vx, 'b-', alpha=0.3, label=r'$\dot x$')
ax_comp.plot(t_plot, vy, 'r-', alpha=0.3, label=r'$\dot y$')
ax_comp.plot(t_plot, ax_vals, 'b--', alpha=0.3, label=r'$\ddot x$')
ax_comp.plot(t_plot, ay, 'r--', alpha=0.3, label=r'$\ddot y$')
ax_comp.set_xlabel(r'$t$ [s]')
ax_comp.set_title('Components')
ax_comp.grid(True)
ax_comp.legend(loc='upper left')

dot_vx, = ax_comp.plot([], [], 'bo', markersize=6)
dot_vy, = ax_comp.plot([], [], 'ro', markersize=6)
dot_ax, = ax_comp.plot([], [], 'b^', markersize=6)
dot_ay, = ax_comp.plot([], [], 'r^', markersize=6)
time_line = ax_comp.axvline(0, color='k', alpha=0.3, linestyle=':')

v_scale = 0.3 / 1.5
a_scale = 0.3 / 2.5
frames = np.linspace(0.1, 5, 200)

def update(ti):
    pos = rr_func(ti).flatten()
    vel = rr_dot_func(ti).flatten()
    acc = rr_ddot_func(ti).flatten()
    vs = vel * v_scale
    acs = acc * a_scale
    dot.set_data([pos[0]], [pos[1]])
    v_arrow.xy = (pos[0] + vs[0], pos[1] + vs[1])
    v_arrow.set_position((pos[0], pos[1]))
    a_arrow.xy = (pos[0] + acs[0], pos[1] + acs[1])
    a_arrow.set_position((pos[0], pos[1]))
    dot_vx.set_data([ti], [vel[0]])
    dot_vy.set_data([ti], [vel[1]])
    dot_ax.set_data([ti], [acc[0]])
    dot_ay.set_data([ti], [acc[1]])
    time_line.set_xdata([ti, ti])
    return dot, v_arrow, a_arrow, dot_vx, dot_vy, dot_ax, dot_ay, time_line

fig.tight_layout()
ani = FuncAnimation(fig, update, frames=frames, blit=False, interval=50)
ani.save('graphics/kinematics_polar_arm.mp4', writer=FFMpegWriter(fps=30))
plt.close(fig)

Example 8: Curvilinear motion

The curvilinear motion of a particle is defined by \(v_x = 50 - 16t\) m/s and \(y = 100 - 4t^2\) m, with \(x(0) = 0\). Plot the path of the particle and determine the velocity and acceleration when \(y = 0\).

Solution

We integrate \(v_x\) with the initial condition \(x(0) = 0\) to obtain \(x(t)\), then assemble the position vector and differentiate to get velocity and acceleration.

import sympy as sp
from mechanicskit import la, ltx
t = sp.symbols('t', real=True, positive=True)
x = sp.Function('x')(t)
y_expr = 100 - 4*t**2

DE = sp.Eq(x.diff(t), 50 - 16*t)
ICs = {x.subs(t, 0): 0}
x_sol = sp.dsolve(DE, ics=ICs).rhs

ltx(r"x(t) &=", x_sol, r"\\ y(t) &=", y_expr, aligned=True)

The particle reaches \(y = 0\) when

t_y0 = sp.solve(sp.Eq(y_expr, 0), t)[0]
ltx(r"y = 0 \;\Rightarrow\; t =", t_y0, r"\text{ s}")

Assembling the position vector and differentiating twice gives the velocity and acceleration.

rr = sp.Matrix([x_sol, y_expr])
rr_dot = rr.diff(t)
rr_ddot = rr.diff(t, 2)

ltx(r"\bm r =", rr, r",\quad \dot{\bm r} =", rr_dot,
    r",\quad \ddot{\bm r} =", rr_ddot, arraystretch=1.5)

Evaluating the velocity at \(t = 5\) s gives the answer. The acceleration is already constant.

v_y0 = rr_dot.subs(t, t_y0)

ltx(r"\dot{\bm r}(5) =", v_y0, r"\text{ m/s},\quad \ddot{\bm r} =", rr_ddot,
    r"\text{ m/s}^2", arraystretch=1.5)

We plot the path with markers at integer seconds, and the velocity and acceleration vectors at \(y = 0\).

from mechanicskit import fplot, arrow

t1 = float(t_y0)
pos = rr.subs(t, t_y0)

fig, ax = plt.subplots(figsize=(7, 7))
fplot(rr[0], rr[1], (t, 0, t1), ax=ax, color='b', linewidth=2)

for ti in range(6):
    p = rr.subs(t, ti)
    ax.plot(float(p[0]), float(p[1]), 'ko', markersize=4)
    ax.annotate(rf'$t={ti}$', (float(p[0]), float(p[1])),
                textcoords='offset points', xytext=(8, -4), fontsize=11)

arrow(pos, rr_dot.subs(t, t_y0), ax=ax, color='b', label=r'$\dot{r}$')
arrow(pos, rr_ddot, ax=ax, color='r', label=r'$\ddot{r}$')

ax.set_aspect('equal')
ax.set(xlabel=r'$x(t)$ [m]', ylabel=r'$y(t)$ [m]', title='Path of the particle')
ax.grid(True)
ax.legend(loc='lower left')
plt.show()

The acceleration is constant because \(v_x\) is linear in \(t\), which gives \(\ddot x = -16\), and \(y\) is quadratic in \(t\), which gives \(\ddot y = -8\).

Finally we animate the motion with the velocity and acceleration vectors attached to the moving particle.

Code

from matplotlib.animation import FuncAnimation
from mechanicskit import to_responsive_html

t1 = float(t_y0)
rr_n = sp.lambdify(t, rr, 'numpy')
rr_dot_n = sp.lambdify(t, rr_dot, 'numpy')
acc = np.array([float(rr_ddot[0]), float(rr_ddot[1])])

fig, ax = plt.subplots(figsize=(7, 7))
fplot(rr[0], rr[1], (t, 0, t1), ax=ax, color='0.7', linewidth=1.5)
ax.set_aspect('equal')
ax.set(xlim=(-25, 85), ylim=(-50, 110),
       xlabel=r'$x$ [m]', ylabel=r'$y$ [m]', title=f'Curvilinear motion, t=0')
ax.grid(True)

dot, = ax.plot([], [], 'ko', markersize=7)
v_arrow = arrow([0, 0], [1, 0], ax=ax, color='b', label=r'$\dot{r}$')
a_arrow = arrow([0, 0], [1, 0], ax=ax, color='r', label=r'$\ddot{r}$')
ax.legend(loc='lower left')

def update(frame):
    ti = t1 * frame / 59
    ax.set_title(f'Curvilinear motion, t={ti:.2f}')
    p = rr_n(ti).flatten()
    v = rr_dot_n(ti).flatten()
    dot.set_data([p[0]], [p[1]])
    v_arrow.set_positions((p[0], p[1]), (p[0]+v[0], p[1]+v[1]))
    a_arrow.set_positions((p[0], p[1]), (p[0]+acc[0], p[1]+acc[1]))
    return dot, v_arrow, a_arrow

anim = FuncAnimation(fig, update, frames=60, interval=50, blit=False)
plt.close(fig)
to_responsive_html(anim)

Once Loop Reflect

Example 9: Lead screw, 3D helical motion

A power screw with lead \(L\) (the axial advancement per revolution) starts from rest and is rotated with an angular speed that increases uniformly with time, \(\dot\theta = k\,t\). The cross arm carries two balls at radial distance \(b\) from the screw axis. Determine the velocity and acceleration of the centre of ball \(A\) when the screw has turned through one complete revolution from rest, using \(k = 2\) rad/s\(^2\), \(b = 20\) mm, and \(L = 3\) mm.

Solution

The motion is a circular helix: ball \(A\) orbits the screw axis at constant radius \(b\) while the assembly descends along the axis at a rate proportional to the angular speed. We integrate the prescribed angular speed to obtain \(\theta(t)\), build the position vector in cylindrical form, and differentiate twice. Evaluating at the time when the screw has completed one full revolution yields the answer.

import sympy as sp
from mechanicskit import la, ltx
t = sp.symbols('t', real=True, positive=True)

k, b, L = 2, 20, 3
theta_f = sp.Function('theta')(t)
DE = sp.Eq(theta_f.diff(t), k*t)
ICs = {theta_f.subs(t, 0): 0}
theta = sp.dsolve(DE, ics=ICs).rhs

ltx(r"\theta(t) =", theta)

The screw has completed one full revolution when \(\theta = 2\pi\), which fixes the time of interest.

t_1 = sp.solve(sp.Eq(theta, 2*sp.pi), t)[0]
ltx(r"t_1 =", t_1, r"\;\text{s} \approx", float(t_1), r"\;\text{s}")

The position vector of the ball centre superposes the circular motion of the cross arm onto the axial descent of the screw,

\[ \bm r(t) = \left[b\cos\theta,\; b\sin\theta,\; - \frac{L\,\theta}{2\pi}\right]^\mathsf T. \]

Differentiating twice with respect to time gives the velocity and acceleration.

rr = sp.Matrix([b*sp.cos(theta), b*sp.sin(theta), -L*theta/(2*sp.pi)])
rr_dot = rr.diff(t)
rr_ddot = rr.diff(t, 2)

ltx(r"\bm r =", rr, r",\;\; \dot{\bm r} =", rr_dot,
    r",\;\; \ddot{\bm r} =", rr_ddot, arraystretch=1.5)

Substituting \(t_1\) delivers the position, velocity, and acceleration at the moment of one full revolution.

p1 = rr.subs(t, t_1)
v1 = sp.simplify(rr_dot.subs(t, t_1))
a1 = sp.simplify(rr_ddot.subs(t, t_1))

ltx(r"\bm r(t_1) =", p1, r",\;\; \dot{\bm r}(t_1) =", v1,
    r",\;\; \ddot{\bm r}(t_1) =", a1, arraystretch=1.5)

Numerically the velocity and acceleration at \(t_1\) are

import numpy as np
v1n = np.array(v1, dtype=float).flatten()
a1n = np.array(a1, dtype=float).flatten()

ltx(r"\dot{\bm r}(t_1) \approx", v1n.reshape(3, 1),
    r"\;\text{mm/s},\quad \ddot{\bm r}(t_1) \approx",
    a1n.reshape(3, 1), r"\;\text{mm/s}^2", arraystretch=1.5)

Two observations stand out. The horizontal velocity component points purely in the \(y\) direction because after a full revolution \(\theta = 2\pi\) places ball \(A\) back on the \(+x\) axis, where the tangent to the circle is along \(\bm j\). The acceleration is dominated by the centripetal term \(b\,\dot\theta^2\) in the \(-x\) direction, which at \(t_1\) equals \(b\,(2t_1)^2 = 8\pi b \approx 503\) mm/s\(^2\) and dwarfs the tangential and axial contributions.

We visualise the helical path from start to one full revolution, with the velocity and acceleration vectors attached at \(t_1\).

import matplotlib.pyplot as plt

rr_n = sp.lambdify(t, rr, 'numpy')
rr_dot_n = sp.lambdify(t, rr_dot, 'numpy')
rr_ddot_n = sp.lambdify(t, rr_ddot, 'numpy')

t_vals = np.linspace(0, float(t_1), 200)
path = np.array([rr_n(ti).flatten() for ti in t_vals])
p1n = rr_n(float(t_1)).flatten()

fig = plt.figure(figsize=(8, 7))
ax = fig.add_subplot(111, projection='3d')
ax.plot(path[:, 0], path[:, 1], path[:, 2], 'b-', linewidth=2, label='path')
ax.plot([b], [0], [0], 'ko', markersize=6, label='start')
ax.plot([p1n[0]], [p1n[1]], [p1n[2]], 'bo', markersize=6,
        label=r'$\theta = 2\pi$')

ax.quiver(*p1n, *(v1n*0.1), color='b', linewidth=2, label='velocity')
ax.quiver(*p1n, *(a1n*0.02), color='r', linewidth=2, label='acceleration')

ax.set(xlabel=r'$x$ [mm]', ylabel=r'$y$ [mm]', zlabel=r'$z$ [mm]')
ax.set_xlim(-25, 25); ax.set_ylim(-25, 25); ax.set_zlim(-25, 25)
ax.view_init(elev=30, azim=48)
ax.legend(loc='upper left')
plt.show()

We extend the simulation to ten full revolutions to expose the helical structure of the trajectory and animate the velocity and acceleration vectors as they sweep around the screw axis.

from matplotlib.animation import FuncAnimation, FFMpegWriter

t_end = float(sp.solve(sp.Eq(theta, 20*sp.pi), t)[0])
t_anim = np.linspace(0, t_end, 600)
path_long = np.array([rr_n(ti).flatten() for ti in t_anim])

fig = plt.figure(figsize=(8, 7))
ax = fig.add_subplot(111, projection='3d')
ax.plot(path_long[:, 0], path_long[:, 1], path_long[:, 2],
        color='0.7', linewidth=1.0)
zmin = float(path_long[:, 2].min()) - 2
ax.set(xlabel=r'$x$ [mm]', ylabel=r'$y$ [mm]', zlabel=r'$z$ [mm]')
ax.set_xlim(-25, 25); ax.set_ylim(-25, 25); ax.set_zlim(zmin, 2)
ax.view_init(elev=25, azim=48)

dot, = ax.plot([], [], [], 'ko', markersize=6)
quivers = {'v': ax.quiver(0, 0, 0, 0, 0, 0, color='b', linewidth=2),
           'a': ax.quiver(0, 0, 0, 0, 0, 0, color='r', linewidth=2)}
title = ax.set_title('')

frames = 300
def update(frame):
    ti = t_end * frame / (frames - 1)
    p = rr_n(ti).flatten()
    v = rr_dot_n(ti).flatten()
    a = rr_ddot_n(ti).flatten()
    dot.set_data([p[0]], [p[1]])
    dot.set_3d_properties([p[2]])
    quivers['v'].remove()
    quivers['a'].remove()
    quivers['v'] = ax.quiver(*p, *(v*0.02), color='b', linewidth=2)
    quivers['a'] = ax.quiver(*p, *(a*0.004), color='r', linewidth=2)
    title.set_text(rf'$t = {ti:.2f}$ s, $\theta = {float(theta.subs(t, ti)):.2f}$ rad')
    return dot, quivers['v'], quivers['a'], title

anim = FuncAnimation(fig, update, frames=frames, interval=50, blit=False)
anim.save('graphics/kinematics_lead_screw.mp4', writer=FFMpegWriter(fps=30))
plt.close(fig)

Example 10: Tractor hoisting, constrained motion

A tractor at \(A\) moves horizontally with forward velocity \(v_A\) and hoists a bale \(B\) via a rope routed over a fixed pulley at \(C\), which sits at height \(h\) above the ground. The rope passes from the tractor up to the top pulley, descends to a movable pulley fastened to the bale, and returns to a fixed anchor at the same height. The bale rises along a vertical line directly below \(C\). Determine the upward velocity \(v_B\) of the bale as a function of \(v_A\) and the tractor position \(x\).

Solution

The kinematic relation between the two velocities follows entirely from geometry and the inextensibility of the rope. Two constraints capture the setup. The right triangle formed by the tractor, the pulley, and the ground gives

\[ x^2 + h^2 = l^2, \]

where \(l\) is the slanted rope segment from \(A\) to \(C\). The total rope length \(L\) is fixed,

\[ L = 2(h - y) + l, \]

since the rope between the top pulley and the bale comprises two parallel segments of length \(h - y\). Differentiating both constraints with respect to time and solving the resulting linear pair for \(v_B\) delivers the answer.

import sympy as sp
from mechanicskit import la, ltx
t = sp.symbols('t', real=True, positive=True)
h, L = sp.symbols('h L', positive=True)
v_A, v_B, l_dot = sp.symbols('v_A v_B l_dot', real=True)
x_t = sp.Function('x')(t)
y_t = sp.Function('y')(t)
l_t = sp.Function('l')(t)

c1 = sp.Eq(x_t**2 + h**2, l_t**2)
c2 = sp.Eq(L, 2*(h - y_t) + l_t)

ltx(r"\text{(geometry)}\quad ", c1, r",\quad \text{(rope)}\quad ", c2)

Differentiating each constraint with respect to \(t\) and replacing \(\dot x\), \(\dot y\), \(\dot l\) with the named velocities \(v_A\), \(v_B\), \(\dot l\) produces a linear system in the rate variables.

subs = {x_t.diff(t): v_A, y_t.diff(t): v_B, l_t.diff(t): l_dot}
c1_dot = sp.Eq(c1.lhs.diff(t), c1.rhs.diff(t)).subs(subs)
c2_dot = sp.Eq(c2.lhs.diff(t), c2.rhs.diff(t)).subs(subs)

ltx(c1_dot, r",\quad", c2_dot)

We solve the system for \(v_B\) and \(\dot l\), then eliminate \(l\) in favour of \(x\) using \(l = \sqrt{x^2 + h^2}\).

x = sp.symbols('x', positive=True)
sol = sp.solve([c1_dot, c2_dot], [v_B, l_dot])
v_B_expr = sp.simplify(sol[v_B].subs(l_t, sp.sqrt(x**2 + h**2)).subs(x_t, x))
y_expr = sp.solve(c2.subs(l_t, sp.sqrt(x**2 + h**2)), y_t)[0]

ltx(r"v_B =", v_B_expr, r",\quad y(x) =", y_expr)

The bale velocity vanishes at \(x = 0\) and approaches \(v_A/2\) as \(x \to \infty\), since the rope segment \(l\) then aligns with the horizontal and the tractor essentially reels in rope at its own speed, with the factor of one half coming from the two-to-one mechanical advantage of the movable pulley.

For the numerical case we take \(h = 5\) m and a rope length \(L = 3h = 15\) m, so that the bale rests on the ground when the tractor sits beneath the pulley. The tractor advances until the bale reaches a height of \(y = 4\) m, which by inversion of \(y(x)\) corresponds to a tractor position \(x = 12\) m.

params = {h: 5, L: 15}
y_x = y_expr.subs(params)
v_B_over_v_A = (v_B_expr/v_A).subs(params)
x_final = sp.solve(sp.Eq(y_x, 4), x)[0]

ltx(r"y(x) =", y_x, r",\quad \frac{v_B}{v_A} =", v_B_over_v_A,
    r",\quad x|_{y=4} =", x_final)

We sketch the geometry at \(x = 4\) m, drawing the two rope segments from the tractor to the top pulley and from the pulley down to the bale.

import numpy as np
import matplotlib.pyplot as plt

y_n = sp.lambdify(x, y_x, 'numpy')
x_i = 4.0
A_i = np.array([x_i, 0.0])
B_i = np.array([0.0, float(y_n(x_i))])
C = np.array([0.0, 5.0])

fig, ax = plt.subplots(figsize=(8, 4))
ax.plot([C[0], A_i[0]], [C[1], A_i[1]], '-ko', markerfacecolor='k', linewidth=2)
ax.plot([C[0], B_i[0]], [C[1], B_i[1]], '-ko', markerfacecolor='k', linewidth=2)
ax.plot([0, 15], [0, 0], 'k-', linewidth=0.8)
ax.annotate(r'$A$', A_i, xytext=(8, 8), textcoords='offset points', fontsize=12)
ax.annotate(r'$B$', B_i, xytext=(8, 8), textcoords='offset points', fontsize=12)
ax.annotate(r'$C$', C,   xytext=(8, 8), textcoords='offset points', fontsize=12)
ax.set_aspect('equal')
ax.set(xlim=(-1, 15), ylim=(-0.5, 6.5),
       xlabel=r'$x$ [m]', ylabel=r'$y$ [m]')
ax.grid(True)
plt.show()

Animating \(x\) from \(0\) to \(12\) m shows how the bale climbs and the rope geometry evolves. The animation is exported to mp4 and embedded below.

from matplotlib.animation import FuncAnimation, FFMpegWriter

fig, ax = plt.subplots(figsize=(8, 4))
ax.plot([0, 15], [0, 0], 'k-', linewidth=0.8)
line_AC, = ax.plot([], [], '-ko', markerfacecolor='k', linewidth=2)
line_BC, = ax.plot([], [], '-ko', markerfacecolor='k', linewidth=2)
ax.plot(*C, 'ko', markerfacecolor='k')
txt_A = ax.text(0, 0, r'$A$', fontsize=12)
txt_B = ax.text(0, 0, r'$B$', fontsize=12)
ax.text(C[0]+0.2, C[1]+0.2, r'$C$', fontsize=12)
ax.set_aspect('equal')
ax.set(xlim=(-1, 15), ylim=(-0.5, 6.5),
       xlabel=r'$x$ [m]', ylabel=r'$y$ [m]')
ax.grid(True)
title = ax.set_title('')

frames = 150
x_path = np.linspace(0, float(x_final), frames)

def update(frame):
    xi = x_path[frame]
    A = np.array([xi, 0.0])
    B = np.array([0.0, float(y_n(xi))])
    line_AC.set_data([C[0], A[0]], [C[1], A[1]])
    line_BC.set_data([C[0], B[0]], [C[1], B[1]])
    txt_A.set_position((A[0]+0.2, A[1]+0.2))
    txt_B.set_position((B[0]+0.2, B[1]+0.2))
    title.set_text(rf'$x = {xi:.2f}$ m, $y = {B[1]:.2f}$ m')
    return line_AC, line_BC, txt_A, txt_B, title

anim = FuncAnimation(fig, update, frames=frames, interval=50, blit=False)
anim.save('graphics/kinematics_tractor_hoist.mp4', writer=FFMpegWriter(fps=30))
plt.close(fig)

Example 11: Crank-slider mechanism

Study the motion of a slider-crank mechanism running at constant angular velocity \(\omega = 25\) rad/s. The crank has length \(r = 125\) mm and the connecting rod has length \(L = 350\) mm. The centre of mass \(G\) of the connecting rod sits \(100\) mm from the crank pin \(B\) along the rod, so \(250\) mm from the slider pin \(A\). Determine the velocity and acceleration of points \(A\), \(B\), and \(G\) as functions of the crank angle \(\theta\), and plot their magnitudes over a full revolution.

Solution

We place the crank pivot \(O\) at the origin and let the slider \(A\) travel along the negative \(x\)-axis. With the crank angle \(\theta\) measured from the horizontal, the crank pin sits at

\[ \bm r_B = r\bigl(-\cos\theta,\; \sin\theta\bigr). \]

The slider pin lies on the \(x\)-axis at horizontal distance \(a + b\) from \(O\), where \(a = r\cos\theta\) is the horizontal projection of the crank and \(b = \sqrt{L^2 - h^2}\) comes from Pythagoras applied to the connecting rod with \(h = r\sin\theta\). Hence

\[ \bm r_A = \bigl(-(a + b),\; 0\bigr). \]

The centre of mass \(G\) partitions the rod with \(BG = 100\) mm out of \(L = 350\) mm, so

\[ \bm r_G = \bm r_B + \frac{BG}{L}\bigl(\bm r_A - \bm r_B\bigr). \]

Since \(\theta = \omega t\) with \(\dot\theta = \omega\) constant and \(\ddot\theta = 0\), time derivatives reduce to derivatives with respect to \(\theta\),

\[ \dot{\bm r} = \omega\,\frac{\mathrm d \bm r}{\mathrm d \theta}, \qquad \ddot{\bm r} = \omega^2\,\frac{\mathrm d^2 \bm r}{\mathrm d \theta^2}, \]

which lets us treat \(\theta\) as the independent variable throughout.

import sympy as sp
from mechanicskit import la, ltx
theta, omega = sp.symbols('theta omega', positive=True)
r, L, BG = sp.Rational(125, 1000), sp.Rational(350, 1000), sp.Rational(100, 1000)
omega_val = 25

a = r*sp.cos(theta)
h = r*sp.sin(theta)
b = sp.sqrt(L**2 - h**2)

rr_B = sp.Matrix([-r*sp.cos(theta), r*sp.sin(theta)])
rr_A = sp.Matrix([-(a + b), 0])
rr_G = sp.simplify(rr_B + (BG/L)*(rr_A - rr_B))

ltx(r"\bm r_A &=", rr_A, r",\\[1em]  \bm r_B &=", rr_B,
    r",\\[1em] \bm r_G &=", rr_G, aligned=True,arraystretch=1.5)

Differentiating each position with respect to \(\theta\) (and multiplying by \(\omega\) or \(\omega^2\) for the time derivatives) gives the velocity and acceleration. The crank pin moves on a circle, so \(|\dot{\bm r}_B| = r\omega\) and \(|\ddot{\bm r}_B| = r\omega^2\) are constants. The slider velocity vanishes twice per revolution, at the dead-centre positions \(\theta = 0\) and \(\theta = \pi\), which forces large slider accelerations near those angles.

def deriv(rr):
    return omega*rr.diff(theta), omega**2*rr.diff(theta, 2)

rr_A_dot, rr_A_ddot = deriv(rr_A)
rr_B_dot, rr_B_ddot = deriv(rr_B)
rr_G_dot, rr_G_ddot = deriv(rr_G)

Substituting \(\omega = 25\) rad/s, we plot the slider velocity \(\dot x_A\) alongside the magnitudes \(|\dot{\bm r}_B|\) and \(|\dot{\bm r}_G|\). The crank pin runs at the constant \(|\dot{\bm r}_B| = r\omega \approx 3.13\) m/s, while the slider velocity oscillates between \(\pm\) a peak just above \(r\omega\) and vanishes at the dead-centre positions \(\theta = 0\) and \(\theta = \pi\).

import numpy as np
import matplotlib.pyplot as plt

vAx_n = sp.lambdify(theta, rr_A_dot[0].subs(omega, omega_val), 'numpy')
vB_mag = sp.sqrt(rr_B_dot.dot(rr_B_dot)).subs(omega, omega_val)
vG_mag = sp.sqrt(rr_G_dot.dot(rr_G_dot)).subs(omega, omega_val)
vB_n = sp.lambdify(theta, vB_mag, 'numpy')
vG_n = sp.lambdify(theta, vG_mag, 'numpy')

th = np.linspace(0, 2*np.pi, 400)
deg = np.degrees(th)
fig, ax = plt.subplots(figsize=(8, 5))
ax.plot(deg, vAx_n(th), color='C0', label=r'$\dot{x}_A$', linewidth=2)
ax.plot(deg, np.full_like(th, float(vB_n(0))), color='C1',
        label=r'$|\dot{r}_B|$', linewidth=2)
ax.plot(deg, vG_n(th), color='C2', label=r'$|\dot{r}_G|$', linewidth=2)
ax.axhline(0, color='k', linewidth=0.5)
ax.set(xlabel=r'$\theta$ [deg]', ylabel='velocity [m/s]',
       title='Crank-slider velocities')
ax.set_xticks(range(0, 361, 45))
ax.grid(True); ax.legend(loc='lower left')
plt.show()

The acceleration plot reveals the asymmetry inherent to a finite connecting rod: the slider acceleration \(\ddot r_{Ax}\) swings between large positive values near top dead centre and smaller negative values around bottom dead centre, while \(|\ddot{\bm r}_B| = r\omega^2 \approx 78.1\) m/s\(^2\) remains constant.

aAx = rr_A_ddot[0].subs(omega, omega_val)
aB_mag = sp.sqrt(rr_B_ddot.dot(rr_B_ddot)).subs(omega, omega_val)
aG_mag = sp.sqrt(rr_G_ddot.dot(rr_G_ddot)).subs(omega, omega_val)

aAx_n = sp.lambdify(theta, aAx, 'numpy')
aB_n = sp.lambdify(theta, aB_mag, 'numpy')
aG_n = sp.lambdify(theta, aG_mag, 'numpy')

fig, ax = plt.subplots(figsize=(8, 5))
ax.plot(np.degrees(th), aAx_n(th), label=r'$\ddot{r}_{Ax}$', linewidth=2)
ax.plot(np.degrees(th), np.full_like(th, aB_n(0)),
        label=r'$|\ddot{r}_B|$', linewidth=2)
ax.plot(np.degrees(th), aG_n(th), label=r'$|\ddot{r}_G|$', linewidth=2)
ax.set(xlabel=r'$\theta$ [deg]', ylabel=r'acceleration [m/s$^2$]',
       title='Crank-slider accelerations')
ax.set_xticks(range(0, 361, 45))
ax.grid(True); ax.legend(loc='lower left')
plt.show()

Finally we animate one full revolution as a synchronised triple view: the mechanism on the left with velocity (blue) and acceleration (red) arrows attached to \(A\), \(B\), and \(G\), the speed plot in the upper right with markers tracking the current angle, and the acceleration plot in the lower right with the same. The animation is rendered to mp4 at 60 fps and embedded below.

from matplotlib.animation import FuncAnimation, FFMpegWriter
from mechanicskit import arrow

rr_A_n = sp.lambdify(theta, rr_A, 'numpy')
rr_B_n = sp.lambdify(theta, rr_B, 'numpy')
rr_G_n = sp.lambdify(theta, rr_G, 'numpy')

vA_v = sp.lambdify(theta, rr_A_dot.subs(omega, omega_val), 'numpy')
vB_v = sp.lambdify(theta, rr_B_dot.subs(omega, omega_val), 'numpy')
vG_v = sp.lambdify(theta, rr_G_dot.subs(omega, omega_val), 'numpy')
aA_v = sp.lambdify(theta, rr_A_ddot.subs(omega, omega_val), 'numpy')
aB_v = sp.lambdify(theta, rr_B_ddot.subs(omega, omega_val), 'numpy')
aG_v = sp.lambdify(theta, rr_G_ddot.subs(omega, omega_val), 'numpy')

th_curve = np.linspace(0, 2*np.pi, 400)
deg_curve = np.degrees(th_curve)
vAx_curve = vAx_n(th_curve)
vB_curve = np.full_like(th_curve, float(vB_n(0)))
vG_curve = vG_n(th_curve)
aA_curve = aAx_n(th_curve)
aB_curve = np.full_like(th_curve, float(aB_n(0)))
aG_curve = aG_n(th_curve)

fig = plt.figure(figsize=(14, 7))
gs = fig.add_gridspec(2, 2, width_ratios=[1.6, 1])
ax_mech = fig.add_subplot(gs[:, 0])
ax_v = fig.add_subplot(gs[0, 1])
ax_a = fig.add_subplot(gs[1, 1])

ax_mech.plot(0, 0, 'ks', markersize=8)
ax_mech.text(0.005, -0.04, r'$O$', fontsize=12)
crank, = ax_mech.plot([], [], '-o', color='C0', markerfacecolor='C0',
                     markersize=8, linewidth=3)
rod, = ax_mech.plot([], [], '-o', color='C1', markerfacecolor='C1',
                   markersize=8, linewidth=3)
slider, = ax_mech.plot([], [], 's', color='k', markerfacecolor='0.7',
                      markersize=18)
gpoint, = ax_mech.plot([], [], 'o', color='C3', markerfacecolor='C3',
                      markersize=10)
txt_A = ax_mech.text(0, 0, r'$A$', fontsize=12)
txt_B = ax_mech.text(0, 0, r'$B$', fontsize=12)
txt_G = ax_mech.text(0, 0, r'$G$', fontsize=12)

vs, as_ = 0.04, 0.001
v_arr = {pt: arrow([0, 0], [1e-3, 0], ax=ax_mech, color='b', linewidth=2)
         for pt in 'ABG'}
a_arr = {pt: arrow([0, 0], [1e-3, 0], ax=ax_mech, color='r', linewidth=2)
         for pt in 'ABG'}

ax_mech.set_aspect('equal')
ax_mech.set(xlim=(-0.6, 0.2), ylim=(-0.22, 0.22),
           xlabel=r'$x$ [m]', ylabel=r'$y$ [m]',
           title=f'$\\theta=0^\\circ$')
ax_mech.grid(True)

ax_v.plot(deg_curve, vAx_curve, color='C0', label=r'$\dot{x}_A$', linewidth=1.5)
ax_v.plot(deg_curve, vB_curve, color='C1', label=r'$|\dot{r}_B|$', linewidth=1.5)
ax_v.plot(deg_curve, vG_curve, color='C2', label=r'$|\dot{r}_G|$', linewidth=1.5)
ax_v.axhline(0, color='k', linewidth=0.5)
mvA, = ax_v.plot([], [], 'o', color='C0', markersize=8)
mvB, = ax_v.plot([], [], 'o', color='C1', markersize=8)
mvG, = ax_v.plot([], [], 'o', color='C2', markersize=8)
ax_v.set(xlabel=r'$\theta$ [deg]', ylabel='velocity [m/s]', xlim=(0, 360))
ax_v.set_xticks(range(0, 361, 90))
ax_v.grid(True); ax_v.legend(loc='lower left')

ax_a.plot(deg_curve, aA_curve, color='C0', label=r'$\ddot{r}_{Ax}$', linewidth=1.5)
ax_a.plot(deg_curve, aB_curve, color='C1', label=r'$|\ddot{r}_B|$', linewidth=1.5)
ax_a.plot(deg_curve, aG_curve, color='C2', label=r'$|\ddot{r}_G|$', linewidth=1.5)
maA, = ax_a.plot([], [], 'o', color='C0', markersize=8)
maB, = ax_a.plot([], [], 'o', color='C1', markersize=8)
maG, = ax_a.plot([], [], 'o', color='C2', markersize=8)
ax_a.set(xlabel=r'$\theta$ [deg]', ylabel=r'accel [m/s$^2$]', xlim=(0, 360))
ax_a.set_xticks(range(0, 361, 90))
ax_a.grid(True); ax_a.legend(loc='lower left')
fig.tight_layout()

frames = 240
def update(frame):
    th_i = 2*np.pi * frame / (frames - 1)
    ax_mech.set_title(f'$\\theta=${int(np.degrees(th_i))}$^\\circ$')
    deg_i = np.degrees(th_i)
    A = rr_A_n(th_i).flatten()
    B = rr_B_n(th_i).flatten()
    G = rr_G_n(th_i).flatten()
    pts = {'A': A, 'B': B, 'G': G}
    vels = {'A': np.array(vA_v(th_i), dtype=float).flatten(),
            'B': np.array(vB_v(th_i), dtype=float).flatten(),
            'G': np.array(vG_v(th_i), dtype=float).flatten()}
    accs = {'A': np.array(aA_v(th_i), dtype=float).flatten(),
            'B': np.array(aB_v(th_i), dtype=float).flatten(),
            'G': np.array(aG_v(th_i), dtype=float).flatten()}

    crank.set_data([0, B[0]], [0, B[1]])
    rod.set_data([B[0], A[0]], [B[1], A[1]])
    slider.set_data([A[0]], [A[1]])
    gpoint.set_data([G[0]], [G[1]])
    txt_A.set_position((A[0]-0.015, A[1]+0.025))
    txt_B.set_position((B[0]+0.005, B[1]+0.020))
    txt_G.set_position((G[0]+0.005, G[1]+0.020))

    for pt in 'ABG':
        p = pts[pt]; v = vels[pt]; a = accs[pt]
        v_arr[pt].set_positions((p[0], p[1]),
                                (p[0]+v[0]*vs, p[1]+v[1]*vs))
        a_arr[pt].set_positions((p[0], p[1]),
                                (p[0]+a[0]*as_, p[1]+a[1]*as_))

    mvA.set_data([deg_i], [vels['A'][0]])
    mvB.set_data([deg_i], [np.linalg.norm(vels['B'])])
    mvG.set_data([deg_i], [np.linalg.norm(vels['G'])])
    maA.set_data([deg_i], [accs['A'][0]])
    maB.set_data([deg_i], [np.linalg.norm(accs['B'])])
    maG.set_data([deg_i], [np.linalg.norm(accs['G'])])
    return ()

anim = FuncAnimation(fig, update, frames=frames, interval=1000/60, blit=False)
anim.save('graphics/kinematics_crank_slider.mp4',
          writer=FFMpegWriter(fps=60))
plt.close(fig)

Three features of the animation deserve comment. The slider velocity and acceleration run out of phase. At the dead-centre positions \(\theta = 0\) and \(\theta = \pi\) the velocity \(\dot x_A\) vanishes, yet the acceleration \(\ddot r_{Ax}\) takes its largest values there. Conversely, where \(\dot x_A\) peaks, \(\ddot r_{Ax}\) has already crossed zero. The relationship is not the rigid quarter-period offset of simple harmonic motion, since the finite connecting rod skews both curves.

The velocity and acceleration at \(G\) are not always perpendicular to one another. The crank pin \(B\) travels on a circle, so its acceleration is purely centripetal and orthogonal to the tangential velocity at every instant. The point \(G\) traces a more complicated curve, and its acceleration carries components both along and across the velocity direction, which the animation makes visible as the blue and red arrows at \(G\) rarely meet at right angles.

Finally, \(\ddot r_{Ax}\) vanishes at the angles where the acceleration vectors at \(B\) and \(G\) stand vertical. The slider can only be accelerated through the horizontal projection of the pin accelerations transmitted by the connecting rod. When that projection cancels, the rod imparts no horizontal impulse to \(A\), and the slider coasts at its peak speed.