Table of contents

Preliminaries

We start with kinematics and the mathematical descriptions of a bodies position, velocity and acceleration. This is the foundation of dynamics, including kinetics.

Assume an universal coordinate system and stick to it! We can choose between cartesian, normal and polar coordinate systems.

⚠ Note
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!

Definition of particle kinematics

If a body has an insignificant rotation compared to the path its center of gravity travels, it is said to be a particle. In later chapters we will make an proper introduction to generalized Newtonian-Euler motion equations, but for now, suffice to say if either the rotational acceleration or the size of the body is small enough, then we can neglect the change of angular momentum

α=0 or Iˉ=0MG=Iˉα=0\alpha =0 \text{ or } \bar{I} = 0 \Rightarrow \sum M_G =\bar{I} \alpha =0

1: The flywheel must be treated as a rigid body; it is being accelerated and de-accelerated, thus it generates a rotational moment around its axis, by design, providing an energy reservoir to the system.

2: This marble can be treated as a particle since it rotational speed does not significantly influences its path.

Vector notation

The position of a bodies center of gravity as a function of time tt can be expressed using the position vector r(t)=[x(t),y(t),z(t)]T\bm r\left(t\right)={\left\lbrack x\left(t\right),y\left(t\right),z\left(t\right)\right\rbrack }^T where each component is a function of the same independent variable tt. The components are called coordinate functions and r(t)\bm r\left(t\right) is known as the parameter curve with the parameter tt. The parameter curve is also known as the curve or path.

Velocity

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

limΔt0ΔrΔt=drdt=:r˙=v \boxed{ \lim_{\Delta t\to 0} \dfrac{\Delta \bm r}{\Delta t}=\dfrac{d\bm r}{dt}=:\dot{\bm r} =\bm v }

As can be seen in the animation, the velocity vector r˙\dot{\bm r} at tt is always tangential to the path.

