Applied Mechanics

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.

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

Vector notation

The position of a bodies center of gravity as a function of time $t$ can be expressed using the position vector $\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 $t$. The components are called coordinate functions and $\bm r\left(t\right)$ is known as the parameter curve with the parameter $t$. The parameter curve is also known as the curve or path.

Velocity

The mean velocity between two positions $\bm r\left(t\right)$ and $\bm r\left(t+\Delta t\right)$ can be described as $\bar{\bm v} =\frac{\Delta \bm 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

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

The dot notation in $\dot{\bm 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{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, $\frac{d\dot{r} }{dt}$. The average acceleration at time $t$ can be formulated as

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

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

solving for the time increment $dt$, we get

from which we get

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\;\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 $\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., $\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 $\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$.

The acceleration (vector) is the time derivative of the velocity (vector) $\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 $\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 recommended way

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

To get between the different quantities position, $\bm r\left(t\right)$, velocity $\dot{\bm r} \left(t\right)$ and acceleration $\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 $\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, $\dot{\bm r}$ on the way down to the final position $\bm r$.

Newtonian physics is formulated as the kinetic equation

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.

• We formulate an ODE which needs to be (at least piece-wise) valid for all $t$ and apply initial conditions at the beginning of the analysis and/or at some other known state as boundary conditions (BC) to fix the unknown constants, where the number of constants depend on the order of the differential equation.

• In this approach, the time $t$ is the central parameter that ties the quantities ($\bm r$, $\dot{\bm r}$ and $\ddot{\bm r}$) together, which is why it often pops out in the solution without being maybe being explicitly asked for. This is the nature of Newtonian mechanics, it is inherently transient, meaning the state of the body is known for every time step. As we shall see in upcoming chapters, there are other methods which can greatly simplify calculations and modeling if the dependency of time is not important, e.g., equation (7). Note however that a complete dynamic analysis requires the solution of the ODE.

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., $\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, $\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.

We shall here derive the most common coordinate systems and their base vectors such that this projection can be done.

• Cartesian, our typical $\textrm{xyz}$-system, here denoted, ${\bm e}_x ,{\;\bm e}_y ,{\;\bm e}_z$, with ${\bm e}_x =\left\lbrack 1,0,0\right\rbrack^T$, ${\bm e}_y =\left\lbrack 0,1,0\right\rbrack ^T$ and ${\bm e}_z =\left\lbrack 0,0,1\right\rbrack ^T$.

• Polar or circular (2D) or cylindrical (3D), ${\bm e}_r ,{\bm e}_\theta ,{\bm e}_z$.

• Natural (tangent, normal and bi-normal), ${\bm e}_t ,{\bm e}_n ,{\bm e}_b$.

Sometimes the bases can be denoted simply $\bm i$, $\bm j$ and $\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 ${\bm e}_t$ is always pointing in the parameter direction, tangential to the path and is defined as the direction of the velocity

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

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,

which defines a right-hand system.

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

Polar coordinates

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

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

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., $\bm e\cdot \bm e=1$. Now we can take the time derivative of this expression $\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 $|\dot{\bm e} |\not= 1$.

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

normalize = @(v)v/sqrt(v(:).'*v(:));
syms t positive
syms theta(t)
r = [cos(theta(t)),sin(theta(t))].'

e_r = simplify(normalize(r))

e_r_dot = simplify(diff(e_r,t))

e_theta = subs(simplify(normalize(e_r_dot)), diff(theta,t), 'theta_dot')

Polar coordinates $\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 $\dot{\bm r}$ and $\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.

We can express the position vector using the distance and base

where $r$ 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

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

where

is recognized as the angular velocity.

Acceleration in Polar Coordinates

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

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

Natural coordinates $\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.

Curvature

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

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=\frac{ds}{dt}$ and using the arc length relation we get

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

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

Thus we arrive at the relation

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

Acceleration in Natural coordinates

and with (36) we get

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 $|\dot{\bm r} \left(t\right)|$ such as

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

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

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 $t$. We note that the Natural coordinates are mostly used for this type of analysis.

The torsion of a curve, $\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

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

clear
syms R t positive
r = R*[cos(t), sin(t), 0].'

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)

tau(r)

Osculating circle

Finally we visualize the osculating circle on its journey accompanying a particle moving along the path $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]

$$\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)].'

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

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

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_{\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., $|{\bm e}_r |=|{\bm e}_{\theta } |=1$. We recall that the dot product is just a matrix multiplication

r_dot_rtheta = simplify([rr_dot.'*e_r, rr_dot.'*e_theta])

r_ddot_rtheta = simplify([rr_ddot.'*e_r, rr_ddot.'*e_theta])

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 $r$ to zero

r_dot = 0; r_ddot = 0;
v = subs(r_dot_rtheta)

a = subs(r_ddot_rtheta)

To summarize planar circular motion: with a constant $r$, we have $\dot{r} =\ddot{r} =0$ and $\bm r\bot \dot{\bm r}$ since $\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{\dot{\theta} }^2 {\bm e}_r$, it is directed in the negative direction of ${\bm 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.

Note that we can express bases using the cross product

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

Another easy way is to explicitly rotate the base ${\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 $B$ as well as the relative position of particle $A$ with respect to $B$, that is ${\bm r}_{A/B}$ (Can be read: A seen from B) we can express the absolute position of $A$ as

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

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

and

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}^{\circ }$ in the north-east direction, although to the passengers it appears that the JAS is moving away at an angle of ${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

e_BA = [-cosd(60); sind(60)]

e_BA = 2x1    
   -0.5000
    0.8660

V_BA = v_BA*e_BA

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

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

[v_B, v_BA] = solve(eqn,[v_B,v_BA])

v_B = vpa(v_B,5)

v_BA = vpa(v_BA,5)

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:

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

Solve for the constants, $c_i$, using the least squares method:

syms t
n = 2;
t.^[0:n]

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]

to get the velocity, here using 4 significant numbers

v = vpa(t.^[0:n]*c,4)

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)

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

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

Example 3 - Rectilinear motion, piecewise

A particle moves along the $x$-axis with an initial velocity $v_x =50\;\mathrm{m}/\mathbf{s}$ at the origin when $t=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 $a_x =-10\;\mathrm{m}/{\mathrm{s}}^2$. Compute the velocity and the $x$-coordinate of the particle for the conditions of $t=8$s and $t=12$s and find the maximum positive $x$-coordinate reached by the particle.

Solution:

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

v = diff(x,t)

IV = [v(4) == 50, x(4) == 50*4]

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)

