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:
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.
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.2.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.
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.2.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.
In the special case of rectilinear motion (motion without curvature), such as described in Figure 6.2.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.
Using scalar notation \(s = |\boldsymbol{r}|\), \(v = |\dot{\boldsymbol{r}}|\), \(a = |\ddot{\boldsymbol{r}}|\) and eliminating \(dt\) from \(v = ds/dt\) and \(a = dv/dt\) gives
\[ \boxed{a \, ds = v \, dv} \tag{6.2.3}\]
This relation is useful when time is not of interest, but it only applies to rectilinear (or scalar) problems and hides the assumptions behind it. We generally prefer the full vector formulation and the differential equation approach described next.
The proposed working method is described in what follows.
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.
The workflow in computational dynamics is summarized below.
To move between position \(\boldsymbol{r}(t)\), velocity \(\dot{\boldsymbol{r}}(t)\) and acceleration \(\ddot{\boldsymbol{r}}(t)\) we differentiate in one direction and integrate in the other. In practice we rarely integrate by hand; instead we solve an ODE whose solution gives \(\boldsymbol{r}(t)\) directly, and then differentiate to obtain the quantities we need. For nonlinear ODEs we typically resort to numerical integration, which generates the velocity \(\dot{\boldsymbol{r}}\) as an intermediate result on the way to the position \(\boldsymbol{r}\).
Newtonian physics is formulated as the second-order differential equation
\[ \sum \boldsymbol{F} = m\ddot{\boldsymbol{r}} \]
but in kinematics we see equally often first-order equations.
We formulate an ODE that must be (at least piecewise) valid for all \(t\) and apply initial conditions at the start of the analysis, or at some other known state, to fix the unknown constants. The number of constants depends on the order of the differential equation. In this approach, time \(t\) is the central parameter tying \(\boldsymbol{r}\), \(\dot{\boldsymbol{r}}\) and \(\ddot{\boldsymbol{r}}\) together, which is why it often appears in the solution even when we did not explicitly ask for it. This is the nature of Newtonian mechanics: it is inherently transient, meaning the state of the body is known at every instant. As we shall see in later chapters, there are methods that bypass the time dependence when it is not of interest (e.g., 6.2.3), but a complete dynamic analysis always requires solving the ODE.
The questions we ask of the model fall into two categories. In the first, the equations are explicitly functions of time, e.g., \(\dot{\boldsymbol{r}}(t)\), so we can evaluate any quantity at a given instant. In the second, we want a quantity expressed as a function of another quantity, for example the velocity at a given position \(\dot{\boldsymbol{r}}(\boldsymbol{r})\), and must solve for the time that is common to both. A useful trick to avoid solving is to simply plot \(\dot{\boldsymbol{r}}(t)\) against \(\boldsymbol{r}(t)\) and read off the answer graphically.
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 can avoid this whole exercise of practicing to compute velocities and acceleration in various bases by just formulating relations in Cartesian coordinates and directly work with the position vector \(\boldsymbol{r}\) which can be easily projected onto any other bases if needed. Let the computer do the tedious mathematics!
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:
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 natural coordinate system is attached to the particle and moves with it along its path. This makes it particularly useful when we are interested in how fast the particle travels along the path and how sharply it turns, rather than its absolute position in space.
The first basis vector, the tangent vector \(\boldsymbol{e}_t\), points in the direction the particle is currently moving. Since the velocity \(\dot{\boldsymbol{r}}\) is always tangent to the path, we obtain the tangent direction by normalizing it,
\[ \boxed{\boldsymbol{e}_t = \frac{\dot{\boldsymbol{r}}}{|\dot{\boldsymbol{r}}|} } \]
The second basis vector, the normal vector \(\boldsymbol{e}_n\), points toward the center of curvature, that is, in the direction the path is currently bending. We obtain it from the time derivative of \(\boldsymbol{e}_t\) which, as shown in Section 17.8, is always perpendicular to \(\boldsymbol{e}_t\) itself,
\[ \boxed{\boldsymbol{e}_n = \frac{\dot{\boldsymbol{e}}_t}{|\dot{\boldsymbol{e}}_t|}} \]
Together, \(\boldsymbol{e}_t\) and \(\boldsymbol{e}_n\) span the plane in which the path locally curves. The third basis vector, the bi-normal \(\boldsymbol{e}_b\), completes the right-handed system by pointing out of that plane,
\[ \boxed{ \boldsymbol{e}_b := \boldsymbol{e}_t \times \boldsymbol{e}_n } \]
For planar problems (2D) the path always curves in the same plane, so \(\boldsymbol{e}_b = \boldsymbol{e}_z\) throughout the motion.
The polar coordinate system is defined relative to a fixed origin. It separates the motion into a radial component (toward or away from the origin) and an angular component (around it), which makes it particularly useful when the radius changes over time or when that change is of interest.
The radial basis vector \(\boldsymbol{e}_r\) points from the origin toward the particle. We obtain it by normalizing the position vector,
\[ \boxed{\boldsymbol{e}_r = \frac{\boldsymbol{r}(t)}{|\boldsymbol{r}(t)|}} \]
The angular basis vector \(\boldsymbol{e}_\theta\) points in the direction of increasing \(\theta\), perpendicular to \(\boldsymbol{e}_r\). As the particle moves, \(\boldsymbol{e}_r\) rotates, and its time derivative (see Section 17.8) gives precisely this perpendicular direction,
\[ \boxed{\boldsymbol{e}_\theta = \frac{\dot{\boldsymbol{e}}_r}{|\dot{\boldsymbol{e}}_r|}} \]
In three dimensions a third basis vector \(\boldsymbol{e}_z\) is added along the axis of rotation, giving the cylindrical coordinate system.
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}\).
Remark 17.1. We verify the result above using \(\boldsymbol{r}(t) = [\cos{\theta(t)}, \sin{\theta(t)}]\).
Be careful when defining angles as degrees! Angles will be functions of time in what follow, meaning implicit derivation will yield expressions where angles are directly multiplied with distances or velocities, this only works if the angles are defined as radians!
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)\[ \bm r =\left[\begin{matrix}\cos{\left(\theta{\left(t \right)} \right)}\\\sin{\left(\theta{\left(t \right)} \right)}\end{matrix}\right] \]
Normalizing gives \(\boldsymbol{e}_r\),
ee_r = sp.trigsimp(rr.normalized())
ltx("\\bm e_r =", ee_r)\[ \bm e_r =\left[\begin{matrix}\cos{\left(\theta{\left(t \right)} \right)}\\\sin{\left(\theta{\left(t \right)} \right)}\end{matrix}\right] \]
We cannot simplify the trigonometric expressions if the angle \(\theta\) is not assumed to be a real valued function.
The time derivative of the unit vector \(\boldsymbol{e}_r\) is
ee_r_dot = ee_r.diff(t)
ltx(r"\dot{\bm e}_r =", ee_r_dot, arraystretch=2.5)\[ {\def\arraystretch{2.5}\dot{\bm e}_r =\left[\begin{matrix}- \sin{\left(\theta{\left(t \right)} \right)} \dfrac{d}{d t} \theta{\left(t \right)}\\\cos{\left(\theta{\left(t \right)} \right)} \dfrac{d}{d t} \theta{\left(t \right)}\end{matrix}\right]} \]
The result is a scaled vector orthogonal to \(\boldsymbol{e}_r\). Note that \(\dot{\theta}\) appears from implicit differentiation, which is why \(|\dot{\boldsymbol{e}}_r| \neq 1\) in general. The factor \(\dot{\theta}\) is the angular velocity \(\omega\), and its sign gives the direction of rotation.
We can therefore write
\[ \boxed{\dot{\boldsymbol{e}}_r = \dot{\theta}\boldsymbol{e}_\theta} \]
and obtain \(\boldsymbol{e}_\theta\) by normalizing \(\dot{\boldsymbol{e}}_r\),
\[ \boxed{\boldsymbol{e}_\theta = \frac{\dot{\boldsymbol{e}}_r}{|\dot{\boldsymbol{e}}_r|}} \]
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)\[ {\def\arraystretch{2.5}{\bm e}_\theta =\left[\begin{matrix}- \dfrac{\dot{\theta} \sin{\left(\theta{\left(t \right)} \right)}}{\left|{\dot{\theta}}\right|}\\\dfrac{\dot{\theta} \cos{\left(\theta{\left(t \right)} \right)}}{\left|{\dot{\theta}}\right|}\end{matrix}\right]} \]
The factor \(\frac{\dot{\theta}}{|\dot{\theta}|}\) is simply the sign of \(\dot{\theta}\), so for \(\dot{\theta} > 0\) we get \(\boldsymbol{e}_\theta = [-\sin\theta, \cos\theta]^\mathsf{T}\). We confirm orthogonality by computing the dot product,
\[ \boldsymbol{e}_r \cdot \boldsymbol{e}_\theta = [\cos\theta, \sin\theta] \cdot [-\sin\theta, \cos\theta] = -\cos\theta\sin\theta + \sin\theta\cos\theta = 0 \]
which verifies that \(\boldsymbol{e}_r \perp \dot{\boldsymbol{e}}_r\) and \(\boldsymbol{e}_r \perp \boldsymbol{e}_\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.
Taking the derivative of the position vector using the product rule gives
\[ \dot{\boldsymbol{r}} = \frac{d}{dt}(r(t)\boldsymbol{e}_r(t)) = \dot{r}\boldsymbol{e}_r + r\dot{\boldsymbol{e}}_r \tag{6.2.4}\]
and with \(\dot{\boldsymbol{e}}_r = \dot{\theta}\boldsymbol{e}_\theta\) from the previous section we have
\[ \boxed{\dot{\boldsymbol{r}} = \dot{r}\boldsymbol{e}_r + r\dot{\theta}\boldsymbol{e}_\theta} \]
The quantity \(\omega = \dot{\theta}\) is the angular velocity, and \(v_\theta = r\dot{\theta}\) is the tangential speed. The first term \(\dot{r}\boldsymbol{e}_r\) captures how fast the particle moves toward or away from the origin, while the second term \(r\dot{\theta}\boldsymbol{e}_\theta\) captures how fast it sweeps around it.
Taking the time derivative of 6.2.4 and applying the product rule to each term gives
\[ \ddot{\boldsymbol{r}} = (\ddot{r}\boldsymbol{e}_r + \dot{r}\dot{\boldsymbol{e}}_r) + (\dot{r}\dot{\theta}\boldsymbol{e}_\theta + r\ddot{\theta}\boldsymbol{e}_\theta + r\dot{\theta}\dot{\boldsymbol{e}}_\theta) \]
We already know that \(\dot{\boldsymbol{e}}_r = \dot{\theta}\boldsymbol{e}_\theta\). By the same argument (differentiating \(\boldsymbol{e}_\theta \cdot \boldsymbol{e}_\theta = 1\)), the time derivative \(\dot{\boldsymbol{e}}_\theta\) must be perpendicular to \(\boldsymbol{e}_\theta\), and direct computation shows \(\dot{\boldsymbol{e}}_\theta = -\dot{\theta}\boldsymbol{e}_r\). Substituting both and collecting terms,
\[ \boxed{\ddot{\boldsymbol{r}} = (\ddot{r} - r\dot{\theta}^2)\boldsymbol{e}_r + (r\ddot{\theta} + 2\dot{r}\dot{\theta})\boldsymbol{e}_\theta} \]
The radial component has two contributions: \(\ddot{r}\) is the acceleration of the particle outward from the origin, while \(r\dot{\theta}^2\) points inward and is the centripetal acceleration that arises from the circular part of the motion. In the angular direction, \(r\ddot{\theta} = r\alpha\) is the angular acceleration and \(2\dot{r}\dot{\theta}\) is the Coriolis term. To see where it comes from, consider a particle moving outward (\(\dot{r} > 0\)) while rotating (\(\dot{\theta} > 0\)). As it moves to a larger radius, it must also increase its tangential speed \(r\dot{\theta}\) just to keep up with the same angular velocity. This “extra” angular acceleration is not due to a change in \(\dot{\theta}\) but due to the radius changing, and it appears twice (factor of 2) because \(\dot{r}\) affects both the transport of existing angular momentum outward and the change of the basis vector \(\boldsymbol{e}_\theta\) itself.
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|} \]
At any point along the path we can fit a circle that locally matches the curvature. The center of this circle is the center of curvature \(C\), and its radius is the radius of curvature \(R\). For an infinitesimal arc angle \(d\beta\), the arc length \(ds\) and the radius are related by
The circle arc length relation \(ds = R\,d\beta\), together with the Pythagorean theorem and trigonometry, forms the foundation of our kinematic modelling.
The speed along the path is \(v = |\dot{\boldsymbol{r}}| = \frac{ds}{dt}\), and using the arc length relation \(ds = R\,d\beta\) we get
\[ \dot{\boldsymbol{r}} = v\,\boldsymbol{e}_t = \frac{ds}{dt}\boldsymbol{e}_t = R\dot{\beta}\boldsymbol{e}_t \]
\[ \boxed{v = R\dot{\beta}} \]
This relation is useful when the angular rate \(\dot{\beta}\) is sought.
As seen in the figure above, a small change in the angle \(d\beta\) causes a corresponding change \(d\boldsymbol{e}_t\) perpendicular to \(\boldsymbol{e}_t\), that is, in the normal direction \(\boldsymbol{e}_n\). This gives
\[ d\boldsymbol{e}_t = \boldsymbol{e}_n\, d\beta \quad \Rightarrow \quad \dot{\boldsymbol{e}}_t = \dot{\beta}\,\boldsymbol{e}_n = \frac{v}{R}\boldsymbol{e}_n \tag{6.2.5}\]
From 6.2.5 we can also read off the radius of curvature and the angular rate,
\[ R = \frac{v}{|\dot{\boldsymbol{e}}_t|}, \qquad \dot{\beta} = \frac{v}{R} \]
The acceleration vector can be split into a tangential and a normal component, \(\ddot{\boldsymbol{r}} = a_t\boldsymbol{e}_t + a_n\boldsymbol{e}_n\). We find these by differentiating the velocity using the product rule,
\[ \boldsymbol{a} = \ddot{\boldsymbol{r}} = \frac{d}{dt}(v\,\boldsymbol{e}_t) = \dot{v}\,\boldsymbol{e}_t + v\,\dot{\boldsymbol{e}}_t \]
and with 6.2.5 we get
\[ \boxed{\ddot{\boldsymbol{r}} = \underbrace{\dot{v}\,\boldsymbol{e}_t}_{a_t} + \underbrace{\frac{v^2}{R}\boldsymbol{e}_n}_{a_n}} \]
The tangential component \(a_t = \dot{v}\) reflects changes in speed, while the normal component \(a_n = v^2/R\) is the centripetal acceleration directed toward the center of curvature. We always have a normal acceleration as long as the speed is nonzero and the path is curved (\(R < \infty\)). This decomposition will be revisited in the kinetics chapter.
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.
In the special case when \(R(t)\) is constant the bases of polar coordinates coincide with the bases of the natural coordinates.
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.
Remark 17.2. We compute the curvature and torsion for a simple circular path.
Let \(\boldsymbol{r}(t) = R [\cos(t), \sin(t), 0]^\mathsf{T}\), a circle of radius \(R\) in the \(xy\)-plane.
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)\[ \bm r =\left[\begin{matrix}R \cos{\left(t \right)}\\R \sin{\left(t \right)}\\0\end{matrix}\right] \]
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)\[ \kappa =\dfrac{1}{R} \]
# 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)\[ \tau =0 \]
The curvature is \(\kappa = 1/R\) as expected for a circle, and the torsion is zero since the path lies entirely in the \(xy\)-plane.
We illustrate the osculating circle on the path \(\boldsymbol{r}(t) = [t,\, \sqrt{t}\sin(2t)]^\mathsf{T}\). The key ingredients are the curvature \(\kappa\), the radius of curvature \(R = 1/\kappa\), the unit normal \(\boldsymbol{e}_n\), and the center of curvature \(\boldsymbol{r}_C = \boldsymbol{r} + R\,\boldsymbol{e}_n\).
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)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)\[ \bm r =\left[\begin{matrix}r{\left(t \right)} \cos{\left(\theta{\left(t \right)} \right)}\\r{\left(t \right)} \sin{\left(\theta{\left(t \right)} \right)}\end{matrix}\right] \]
rr_dot = rr.diff(t)
ltx(r"\dot{\bm r} =", rr_dot, arraystretch=2.5)\[ {\def\arraystretch{2.5}\dot{\bm r} =\left[\begin{matrix}- r{\left(t \right)} \sin{\left(\theta{\left(t \right)} \right)} \dfrac{d}{d t} \theta{\left(t \right)} + \cos{\left(\theta{\left(t \right)} \right)} \dfrac{d}{d t} r{\left(t \right)}\\r{\left(t \right)} \cos{\left(\theta{\left(t \right)} \right)} \dfrac{d}{d t} \theta{\left(t \right)} + \sin{\left(\theta{\left(t \right)} \right)} \dfrac{d}{d t} r{\left(t \right)}\end{matrix}\right]} \]
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)\[ {\def\arraystretch{1.5}\dot{\bm r} =\left[\begin{matrix}- \dot{\theta} r{\left(t \right)} \sin{\left(\theta{\left(t \right)} \right)} + \dot{r} \cos{\left(\theta{\left(t \right)} \right)}\\\dot{\theta} r{\left(t \right)} \cos{\left(\theta{\left(t \right)} \right)} + \dot{r} \sin{\left(\theta{\left(t \right)} \right)}\end{matrix}\right]} \]
The acceleration is
rr_ddot = rr.diff(t, 2).subs(dot_subs)
ltx(r"\ddot{\bm r} =", rr_ddot, arraystretch=1.5)\[ {\def\arraystretch{1.5}\ddot{\bm r} =\left[\begin{matrix}- \ddot{\theta} r{\left(t \right)} \sin{\left(\theta{\left(t \right)} \right)} + \ddot{r} \cos{\left(\theta{\left(t \right)} \right)} - \dot{\theta}^{2} r{\left(t \right)} \cos{\left(\theta{\left(t \right)} \right)} - 2 \dot{\theta} \dot{r} \sin{\left(\theta{\left(t \right)} \right)}\\\ddot{\theta} r{\left(t \right)} \cos{\left(\theta{\left(t \right)} \right)} + \ddot{r} \sin{\left(\theta{\left(t \right)} \right)} - \dot{\theta}^{2} r{\left(t \right)} \sin{\left(\theta{\left(t \right)} \right)} + 2 \dot{\theta} \dot{r} \cos{\left(\theta{\left(t \right)} \right)}\end{matrix}\right]} \]
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)\[ {\def\arraystretch{1.5}\dot{\bm r}_{r\theta} =\left[\begin{matrix}\dot{r}\\\dot{\theta} r{\left(t \right)}\end{matrix}\right]} \]
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)\[ {\def\arraystretch{1.5}\ddot{\bm r}_{r\theta} =\left[\begin{matrix}\ddot{r} - \dot{\theta}^{2} r{\left(t \right)}\\\ddot{\theta} r{\left(t \right)} + 2 \dot{\theta} \dot{r}\end{matrix}\right]} \]
We recover the same expressions derived by hand in Section 17.9.
For constant radius \(r(t) = R\) we set \(\dot{r} = \ddot{r} = 0\) and get the well-known expressions for circular motion.
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)\[ {\def\arraystretch{1.5}\dot{\bm r}_{r\theta} =\left[\begin{matrix}0\\\dot{\theta} r{\left(t \right)}\end{matrix}\right]} \]
ltx(r"\ddot{\bm r}_{r\theta} =", rr_ddot_polar.subs(circ_subs), arraystretch=1.5)\[ {\def\arraystretch{1.5}\ddot{\bm r}_{r\theta} =\left[\begin{matrix}- \dot{\theta}^{2} r{\left(t \right)}\\\ddot{\theta} r{\left(t \right)}\end{matrix}\right]} \]
The first component of the velocity is the radial velocity (zero for circular motion) and the second is the tangential velocity. For the acceleration, the first term is the centripetal acceleration, always pointing toward the center, and the second is the tangential acceleration, which vanishes when the angular velocity is constant.
To summarize planar circular motion: with a constant \(r\), we have \(\dot{r} = \ddot{r} = 0\) and \(\boldsymbol{r} \perp \dot{\boldsymbol{r}}\) since \(\frac{d}{dt}(r^2) = 0 = \frac{d}{dt}(\boldsymbol{r} \cdot \boldsymbol{r}) = \dot{\boldsymbol{r}} \cdot \boldsymbol{r} + \boldsymbol{r} \cdot \dot{\boldsymbol{r}} \Leftrightarrow \boldsymbol{r} \cdot \dot{\boldsymbol{r}} = 0\). Note especially the centripetal acceleration \(-r\dot{\theta}^2\boldsymbol{e}_r\), it is directed in the negative direction of \(\boldsymbol{e}_r\), towards the center of the circle.
If the rotational velocity is constant, i.e., if \(\omega = \dot{\theta} = \text{constant}\) or \(\alpha = \ddot{\theta} = 0\) then we only have centripetal acceleration, caused by a curvilinear motion and perpendicular to the tangent direction, directed inwards, towards the center.
\[ \begin{cases} \boldsymbol{r}(t) = r\,\boldsymbol{e}_r \\ \boldsymbol{v}(t) = r\dot{\theta}\boldsymbol{e}_\theta \\ \boldsymbol{a}(t) = -r\dot{\theta}^2\boldsymbol{e}_r + r\ddot{\theta}\boldsymbol{e}_\theta \end{cases} \]
especially when \(\omega = \dot{\theta} = \text{constant}\)
\[ \begin{cases} \boldsymbol{r}(t) = r\,\boldsymbol{e}_r \\ \boldsymbol{v}(t) = r\omega\boldsymbol{e}_\theta \\ \boldsymbol{a}(t) = -r\omega^2\boldsymbol{e}_r = \{v = |\boldsymbol{v}| = r\omega\} = -\frac{v^2}{r}\boldsymbol{e}_r \end{cases} \]
Note that we can also express \(\boldsymbol{e}_\theta\) using the cross product, \(\boldsymbol{e}_\theta = \boldsymbol{e}_z \times \boldsymbol{e}_r\), which is convenient since circular motion is common and we often have angular velocity as input data. This gives
\[ \boldsymbol{v}(t) = \dot{\theta}\,\boldsymbol{e}_z \times \boldsymbol{r} = r\dot{\theta}\,\boldsymbol{e}_\theta = r\omega\,\boldsymbol{e}_\theta \]
The advantage over explicitly rotating \(\boldsymbol{e}_r\) by 90 degrees is that the cross product automatically gives the correct direction of rotation from the sign of \(\omega\).
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} \]
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
the time for the particle to reach a velocity of 72 m/s from its initial condition at \(t = 0\),
the acceleration when \(v = 30\) m/s, and
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)\[ s(t) =2 t^{3} - 24 t + 6\\v(t) =6 t^{2} - 24\\ a(t) =12 t \]
(a) Find \(t\) when \(v = 72\) m/s.
t1 = sp.solve(sp.Eq(v, 72), t)
ltx(r"v(t) = 72 \Rightarrow t =", t1)\[ v(t) = 72 \Rightarrow t =\begin{bmatrix}4\end{bmatrix} \]
(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")\[ v(t) = 30 \Rightarrow t =3\\ a(t_2) =36\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}")\[ s(4) - s(1) =38--16=54\text{ m} \]
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)\[ {\def\arraystretch{2.5}\left[\begin{matrix}\dfrac{\sqrt{2} v_{B}}{2}\\\dfrac{\sqrt{2} v_{B}}{2}\end{matrix}\right] = \left[\begin{matrix}800 - \dfrac{v_{B/A}}{2}\\\dfrac{\sqrt{3} v_{B/A}}{2}\end{matrix}\right]} \]
sol = sp.solve(eq, (v_B, v_BA))
sol[v_B] | la.evalf(3)\[ {\def\arraystretch{1.0}717.0} \]
sol[v_BA].evalf(3) | la\[ {\def\arraystretch{1.0}586.0} \]
The Gripen’s true speed is 717 km/h, and it moves away from the passengers at 586 km/h.
A striking consequence of relative motion is observed in aviation: when two aircraft are on a collision course, the other aircraft appears stationary in the pilot’s field of view. The relative velocity vector points directly at the observer, so there is no angular motion, only a growing silhouette.
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)\[ \bm X =\begin{bmatrix}1.00 & 1.15 & 1.32 \\ 1.00 & 2.06 & 4.23 \\ 1.00 & 3.15 & 9.91 \\ 1.00 & 4.04 & 16.36 \\ 1.00 & 5.02 & 25.17 \\ 1.00 & 6.26 & 39.13 \\ \vdots & \vdots & \vdots \\ 1.00 & 10.18 & 103.53\end{bmatrix}_{10 \times 3} \]
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)\[ {\def\arraystretch{1.0}\bm c =\begin{bmatrix}1.06 \\ 7.75 \\ -0.31\end{bmatrix}_{3 \times 1}} \]
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)\[ {\def\arraystretch{1.0}\begin{bmatrix} 1.06 \\ 7.75 \\ -0.31 \end{bmatrix}_{3 \times 1}} \]
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)\[ v(t) =- 0.31 t^{2} + 7.75 t + 1.06 \]
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)\[ a(t) =7.75 - 0.62 t \]
# 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_{tot} = \int_0^{10.175} v(t)\, dt =303.4 \]
s = sp.integrate(v)
ltx(r"s(t) =", s, precision=3)\[ s(t) =- 0.103 t^{3} + 3.88 t^{2} + 1.06 t \]
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()
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\[ {\def\arraystretch{1.0}\dfrac{d^{2}}{d t^{2}} x{\left(t \right)} = - k^{2} x{\left(t \right)}} \]
x_sol = sp.simplify(sp.dsolve(DE, ics=ICs).rhs)
ltx(r"x(t) =", x_sol)\[ x(t) =\dfrac{v_{0} \sin{\left(k t \right)}}{k} \]
v_sol = x_sol.diff(t)
a_sol = v_sol.diff(t)
ltx(r"v(t) =", v_sol)\[ v(t) =v_{0} \cos{\left(k t \right)} \]
ltx(r"a(t) =", a_sol)\[ a(t) =- k v_{0} \sin{\left(k t \right)} \]
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)\[ x(t) =\dfrac{\sin{\left(4 t \right)}}{2},\quad v(t) =2 \cos{\left(4 t \right)},\quad a(t) =- 8 \sin{\left(4 t \right)} \]
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)\[ \int_0^{2\pi} v\, dt =0,\qquad \int_0^{2\pi} |v|\, dt =2 \]
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()
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(t) =\begin{cases} t^{2} & \text{for}\: t < 10 \\20 t - 100 & \text{otherwise} \end{cases} \]
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)\[ \dot x(t) =\begin{cases} 2 t & \text{for}\: t < 10 \\20 & \text{otherwise} \end{cases},\quad\ddot x(t) =\begin{cases} 2 & \text{for}\: t < 10 \\0 & \text{otherwise} \end{cases} \]
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()
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.
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)\[ x(t) =80 t - 800 + 800 e^{- \dfrac{t}{10}},\quad v(t) =80 - 80 e^{- \dfrac{t}{10}} \]
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()
How long to reach top speed? Since \(v(t) = 80(1 - e^{-t/10})\) approaches 80 m/s asymptotically, we never truly reach it. In practice we can ask when \(v = 0.99 \cdot v_1\).
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}")\[ x(t_{99\%}) =2892.0\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}")\[ x_{0-100} =63.44\text{ m} \]
We approximate the derivatives using finite differences,
\[ \ddot{x} \approx \frac{\dot{x}_{i+1} - \dot{x}_i}{\Delta t} = k(v_1 - \dot{x}_i) \quad \Rightarrow \quad \dot{x}_{i+1} = \dot{x}_i + k(v_1 - \dot{x}_i)\Delta t \]
\[ \dot{x} \approx \frac{x_{i+1} - x_i}{\Delta t} \quad \Rightarrow \quad x_{i+1} = x_i + \dot{x}_{i+1} \Delta t \]
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] + dtx_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()
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)\[ \theta(t) =\dfrac{t^{3}}{50} + \dfrac{t}{5},\qquad r(t) =\dfrac{t^{2}}{25} + \dfrac{1}{5} \]
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)\[ {\def\arraystretch{2.0}\dot{\bm r}_{r\theta} =\left[\begin{matrix}\dfrac{2 t}{25}\\\dfrac{3 t^{4}}{1250} + \dfrac{t^{2}}{50} + \dfrac{1}{25}\end{matrix}\right],\quad \ddot{\bm r}_{r\theta} =\left[\begin{matrix}- \dfrac{9 t^{6}}{62500} - \dfrac{21 t^{4}}{12500} - \dfrac{4 t^{2}}{625} + \dfrac{9}{125}\\\dfrac{t \left(9 t^{2} + 35\right)}{625}\end{matrix}\right]} \]
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)