Buy this book at Amazon.com

This HTML version of the book is provided as a convenience, but some math equations are not translated correctly. The PDF version is more reliable.

Chapter 11  Second-order Systems

11.1  Spatial vectors

Before we can start modelling second-order systems, we need to discuss vectors. The word “vector” means different things to different people. In MATLAB, a vector is a matrix that has either one row or one column. So far we have used MATLAB vectors to represent

sequences:

A sequence is a set of values identified by integer indices; it is natural to store the elements of the sequence as elements of a MATLAB vector.

state vectors:

A state vector is a set of values that describes the state of a physical system. When you call ode45, you give it initial conditions in a state vector. Then when ode45 calls your rate function, it gives you a state vector.

discrete maps:

If you have two vectors with the same length, you can think of them as a mapping from the elements of one vector to the corresponding elements of the other. For example, in Section 9.5, the results from ode45 are vectors, T and Y, that represent a mapping from the time values in T to the population values in Y.

In this chapter we will see another use of MATLAB vectors: representing spatial vectors. A spatial vector is a geometric object that can be used to represent a multidimensional physical quantity like position, velocity, acceleration, or force.

These quantities cannot be described with a single number because they contain multiple components (or equivalently—as we’ll see shortly —because they have both a magnitude and a direction). For example, in three-dimensional Cartesian coordinates, it takes three numbers to specify a position in space; these numbers are usually called x, y, and z coordinates. Thus we can write a position vector r(t), representing an object’s location at time t, as

As another example, in two-dimensional polar coordinates, it takes two numbers to specify a position in space: a radial distance ρ and an angle θ. Here we can write a two-dimensional velocity vector v(t), representing an object’s velocity at time t, as

Notice that both vectors r(t) and v(t) are presented in bold, indicating that they are vectors.

An equivalent way of describing a spatial vector is with a magnitude¹ and a direction. If you imagine the vector r(t) from equation (1) as an arrow starting from the origin and terminating at the point (x(t),y(t),z(t)), then the magnitude of r(t) is how long the arrow is and the direction of r(t) is where the arrow points.

The magnitude of r(t) is usually written as

but when there’s little room for confusion, it’s often more convenient to write

Exactly how you calculate the magnitude of a vector from its components depends on the coordinate system.

The direction of r(t) is conveniently described by

where r(t) is a unit vector, as indicated by its hat. Notice this equation (5) contains the product of a vector r(t) and a scalar 1/| r(t) |; this operation is called scalar multiplication, wherein the scalar multiplies (or scales) the magnitude of the vector without changing its direction.² For unit vectors we always have

which is why unit vectors are so convenient for representing directions.

11.2  Newtonian motion

Newton’s second law of motion is often written like this:

where F is the net force acting on an object, m is the mass of the object, and a is the acceleration of the object. In the simple case of one-dimensional motion, F and a are scalars, letting us write

but in general they are vectors. If F and a vary in time, then they can be thought of as functions that return vectors; that is, F is a function and the result of evaluating F(t) is a vector that describes the net force at time t. So a more explicit way to write Newton’s second law is

The arrangement of this equation (9) suggests that if you know m and a you can compute force, which is true, but in most physical simulations it is the other way around: based on a physical model, you know F and m, and you compute a.

So if we know acceleration, a, as a function of time, how do we find the position of the object, r? Well, we know that acceleration is the second derivative of position, so we can write a differential equation

where a and r are functions of time that return vectors, and d²r/dt² is the second time derivative of r.

Because equation (10) includes a second derivative, it is a second-order ODE. ode45 can’t solve the equation this form, but by introducing a new variable, v, for velocity, we can rewrite it as a system of first-order ODEs:

Equation (11) says that the first derivative of r is v; equation (12) says that the first derivative of v is a.

11.3  Freefall

Let’s start with a simple example: an object in freefall in a vacuum (where there’s no air resistance). Near the surface of the earth, the acceleration of gravity is g = −9.8 m/s², where the minus sign indicates that gravity pulls down.

If the object falls straight down, we can describe its position with a scalar value y representing altitude. So this will be a one-dimensional problem, at least for now.

Here is a rate function we can use with ode45 to solve this problem:

function res = freefall(t, X)
    y = X(1);      % the first element is position
    v = X(2);      % the second element is velocity

    dydt = v;
    dvdt = acceleration(t, y, v);

    res = [dydt; dvdt];    % pack the results in a column vector
end

function res = acceleration(t, y, v)
    g = 9.8;      % gravity in m/s^2
    res = -g;
end

The first function, freefall, is the rate function. It gets t and X as input variables, where the elements of X are understood to be precisely the position and velocity of the object. It returns a column vector that contains dydt and dvdt, which are velocity and acceleration, respectively.