v(t) = diff(x,t)

v(8)

v(12)

t_max = solve(v==0, t)

x(t_max)

v(t) = piecewise(0 <= t < 4, 50, ...
                 4 < t,      v)

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 $x$-axis and has an initial velocity $v_0 =2$m/s in the s-direction as it crosses the mid-position where $s=0$ and $t=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=-k^2 s$, where $k=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

v = diff(x,t)

BV = [v(0) == v0, x(0) == 0]

x(t) = dsolve([DE,BV])

x = simplify(x)

v = diff(x,t)

a = diff(v,t)

x = subs(x,[v0, k],[2,4])

v = subs(v,[v0, k],[2,4])

a = subs(a,[v0, k],[2,4])

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\pi$is 8 m. But the signed area is zero, since the sliders position after $2\pi$ seconds is $x=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\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\in \left\lbrack 0,30\right\rbrack \mathrm{s}$!

syms t
x = piecewise(0 <= t < 10, t^2,...
              t >= 10, 20*t-100)

Now we just differentiate

x_dot = diff(x,t)

x_ddot = diff(x,t,2)

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

• Draw the distance and velocity curve

• How long does is take to reach top speed and how far has the car traveled?

• What is the 0-100 km/h time and distance?

Solution:

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

clear
syms x(t)
x_dot = diff(x,t)

v_1 = 80; k = 0.1;
DE = diff(x,t,2) == k*(v_1-x_dot)

IV = [x_dot(0)==0, x(0)==0]

x = dsolve(DE,IV)

v = diff(x,t)

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)

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_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 $v_x =50-16t$ and $y=100-4t^2$, where $v_x$ is in meter per second, $y$ is in meters and $t$ is in seconds. It is also known that $x=0$ when $t=0$. Plot the path of the particle and determine its velocity and acceleration when the position $y=0$ is reached.

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

y = 100 - 4*t^2

DE = diff(x,t) == 50 - 16*t

BV = x(0) == 0

x = dsolve(DE,BV)

t_y0 = solve(y==0,t)

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].'

r_dot = diff(r,t)

r_ddot = diff(r,t,2)

Velocity when $y=0$

r_dot_y0 = subs(r_dot, t, t_y0)