The dot notation in r˙\dot{\bm r} or in r¨¨¨¨˙=d9dt9r\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, r˙|\dot{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 given by the change of velocity over time, dr˙dt\frac{d\dot{r} }{dt}. The average acceleration at time tt can be formulated as

limΔt0Δr˙Δt=dr˙dt=d2rdt2=:r¨=a \boxed{ \lim_{\Delta t\to 0} \dfrac{\Delta \dot{\bm r} }{\Delta t}=\dfrac{d\dot{\bm r} }{dt} = \dfrac{d^2 \bm r}{d t^2}=:\ddot{\bm r} =\bm a }

From the animation we can see that the acceleration is pointing inwards, towards the center of curvature. Also, we see that the acceleration in general differs from the direction of the position vector. The acceleration is in fact orthogonal to the velocity as we shall see in the following sections.

In the following sections, we will analyze curvilinear motion in more detail, but next we'll discuss the special case of no curvature...

Rectilinear motion

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

we sometimes use the (simple) scalar notation

s=r,v=r˙  and  a=r¨s=|\bm r|,v=|\dot{\bm r} |\;\textrm{and}\; a=|\ddot{\bm r} |

We can express these scalar values using the Leibniz notation for time derivative

v=dsdt  and  a=dvdtv=\dfrac{ds}{dt}\;\textrm{and}\;a=\dfrac{dv}{dt}

solving for the time increment dtdt, we get

dt=dsv=dvadt=\dfrac{ds}{v}=\dfrac{dv}{a}

from which we get

a  ds=v  dv \boxed{ a\;ds=v\;dv }

Which is convenient to use in special cases where we are not interested in time explicitly, but should be used with care and really be avoided all together, instead the proposed working method is described in what follows.

Differential equations

Working purely with the scalar notation from the previous section is troublesome. A strong recommendation is to avoid the use of derived formulas often found in classical books in mechanics and formularies, since it is often more tedious to figure out the assumptions under which these formulas are valid than to create a model using the differential equations directly. Many ready-to-go formulas assume the acceleration to be constant. We shall also refrain from using the notations s,v  and  as,v\;\textrm{and}\;a for varying position, velocity and acceleration since it is harder to let go of the "thinking by formularies" and move towards "Concept Based Modeling" with "Computational Thinking". Instead use r,r˙  and  r¨\bm r,\dot{\bm r} \;\textrm{and}\;\ddot{\bm r} along with knowledge of differential equations!

Remember that the derivative and integral operator works naturally on vectors i.e., r˙(t)=[x˙(t),y˙(t),z˙(t)]T\dot{\bm r} \left(t\right)={\left\lbrack \dot{x} \left(t\right),\dot{y} \left(t\right),\dot{z} \left(t\right)\right\rbrack }^T and r(t)  dt=[x(t)  dt,y(t)  dt,z(t)  dt]T\int \bm r\left(t\right)\;dt={\left\lbrack \int x\left(t\right)\;dt,\int y\left(t\right)\;dt,\int z\left(t\right)\;dt\right\rbrack }^T.

⚠ Note
The velocity (vector) is the time derivative of the position (vector), r˙=d  rdt\dot{\bm r} =\frac{d\;\bm r}{dt} with respect to time, tt. The magnitude r˙|\dot{\bm r} | is known as speed (scalar value). The direction is tangent to the trajectory path.

The acceleration (vector) is the time derivative of the velocity (vector) r¨=d  r˙d  t=ddt(d  rdt)=d2rdt2\ddot{\bm r} =\frac{d\;\dot{r} }{d\;t}=\frac{d}{dt}\left(\frac{d\;\bm r}{dt}\right)=\frac{d^2 \bm r}{ dt^2 }, or the second time derivative of the position (vector).

In what follows we shall rely more on working with the vector form of r\bm r and formulate the differentials as ordinary differential equations which we can easy solve using either a symbolic differential equation solve, e.g., dsolve() or numerically by utilizing e.g., the Euler method.

The modern workflow in computational dynamics can be summarized in the figure below

To get between the different quantities position, r(t)\bm r\left(t\right), velocity r˙(t)\dot{\bm r} \left(t\right) and acceleration r¨(t)\ddot{\bm r} \left(t\right) we take the derivative in one direction and integrate in the other. In practice one does not typically manually integrate expressions, one typically solves an ordinary differential equation (ODE) instead, which takes us all the way down to r(t)\bm r\left(t\right) which is the solution to an ODE. Then we can take the derivative to determine the needed quantities. This is typical for a symbolic solution. Non-linear ODEs, on the other hand, are typically solved using various numerical algorithms which will first generate velocity, r˙\dot{\bm r} on the way down to the final position r\bm r.

Newtonian physics is formulated as the kinetic equation

F=mr¨ \sum \bm F = m \ddot{\bm r}

which is a second order ODE, but in kinematics we se equally often differential equations of first order. Regardless, it is very convenient to solve any equation using matlab.

The computational approach to dynamics is to formulate an initial value problem (IVP) and solve it using a differential equation solver.

The inquiries to the model can be essentially divided into two categories:

  1. Explicit time dependent, the equations are explicitly functions of a time variable, e.g., r˙(t)\dot{\bm r} \left(t\right), such that we can evaluate the quantities for a given time.

  2. Implicit time dependent, where we want to know a quantity expressed as a function of another quantity and thus we need to solve an equation to get the time which is common for both quantities. E.g., find the velocity when the position is a given value, r˙(r)\dot{r}(r).

Coordinate systems

There are several types of bases (or coordinate systems) on which we can formulate kinematic expressions and work with derivatives, traditionally this is done to simplify calculations, done by hand, and the resulting symbolic expressions.

⚠ Note
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 r\bm 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 base vectors such that this projection can be done.

Sometimes the bases can be denoted simply i\bm i, j\bm j and k\bm k, but throughout this treatment we will try to be explicit and denote en basis vectors and unit vectors with a specific sub-index for clarity.

Natural coordinates

The tangent vector et{\bm e}_t is always pointing in the parameter direction, tangential to the path and is defined as the direction of the velocity

et=r˙r˙ \boxed{ {\bm e}_t =\frac{\dot{\bm r} }{|\dot{\bm r} |} }

The normal vector en{\bm e}_n is perpendicular to et{\bm e}_t and always points inwards to the center of curvature

en=et˙et˙ \boxed{\bm{e}_{n}=\dfrac{\dot{\bm{e}_{t}}}{|\dot{\bm{e}_{t}}|}}

here, suddenly a time derivate of a unit vector appears, see this section for an explanation.

Finally the lesser known bi-normal vector is defined as the cross product between the tangent and normal vectors,

eb:=et×en \boxed{ {\bm e}_b :={\bm e}_t \times {\bm e}_n }

which defines a right-hand system.

Note that for planar problems (2D) eb=ez{\bm e}_b ={\bm e}_z.

Polar coordinates

er{\bm e}_r has the same direction as r(t)\bm r\left(t\right) and is given by

er=r(t)r(t) \boxed{ {\bm e}_r =\frac{\bm r\left(t\right)}{|\bm r\left(t\right)|} }

eθ{\bm e}_{\theta } is given by

eθ=e˙re˙r \boxed{ {\bm e}_{\theta } =\frac{ {\dot{\bm e} }_r }{|{\dot{\bm e} }_r |} }

The time derivative of a unit vector

In order to ensure orthogonality of two vectors, we can take their dot product, if the resulting scalar is zero, they must be orthogonal. If two vectors are parallel then their dot product must be one, i.e., ee=1\bm e\cdot \bm e=1. Now we can take the time derivative of this expression ddt(ee=1)=product  ruleddtee+eddte=0=e˙e+ee˙=0e˙e=0ee˙\frac{d}{dt}\left(\bm e\cdot \bm e=1\right)\overset{\textrm{product}\;\textrm{rule}}{=} \frac{d}{dt}\bm e\cdot \bm e+\bm e\cdot \frac{d}{dt}\bm e=0=\dot{\bm e} \cdot \bm e+\bm e\cdot \dot{\bm e} =0\Leftrightarrow \dot{\bm e} \cdot \bm e=0\Rightarrow \bm e\bot \dot{\bm e}

Also note that typically e˙1|\dot{\bm e} |\not= 1.

So at the limit Δt0\Delta t\to 0 we clearly see that ee˙\bm e\bot \dot{\bm e}.

Let r(t)=[cosθ(t),sinθ(t)]\bm r\left(t\right)=\left\lbrack \cos \theta \left(t\right),\sin \theta \left(t\right)\right\rbrack

normalize = @(v)v/sqrt(v(:).'*v(:));
syms t positive
syms theta(t)
r = [cos(theta(t)),sin(theta(t))].'
(cos(θ(t))sin(θ(t))) \left(\begin{array}{c} \cos \left(\theta \left(t\right)\right)\\ \sin \left(\theta \left(t\right)\right) \end{array}\right)

We compute the normalized vector er{\bm e}_r

e_r = simplify(normalize(r))
(cos(θ(t))sin(θ(t))) \left(\begin{array}{c} \cos \left(\theta \left(t\right)\right)\\ \sin \left(\theta \left(t\right)\right) \end{array}\right)

The time derivative of er{\bm e}_r is then given by er˙\dot{ {\bm e}_{r} }

e_r_dot = simplify(diff(e_r,t))
(sin(θ(t))t  θ(t)cos(θ(t))t  θ(t)) \def\arraystretch{2.5} \left(\begin{array}{c} -\sin \left(\theta \left(t\right)\right)\,\dfrac{\partial }{\partial t}\;\theta \left(t\right)\\ \cos \left(\theta \left(t\right)\right)\,\dfrac{\partial }{\partial t}\;\theta \left(t\right) \end{array}\right)

What we computed above is a scaled direction which is perpendicular to the direction vector er{\bm e}_r. Note that the time derivative of θ\theta pops out from the chain rule, which is why e˙1|\dot{e} |\not= 1 necessarily since θ˙1|\dot{\theta} |\not= 1 necessarily. Note that the direction of rotation is given by the sign of θ˙\dot{\theta} i.e., {-1,1}.

The above we re-write in Newton notation and define the basis vector as eθ\bm e_\theta

e˙r=θ˙eθ \boxed{ \dot{\bm e}_r = \dot \theta \bm e_{\theta} }

Thus we can get the perpendicular vector, eθ\bm e_\theta, by normalizing e˙\dot{e}, i.e.,

eθ=e˙re˙r \boxed{ \bm e_{\theta } =\frac{\dot{\bm e}_r }{| \dot{\bm e}_r |} }
e_theta = subs(simplify(normalize(e_r_dot)), diff(theta,t), 'theta_dot')
(θ˙sin(θ(t))θ˙2θ˙cos(θ(t))θ˙2) \def\arraystretch{2.5} \left(\begin{array}{c} -\dfrac{\dot{\theta} \,\sin \left(\theta \left(t\right)\right)}{\sqrt{ {\dot{\theta} }^2 } }\\ \dfrac{\dot{\theta} \,\cos \left(\theta \left(t\right)\right)}{\sqrt{ {\dot{\theta} }^2 } } \end{array}\right)

Here we can see that θ˙θ˙2\frac{\dot{\theta} }{\sqrt{ {\dot{\theta} }^2 }} only returns the sign of θ˙\dot{\theta}. Thus in this example eθ=[sinθ,cosθ]Te_{\theta } ={\left\lbrack -\sin \theta ,\cos \theta \right\rbrack }^T and we get

ereθ=[cosθ,sinθ][sinθ,cosθ]=cosθsinθ+sinθcosθ=0 \bm e_{r} \cdot \bm e_{\theta} = [\cos \theta, \sin\theta] \cdot [-\sin \theta, \cos\theta] = -\cos\theta\sin\theta + \sin\theta\cos\theta = 0

Thus we have shown that

ere˙r,and ereθ \boxed{\bm{e}_{r}\perp\dot{\bm{e}}_{r},\text{and }\bm{e}_{r}\perp\bm{e}_{\theta}}

Polar coordinates (rθ)\left(r-\theta \right)

We shall here derive the kinematic equations explicitly in polar coordinates. Note again that in practice, working with a symbolic math manipulator like matlab, we can just formulate r˙\dot{\bm r} and r¨\ddot{\bm r} and project onto the base. There is no need to learn the derived expressions below by heart in a modern workflow, instead we explore the connections between these entities and their physical meaning.

We begin by stating the earlier established bases.

er=rr{\bm e}_r =\dfrac{\bm r}{|\bm r|} eθ=e˙re˙r{\bm e}_{\theta } =\dfrac{\dot{\bm e}_r }{|{\dot{\bm e} }_r |}

We can express the position vector using the distance and base

r=rer,\bm r=r{\bm e}_r,

where rr is the scalar length of the position vector.

Velocity in Polar Coordinates

Taking the derivative of the position vector we get, using the product rule

r˙=ddt(r(t)e(t))=r˙er+re˙r\dot{\bm r} =\dfrac{d}{dt}\left(r(t) \bm e(t)\right)=\dot{r} {\bm e}_r +r{\dot{\bm e} }_r

and with e˙r=θ˙eθ{\dot{\bm e} }_r =\dot{\theta} {\bm e}_{\theta } from (17) we have

r˙=r˙er+rθ˙eθ\boxed{ \dot{\bm r} =\dot{r} {\bm e}_r +r\dot{\theta} {\bm e}_{\theta } }

where

ω=rθ˙\boxed{ \omega = r\dot{\theta} }

is recognized as the angular velocity.

Acceleration in Polar Coordinates

Taking the time derivative of the expression for velocity above we get

r¨=ddt(r˙)=(r¨er+r˙e˙r)+(r˙θ˙eθ+rθ¨eθ+rθ˙e˙θ)\ddot{\bm r} =\dfrac{d}{dt}\left(\dot{\bm r} \right)=\left(\ddot{r} {\bm e}_r +\dot{r} {\dot{\bm e} }_r \right)+\left(\dot{r} \dot{\theta} {\bm e}_{\theta } +r\ddot{\theta} {\bm e}_{\theta } +r\dot{\theta} {\dot{\bm e} }_{\theta } \right)

and with e˙r=θ˙eθ{\dot{e} }_r =\dot{\theta} {\bm e}_{\theta } and e˙θ=θ˙er{\dot{\bm e} }_{\theta } =-\dot{\theta} {\bm e}_r we have the expression of the acceleration vector expressed in polar coordinates

r¨=(r¨rθ˙2)er+(rθ¨+2r˙θ˙)eθ \boxed{ \ddot{\bm r} =\left(\ddot{r} -r{\dot{\theta} }^2 \right){\bm e}_r +\left(r\ddot{\theta} +2\dot{r} \dot{\theta} \right){\bm e}_{\theta } }

Natural coordinates (tn)\left(t-n\right)

Arguably more important are the kinematic relation expressed in Natural coordinates and even more useful are the acceleration terms which are split into normal and tangential directions. Here we shall derive these expressions explicitly, but note that we can always get these by just projecting the total velocity or acceleration vector onto any bases. Thus, there is no need to know these formulas or use them only for the sake of computing the quantities in one of these special bases.

We begin by stating the earlier established bases.

et=r˙(t)r˙(t){\bm e}_t =\dfrac{\dot{\bm r} \left(t\right)}{|\dot{\bm r} \left(t\right)|} en=e˙te˙t {\bm e}_n =\dfrac{ \dot{\bm e}_t }{|{\dot{\bm e} }_t |}

Curvature

Next, we introduce the center of curvature, CC, at some point r(t)\bm r\left(t\right) and the corresponding curvature radius RR. We then can fit a circle arc length in the neighborhood of r(t)\bm r\left(t\right) such that we get a constant radius RR for an infinitesimal circle arc angle dβd\beta . This way we can get a relation between the section distance dsds and the angle

ds=R  dβds=R\;d\beta

This circle arc length relation, along with the Pythagorean theorem and trigonometry is a corner stone in dynamics and the main building blocks in our modelling toolbox.

Velocity in Natural coordinates

We remember the definition of velocity as v=dsdtv=\frac{ds}{dt} and using the arc length relation we get

r˙=r˙et=dsdt  et=Rβ˙et\dot{\bm r} =\dot{r} {\bm e}_t =\dfrac{ds}{dt}\;{\bm e}_t = R \dot{\beta} {\bm e}_t r˙=Rβ˙ \boxed{ \dot{r} = R \dot{\beta} }

This relation is useful for when the angular rate β˙\dot{\beta} is sought.

For a small change detd\bm e_t by dβd\beta we get the perpendicular direction det{d\bm e}_t which has the same direction as en{\bm e}_{\bm n}, we can thus formulate the relation

det=endβdetdt=dβdtenet˙=β˙en=r˙Rend\bm e_t = \bm e_n d\beta \Rightarrow \dfrac{d\bm e_t }{dt}=\dfrac{d\beta }{dt}{\bm e}_n \Rightarrow \dot{\bm e_t } =\dot{\beta} {\bm e}_n =\dfrac{\dot{r} }{R}{\bm e}_n

Thus we arrive at the relation

et˙=r˙Ren \dot{ {\bm e}_t } =\dfrac{\dot{r} }{R}{\bm e}_n R=r˙et˙R=\dfrac{\dot{r} }{|\dot{ {\bm e}_t } |}

Since r˙=Rβ˙|\dot{\bm r} |=R\dot{\beta} we get

β˙=r˙R=r˙2et˙\dot{\beta} =\dfrac{\dot{r} }{R}=\dfrac{ {\dot{r} }^2 }{|\dot{ {\bm e}_t } |}

Acceleration in Natural coordinates

a=r¨=ddt(r˙et)=dr˙dtet+r˙d  etdt\bm a=\ddot{\bm r} =\dfrac{d}{dt}\left(\dot{r} {\bm e}_t \right)=\dfrac{d\dot{r} }{dt}{\bm e}_t +\dot{r} \dfrac{ {d\;\bm e}_t }{dt} r¨=r¨et+r˙e˙t\ddot{\bm r} =\ddot{r} {\bm e}_t +\dot{r} {\dot{\bm e} }_t

and with (36) we get

r¨=r˙2Renan+r¨etat\ddot{\bm r} =\underset{ {\bm a}_n }{\underbrace{\dfrac{ {\dot{r} }^2 }{R}{\bm e}_n } } +\underset{ {\bm a}_t }{\underbrace{\ddot{r} {\bm e}_t } }

From this expression we note that we always have an acceleration as long as there is non-zero velocity and a curved path.

This split into the normal and tangential directions will be revisited and used once we get to the kinetics chapter.

A tiny view on differential geometry

The theory on curves is extensive and belongs to the field of mathematics called differential geometry. Here we shall only introduce som very basic concepts, such as the curve length, for which there might not exist a closed form solution, it needs to be integrated over all small contributions r˙(t)|\dot{\bm r} \left(t\right)| such as

s(t)=0tr˙(t)dts\left(t\right)=\int_0^t |\dot{\bm r} \left(t\right)|dt

This can be useful if the derivatives are needed to be expressed as dsds instead of dtdt which can lead to simpler expressions, but since ss is a function of tt it does not make the practical computations easier when working computer based.

Another important concept is the curvature κ(t)=1R(t)\kappa \left(t\right)=\frac{1}{R\left(t\right)}, which as we can see tends to zero as the curvature radius R(t)R\left(t\right) tends to infinity, which makes sense for a flat path. The curvature center rC(t){\bm r}_C \left(t\right) is given by

rC(t)=r(t)+R(t)en(t){\bm r}_C \left(t\right)=\bm r\left(t\right)+R\left(t\right){\bm e}_n \left(t\right)

This center-point lies in the center of what is known as the osculating circle, named by Leibniz (Circulus Osculans), it is the best fitting circle to the curve at the time tt. We note that the Natural coordinates are mostly used for this type of analysis.

⚠ Note
In the special case when R(t)R\left(t\right) is constant the bases of polar coordinates coincide with the bases of the natural coordinates.

The torsion of a curve, τ(t)\tau \left(t\right) is a measure of how much a curve is being bent into the third dimension (compared to the curvature). Think of it like the pitch of a screw. Both the expressions for curvature and torsion are tedious to derive, hence we just state them here

κ(t)=r˙(t)×r¨(t)r˙(t)3\kappa \left(t\right)=\dfrac{|\dot{\bm r} \left(t\right) \times \ddot{\bm r} \left(t\right)|}{|\dot{\bm r} \left(t\right)|^3 } τ(t)=(r˙(t)×r¨(t))r¨˙(t)r˙(t)×r¨(t)2\tau \left(t\right)=-\dfrac{\left(\dot{\bm r} \left(t\right) \times \ddot{\bm r} \left(t\right)\right)\cdot \dot{\ddot{\bm r}} \left(t\right)}{|\dot{\bm r} \left(t\right)\times \ddot{\bm r} \left(t\right)|^2 }

where r¨˙\dot{\ddot{\bm r}} denotes the third order derivative of time, also known as jerk or jolt, for higher order time derivatives see this list.

Let us work out the curvature and torsion on a simple example, let

r(t)=R  [cos(t),sin(t),0]Tr\left(t\right)=R\;{\left\lbrack \cos \left(t\right),\sin \left(t\right),0\right\rbrack }^T
clear
syms R t positive
r = R*[cos(t), sin(t), 0].'
(Rcos(t)Rsin(t)0) \left(\begin{array}{c} R\,\cos \left(t\right)\\ R\,\sin \left(t\right)\\ 0 \end{array}\right)
r_dot = @(r) diff(r,t);
r_ddot = @(r) diff(r,t,2);
r_dddot = @(r) diff(r,t,3);

kappa = @(r) simplify( norm(cross(r_dot(r),r_ddot(r))) / norm(r_dot(r))^3 );
tau = @(r) - dot(cross(r_dot(r),r_ddot(r)), r_dddot(r)) / norm( cross(r_dot(r),r_ddot(r)) )^2;
kappa(r)
1R \dfrac{1}{R}
tau(r)
00

The results seem reasonable, torsion is zero since the path is zero in the zz-direction.

Osculating circle

Finally we visualize the osculating circle on its journey accompanying a particle moving along the path x(t)=t,  y(t)=tsin(2t)x\left(t\right)=t,\;y\left(t\right)=\sqrt{t}\sin \left(2t\right)

clear
syms t positive
y = 1*sin(2*t)*t^0.5;
y = matlabFunction(y);

r = [t,sym(y),0]
(ttsin(2t)0) \left(\begin{array}{ccc} t & \sqrt{t}\,\sin \left(2\,t\right) & 0 \end{array}\right)
r_dot = diff(r,t);
r_ddot = diff(r,t,2);
e_t = matlabFunction( simplify(r_dot/norm(r_dot)) );
e_n = matlabFunction( diff(e_t, t)/norm(diff(e_t, t)) );
e_t_dot = diff(r_dot/norm(r_dot),t);
e_t_dot_norm = matlabFunction(simplify(norm(e_t_dot)) );
r_dot = matlabFunction(r_dot);
r_ddot = matlabFunction(r_ddot);

figure
fplot(y,[0,10], 'Linewidth',2)
axis equal
xlabel('$x$','Interpreter','latex')
ylabel('$y$','Interpreter','latex')
c_bg = [33,33,33]/255;
set(gcf,'Color',c_bg)
set(gca,'FontSize',14,'FontName','Times','Color',c_bg,'XColor','w','YColor','w','ZColor','w')
hold on
box off

t = 1.2;
x = t;

e_t_N = e_t(t);
e_n_N = e_n(t);
c_green = [0, 176, 80]/255;
c_blue = [0,176,240]/255;

v = r_dot(t)
v = 1x3    
    1.0000   -1.3072         0
kappa = norm(cross(r_dot(t),r_ddot(t))) / norm(r_dot(t))^3
kappa = 0.9946
R = double( norm(v)/e_t_dot_norm(t) )
R = 1.0054
R = 1/kappa
R = 1.0054
r_N = [t, y(t),0];

r_C_N = r_N + e_n_N*R;

hc = plot(r_C_N(1), r_C_N(2),'ro','MarkerFaceColor','r','MarkerSize',3);

tt = linspace(0,2*pi,100)';
circle = @(R,r)R*[cos(tt),sin(tt)]+[r(1), r(2)]
circle = 
    @(R,r)R*[cos(tt),sin(tt)]+[r(1),r(2)]
C = circle(R,r_C_N);
hcircle = plot(C(:,1), C(:,2),'r-');
hp = plot(x,y(t),'ko','MarkerFaceColor','k','MarkerSize',8);
het = quiver(x,y(t), e_t_N(1), e_t_N(2), 1, 'Color',c_green, 'Linewidth',2);
hen = quiver(x,y(t), e_n_N(1), e_n_N(2), 1, 'Color',c_blue, 'Linewidth',2);
axis([0,10,-3,3])
% 
% for t = linspace(0,10,20)
%     x = t;
%     e_t_N = e_t(t);
%     e_n_N = e_n(t);
%     kappa = norm(cross(r_dot(t),r_ddot(t))) / norm(r_dot(t))^3;
%     R = 1/kappa;
%     r_N = [t, y(t),0];
%     r_C_N = r_N + e_n_N*R;
%     C = circle(R,r_C_N);
%     set(hp, 'XData', x, 'Ydata', y(t))
%     set(het,'XData', x, 'Ydata', y(t),'UData',e_t_N(1),'VData',e_t_N(2))
%     set(hen,'XData', x, 'Ydata', y(t),'UData',e_n_N(1),'VData',e_n_N(2))
%     set(hc, 'XData', r_C_N(1), 'Ydata', r_C_N(2))
%     set(hcircle, 'XData', C(:,1), 'Ydata', C(:,2))
%     axis([0,10,-3,3])
%     drawnow
% end

Circular motion

A special case of curvilinear motion is planar circular motion.

Let us study the case of a particle in polar coordinates using the arm below

We begin by expressing the particles position

clear

syms theta(t) r(t) 

rr = r*[cos(theta), sin(theta)].'
(cos(θ(t))r(t)sin(θ(t))r(t)) \left(\begin{array}{c} \cos \left(\theta \left(t\right)\right)\,r\left(t\right)\\ \sin \left(\theta \left(t\right)\right)\,r\left(t\right) \end{array}\right)

Then we let matlab compute the derivatives

rr_dot = diff(rr,t);
syms r_dot theta_dot r_ddot theta_ddot
rr_dot = subs(rr_dot, [diff(r,t), diff(theta,t)], [r_dot,theta_dot])
(r˙cos(θ(t))θ˙sin(θ(t))r(t)r˙sin(θ(t))+θ˙cos(θ(t))r(t)) \left(\begin{array}{c} \dot{r} \,\cos \left(\theta \left(t\right)\right)-\dot{\theta} \,\sin \left(\theta \left(t\right)\right)\,r\left(t\right)\\ \dot{r} \,\sin \left(\theta \left(t\right)\right)+\dot{\theta} \,\cos \left(\theta \left(t\right)\right)\,r\left(t\right) \end{array}\right)
rr_ddot = diff(rr,t,2);
rr_ddot = subs(rr_ddot, [diff(r,t,2), diff(theta,t,2)], ...
                        [r_ddot, theta_ddot]);
rr_ddot = subs(rr_ddot, [diff(r,t), diff(theta,t)], ...
                        [r_dot, theta_dot])
(cos(θ(t))r(t)θ˙22r˙sin(θ(t))θ˙+r¨cos(θ(t))θ¨sin(θ(t))r(t)sin(θ(t))r(t)θ˙2+2r˙cos(θ(t))θ˙+r¨sin(θ(t))+θ¨cos(θ(t))r(t)) \left(\begin{array}{c} -\cos \left(\theta \left(t\right)\right)\,r\left(t\right)\,{\dot{\theta} }^2 -2\,\dot{r} \,\sin \left(\theta \left(t\right)\right)\,\dot{\theta} +\ddot{r} \,\cos \left(\theta \left(t\right)\right)-\ddot{\theta} \,\sin \left(\theta \left(t\right)\right)\,r\left(t\right)\\ -\sin \left(\theta \left(t\right)\right)\,r\left(t\right)\,{\dot{\theta} }^2 +2\,\dot{r} \,\cos \left(\theta \left(t\right)\right)\,\dot{\theta} +\ddot{r} \,\sin \left(\theta \left(t\right)\right)+\ddot{\theta} \,\cos \left(\theta \left(t\right)\right)\,r\left(t\right) \end{array}\right)

In the above expressions, we made the derivatives easier to read by converting to Newtonian notation. Clearly, Matlab can handle tedious derivation.

Now, we define the polar bases

e_r = [cos(theta), sin(theta)].';
e_theta = [-sin(theta), cos(theta)].';

where eθe_{\theta } is orthogonal according to the right-hand rule (rotate 90 degrees counter clockwise).

Now just project the Cartesian vectors onto the polar bases. We can get away with the "low budget" projection since the bases have unit length, i.e., er=eθ=1|{\bm e}_r |=|{\bm e}_{\theta } |=1. We recall that the dot product is just a matrix multiplication

[a1,a2][b1b2]=a1b1+a2b2\left\lbrack a_1 ,a_2 \right\rbrack \left\lbrack \begin{array}{c} b_1 \\ b_2 \end{array}\right\rbrack =a_1 b_1 +a_2 b_2
r_dot_rtheta = simplify([rr_dot.'*e_r, rr_dot.'*e_theta])
(r˙θ˙r(t)) \left(\begin{array}{cc} \dot{r} & \dot{\theta} \,r\left(t\right) \end{array}\right)
r_ddot_rtheta = simplify([rr_ddot.'*e_r, rr_ddot.'*e_theta])
(r¨θ˙2r(t)2r˙θ˙+θ¨r(t)) \left(\begin{array}{cc} \ddot{r} -{\dot{\theta} }^2 \,r\left(t\right) & 2\,\dot{r} \,\dot{\theta} +\ddot{\theta} \,r\left(t\right) \end{array}\right)

We can see that we end up with the same expressions as we have tediously derived using the classical approach in the section on polar coordinates above.

For the special case of circular motion, we set the derivatives of rr to zero

r_dot = 0; r_ddot = 0;
v = subs(r_dot_rtheta)
(0θ˙r(t)) \left(\begin{array}{cc} 0 & \dot{\theta} \,r\left(t\right) \end{array}\right)
a = subs(r_ddot_rtheta)
(θ˙2r(t)θ¨r(t)) \left(\begin{array}{cc} -{\dot{\theta} }^2 \,r\left(t\right) & \ddot{\theta} \,r\left(t\right) \end{array}\right)

To summarize planar circular motion: with a constant rr, we have r˙=r¨=0\dot{r} =\ddot{r} =0 and rr˙\bm r\bot \dot{\bm r} since ddt(r2)=0=ddt(rr)=r˙r+rr˙rr˙=0\frac{d}{dt}\left(r^2 \right)=0=\frac{d}{dt}\left(\bm r\cdot \bm r\right)=\dot{\bm r} \cdot \bm r+\bm r\cdot \dot{\bm r} \Leftrightarrow \bm r\cdot \dot{\bm r} =0. Note especially the centripetal acceleration rθ˙2er-r{\dot{\theta} }^2 {\bm e}_r, it is directed in the negative direction of er{\bm e}_r, towards the center of the circle.

If the rotational velocity is constant, i.e., if ω=θ˙=constant\omega =\dot{\theta} = \text{constant} or α=θ¨=0\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.

{r(t)=r  erv(t)=rθ˙eθa(t)=rθ˙2er+rθ¨eθ \left\lbrace \begin{array}{ll} \bm r\left(t\right)=r\;{\bm e}_r & \\ \bm v\left(t\right)=r\dot{\theta} {\bm e}_{\theta } & \\ \bm a\left(t\right)=-r{\dot{\theta} }^2 {\bm e}_r +r\ddot{\theta} {\bm e}_{\theta } & \end{array}\right. especially when   ω=θ˙=constant  {r(t)=r  erv(t)=rωeθa(t)=rω2er={v=v=rω}=v2rer \text{especially when } \; \omega =\dot{\theta} =\textrm{constant}\;\left\lbrace \begin{array}{ll} \bm r\left(t\right)=r\;{\bm e}_r & \\ \bm v\left(t\right)=r\omega {\bm e}_{\theta } & \\ \bm a\left(t\right)=-r\omega^2 {\bm e}_r =\left\lbrace v=|\bm v|=r\omega \right\rbrace =-\dfrac{v^2 }{r}{\bm e}_r & \end{array}\right.

Note that we can express bases using the cross product

eθ=ez×er{\bm e}_{\theta } ={\bm e}_z \times {\bm e}_r

which is very convenient since circular motion is common and we often have angular velocity as input data, thus we can formulate

v(t)=θ˙ez×r=rθ˙eθ=rωeθ\bm v\left(t\right)=\dot{\theta} {\bm e}_z \times \bm r=r\dot{\theta} {\bm e}_{\theta } =r\omega {\bm e}_{\theta }

Another easy way is to explicitly rotate the base er{\bm e}_r by 90 degrees into the direction of curvature, the drawback is we need to know which direction to rotate into, which is why the previous method is preferred, since the direction of rotation is determined by the sign of the angular velocity ω\omega.

Relative Motion

Many times we can express position and motion (velocities and acceleration) relative to some other point rather than the origin. This is done such that the modeling becomes significantly easier and less error prone. Besides, in real world applications, engineers need to make measurements and many times it is more accurate to make relative measurements rather than absolute. This approach is called relative motion analysis and will be used extensively in rigid body dynamics.

If we know the position of particle BB as well as the relative position of particle AA with respect to BB, that is rA/B{\bm r}_{A/B} (Can be read: A seen from B) we can express the absolute position of AA as

rA=rB+rA/B{\bm r}_A ={\bm r}_B +{\bm r}_{A/B}

Again it might be more convenient to express rA/B{\bm r}_{A/B} and rB{\bm r}_B than rA{\bm r}_A directly.

If the position vectors are expressed as variables of time, we can differentiate them with respect to time and get

ddt(rA=rB+rA/B)=r˙A=r˙B+r˙A/B      or      vA=vB+vA/B\dfrac{d}{dt}\left({\bm r}_A ={\bm r}_B +{\bm r}_{A/B} \right)={\dot{\bm r} }_A ={\dot{\bm r} }_B +{\dot{\bm r} }_{A/B} \;\;\;\textrm{or}\;\;\;{\bm v}_A ={\bm v}_B +{\bm v}_{A/B}

and

r¨A=r¨B+r¨A/B      or      aA=aB+aA/B{\ddot{\bm r} }_A ={\ddot{\bm r} }_B +{\ddot{\bm r} }_{A/B} \;\;\;\textrm{or}\;\;\;{\bm a}_A ={\bm a}_B +{\bm a}_{A/B}

Example

Passengers in the A380 flying east at a speed of 800 km/h observe a stunning JAS 39 Gripen passing underneath with a heading of 45{45}^{\circ } in the north-east direction, although to the passengers it appears that the JAS is moving away at an angle of 60{60}^{\circ } as shown in the figure. determine the true velocity of the JAS 39 Gripen.

Solution:

clear
syms v_B v_BA

v_A = 800
v_A = 800
V_A = v_A*[1;0]
V_A = 2x1    
   800
     0
e_B = [cosd(45); sind(45)]
e_B = 2x1    
    0.7071
    0.7071
V_B = v_B*e_B
(2vB22vB2) \left(\begin{array}{c} \dfrac{\sqrt{2}\,v_B }{2}\\ \dfrac{\sqrt{2}\,v_B }{2} \end{array}\right)
e_BA = [-cosd(60); sind(60)]
e_BA = 2x1    
   -0.5000
    0.8660
V_BA = v_BA*e_BA
(vBA23vBA2) \left(\begin{array}{c} -\dfrac{v_{\textrm{BA} } }{2}\\ \dfrac{\sqrt{3}\,v_{\textrm{BA} } }{2} \end{array}\right)

Now we can formulate the equation, which is good to write out, that way we can see the number of variables and the number of equations.

eqn = V_B == V_A + V_BA
(2vB2=800vBA22vB2=3vBA2) \left(\begin{array}{c} \dfrac{\sqrt{2}\,v_B }{2}=800-\dfrac{v_{\textrm{BA} } }{2}\\ \dfrac{\sqrt{2}\,v_B }{2}=\dfrac{\sqrt{3}\,v_{\textrm{BA} } }{2} \end{array}\right)

We then solve for the unknowns and convert to numerical values.

[v_B, v_BA] = solve(eqn,[v_B,v_BA])
16003+1 \dfrac{1600}{\sqrt{3}+1}
v_B = vpa(v_B,5)
717.26 717.26
v_BA = vpa(v_BA,5)
585.64 585.64

To conclude, the JAS is moving away from the passengers at 585 km/h and its absolute speed is 717 km/h.

Example 2 - 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 curve as well as the acceleration curve.

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:

Fit a quadratic function to the data, in general we have:

vi=c0+c1t+c2t2+cntn  i[1,m]v_i =c_0 +c_1 t+c_2 t^2 +\cdots c_n t^n ~~\forall i\in [1,m]

where m=10m=10 in our case. On matrix form this becomes:

[1t1t12t1n1t2t22t2n1tmtm2tmn]X[c0c1cm]c=[v1v2vm]v\underset{\mathbf{X}}{\underbrace{\left\lbrack \begin{array}{ccccc} 1 & t_1 & t_1^2 & \cdots & t_1^n \\ 1 & t_2 & t_2^2 & \cdots & t_2^n \\ \vdots & \vdots & \vdots & \cdots & \vdots \\ 1 & t_m & t_m^2 & \cdots & t_m^n \end{array}\right\rbrack } } \underset{\mathbf{c}}{\underbrace{\left\lbrack \begin{array}{c} c_0 \\ c_1 \\ \vdots \\ c_m \end{array}\right\rbrack } } =\underset{\mathbf{v}}{\underbrace{\left\lbrack \begin{array}{c} v_1 \\ v_2 \\ \vdots \\ v_m \end{array}\right\rbrack } }

Solve for the constants, cic_i, using the least squares method:

c=(XTX)1XTv\mathbf{c}={\left({\mathbf{X}}^T \mathbf{X}\right)}^{-1} {\mathbf{X}}^T \mathbf{v}
syms t
n = 2;
t.^[0:n]
(1tt2) \left(\begin{array}{ccc} 1 & t & t^2 \end{array}\right)
X = @(t)t.^[0:n];
X = X(T)
X = 10x3    
    1.0000    1.1486    1.3193
    1.0000    2.0569    4.2308
    1.0000    3.1485    9.9131
    1.0000    4.0443   16.3564
    1.0000    5.0165   25.1653
    1.0000    6.2552   39.1275
    1.0000    7.1682   51.3831
    1.0000    8.2789   68.5402
    1.0000    9.2090   84.8057
    1.0000   10.1750  103.5306
c = (X'*X)\X'*V
c = 3x1    
    1.0555
    7.7545
   -0.3097
c = X\V
c = 3x1    
    1.0555
    7.7545
   -0.3097

The constants can now be multiplied with

t.^[0:n]
(1tt2) \left(\begin{array}{ccc} 1 & t & t^2 \end{array}\right)

to get the velocity, here using 4 significant numbers

v = vpa(t.^[0:n]*c,4)
0.3097t2+7.755t+1.055 -0.3097\,t^2 +7.755\,t+1.055

Let us plot the data and the quadratic function

figure; hold on
plot(T,V,'ro')
fplot(v,[0,max(T)],'b-')
xlabel('Time, t [s]')
ylabel('Velocity, v [m/s]')
grid on
tn = linspace(0,max(T),100);
vn = double(subs(v,t,tn));
area(tn, vn, 'FaceColor','b','FaceAlpha', 0.2)

The fit looks quite good, let us look at other methods of applying the least squares method. In Matlab the \ operator is really only short for linsolve(), which is a powerful function, and for matrix system which are not square, linsolve() automatically computes the least squares solution. Let's try it out

X\V
ans = 3x1    
    1.0555
    7.7545
   -0.3097

In addition to this, there are specific functions for fitting polynomials to data, called polyfit():

[c, S] = polyfit(T,V,2)
c = 1x3    
   -0.3097    7.7545    1.0555

S = 
   R: [3x3 double]
   df: 7
   normr: 1.8563

Now back to the problem at hand, first the acceleration, we just integrate the velocity

a = vpa(diff(v,t),4)
7.7550.6195t 7.755-0.6195\,t

Next we plot the acceleration

figure;
fplot(a,[0,max(T)],'m-')
xlabel('Time, t [s]')
ylabel('Acceleration, a [m/s^2]')
grid on

We note that the car is linearly de-accelerating.

Finally, the total distance, which is given by integrating the velocity over the complete time range, i.e., we compute the area under the velocity curve

s = int(v,[0,max(T)])
303.3946047149794068076289477176 303.3946047149794068076289477176

The total distance traveled at the last data-point is 303m.

Example 3 - Rectilinear motion, piecewise

A particle moves along the xx-axis with an initial velocity vx=50  m/sv_x =50\;\mathrm{m}/\mathbf{s} at the origin when t=0t=0. For the first 4 seconds it has no acceleration, and thereafter it is acted on by a retarding force which gives it a constant acceleration ax=10  m/s2a_x =-10\;\mathrm{m}/{\mathrm{s}}^2. Compute the velocity and the xx-coordinate of the particle for the conditions of t=8t=8s and t=12t=12s and find the maximum positive xx-coordinate reached by the particle.

Solution:

given:{axt4s=x¨=10vxt<4s=x˙=50sxt=4s=x=vxt=504=200m\textrm{given}:\left\lbrace \begin{array}{ll} a_x |_{t\ge 4\mathrm{s}} =\ddot{x} =-10 & \\ v_x |_{t<4\mathrm{s}} =\dot{x} =50 & \\ s_x |_{t=4\mathrm{s}} =x=v_x t=50\cdot 4=200\mathrm{m} & \end{array}\right.

Here we have written the given data as an ODE with initial conditions, now we can solve

clear
syms x(t)
DE = diff(x,t,2) == -10
2t2  x(t)=10 \frac{\partial^2 }{\partial t^2 }\;x\left(t\right)=-10
v = diff(x,t)
t  x(t) \frac{\partial }{\partial t}\;x\left(t\right)
IV = [v(4) == 50, x(4) == 50*4]
(((t  x(t))t=4)=50x(4)=200) \left(\begin{array}{cc} {\left({ {\left(\frac{\partial }{\partial t}\;x\left(t\right)\right)}\left|\right.}_{t=4} \right)}=50 & x\left(4\right)=200 \end{array}\right)

We first solve the ODE to get the position, which we then can differentiate w.r.t. time and answer the inquiries

x(t) = dsolve(DE,IV)
t(5t90)80 -t\,{\left(5\,t-90\right)}-80
v(t) = diff(x,t)
9010t 90-10\,t
v(8)
10 10
v(12)
30 -30
t_max = solve(v==0, t)
9 9
x(t_max)
325 325
v(t) = piecewise(0 <= t < 4, 50, ...
                 4 < t,      v)
{50if t[0,4)9010tif 4<t \begin{cases} 50 & \text{if } t \in \left[0,4\right)\\ 90-10t & \text{if } 4 < t \end{cases}

Finally a graph over the situation

figure
hold on
fplot(50*t, [0,4], 'b', 'Linewidth',2)
fplot(x, [4,12],'r', 'Linewidth',2)
xlabel('Time $t$ [s]','Interpreter','Latex')
ylabel('Distance $s(t)$ [m]','Interpreter','Latex')

The area under the graph

tn1 = linspace(0,4,100);
tn2 = linspace(4,12,100);
xn = matlabFunction(x);
area(tn1, 50*tn1, 'Facecolor','b','FaceAlpha',0.2)
area(tn2, xn(tn2), 'Facecolor','r','FaceAlpha',0.2)

And the velocity graph

figure; hold on
fplot(v, [0,12], 'b', 'Linewidth',2)
xlabel('Time $t$ [s]','Interpreter','Latex')
ylabel('Velocity $v(t)$ [m/s]','Interpreter','Latex')
axis([0,12,-30,55])

Example 4 - Oscillating slider

A spring mounted slider is being moved along the xx-axis and has an initial velocity v0=2v_0 =2m/s in the s-direction as it crosses the mid-position where s=0s=0 and t=0t=0. The two springs together exert a retarding force to the motion of the slider, which gives it an acceleration proportional to the displacement but oppositely directed and equal to a=k2sa=-k^2 s, where k=4k=4. Determine the displacement, velocity and acceleration of the slider and plot.

syms x(t) k

syms v0
DE = diff(x,t,2) == -k^2*x
2t2  x(t)=k2x(t) \frac{\partial^2 }{\partial t^2 }\;x\left(t\right)=-k^2 \,x\left(t\right)
v = diff(x,t)
t  x(t) \frac{\partial }{\partial t}\;x\left(t\right)
BV = [v(0) == v0, x(0) == 0]
(((t  x(t))t=0)=v0x(0)=0) \left(\begin{array}{cc} {\left({ {\left(\frac{\partial }{\partial t}\;x\left(t\right)\right)}\left|\right.}_{t=0} \right)}=v_0 & x\left(0\right)=0 \end{array}\right)
x(t) = dsolve([DE,BV])
v0ekti(e2kti1)i2k -\frac{v_0 \,{\mathrm{e} }^{-k\,t\,\mathrm{i} } \,{\left({\mathrm{e} }^{2\,k\,t\,\mathrm{i} } -1\right)}\,\mathrm{i} }{2\,k}
x = simplify(x)
v0sin(kt)k \frac{v_0 \,\sin \left(k\,t\right)}{k}
v = diff(x,t)
v0cos(kt) v_0 \,\cos \left(k\,t\right)
a = diff(v,t)
kv0sin(kt) -k\,v_0 \,\sin \left(k\,t\right)
x = subs(x,[v0, k],[2,4])
sin(4t)2 \frac{\sin \left(4\,t\right)}{2}
v = subs(v,[v0, k],[2,4])
2cos(4t) 2\,\cos \left(4\,t\right)
a = subs(a,[v0, k],[2,4])
8sin(4t) -8\,\sin \left(4\,t\right)
t1 = 2*pi;
figure; hold on;
fplot(x,[0,t1],'Linewidth',2,'DisplayName','Position [m]')
fplot(v,[0,t1],'Linewidth',2,'DisplayName','Velocity [m/s]')
fplot(a,[0,t1],'Linewidth',2,'DisplayName','Acceleration [m/s^2]')
set(gca,'XAxisLocation','origin')
legend('Show')
xlabel('Time $t$ [s]','Interpreter','Latex')
tn = linspace(0,t1,100);
vn = matlabFunction(v);
figure; hold on
fplot(v, [0,t1], 'Linewidth',2,'DisplayName','Velocity [m/s]')
area(tn, vn(tn), 'FaceColor','b','FaceAlpha',0.2,'DisplayName','Position [m]')
set(gca,'XAxisLocation','origin')
legend('Show')
xlabel('Time $t$ [s]','Interpreter','Latex')

The total area is the distance the slider has traveled, which for the time t=2πt=2\piis 8 m. But the signed area is zero, since the sliders position after 2π2\pi seconds is x=0  mx=0\;\mathrm{m}!

double(int([v, abs(v)],0,t1))
ans = 1x2    
     0     8

Finally a graph over the phase space

figure; hold on;
fplot(x,v,[0,t1],'Linewidth',2)
set(gca,'XAxisLocation','origin')
% axis equal
xlabel('Distance $x(t)$ [m]','Interpreter','Latex')
ylabel('Velocity $v(t)$ [m/s]','Interpreter','Latex')
set(gca,'FontSize',14,'FontName','Times')
axis([-1,1,-2.1,2.1])

Example 5 - Piecewise functions

We have noticed the position of a particular nice motorcycle x(t)={t20t<1020t100t10x\left(t\right)=\left\lbrace \begin{array}{ll} t^2 & 0\le t<10\\ 20t-100 & t\ge 10 \end{array}\right.

Plot its velocity and acceleration curves for t[0,30]st\in \left\lbrack 0,30\right\rbrack \mathrm{s}!

syms t
x = piecewise(0 <= t < 10, t^2,...
              t >= 10, 20*t-100)
{t2  if    t[0,10)20t100  if    10t \left\lbrace \begin{array}{cl} t^2 & \;\textrm{if}\;\;t\in \left\lbrack 0,10\right)\\ 20\,t-100 & \;\textrm{if}\;\;10\le t \end{array}\right.

Now we just differentiate

x_dot = diff(x,t)
{2t  if    t(0,10)20  if    10<t \left\lbrace \begin{array}{cl} 2\,t & \;\textrm{if}\;\; t \in \left(0,10\right)\\ 20 & \;\textrm{if}\;\;10 < t \end{array}\right.
x_ddot = diff(x,t,2)
{2  if    t(0,10)0  if    10<t \left\lbrace \begin{array}{cl} 2 & \;\textrm{if}\;\;t\in \left(0,10\right)\\ 0 & \;\textrm{if}\;\;10< t \end{array}\right.
figure;
fplot(x,[0,30],'b','Linewidth',2)
xlabel('Time [s]')
ylabel('Distance [m]')
title('Distance')
figure;
fplot(x_dot,[0,30],'r','Linewidth',2)
xlabel('Time [s]')
ylabel('Velocity [m]')
title('Velocity')
figure;
fplot(x_ddot,[0,30],'m','Linewidth',2)
xlabel('Time [s]')
ylabel('Acceleration [m]')
title('Acceleration')

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  s1k=0\ldotp 1\;s^{-1} with respect to the difference in wanted top speed v1=80  m/sv_1 =80\;\mathrm{m}/\mathrm{s} and the current speed vv. Assume the car starts from rest and the driver puts the pedal to the metal.

Solution:

From the problem description we create the model using the ODE:

ODE={x¨(t)=k(v1x˙(t))x˙(0)=0x(0)=0\textrm{ODE}=\left\lbrace \begin{array}{ll} \ddot{x} \left(t\right)=k\left(v_1 -\dot{x} \left(t\right)\right) & \\ \dot{x} \left(0\right)=0 & \\ x\left(0\right)=0 & \end{array}\right.
clear
syms x(t)
x_dot = diff(x,t)
t  x(t) \dfrac{\partial }{\partial t}\;x\left(t\right)
v_1 = 80; k = 0.1;
DE = diff(x,t,2) == k*(v_1-x_dot)
2t2  x(t)=8t  x(t)10 \dfrac{\partial^2 }{\partial t^2 }\;x\left(t\right)=8-\dfrac{\dfrac{\partial }{\partial t}\;x\left(t\right)}{10}
IV = [x_dot(0)==0, x(0)==0]
(((t  x(t))t=0)=0x(0)=0) \left(\begin{array}{cc} {\left({ {\left(\dfrac{\partial }{\partial t}\;x\left(t\right)\right)}\left|\right.}_{t=0} \right)}=0 & x\left(0\right)=0 \end{array}\right)
x = dsolve(DE,IV)
80t+800et10800 80\,t+800\,{\mathrm{e} }^{-\dfrac{t}{10} } -800
v = diff(x,t)
8080et10 80-80\,{\mathrm{e} }^{-\dfrac{t}{10} }
figure; 
fplot(v*3.6, [0,70], 'Linewidth',2)
xlabel('Time [s]')
ylabel('Velocity [m/s]')
axis([0,70,0,300])
hold on
fplot(100, [0,70],'b--')
vn = matlabFunction(v);

v_max = limit(v*3.6, t, inf)
288 288
vn(70)*3.6
ans = 287.7374
s1 = double(int(v,0,70))
s1 = 4.8007e+03

The car reaches a maximum theoretical speed of 288 km/h which takes around 70 seconds and the car has then traveled around 4.8 km.

syms t positive
t100 = double(solve(v*3.6==100,t))
t100 = 4.2652
double(int(v,0,t100))
ans = 63.4370

The car reaches 100 km/h after 4.27 seconds, not great, not terrible either, by then the car has traveled around 63 meters.

The solution using the numerical Euler forward method

x¨=d2xdt2=dx˙dtΔx˙Δt=x˙i+1x˙iΔt=0.1(80x˙i)x˙i+1=x˙i+0.1(80x˙i)Δtx˙=dxdtΔxΔt=xi+1xiΔt=x˙i+1xi+1=xi+x˙i+1Δt \def\arraystretch{2.5} \begin{array}{l} \ddot{x} =\dfrac{d^2 x}{ {dt}^2 }=\dfrac{d\dot{x} }{dt}\approx \dfrac{\Delta \dot{x} }{\Delta t}=\dfrac{ {\dot{x} }_{i+1} -{\dot{x} }_i }{\Delta t}=0\ldotp 1\left(80-{\dot{x} }_i \right)\Rightarrow {\dot{x} }_{i+1} ={\dot{x} }_i +0\ldotp 1\left(80-{\dot{x} }_i \right)\Delta t\\ \dot{x} =\dfrac{dx}{dt}\approx \dfrac{\Delta x}{\Delta t}=\dfrac{x_{i+1} -x_i }{\Delta t}={\dot{x} }_{i+1} \Rightarrow x_{i+1} =x_i +{\dot{x} }_{i+1} \Delta t \end{array}
x_i = 0;
x_dot_i = 0;
ti = 0;
dt = 0.1;

t1 = 70;
n = length(0:dt:t1)
n = 701
x = NaN(n,1);
t =(0:dt:t1)';
x_dot = x;

fprintf('i | t     | x_dot    |  x_i')
fprintf('------------------------------')
i = 1;
while ti <= t1
    ti = ti + dt;
    x_dot_i = x_dot_i + 0.1*(80-x_dot_i)*dt;
    x_i = x_i + x_dot_i*dt;
    fprintf('%d | %0.2f  | %0.4f   |  %0.4f  \n',i,ti,x_dot_i, x_i)

    t(i) = ti;
    x(i) = x_i;
    x_dot(i) = x_dot_i;

    i = i+1;
end
i | t     | x_dot    |  x_i
------------------------------
1 | 0.10  | 0.8000   |  0.0800  
2 | 0.20  | 1.5920   |  0.2392  
3 | 0.30  | 2.3761   |  0.4768  
4 | 0.40  | 3.1523   |  0.7920  
5 | 0.50  | 3.9208   |  1.1841  
6 | 0.60  | 4.6816   |  1.6523  
7 | 0.70  | 5.4348   |  2.1958  
...
698 | 69.80  | 79.9281   |  4792.7114  
699 | 69.90  | 79.9289   |  4800.7042  
700 | 70.00  | 79.9296   |  4808.6972
figure;
plot(t,x,'b-','Linewidth',2)
xlabel('Time [s]')
ylabel('Distance [m]')
title('Numerical solution')
figure;
plot(t,x_dot,'r-','Linewidth',2)
xlabel('Time [s]')
ylabel('Velocity [m/s]')
title('Numerical solution')

Example 7 - Curvilinear motion

The curvilinear motion of a particle is defined by vx=5016tv_x =50-16t and y=1004t2y=100-4t^2, where vxv_x is in meter per second, yy is in meters and tt is in seconds. It is also known that x=0x=0 when t=0t=0. Plot the path of the particle and determine its velocity and acceleration when the position y=0y=0 is reached.

clear
syms x(t) y(t)
syms t positive

y = 100 - 4*t^2
1004t2 100-4\,t^2
DE = diff(x,t) == 50 - 16*t
t  x(t)=5016t \dfrac{\partial }{\partial t}\;x\left(t\right)=50-16\,t
BV = x(0) == 0
x(0)=0 x\left(0\right)=0
x = dsolve(DE,BV)
2t(4t25) -2\,t\,{\left(4\,t-25\right)}
t_y0 = solve(y==0,t)
5 5
figure
fplot(x,y,[0,double(t_y0)],'Linewidth',2)
grid on
xlabel('$x(t)$ [m]','Interpreter','latex');
ylabel('$y(t)$ [m]','Interpreter','latex');
axis equal
title("The path of the particle")

Now we can compute the velocity and acceleration, first simply in Cartesian coordinates

r = [x,y].'
(2t(4t25)1004t2) \left(\begin{array}{c} -2\,t\,{\left(4\,t-25\right)}\\ 100-4\,t^2 \end{array}\right)
r_dot = diff(r,t)
(5016t8t) \left(\begin{array}{c} 50-16\,t\\ -8\,t \end{array}\right)
r_ddot = diff(r,t,2)
(168) \left(\begin{array}{c} -16\\ -8 \end{array}\right)

Velocity when y=0y=0

r_dot_y0 = subs(r_dot, t, t_y0)
(3040) \left(\begin{array}{c} -30\\ -40 \end{array}\right)

t1 = double(t_y0);
figure; hold on
fplot(x,y,[0,t1],'Linewidth',2)
grid on
xlabel('$x(t)$ [m]','Interpreter','latex');
ylabel('$y(t)$ [m]','Interpreter','latex');
axis equal
title("The path of the particle")

xn = matlabFunction(x);
yn = matlabFunction(y);

v = matlabFunction(r_dot);
a = double(r_ddot)
a = 2x1    
   -16
    -8
t = 0;
vi = v(t);
xi = xn(t);
yi = yn(t);

plot(xn(0),yn(0),'ko','MarkerFacecolor','k','MarkerSize',4)
text(xn(0),yn(0)+5,'$t=0$','Interpreter','latex')

plot(xn(1),yn(1),'ko','MarkerFacecolor','k','MarkerSize',4)
text(xn(1)-10,yn(1)-5,'$t=1$','Interpreter','latex')

plot(xn(2),yn(2),'ko','MarkerFacecolor','k','MarkerSize',4)
text(xn(2)-15,yn(2)-5,'$t=2$','Interpreter','latex')

plot(xn(3),yn(3),'ko','MarkerFacecolor','k','MarkerSize',4)
text(xn(3)-18,yn(3)-1,'$t=3$','Interpreter','latex')

plot(xn(4),yn(4),'ko','MarkerFacecolor','k','MarkerSize',4)
text(xn(4)-18,yn(4)-1,'$t=4$','Interpreter','latex')

plot(xn(5),yn(5),'ko','MarkerFacecolor','k','MarkerSize',4)
text(xn(5)+5,yn(5)-1,'$t=5$','Interpreter','latex')

hqv = quiver(xi, yi, vi(1), vi(2), 1,'b');
hqa = quiver(xi, yi, a(1), a(2), 1,'r');
hp = plot(xi,yi,'ko','MarkerFacecolor','k');
axis([-10,100,-45,110])
% for t = linspace(0,t1,100)
%     vi = v(t);
%     xi = xn(t);
%     yi = yn(t);
%     set(hp,'XData',xi,'Ydata',yi)
%     set(hqv,'XData',xi,'Ydata',yi,'UData',vi(1),'VData',vi(2))
%     set(hqa,'XData',xi,'Ydata',yi,'UData',a(1),'VData',a(2))
%     drawnow
% end

We note that the acceleration is constant in this perticular examples, which is due to the velocity constraint in the x-direction.

Example 8 - Curvelinear motion, polar coordinates

The motion of a particle, BB, confined to the radial screw of an arm is to be studied. Let r(t)r\left(t\right) and θ(t)\theta \left(t\right) and analyse the motion of the particle in Cartesian and polar coordinates!

We begin by expressing the particles position

clear

syms theta(t) r(t) 

rr = r*[cos(theta), sin(theta)].'
(cos(θ(t))r(t)sin(θ(t))r(t)) \left(\begin{array}{c} \cos \left(\theta \left(t\right)\right)\,r\left(t\right)\\ \sin \left(\theta \left(t\right)\right)\,r\left(t\right) \end{array}\right)

Then we let matlab compute the derivatives

rr_dot = diff(rr,t);
syms r_dot theta_dot r_ddot theta_ddot
rr_dot = subs(rr_dot, [diff(r,t), diff(theta,t)], [r_dot,theta_dot])
(r˙cos(θ(t))θ˙sin(θ(t))r(t)r˙sin(θ(t))+θ˙cos(θ(t))r(t)) \left(\begin{array}{c} \dot{r} \,\cos \left(\theta \left(t\right)\right)-\dot{\theta} \,\sin \left(\theta \left(t\right)\right)\,r\left(t\right)\\ \dot{r} \,\sin \left(\theta \left(t\right)\right)+\dot{\theta} \,\cos \left(\theta \left(t\right)\right)\,r\left(t\right) \end{array}\right)
rr_ddot = diff(rr,t,2);
rr_ddot = subs(rr_ddot, [diff(r,t,2), diff(theta,t,2)], ...
                        [r_ddot, theta_ddot]);
rr_ddot = subs(rr_ddot, [diff(r,t), diff(theta,t)], ...
                        [r_dot, theta_dot])
r¨(t)=(cos(θ(t))r(t)θ˙22r˙sin(θ(t))θ˙+r¨cos(θ(t))θ¨sin(θ(t))r(t)sin(θ(t))r(t)θ˙2+2r˙cos(θ(t))θ˙+r¨sin(θ(t))+θ¨cos(θ(t))r(t)) \ddot{\bm r}(t) = \left(\begin{array}{c} -\cos \left(\theta \left(t\right)\right)\,r\left(t\right)\,{\dot{\theta} }^2 -2\,\dot{r} \,\sin \left(\theta \left(t\right)\right)\,\dot{\theta} +\ddot{r} \,\cos \left(\theta \left(t\right)\right)-\ddot{\theta} \,\sin \left(\theta \left(t\right)\right)\,r\left(t\right)\\ -\sin \left(\theta \left(t\right)\right)\,r\left(t\right)\,{\dot{\theta} }^2 +2\,\dot{r} \,\cos \left(\theta \left(t\right)\right)\,\dot{\theta} +\ddot{r} \,\sin \left(\theta \left(t\right)\right)+\ddot{\theta} \,\cos \left(\theta \left(t\right)\right)\,r\left(t\right) \end{array}\right)

In the above expressions, we made the derivatives easier to read by converting to Newtonian notation. Clearly, Matlab can handle tedious derivation.

Now, we define the polar bases

e_r = [cos(theta), sin(theta)].';
e_theta = [-sin(theta), cos(theta)].';

where eθe_{\theta } is orthogonal according to the right-hand rule (rotate 90 degrees counter clockwise).

Now just project the Cartesial vectors onto the polar bases. We can get away with the "low budget" projection since the bases have unit length, i.e., er=eθ=1|{\bm e}_r |=|{\bm e}_{\theta } |=1. We recall that the dot product is just a matrix multiplication

[a1,a2][b1b2]=a1b1+a2b2\left\lbrack a_1 ,a_2 \right\rbrack \left\lbrack \begin{array}{c} b_1 \\ b_2 \end{array}\right\rbrack =a_1 b_1 +a_2 b_2
r_dot_rtheta = simplify([rr_dot.'*e_r, rr_dot.'*e_theta])
(r˙θ˙r(t)) \left(\begin{array}{cc} \dot{r} & \dot{\theta} \,r\left(t\right) \end{array}\right)
r_ddot_rtheta = simplify([rr_ddot.'*e_r, rr_ddot.'*e_theta])
(r¨θ˙2r(t)2r˙θ˙+θ¨r(t)) \left(\begin{array}{cc} \ddot{r} -{\dot{\theta} }^2 \,r\left(t\right) & 2\,\dot{r} \,\dot{\theta} +\ddot{\theta} \,r\left(t\right) \end{array}\right)

We can see that we end up with the same expressions as we have tediously derived using the classical approach in the sections above.

Now, let θ(t)=0.2t+0.02t3\theta \left(t\right)=0\ldotp 2t+0\ldotp 02t^3 and r(t)=0.2+0.04t2r\left(t\right)=0\ldotp 2+0\ldotp 04t^2 and create plots over the situation.

syms t
theta = 0.2*t+0.02*t^3;
r = 0.2 + 0.04*t^2;
rr = r*[cos(theta), sin(theta)].';
x = rr(1)
cos(t350+t5)(t225+15) \cos \left(\dfrac{t^3 }{50}+\dfrac{t}{5}\right)\,{\left(\dfrac{t^2 }{25}+\dfrac{1}{5}\right)}
y = rr(2)
sin(t350+t5)(t225+15) \sin \left(\dfrac{t^3 }{50}+\dfrac{t}{5}\right)\,{\left(\dfrac{t^2 }{25}+\dfrac{1}{5}\right)}
xn = matlabFunction(x);
yn = matlabFunction(y);

figure
fplot(rr(1), rr(2), [0,5], 'LineWidth', 2)
axis equal
grid on
xlabel('$x(t)$ [m]','Interpreter','latex');
ylabel('$y(t)$ [m]','Interpreter','latex');
title("The path of the particle")
axis([-1.3,0.5,-0.5,0.8])
text(xn(1),yn(1),'$t=1$','Interpreter','latex','Backgroundcolor','w')
text(xn(2),yn(2),'$t=2$','Interpreter','latex','Backgroundcolor','w')
text(xn(3),yn(3),'$t=3$','Interpreter','latex','Backgroundcolor','w')
text(xn(4),yn(4),'$t=4$','Interpreter','latex','Backgroundcolor','w')
text(xn(5),yn(5),'$t=5$','Interpreter','latex','Backgroundcolor','w')
figure; hold on;
fplot(diff(x,t),[0,5],'LineWidth', 2,'DisplayName','$\dot x$')
fplot(diff(y,t),[0,5],'LineWidth', 2,'DisplayName','$\dot y$')
fplot(diff(x,t,2),[0,5],'LineWidth', 2,'DisplayName','$\ddot x$')
fplot(diff(y,t,2),[0,5],'LineWidth', 2,'DisplayName','$\ddot y$')
legend('show','Interpreter','latex')

syms t
xy_dot = matlabFunction( simplify(diff([x,y],t)) );
xy_ddot = matlabFunction( simplify(diff([x,y],t,2)) );

figure; hold on
fplot(rr(1), rr(2), [0,5], 'LineWidth', 2)
axis equal
grid on
xlabel('$x(t)$ [m]','Interpreter','latex');
ylabel('$y(t)$ [m]','Interpreter','latex');
title("The path of the particle")
axis([-1.5,0.7,-0.5,1.3])

t = 1;
vi = xy_dot(t);
ai = xy_ddot(t);
xi = xn(t);
yi = yn(t);
hp = plot([0,xn(t)],[0,yn(t)],'k-');
hv = quiver(xi, yi, vi(1), vi(2),0.4,'b');
ha = quiver(xi, yi, ai(1), ai(2),0.4,'r');
for t = linspace(0,5,200)
    vi = xy_dot(t);
    ai = xy_ddot(t);
    xi = xn(t);
    yi = yn(t);
    set(hp,'XData',xi,'Ydata',yi)
    set(hv,'XData',xi,'Ydata',yi,'UData',vi(1),'VData',vi(2))
    set(ha,'XData',xi,'Ydata',yi,'UData',ai(1),'VData',ai(2))
    drawnow
end

Example 10 - Radar tracking, 3D, cylindrical and spherical coordinates

An aircraft PP takes off at AA with a velocity v0v_0 of 250 km/h and climbs in the vertical yzy^{\prime } -z^{\prime } plane at the constant 15{15}^{\circ } angle with an acceleration along its flight path of 0.8  m/s20\ldotp 8\;\mathrm{m}/{\mathrm{s}}^2. The flight progress is monitored by radar at point OO.

clear
syms t positive
syms y(t) z(t)

v0 = 250/3.6;
alpha = 15*pi/180;
a = 0.8;

OA = [3000; 0; 0];
eAP = [0; cos(alpha); sin(alpha)];
V0 = v0*eAP;
Acc = a*eAP;

x = 3000;

yp = diff(y,1);
ypp = diff(y,2);
DEy = ypp == Acc(2);
BVy = [yp(0) == V0(2); y(0) == 0];
y = vpa(dsolve(DEy,BVy), 4)
y(t)=0.3864t2+67.08t y(t) = 0.3864\,t^2 +67.08\,t
zp = diff(z,1);
zpp = diff(z,2);
DEz = zpp == Acc(3);
BVz = [zp(0) == V0(3); z(0) == 0];
z = vpa(dsolve(DEz,BVz), 4)
z(t)=0.1035t2+17.97t z(t) = 0.1035\,t^2 +17.97\,t
R = vpa([x; y; z],4)
(3000.00.3864t2+67.08t0.1035t2+17.97t) \left(\begin{array}{c} 3000.0\\ 0.3864\,t^2 +67.08\,t\\ 0.1035\,t^2 +17.97\,t \end{array}\right)
Rp = vpa(diff(R,1),4)
(00.7727t+67.080.2071t+17.97) \left(\begin{array}{c} 0\\ 0.7727\,t+67.08\\ 0.2071\,t+17.97 \end{array}\right)
Rpp = vpa(diff(R,2),4)
(00.77270.2071) \left(\begin{array}{c} 0\\ 0.7727\\ 0.2071 \end{array}\right)

The situation graphically

figure; hold on; grid on; axis equal
plot3(subs(R(1),t,0),subs(R(2),t,0),subs(R(3),t,0),'ob', 'Displayname',"start")
fplot3(R(1),R(2),R(3),[0,60],'-b', 'Displayname',"path")
plot3([0,4000], [0,0], [0,0], '--k', 'HandleVisibility','off')
plot3([0,0], [0,4000], [0,0], '--k', 'HandleVisibility','off')
plot3([0,0], [0,0], [0,4000], '--k', 'HandleVisibility','off')
plot3(subs(R(1),t,60),subs(R(2),t,60),subs(R(3),t,60),'or', 'Displayname',"$t=60$s")
plot3(0,0,0,'ok', 'HandleVisibility','off')
quiver3(subs(R(1),t,60),subs(R(2),t,60),subs(R(3),t,60),subs(Rp(1),t,60),subs(Rp(2),t,60),subs(Rp(3),t,60),15,'k','LineWidth', 1, 'Displayname',"velocity")
quiver3(subs(R(1),t,60),subs(R(2),t,60),subs(R(3),t,60),0,subs(Rp(2),t,60),0,15,'g','LineWidth', 1, 'Displayname',"y-velocity")
quiver3(subs(R(1),t,60),subs(R(2),t,60),subs(R(3),t,60),0,0,subs(Rp(3),t,60),15,'b','LineWidth', 1, 'Displayname',"z-velocity")
view([43 23])
xlabel('x')
ylabel('y')
zlabel('z')
legend("show","location","northwest",Interpreter="latex")
title('Velocity in cartesian coordinates')

Cylindrical coordinates

OP_xy = [x; y; 0];
eOP_xy = OP_xy/(sqrt(x^2 + y^2));
r_cyl = sqrt(x^2 + y^2);
theta = atan(y/x);
rp_cyl = diff(r_cyl);
thetap = diff(theta);
Rp_cyl = rp_cyl * eOP_xy;
Theta_p = thetap * [0;0;1];
Tan_p = cross(Theta_p,r_cyl*eOP_xy);

rp_cyl_60 = vpa(subs(rp_cyl,t,60),5)
rp_cyl_60 =  99.234
thetap_60 = vpa(subs(thetap,t,60),5)
thetap_60 =  0.0088791
zp_60 = vpa(subs(Rp(3),t,60),5)
zp_60 =  30.397
figure
hold on
grid on
axis equal
plot3(subs(R(1),t,0),subs(R(2),t,0),subs(R(3),t,0),'ob', 'Displayname',"start")
fplot3(R(1),R(2),R(3),[0,60],'-b', 'Displayname',"path")
plot3([0,4000], [0,0], [0,0], '--k', 'HandleVisibility','off')
plot3([0,0], [0,4000], [0,0], '--k', 'HandleVisibility','off')
plot3([0,0], [0,0], [0,4000], '--k', 'HandleVisibility','off')
plot3(subs(R(1),t,60),subs(R(2),t,60),subs(R(3),t,60),'or', 'Displayname',"$t=60$s")
plot3(0,0,0,'ok', 'HandleVisibility','off')
quiver3(subs(R(1),t,60),subs(R(2),t,60),subs(R(3),t,60),subs(Rp(1),t,60),subs(Rp(2),t,60),subs(Rp(3),t,60),15,'k','LineWidth', 1, 'Displayname',"$V$")
quiver3(subs(R(1),t,60),subs(R(2),t,60),subs(R(3),t,60),subs(Rp_cyl(1),t,60),subs(Rp_cyl(2),t,60),subs(Rp_cyl(3),t,60),15,'r','LineWidth', 1, 'Displayname',"$V_r$")
quiver3(subs(R(1),t,60),subs(R(2),t,60),subs(R(3),t,60),subs(Tan_p(1),t,60),subs(Tan_p(2),t,60),subs(Tan_p(3),t,60),15,'g','LineWidth', 1, 'Displayname',"$V_\theta$")
quiver3(subs(R(1),t,60),subs(R(2),t,60),subs(R(3),t,60),0,0,subs(Rp(3),t,60),15,'b','LineWidth', 1, 'Displayname',"$V_z$")
view([43 23])
xlabel('x');ylabel('y');zlabel('z')
legend("show","location","northwest",Interpreter="latex")
title('Velocity in cylindrical coordinates')

Spherical coordinates

r_sph = sqrt(x^2 + y^2 + z^2);
OP = [x; y; z];
eOP = OP/r_sph;
phi = atan(z/r_sph);
rp_sph = diff(r_sph);
phip = diff(phi);
Rp_sph = rp_sph * eOP;
Phi_p = phip * cross(eOP,[0;0;1]);
Phi_tan_p = cross(Phi_p,r_sph*eOP);
Tan_p = cross(Theta_p,r_sph*eOP);

rp_sph_60 = vpa(subs(rp_sph,t,60),5)
103.55 103.55
thetap_60 = vpa(subs(thetap,t,60),5)
0.0088791 0.0088791
phip_60 = vpa(subs(phip,t,60),5)
0.0010113 0.0010113
figure; hold on; grid on; axis equal

plot3(subs(R(1),t,0),subs(R(2),t,0),subs(R(3),t,0),'ob', 'Displayname',"start")
fplot3(R(1),R(2),R(3),[0,60],'-b', 'Displayname',"path")
plot3([0,4000], [0,0], [0,0], '--k', 'HandleVisibility','off')
plot3([0,0], [0,4000], [0,0], '--k', 'HandleVisibility','off')
plot3([0,0], [0,0], [0,4000], '--k', 'HandleVisibility','off')
plot3(subs(R(1),t,60),subs(R(2),t,60),subs(R(3),t,60),'or', 'Displayname',"$t=60$s")
plot3(0,0,0,'ok', 'HandleVisibility','off')
quiver3(subs(R(1),t,60),subs(R(2),t,60),subs(R(3),t,60),subs(Rp(1),t,60),subs(Rp(2),t,60),subs(Rp(3),t,60),15,'k','LineWidth', 1, 'Displayname',"$V$")
quiver3(subs(R(1),t,60),subs(R(2),t,60),subs(R(3),t,60),subs(Rp_sph(1),t,60),subs(Rp_sph(2),t,60),subs(Rp_sph(3),t,60),15,'r','LineWidth', 1, 'Displayname',"$V_r$")
quiver3(subs(R(1),t,60),subs(R(2),t,60),subs(R(3),t,60),subs(Tan_p(1),t,60),subs(Tan_p(2),t,60),subs(Tan_p(3),t,60),15,'g','LineWidth', 1, 'Displayname',"$V_\theta$")
quiver3(subs(R(1),t,60),subs(R(2),t,60),subs(R(3),t,60),subs(Phi_tan_p(1),t,60),subs(Phi_tan_p(2),t,60),subs(Phi_tan_p(3),t,60),100,'b','LineWidth', 1, 'Displayname',"$V_\phi$")
view([43 23])
xlabel('x'); ylabel('y'); zlabel('z')
legend("show","location","northwest",Interpreter="latex")
title('Velocity in spherical coordinates')

Example 11 - Lead screw, 3D helical motion

The power screw starts from rest and is given a rotational speed θ˙\dot{\theta} which increases uniformly with time tt according to θ˙=k  t\dot{\theta} =k\;t, where kk is constant. Determine the velocity and acceleration of the center of ball AA when the screw has turned through one complete revolution from rest. The lead of the screw (advancement per revolution) is LL.

clear
syms t positive
syms theta(t) r(t)

Let k=2k=2, r0=20r_0 =20 and L=3L=3.

k = 2;
r0 = 20;
L = 3;
theta_p = diff(theta)
t  θ(t) \dfrac{\partial }{\partial t}\;\theta \left(t\right)
DE = theta_p == k*t
t  θ(t)=2t \dfrac{\partial }{\partial t}\;\theta \left(t\right)=2\,t
BV = theta(0)==0
θ(0)=0 \theta \left(0\right)=0
theta = dsolve(DE,BV)
t2 t^2
t1 = double( max(solve(theta==2*pi,t)) )
t1 = 2.5066
x = cos(theta)*r0
20cos(t2) 20\,\cos \left(t^2 \right)
y = sin(theta)*r0
20sin(t2) 20\,\sin \left(t^2 \right)
z = -theta*L/(2*pi)
3t22π -\dfrac{3\,t^2 }{2\,\pi }
Rp = diff([x;y;z])
(40tsin(t2)40tcos(t2)3tπ) \left(\begin{array}{c} -40\,t\,\sin \left(t^2 \right)\\ 40\,t\,\cos \left(t^2 \right)\\ -\dfrac{3\,t}{\pi } \end{array}\right)
Rpp = diff([x;y;z],2)
(40sin(t2)80t2cos(t2)40cos(t2)80t2sin(t2)3π) \left(\begin{array}{c} -40\,\sin \left(t^2 \right)-80\,t^2 \,\cos \left(t^2 \right)\\ 40\,\cos \left(t^2 \right)-80\,t^2 \,\sin \left(t^2 \right)\\ -\dfrac{3}{\pi } \end{array}\right)
Rn = matlabFunction([x;y;z]);
Rpn = matlabFunction(Rp);
Rppn = matlabFunction(Rpp);

Ri = Rn(t1)
Ri = 3x1    
   20.0000
    0.0000
   -3.0000
Rpi = Rpn(t1)
Rpi = 3x1    
   -0.0000
  100.2651
   -2.3937
Rppi = Rppn(t1)
Rppi = 3x1    
 -502.6548
   40.0000
   -0.9549
figure; hold on; grid on; axis equal
fplot3(x,y,z,[0,4], 'Displayname',"path")
plot3(r0,0,0,'ok', 'Displayname',"starting point")
hp = plot3(Ri(1),Ri(2),Ri(3),'ob', 'Displayname',"$\theta=2\pi$");
hqv = quiver3(Ri(1),Ri(2),Ri(3),Rpi(1),Rpi(2),Rpi(3),0.1,'b','LineWidth', 2, 'Displayname',"velocity");
hqa = quiver3(Ri(1),Ri(2),Ri(3),Rppi(1),Rppi(2),Rppi(3),0.02,'r','LineWidth', 2, 'Displayname',"acceleration");
view([48.4 31.5])
xlabel('x'); ylabel('y'); zlabel('z')
legend("show","location","northwest",Interpreter="latex")
axis([-25,25,-25,25,-25,25])
for ti = linspace(0,4,100)
     Ri = Rn(ti);
     Rpi = Rpn(ti);
     Rppi = Rppn(ti);
     set(hp,'XData',Ri(1), 'YData', Ri(2), 'Zdata', Ri(3))
     set(hqv,'XData',Ri(1), 'YData', Ri(2), 'Zdata', Ri(3), 'UData', Rpi(1), 'VData', Rpi(2), 'WData', Rpi(3))
     set(hqa,'XData',Ri(1), 'YData', Ri(2), 'Zdata', Ri(3), 'UData', Rppi(1), 'VData', Rppi(2), 'WData', Rppi(3))
     drawnow
 end

Example 12 - Tractor hoisting, constrained motion

A tractor is used to hoist the bale BB using a pulley arrangement as shown. If the tractor at AA has a forward velocity vxv_x, determine an expression for the upward velocity vBv_B

Solution:

This is a so called constrained motion. We can tackle this problem by finding the geometric relation between AA and BB. We take note of a triangle and thus utilize the Pythagorean theorem to formulate

x(t)2+h2=l(t)2{x\left(t\right)}^2 +h^2 ={l\left(t\right)}^2

Also, we note the geometry of the rope, which has a constant length LL.

L=2(hy(t))+l(t)L=2\left(h-y\left(t\right)\right)+l\left(t\right)
clear
syms l(t) h x(t) y(t) L

ekv = [x^2 + h^2 == l^2
       L == 2*(h-y) + l];
ekv = [ekv; diff(ekv,t)];

Cleaning up the output and simplifying

syms l_dot v_A v_B positive
ekv = subs(ekv, [diff(x,t), diff(y,t), diff(l,t)], [v_A, v_B, l_dot]);

ekv = subs(ekv,[x,y,l],[sym('x'), sym('y'), sym('l')])
(h2+x2=l2L=2h+l2y2vAx=2ll˙0=l˙2vB) \left(\begin{array}{c} h^2 +x^2 =l^2 \\ L=2\,h+l-2\,y\\ 2\,v_A \,x=2\,l\,\dot{l} \\ 0=\dot{l} -2\,v_B \end{array}\right)
syms y l x positive
[y, v_B, l, l_dot] = solve(ekv,[y, v_B, l, l_dot])
Warning: Solutions are only valid under certain conditions. To include parameters and conditions in the solution, specify the 'ReturnConditions' value as 'true'.
y =
hL2+1h2+x2(h2+x2)2 h-\dfrac{L}{2}+\dfrac{\sqrt{\dfrac{1}{h^2 +x^2 } }\,{\left(h^2 +x^2 \right)} }{2}
v_B =
vAx1h2+x22 \dfrac{v_A \,x\,\sqrt{\dfrac{1}{h^2 +x^2 } } }{2}
l =
1h2+x2(h2+x2) \sqrt{\dfrac{1}{h^2 +x^2 } }\,{\left(h^2 +x^2 \right)}
l_dot =
vAx1h2+x2 v_A \,x\,\sqrt{\dfrac{1}{h^2 +x^2 } }

Animation

h = 5;
L = 3*h;
r_A = [x;0]
(x0) \left(\begin{array}{c} x\\ 0 \end{array}\right)
r_B = simplify(subs([0;y]))
(0x2+25252) \left(\begin{array}{c} 0\\ \dfrac{\sqrt{x^2 +25} }{2}-\dfrac{5}{2} \end{array}\right)
r_C = [0;h]
r_C = 2x1    
     0
     5
r_A_n = matlabFunction(r_A);
r_B_n = matlabFunction(r_B);

xi = 4;
r_A_i = r_A_n(xi)
r_A_i = 2x1    
     4
     0
r_B_i = r_B_n(xi)
r_B_i = 2x1    
         0
    0.7016
x_at_y4 = double(solve(r_B(2)==4,x))
x_at_y4 = 12
figure; hold on; axis equal;
hAC = plot([r_C(1), r_A_i(1)], [r_C(2), r_A_i(2)], '-ko', ...
    'MarkerFaceColor','k','LineWidth',2);
hBC = plot([r_C(1), r_B_i(1)], [r_C(2), r_B_i(2)], '-ko', ...
    'MarkerFaceColor','k','LineWidth',2);
htA = text(r_A_i(1)+0.1,r_A_i(2)+0.2,'$A$','interpreter','latex');
text(r_C(1)+0.1,r_C(2)+0.2,'$C$','interpreter','latex');
htB = text(r_B_i(1)+0.1,r_B_i(2)+0.2,'$B$','interpreter','latex');
axis([-1,15,0,5])

for xi = linspace(0,x_at_y4,50)
    r_A_i = r_A_n(xi);
    r_B_i = r_B_n(xi);
    set(hAC,'XData',[r_C(1), r_A_i(1)],'YData',[r_C(2), r_A_i(2)])
    set(hBC,'XData',[r_C(1), r_B_i(1)],'YData',[r_C(2), r_B_i(2)])
    set(htA,'Position',[r_A_i(1)+0.1,r_A_i(2)+0.2])
    set(htB,'Position',[r_B_i(1)+0.1,r_B_i(2)+0.2])
    drawnow
end

Example 13 - Crank slider mechanism

Study the motion of the slider-crank mechanism. Plot the motion of A,GA,G and BB as well as their velocities and accelerations as a function of θ\theta. Let the angular velocity ω=25\omega =25 rad/s.

Solution:

Typical example, we can tackle this using implicit derivation and let Matlab do the heavy lifting. Model the geometry relations first. Let LL denote the distance AB\textrm{AB}. Then we utilize Pythagoras' theorem to formulate the point AA, using vectors. Let h  h\;denote the height over the x-axis to point BB.

h=r  sinθ(t)  and  a=r  cosθ(t)  as  well  as  b=L2h2h=r\;\sin \theta \left(t\right)\;\textrm{and}\;a=r\;\cos \theta \left(t\right)\;\textrm{as}\;\textrm{well}\;\textrm{as}\;b=\sqrt{L^2 -h^2 }
OA=a+b\textrm{OA}=a+b rA=[a+b,0]{\bm r}_A =\left\lbrack a+b,0\right\rbrack rb=r[cosθ,sinθ]{\bm r}_b =r\left\lbrack \cos \theta ,\sin \theta \right\rbrack rBA=rArB{\bm r}_{\textrm{BA} } ={\bm r}_A -{\bm r}_B rG=rB+100100+250rBA{\bm r}_G ={\bm r}_B +\dfrac{100}{100+250}{\bm r}_{\textrm{BA} }
clear
syms theta(t)
syms L r positive

h = r*sin(theta);
a = r*cos(theta);
b = sqrt(L^2 - h^2);

r_A = [-(a+b); 0]
(rcos(θ(t))L2r2sin(θ(t))20) \left(\begin{array}{c} -r\,\cos \left(\theta \left(t\right)\right)-\sqrt{L^2 -r^2 \,{\sin \left(\theta \left(t\right)\right)}^2 }\\ 0 \end{array}\right)
r_B = r*[-cos(theta); sin(theta)]
(rcos(θ(t))rsin(θ(t))) \left(\begin{array}{c} -r\,\cos \left(\theta \left(t\right)\right)\\ r\,\sin \left(\theta \left(t\right)\right) \end{array}\right)
r_BA = r_A - r_B;
r_G = r_B + 100/(100+250)*r_BA
(rcos(θ(t))2L2r2sin(θ(t))275rsin(θ(t))7) \left(\begin{array}{c} -r\,\cos \left(\theta \left(t\right)\right)-\dfrac{2\,\sqrt{L^2 -r^2 \,{\sin \left(\theta \left(t\right)\right)}^2 } }{7}\\ \dfrac{5\,r\,\sin \left(\theta \left(t\right)\right)}{7} \end{array}\right)
syms theta_dot theta_ddot positive
r_A_dot = subs(diff(r_A,t), [diff(theta,t), theta], [theta_dot, sym('theta', 'positive')]);
r_A_ddot = subs(diff(r_A,t,2), [diff(theta,t), diff(theta,t,2), theta], ...
                               [theta_dot, theta_ddot, sym('theta', 'positive')]);
r_A_dot = formula(r_A_dot)
(rθ˙sin(θ)+r2θ˙cos(θ)sin(θ)L2r2sin(θ)20) \left(\begin{array}{c} r\,\dot{\theta} \,\sin \left(\theta \right)+\dfrac{r^2 \,\dot{\theta} \,\cos \left(\theta \right)\,\sin \left(\theta \right)}{\sqrt{L^2 -r^2 \,{\sin \left(\theta \right)}^2 } }\\ 0 \end{array}\right)
r_A_ddot = formula(r_A_ddot)
(rθ˙2cos(θ)+r2θ˙2cos(θ)2σ1r2θ˙2sin(θ)2σ1+r4θ˙2cos(θ)2sin(θ)2σ13/20)where    σ1=L2r2sin(θ)2 \begin{array}{l} \left(\begin{array}{c} r\,{\dot{\theta} }^2 \,\cos \left(\theta \right)+\dfrac{r^2 \,{\dot{\theta} }^2 \,{\cos \left(\theta \right)}^2 }{\sqrt{\sigma_1 } }-\dfrac{r^2 \,{\dot{\theta} }^2 \,{\sin \left(\theta \right)}^2 }{\sqrt{\sigma_1 } }+\dfrac{r^4 \,{\dot{\theta} }^2 \,{\cos \left(\theta \right)}^2 \,{\sin \left(\theta \right)}^2 }{ {\sigma_1 }^{3/2} }\\ 0 \end{array}\right)\\ \mathrm{}\\ \textrm{where}\\ \mathrm{}\\ \;\;\sigma_1 =L^2 -r^2 \,{\sin \left(\theta \right)}^2 \end{array}
r_B_dot = subs(diff(r_B,t), [diff(theta,t), theta], [theta_dot, sym('theta','positive')]);
r_B_dot = formula(r_B_dot)
(rθ˙sin(θ)rθ˙cos(θ)) \left(\begin{array}{c} r\,\dot{\theta} \,\sin \left(\theta \right)\\ r\,\dot{\theta} \,\cos \left(\theta \right) \end{array}\right)
r_B_dot_m = simplify(norm(r_B_dot))
rθ˙ r\,\dot{\theta}
r_B_ddot = subs(diff(r_B,t,2), [diff(theta,t), diff(theta,t,2), theta], ...
                               [theta_dot, theta_ddot, sym('theta', 'positive')]);
r_B_ddot = formula(r_B_ddot)
(rθ˙2cos(θ)rθ˙2sin(θ)) \left(\begin{array}{c} r\,{\dot{\theta} }^2 \,\cos \left(\theta \right)\\ -r\,{\dot{\theta} }^2 \,\sin \left(\theta \right) \end{array}\right)
r_B_ddot_m = simplify(norm(r_B_ddot))
rθ˙2 r\,{\dot{\theta} }^2
r_G_dot = subs(diff(r_G,t), [diff(theta,t), theta], ...
                            [theta_dot, sym('theta','positive')]);
r_G_dot = formula(r_G_dot)
(rθ˙sin(θ)+2r2θ˙cos(θ)sin(θ)7L2r2sin(θ)25rθ˙cos(θ)7) \left(\begin{array}{c} r\,\dot{\theta} \,\sin \left(\theta \right)+\dfrac{2\,r^2 \,\dot{\theta} \,\cos \left(\theta \right)\,\sin \left(\theta \right)}{7\,\sqrt{L^2 -r^2 \,{\sin \left(\theta \right)}^2 } }\\ \dfrac{5\,r\,\dot{\theta} \,\cos \left(\theta \right)}{7} \end{array}\right)
r_G_ddot = subs(diff(r_G,t,2), [diff(theta,t), diff(theta,t,2), theta], ...
                               [theta_dot, theta_ddot, sym('theta','positive')]);
r_G_ddot = formula(r_G_ddot)
(rθ˙2cos(θ)+2r2θ˙2cos(θ)27σ12r2θ˙2sin(θ)27σ1+2r4θ˙2cos(θ)2sin(θ)27σ13/25rθ˙2sin(θ)7)where    σ1=L2r2sin(θ)2 \begin{array}{l} \left(\begin{array}{c} r\,{\dot{\theta} }^2 \,\cos \left(\theta \right)+\dfrac{2\,r^2 \,{\dot{\theta} }^2 \,{\cos \left(\theta \right)}^2 }{7\,\sqrt{\sigma_1 } }-\dfrac{2\,r^2 \,{\dot{\theta} }^2 \,{\sin \left(\theta \right)}^2 }{7\,\sqrt{\sigma_1 } }+\dfrac{2\,r^4 \,{\dot{\theta} }^2 \,{\cos \left(\theta \right)}^2 \,{\sin \left(\theta \right)}^2 }{7\,{\sigma_1 }^{3/2} }\\ -\dfrac{5\,r\,{\dot{\theta} }^2 \,\sin \left(\theta \right)}{7} \end{array}\right)\\ \mathrm{}\\ \textrm{where}\\ \mathrm{}\\ \;\;\sigma_1 =L^2 -r^2 \,{\sin \left(\theta \right)}^2 \end{array}

Substitute with some numerical data

symData = [L, r, theta_dot, theta_ddot, theta];
numData = [(250+100)/1000, 125/1000, 25*2*pi, 0, sym('theta','positive')];
L = (250+100)/1000; r = 125/1000; theta_dot = 25; theta_ddot = 0;

r_A = subs(r_A);
r_A = subs(r_A,theta,sym('theta','positive'));
r_A = formula(r_A);

r_B = subs(r_B);
r_B = subs(r_B,theta,sym('theta','positive'));
r_B = formula(r_B);

r_G = subs(r_G);
r_G = subs(r_G,theta,sym('theta','positive'));
r_G = formula(r_G)
(cos(θ)8249400sin(θ)26475sin(θ)56) \left(\begin{array}{c} -\dfrac{\cos \left(\theta \right)}{8}-\dfrac{2\,\sqrt{\dfrac{49}{400}-\dfrac{ {\sin \left(\theta \right)}^2 }{64} } }{7}\\ \dfrac{5\,\sin \left(\theta \right)}{56} \end{array}\right)
theta = sym('theta','positive');

r_B_dot = subs(r_B_dot);
r_B_ddot = subs(r_B_ddot);

r_A_dot = simplify(subs(r_A_dot));
r_A_ddot = simplify(subs(r_A_ddot));

r_G_dot = simplify(subs(r_G_dot));
r_G_ddot = simplify(subs(r_G_ddot));

r_A_dot_m = norm(r_A_dot);
r_B_dot_m = norm(r_B_dot);
r_G_dot_m = norm(r_G_dot);

r_A_ddot_m = simplify(r_A_ddot(1))
125(8550cos(θ)2+625cos(θ)4+5cos(θ)(25cos(θ)2+171)3/24275)8(25cos(θ)2+171)3/2 \dfrac{125\,{\left(8550\,{\cos \left(\theta \right)}^2 +625\,{\cos \left(\theta \right)}^4 +5\,\cos \left(\theta \right)\,{ {\left(25\,{\cos \left(\theta \right)}^2 +171\right)} }^{3/2} -4275\right)} }{8\,{ {\left(25\,{\cos \left(\theta \right)}^2 +171\right)} }^{3/2} }
r_B_ddot_m = simplify(norm(r_B_ddot))
6258 \dfrac{625}{8}
r_G_ddot_m = simplify(norm(r_G_ddot))
390625(3420cos(θ)2+250cos(θ)4+7cos(θ)(25cos(θ)2+171)3/21710)2(25cos(θ)2+171)3+9765625sin(θ)256 \dfrac{\sqrt{\dfrac{390625\,{ {\left(3420\,{\cos \left(\theta \right)}^2 +250\,{\cos \left(\theta \right)}^4 +7\,\cos \left(\theta \right)\,{ {\left(25\,{\cos \left(\theta \right)}^2 +171\right)} }^{3/2} -1710\right)} }^2 }{ { {\left(25\,{\cos \left(\theta \right)}^2 +171\right)} }^3 }+9765625\,{\sin \left(\theta \right)}^2 } }{56}
figure; hold on; grid on
fplot(r_A_dot_m,[0,2*pi],'LineWidth',2,'Displayname','$|\dot \mathbf{r}_A|$')
fplot(r_B_dot_m,[0,2*pi],'LineWidth',2,'Displayname','$|\dot \mathbf{r}_B|$')
fplot(r_G_dot_m,[0,2*pi],'LineWidth',2,'Displayname','$|\dot \mathbf{r}_G|$')
legend('show',interpreter='latex')
tickRangeDegrees = 0:45:360;
xticks(tickRangeDegrees*pi/180);
xticklabels(cellstr(num2str(tickRangeDegrees')));
xlabel('$\theta$ [degrees]',interpreter='latex')
ylabel('Speed [m/s]',interpreter='latex')
title('Crank-slider speed')
figure; hold on; grid on
fplot(r_A_ddot_m,[0,2*pi],'LineWidth',2,'Displayname','$\ddot r_{Ax}$')
fplot(r_B_ddot_m,[0,2*pi],'LineWidth',2,'Displayname','$|\ddot \mathbf{r}_B|$')
fplot(r_G_ddot_m,[0,2*pi],'LineWidth',2,'Displayname','$|\ddot \mathbf{r}_G|$')
legend('show',interpreter='latex')
tickRangeDegrees = 0:45:360;
xticks(tickRangeDegrees*pi/180);
xticklabels(cellstr(num2str(tickRangeDegrees')));
xlabel('$\theta$ [degrees]',interpreter='latex')
ylabel('Acceleration [$\mathrm{m/s}^2$]',interpreter='latex')
title('Crank-slider acceleration')

% r_A_n = matlabFunction(r_A);
% r_B_n = matlabFunction(r_B);
% r_G_n = matlabFunction(r_G);
% 
% r_A_dot_n = matlabFunction(r_A_dot);
% r_B_dot_n = matlabFunction(r_B_dot);
% r_G_dot_n = matlabFunction(r_G_dot);
% 
% r_A_ddot_n = matlabFunction(r_A_ddot);
% r_B_ddot_n = matlabFunction(r_B_ddot);
% r_G_ddot_n = matlabFunction(r_G_ddot);
% 
% r_A_dot_m_n = matlabFunction(r_A_dot_m);
% r_B_dot_m_n = matlabFunction(r_B_dot_m);
% r_G_dot_m_n = matlabFunction(r_G_dot_m);
% 
% r_A_ddot_m_n = matlabFunction(r_A_ddot_m);
% r_B_ddot_m_n = matlabFunction(r_B_ddot_m);
% r_G_ddot_m_n = matlabFunction(r_G_ddot_m);
% 
% theta = 30*pi/180;
% 
% r_A_i = r_A_n(theta)*1000;
% r_B_i = r_B_n(theta)*1000;
% r_G_i = r_G_n(theta)*1000;
% 
% r_A_dot_i = r_A_dot_n(theta);
% r_B_dot_i = r_B_dot_n(theta);
% r_G_dot_i = r_G_dot_n(theta);
% 
% r_A_ddot_i = r_A_ddot_n(theta);
% r_B_ddot_i = r_B_ddot_n(theta);
% r_G_ddot_i = r_G_ddot_n(theta);
% 
% xmin = -500; xmax = 150;
% ymin = -150; ymax = 150;
% 
% close all
% figure;
% set(gcf,'Visible','on','Color','w','Position',[680   242   872   636])
% subplot(2,2,[1,2])
% hold on; axis equal;
% hB = plot([0,r_B_i(1)],[0,r_B_i(2)],'-ko','MarkerFaceColor','k');
% text(5,-5,'$O$',interpreter='latex')
% plot([xmin,xmax],[0,0],':k')
% 
% hAB = plot([r_A_i(1),r_B_i(1)],[r_A_i(2),r_B_i(2)],'-ko','MarkerFaceColor','k');
% 
% 
% 
% scale_v = 20;
% scale_a = 1;
% hqBv = quiver(r_B_i(1), r_B_i(2), r_B_dot_i(1), r_B_dot_i(2), scale_v, 'b', 'LineWidth',2);
% hqBa = quiver(r_B_i(1), r_B_i(2), r_B_ddot_i(1), r_B_ddot_i(2), scale_a, 'r', 'LineWidth',2);
% hqAv = quiver(r_A_i(1), r_A_i(2), r_A_dot_i(1), r_A_dot_i(2), scale_v, 'b', 'LineWidth',2);
% hqAa = quiver(r_A_i(1), r_A_i(2), r_A_ddot_i(1), r_A_ddot_i(2), scale_a, 'r', 'LineWidth',2);
% hqGv = quiver(r_G_i(1), r_G_i(2), r_G_dot_i(1), r_G_dot_i(2), scale_v, 'b', 'LineWidth',2);
% hqGa = quiver(r_G_i(1), r_G_i(2), r_G_ddot_i(1), r_G_ddot_i(2), scale_a, 'r', 'LineWidth',2);
% 
% htG = text(r_G_i(1)+5,r_G_i(2)-5, '$G$', interpreter='latex');
% hG = plot(r_G_i(1),r_G_i(2),'ko','MarkerFaceColor','k');
% htA = text(r_A_i(1)+5,r_A_i(2)-5, '$A$', interpreter='latex');
% htB = text(r_B_i(1)+5,r_B_i(2)-5, '$B$', interpreter='latex');
% 
% axis([xmin, xmax, ymin, ymax])
% xlabel('$x$', interpreter='latex'); ylabel('$y$', interpreter='latex')
% 
% htitle = title(['$\theta=$',num2str(theta*180/pi),'$^\circ$'], interpreter='latex');
% 
% subplot(2,2,[3]); hold on; grid on
% fplot(r_A_dot_m,[0,2*pi],'LineWidth',2,'Displayname','$|\dot \mathbf{r}_A|$')
% fplot(r_B_dot_m,[0,2*pi],'LineWidth',2,'Displayname','$|\dot \mathbf{r}_B|$')
% fplot(r_G_dot_m,[0,2*pi],'LineWidth',2,'Displayname','$|\dot \mathbf{r}_G|$')
% legend('show','interpreter','latex','Position',[0.17282 0.44296 0.21282 0.02450],...
%     'Orientation','horizontal')
% tickRangeDegrees = 0:45:360;
% xticks(tickRangeDegrees*pi/180);
% xticklabels(cellstr(num2str(tickRangeDegrees')));
% xlabel('$\theta$ [degrees]',interpreter='latex')
% ylabel('Speed [m/s]',interpreter='latex')
% 
% c1 = [0 0.4470 0.7410];
% c2 = [0.8500 0.3250 0.0980];
% c3 = [0.9290 0.6940 0.1250];
% hpAv = plot(theta, r_A_dot_m_n(theta),'o','color',c1,'MarkerFaceColor',c1,'HandleVisibility','off');
% hpBv = plot(theta, r_B_dot_m_n(theta),'o','color',c2,'MarkerFaceColor',c2,'HandleVisibility','off');
% hpGv = plot(theta, r_G_dot_m_n(theta),'o','color',c3,'MarkerFaceColor',c3,'HandleVisibility','off');
% 
% subplot(2,2,[4]); hold on; grid on
% fplot(r_A_ddot_m,[0,2*pi],'LineWidth',2,'Displayname','$\ddot r_{Ax}$')
% fplot(r_B_ddot_m,[0,2*pi],'LineWidth',2,'Displayname','$|\ddot \mathbf{r}_B|$')
% fplot(r_G_ddot_m,[0,2*pi],'LineWidth',2,'Displayname','$|\ddot \mathbf{r}_G|$')
% legend('show','interpreter','latex','Position',[0.60686 0.44296 0.21282 0.02450],...
%     'Orientation','horizontal')
% tickRangeDegrees = 0:45:360;
% xticks(tickRangeDegrees*pi/180);
% xticklabels(cellstr(num2str(tickRangeDegrees')));
% xlabel('$\theta$ [degrees]',interpreter='latex')
% ylabel('Acceleration [$\mathrm{m/s}^2$]',interpreter='latex')
% hpAa = plot(theta, r_A_ddot_m_n(theta),'o','color',c1,'MarkerFaceColor',c1,'HandleVisibility','off');
% hpBa = plot(theta, r_B_ddot_m_n(),'o','color',c2,'MarkerFaceColor',c2,'HandleVisibility','off');
% hpGa = plot(theta, r_G_ddot_m_n(theta),'o','color',c3,'MarkerFaceColor',c3,'HandleVisibility','off');
% 
% 
% for theta = linspace(0,2*pi, 100)
%     if theta == 0
%         exportgraphics(gcf,"crank-slider.gif","Append",false)
%     end
%     r_A_i = r_A_n(theta)*1000;
%     r_B_i = r_B_n(theta)*1000;
%     r_G_i = r_G_n(theta)*1000;
% 
%     r_A_dot_i = r_A_dot_n(theta);
%     r_B_dot_i = r_B_dot_n(theta);
%     r_G_dot_i = r_G_dot_n(theta);
% 
%     r_A_ddot_i = r_A_ddot_n(theta);
%     r_B_ddot_i = r_B_ddot_n(theta);
%     r_G_ddot_i = r_G_ddot_n(theta);
% 
%     set(hAB,'XData', [r_A_i(1),r_B_i(1)], 'Ydata', [r_A_i(2),r_B_i(2)])
%     set(hB,'XData', [0,r_B_i(1)], 'Ydata', [0,r_B_i(2)])
%     set(hG, 'XData', r_G_i(1), 'YData', r_G_i(2))
%     set(htG,'Position',[r_G_i(1)+5,r_G_i(2)-5])
%     set(htA,'Position',[r_A_i(1)+5,r_A_i(2)-5])
%     set(htB,'Position',[r_B_i(1)+5,r_B_i(2)-5])
%     set(htitle,'String', ['$\theta=$',sprintf('%0.0f',theta*180/pi),'$^\circ$'])
%     set(hqBv,'XData',r_B_i(1),'Ydata',r_B_i(2),'Udata',r_B_dot_i(1),'VData',r_B_dot_i(2))
%     set(hqBa,'XData',r_B_i(1),'Ydata',r_B_i(2),'Udata',r_B_ddot_i(1),'VData',r_B_ddot_i(2))
%     set(hqAv,'XData',r_A_i(1),'Ydata',r_A_i(2),'Udata',r_A_dot_i(1),'VData',r_A_dot_i(2))
%     set(hqAa,'XData',r_A_i(1),'Ydata',r_A_i(2),'Udata',r_A_ddot_i(1),'VData',r_A_ddot_i(2))
%     set(hqGv,'XData',r_G_i(1),'Ydata',r_G_i(2),'Udata',r_G_dot_i(1),'VData',r_G_dot_i(2))
%     set(hqGa,'XData',r_G_i(1),'Ydata',r_G_i(2),'Udata',r_G_ddot_i(1),'VData',r_G_ddot_i(2))
% 
%     set(hpAv,'XData',theta,'Ydata',r_A_dot_m_n(theta))
%     set(hpBv,'XData',theta,'Ydata',r_B_dot_m_n(theta))
%     set(hpGv,'XData',theta,'Ydata',r_G_dot_m_n(theta))
%     set(hpAa,'XData',theta,'Ydata',r_A_ddot_m_n(theta))
%     set(hpBa,'XData',theta,'Ydata',r_B_ddot_m_n())
%     set(hpGa,'XData',theta,'Ydata',r_G_ddot_m_n(theta))
% 
%     exportgraphics(gcf,"crank-slider.gif","Append",true)
%     drawnow
% end