We’re given the velocity as the second element of X, so freefall simply assigns this value to dydt. freefall obtains acceleration by calling the second function, acceleration, which takes time, position, and velocity as input variables. In this example, the net acceleration is a constant, so the input variables are unnecessary, but in Section 11.4 we’ll need them.

Here’s how to run ode45 with this rate function:

>> ode45(@freefall, [0, 30], [4000, 0])

As always, the first argument is the function handle, the second is the time interval (30 seconds) and the third is the initial condition: in this case, the initial altitude is 4000 meters and the initial velocity is 0. You can think of our object as a skydiver jumping out of an airplane at about 12,000 feet.

Here’s what the result looks like:

The bottom line shows velocity starting at zero and dropping linearly. The top line shows position starting at 4000 m and dropping parabolically (but remember that this parabola is a function of time, not a ballistic trajectory).

Notice that ode45 doesn’t know where the ground is, so the skydiver keeps going through zero into negative altitude. We will address this issue in Section 11.9.

11.4  Air resistance

To make this simulation more realistic, we can add air resistance. For large objects moving quickly through air, the force due to air resistance, called “drag,” is proportional to v² and directed oppositely to v:

where b is a drag constant that depends on the density of air, the cross-sectional area of the object and the surface properties of the object.³For this problem, let’s say that b = 0.2.

To convert from force to acceleration, we have to know mass, so let’s say that the skydiver (with equipment) weighs 75 kg.

Here’s a version of acceleration that takes air resistance into account (note that you don’t have to make any changes in the freefall function):

function res = acceleration(t, y, v)
    a_g = -9.8;              % acceleration due to gravity in m/s^2
    b = 0.2;                 % drag constant
    m = 75;                  % mass in kg
    F_d = b * v^2;           % drag force in N
    a_d = F_d / m;           % drag acceleration in m/s^2
    res = a_g + a_d;         % total acceleration
end

The sign of the drag force (and “drag acceleration”) is positive here because the object is falling (i.e., the velocity is in the −y direction). The output variable is simply the sum of the gravity and drag accelerations. Be careful when you are working with forces and accelerations; make sure you only add forces to forces or accelerations to accelerations. In my code, I use comments to remind myself what units the values are in. That helps me avoid nonsense like adding forces to accelerations.

Here’s what the result looks like with air resistance:

Big difference! With air resistance, velocity increases until the drag acceleration equals g; after that, velocity is a constant, known as “terminal velocity,” and position decreases linearly (and much more slowly than it would in a vacuum). To examine the results more closely, we can assign them to variables:

>> [T, M] = ode45(@freefall, [0, 30], [4000, 0]);

And then read the terminal position and velocity:

>> M(end,1)

ans = 2.4412e+03          % altitude in meters

>> M(end,2)

ans = -60.6143            % velocity in m/s

11.5  Two dimensions

So far we have used ode45 for a system of first-order equations and for a single second-order equation. The next logical step is a system of second-order equations, and the next logical example is a projectile. A “projectile” is an object propelled through space, usually toward, and often to the detriment of, a target.

If a projectile stays in a plane, we can think of the system as two-dimensional, with x representing the horizontal distance traveled and y representing the height or altitude. So now instead of a skydiver, think of a circus performer being fired out of a cannon.

According to Wikipedia, the record distance for a human cannonball is 59.05 meters (almost 194 feet).⁴

Here is a general framework for computing the trajectory of a projectile in two dimensions using ode45:

function res = projectile(t, W)
    R = W(1:2);
    V = W(3:4);

    dRdt = V;
    dVdt = acceleration(t, R, V);

    res = [dRdt; dVdt];
end

function res = acceleration(t, R, V)
    g = 9.8;             % gravity in m/s^2
    res = [0; -g];
end

The second argument of the rate function is understood to be a vector, W, with four elements. The first two are assigned to R, which represents position; the last two are assigned to V, which represents velocity. Both R and V have two elements representing the x and y components of the projectile’s motion.

The result from acceleration is also a vector; ignoring air resistance (for now), the acceleration in the x direction is 0; in the y direction it’s −g. Other than that, this code is similar to what we saw in Section 11.3.

If we launch the human projectile from an initial height of 3 meters, with velocities 40 m/s and 30 m/s in the x and y directions, the ode45 call looks like this:

ode45(@projectile, [0,10], [0, 3, 40, 30]);

And the result looks like this:

You might have to think a little to figure out which line is which. It looks like the flight time is about 6 seconds.

11.6  Spacial vectors revisited

