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

 r(t) = (x(t), y(t), z(t))     (1)

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

 v(t) = (vρ(t),vθ(t))     (2)

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 magnitude1 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

 | r(t) |     (3)

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

 r(t) = | r(t) |     (4)

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

r(t) =
 r(t) | r(t) |
(5)

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.2 For unit vectors we always have

 | r(t) | = 1     (6)

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:

 F = m a     (7)

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

 F = m a     (8)

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

 ∀ t: F(t) = m a(t)     (9)

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

 d2r dt2
= a     (10)

where a and r are functions of time that return vectors, and d2r/dt2 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:

 dr dt
= v     (11)

 dv dt
= a     (12)

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/s2, 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 v2 and directed oppositely to v:

 Fd = − b v2 v     (13)

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.3For 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
```
Exercise 1   Increase the mass of the skydiver, and confirm that terminal velocity increases. This relationship is the source of the intuition that heavy objects fall faster; in air, they do!
Exercise 2   The terminal velocity of a 75 kg skydiver is about 60 m/s, which is about 130 mph. If you hit the ground at that speed, you would almost certainly be killed. That’s where parachutes come in.

Modify acceleration so that after 30 seconds of free-fall the skydiver deploys a parachute, which (almost) instantly increases the drag constant to 2.7. Hint: see Section 11.9.

What is the terminal velocity now? How long (after deployment) does it take to reach the ground?

## 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).4

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.

Exercise 3   Extract the x and y components of position, plot the trajectory of the projectile, and estimate the horizontal distance traveled.

## 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, vx and vy, 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)):

Fd = −
 1 2
ρ v2 Cd A   v     (14)

where ρ, Cd, and A are scalars representing fluid density (1.3   kg/m3 for air at sea level), a dimensionless drag coefficient (0.3 is reasonable for a baseball), and cross-sectional area (0.0042   m2 for a baseball).5

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

u | =
 ux2 + uy2
(15)

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 Fd we just need to recall our values for ρ, Cd, 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/s2.

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.

Exercise 4   Suppose the baseball lands in a lake, and breaks the water’s surface with velocity
 v0 = (25, −25)    m/s     (16)

How long will it take to resurface?

Consider the baseball to be fully submerged when it enters the water. You will need to account for the gravitational force and the drag force (given in equation (14)). You will also need to consider the force due to “bouyancy” which is given by Archimedes’ principle:

 Fb = ρ V g   ĵ     (17)

where ρ is the density of the fluid (for water this is 1000   kg / m3), V is the volume of the fluid displaced (equal to the volume of the baseball—take this to be 0.000213   m3), and g the familiar gravitational constant; ĵ is the unit vector that points straight up:

 ĵ = (0, 1)     (18)

## 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/m3)—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 resource6).

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);
```
Exercise 5   How would you modify event to stop when the height of the projectile falls through 3m? If the human projectile falls on a net at a height of 3m, what horizontal distance will she have travelled when she hits the net?

## 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

Exercise 6

In the absense of wind the flight of a baseball is governed by three forces: gravity, drag due to air resistance, and Magnus force due to spin. If we ignore Magnus force, the path of the baseball stays in a plane, so we can model it as a projectile in two dimensions.

As presented in Section 11.6, a simple model for the drag force on a baseball is:

Fd = −
 1 2
ρ v2 Cd A · v     (19)

where Fd is a vector that represents the force on the baseball due to drag, ρ is the density of air (1.3 kg/m3 at sea level), Cd is the drag coefficient (0.3 is a reasonable choice), A is the baseball’s cross sectional area (0.0042 m2), and v is the baseball’s velocity. The baseball has a mass of 0.145 kg.

• Write a function that takes as input variables the initial velocity, the initial height, and the launch angle of the baseball, uses ode45 to compute the trajectory, and returns the range (horizontal distance in flight until the ball impacts the ground) as an output variable.
• Use your implemented function to display the horizontal distance profile with initial velocities from 30m/s to 50m/s at a launch angle of 45 degrees and an initial height of 1 meter. Describe the relation between the horizontal distance and the initial velocity.
• Fix the initial velocity at 48m/s and initial height at 1 meter, and display the horizontal distance profile with launch angles from 20 to 60 degrees. Describe the relation between the horizontal distance and the launch angle. Approximate the optimal launch angle for this problem (for present purposes, the launch angle that gives the furthest horizontal distance to the ball impacting the ground). You can alter the plotted launch angles if needed or desired.

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
3
Notice that equation (13) makes use of the convention v = | v |. See Section 11.1
4
5
6

#### Are you using one of our books in a class?

We'd like to know about it. Please consider filling out this short survey.