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
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 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
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.2 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 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:
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:
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)):
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
>> 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
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:
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:
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:
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?
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:
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.
Are you using one of our books in a class?We'd like to know about it. Please consider filling out this short survey.