Let’s talk about spacial vectors again, this time by way of example. Suppose that you are given the velocity of a baseball in the form of a MATLAB vector with two elements, v_x and v_y, which are the components of velocity in the x and y directions.

>> V = [30; 40];       % velocity in m/s

Also suppose you are asked to compute the total acceleration of the ball due to drag and gravity. Wikipedia specifies the drag force as follows (cf. equation (13)):

where ρ, C_d, and A are scalars representing fluid density (1.3   kg/m³ for air at sea level), a dimensionless drag coefficient (0.3 is reasonable for a baseball), and cross-sectional area (0.0042   m² for a baseball).⁵

To use equation (14), we need both the magnitude and the direction of V. The magnitude of a spacial vector u in two-dimensional Cartesian coordinates is

You could compute the magnitude using equation (15) with hypotenuse from Section 6.5, or you could use the MATLAB function norm (norm is another name for magnitude); let’s use norm:

>> v = norm(V)

v = 50

To find the unit vector Vhat corresponding to V, we simply need to divide V by its magnitude:

>> Vhat = V / v

Vhat =  0.6  0.8

Let’s confirm that the norm of Vhat is unity:

>> norm(Vhat)

ans = 1

Good. Now, to compute F_d we just need to recall our values for ρ, C_d, and A; and then follow equation (14):

>> rho = 1.3;          % kg/m^3
>> Cd = 0.3;           % dimensionless
>> A = 0.0042;         % m^2

>> Fd = - 1/2 * rho * v^2 * Cd * A * Vhat  % N

Fd =  -1.2  -1.6

Remembering Newton’s second law (and taking the mass m of a baseball to be 0.145   kg), we can compute “drag acceleration”:

>> m = 0.145;          % kg
>> Ad = Fd / m         % m/s^2

Ad =  -8.5  -11.3

To represent acceleration due to gravity, we create a vector with two components:

>> Ag = [0; -9.8];

The x component of acceleration due to gravity is 0; the y component is −9.8   m/s².

Finally, we compute the total acceleration by summing Ad and Ag:

>> At = Ad + Ag

At =  -8.5  -1.5

One nice thing about this computation is that we didn’t have to think much about the components of the vectors. By treating spatial vectors as basic quantities, we can express complicated computations concisely.

11.7  Putting it all together

Let’s return to our two-dimensional simulation of the human cannonball from Section 11.5. Our rate function, projectile, looked like this:

function res = projectile(t, W)
    R = W(1:2);
    V = W(3:4);

    dRdt = V;
    dVdt = acceleration(t, R, V);

    res = [dRdt; dVdt];
end

Suppose we want to model the effect of drag on our cannonball’s flight. Luckily, we already did most of the work for this in Section 11.6. We know that we can represent the acceleration due to gravity with the vector

>> Ag = [0; -9.8];

and we know that we can represent the acceleration due to drag with the vector

>> Fd = - 1/2 * rho * v^2 * Cd * A * Vhat;
>> Ad = Fd / m;

We also know how to compute v and Vhat given V. But what constants should we use? Sure, rho is fairly obvious—we can take rho to be the density of air at sea level (1.3   kg/m³)—but what about Cd, A, and m? I want you to figure out these constants: research is as much a part of physical modelling as math and programming are—and it’s the part that’s most fun.

Visualize a human cannonball. If you can’t, then look one up on YouTube. Look up drag coefficients to help estimate the drag coefficient Cd (the “drag coefficient” Wikipedia page is a good resource⁶).

Now you need to decide who to fire out of the cannon. What’s their waist size? What’s their weight? Use these to estimate A and m.

Once you decide on these constants, you can use the following acceleration function to test your model (together with projectile above):

function res = acceleration(t, R, V)
    g = 9.8;             % gravity in m/s^2
    rho = 1.3;           % air density in kg/m^3
    %Cd = ?
    %A = ?
    %m = ?

    Ag = [0; -g];

    v = norm(V);
    Vhat = V/v;
    Fd = - 1/2 * rho * v^2 * Cd * A * Vhat;
    Ad = Fd / m;

    res = Ag + Ad;
end

Use the following to run the simulation:

>> T = [0,10];            % time-span; change this if needed
>> W0 = [0, 3, 40, 30];   % initial position and velocity
                          % think: is this reasonable?

>> ode45(@projectile, T, W0)

How did it go? Were your results reasonable? Why or why not? Consider our initial discussion of modeling from Section 8.2 . Is this model a good choice to simulate a human cannonball?

11.8  What could go wrong?

What could go wrong? Well, vertcat for one. To explain what that means, I’ll start with concatenation, which is the operation of joining two matrices into a larger matrix. “Vertical concatenation” joins the matrices by stacking them on top of each other; “horizontal concatenation” lays them side by side.

Here’s an example of horizontal concatenation with row vectors:

>> x = 1:3

x = 1     2     3

>> y = 4:5

y = 4     5

>> z = [x, y]

z = 1     2     3     4     5

Inside brackets, the comma operator performs horizontal concatenation. The vertical concatenation operator is the semi-colon. Here is an example with matrices:

>> X = zeros(2,3)

X =  0     0     0
     0     0     0

>> Y = ones(2,3)

Y =  1     1     1
     1     1     1

>> Z = [X; Y]

Z =  0     0     0
     0     0     0
     1     1     1
     1     1     1

These operations only work if the matrices are the same size along the dimension where they are glued together. If not, you get:

>> a = 1:3

a = 1     2     3

>> b = a'

b =  1
     2
     3

>> c = [a, b]
Error using horzcat
Dimensions of matrices being concatenated are not
consistent.

>> c = [a; b]
Error using vertcat
Dimensions of matrices being concatenated are not
consistent.

In this example, a is a row vector and b is a column vector, so they can’t be concatenated in either direction.

Reading the error messages, you probably guessed that horzcat is the function that performs horizontal concatenation, and likewise with vertcat and vertical concatenation.

These operations are relevant to our rate function projectile from Section 11.5 because of the last line, which packs dRdt and dVdt into the output variable:

function res = projectile(t, W)
    R = W(1:2);
    V = W(3:4);

    dRdt = V;
    dVdt = acceleration(t, R, V);

    res = [dRdt; dVdt];
end

As long as both dRdt and dVdt are column vectors, the semi-colon performs vertical concatenation, and the result is a column vector with four elements. But if either of them is a row vector, that’s trouble.

ode45 expects the result from projectile to be a column vector, so if you are working with ode45, it is probably a good idea to make everything a column vector.

In general, if you run into problems with horzcat and vertcat, use size to display the dimensions of the operands, and make sure you are clear on which way your vectors go.

11.9  ODE Events

Normally when you call ode45 you have to specify a start time and an end time. But in many cases, you don’t know ahead of time when the simulation should end (this was an issue in Section 11.3’s freefall simulation). Fortunately MATLAB provides a mechanism for dealing with this problem. The bad news is that it is a little awkward. Here’s how it works:

1. Before calling ode45 you use odeset to create an object called options that contains values that control how ode45 works:

   options = odeset('Events', @event);
   

   In this case, the name of the option is Events and the value is a function handle. When ode45 runs, it will invoke event after each timestep. You can call this function anything you want, but a name like event is conventional.

2. The function you provide has to take the same input variables as your rate function. For example, here is an event function that would work with projectile from Section 11.5

   function [value,isterminal,direction] = event(t,X)
       value = X(2);        % Extract the current height
       isterminal = 1;      % Stop if height crosses zero
       direction = -1;      % if the height is decreasing
   end
   

   event returns three output variables:

   The first output variable, value, determines whether an event can occur. In this case value gets the second element of X, which is understood to be the height of the projectile. An “event” can occur when value passes through 0, but whether it does occur depends on the third output variable, direction.

   The second output variable, isterminal, determines what happens when the event occurs. If isterminal=1, the event is “terminal” and the simulation stops. If isterminal=0, the simulation continues, but ode45 does some additional work to make sure that the solution in the vicinity of the event is accurate, and that one of the estimated values in the result is at the time of the event.

   The third output variable, direction, determines whether an event occurs when value is increasing (direction=1), decreasing (direction=-1, or both direction=0.

3. When you call ode45, you pass options as a fourth argument:

   ode45(@projectile, [0,10], [0, 3, 40, 30], options);
   

11.10  Glossary

spatial vector:

A value that represents a multidimensional physical quantity like position, velocity, acceleration or force.

unit vector:

A vector with norm 1, used to indicate direction.

norm:

The magnitude of a vector. Sometimes called “length,” but not to be confused with the number of elements in a MATLAB vector.

concatenation:

The operation of joining two matrices end-to-end to form a new matrix.

11.11  Exercises

1

Magnitude is also called “length” but I will avoid that term because it gets confused with the length function, which returns the number of elements in a MATLAB vector.

2

See https://en.wikipedia.org/wiki/Scalar_multiplication.

3

Notice that equation (13) makes use of the convention v = | v |. See Section 11.1

4

https://en.wikipedia.org/wiki/Human_cannonball

5

https://en.wikipedia.org/wiki/Drag_(physics).

6

https://en.wikipedia.org/wiki/Drag_coefficient

Are you using one of our books in a class?

Think Bayes

Think Python 2e

Think Stats 2e

Think Complexity