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 12 Vectors as vectors12.1 What’s a vector?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 value that represents a multidimensional physical quantity like position, velocity, acceleration or force^{1}. These quantities cannot be described with a single number because they contain multiple components. For example, in a 3dimensional Cartesian coordinate space, it takes three numbers to specify a position in space; they are usually called x, y and z coordinates. As another example, in 2dimensional polar coordinates, you can specify a velocity with two numbers, a magnitude and an angle, often called r and θ. It is convenient to represent spatial vectors using MATLAB vectors because MATLAB knows how to perform most of the vector operations you need for physical modeling. For 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] And suppose you are asked to compute the total acceleration of the ball due to drag and gravity. In math notation, the force due to drag is
where V is a spatial vector representing velocity, v is the magnitude of the velocity (sometimes called “speed”), and V is a unit vector in the direction of the velocity vector. The other terms, ρ, A and C_{d}, are scalars. The magnitude of a vector is the square root of the sum of the squares of the elements. You could compute it with hypotenuse from Section 5.5, or you could use the MATLAB function norm (norm is another name^{2} for the magnitude of a vector): >> v = norm(V) v = 50 V is a unit vector, which means it should have norm 1, and it should point in the same direction as V. The simplest way to compute it is to divide V by its own norm. >> Vhat = V / v Vhat = 0.6 0.8 Then we can confirm that the norm of V is 1: >> norm(Vhat) ans = 1 To compute F_{d} we just multiply the scalar terms by V. Fd =  1/2 * C * rho * A * v^2 * Vhat Similarly, we can compute acceleration by dividing the vector F_{d} by the scalar m. Ad = Fd / m To represent the acceleration of gravity, we create a vector with two components: Ag = [0; 9.8] The x component of gravity is 0; the y component is −9.8 m/s^{2}. Finally we compute total acceleration by adding vector quantities: A = Ag + Ad; 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 complex computations concisely. 12.2 Dot and cross productsMultiplying a vector by a scalar is a straightforward operation; so is adding two vectors. But multiplying two vectors is more subtle. It turns out that there are two vector operations that resemble multiplication: dot product and cross product. The dot product of vectors A and B is a scalar:
where a is the magnitude of A, b is the magnitude of B, and θ is the angle between the vectors. We already know how to compute magnitudes, and you could probably figure out how to compute θ, but you don’t have to. MATLAB provides a function, dot, that computes dot products. d = dot(A, B) dot works in any number of dimensions, as long as A and B have the same number of elements. If one of the operands is a unit vector, you can use the dot product to compute the component of a vector A that is in the direction of a unit vector, î: s = dot(A, ihat) In this example, s is the scalar projection of A onto î. The vector projection is the vector that has magnitude s in the direction of î: V = dot(A, ihat) * ihat The cross product of vectors A and B is a vector whose direction is perpendicular to A and B and whose magnitude is
where (again) a is the magnitude of A, b is the magnitude of B, and θ is the angle between the vectors. MATLAB provides a function, cross, that computes cross products. C = cross(A, B) cross only works for 3dimensional vectors; the result is a 3dimensional vector. A common use of cross is to compute torques. If you represent a moment arm R and a force F as 3dimensional vectors, then the torque is just Tau = cross(R, F) If the components of R are in meters and the components of F are in Newtons, then the torques in Tau are in Newtonmeters. 12.3 Celestial mechanicsModeling celestial mechanics is a good opportunity to compute with spatial vectors. Imagine a star with mass m_{1} at a point in space described by the vector P_{1}, and a planet with mass m_{2} at point P_{2}. The magnitude of the gravitational force^{3} between them is
where r is the distance between them and G is the universal gravitational constant, which is about 6.67 × 10^{−11} N m^{2} / kg^{2}. Remember that this is the appropriate value of G only if the masses are in kilograms, distances in meters, and forces in Newtons. The direction of the force on the star at P_{1} is in the direction toward P_{2}. We can compute relative direction by subtracting vectors; if we compute R = P2  P1, then the direction of R is from P1 to P2. The distance between the planet and star is the length of R: r = norm(R) The direction of the force on the star is R: rhat = R / r Exercise 1
Write a sequence of MATLAB statements that computes F12, a vector
that represents the force on the star due to the planet, and F21, the force on the planet due to the star.
Exercise 2
Encapsulate these statements in a function named gravity_force_func that takes P1, m1, P2, and
m2 as input variables and returns F12.
Exercise 3
Write a simulation of the orbit of Jupiter around the Sun. The mass
of the Sun is about 2.0 × 10^{30} kg. You can get the mass of
Jupiter, its distance from the Sun and orbital velocity from
http://en.wikipedia.org/wiki/Jupiter. Confirm that it takes
about 4332 days for Jupiter to orbit the Sun.
12.4 AnimationAnimation is a useful tool for checking the results of a physical model. If something is wrong, animation can make it obvious. There are two ways to do animation in MATLAB. One is to use getframe to capture a series of images and movie to play them back. The more informal way is to draw a series of plots. Here is an example I wrote for Exercise 12.3: function animate_func(T,M) X1 = M(:,1); Y1 = M(:,2); X2 = M(:,3); Y2 = M(:,4); minmax = [min([X1;X2]), max([X1;X2]), min([Y1;Y2]), max([Y1;Y2])]; for i=1:length(T) clf; axis(minmax); hold on; draw_func(X1(i), Y1(i), X2(i), Y2(i)); drawnow; end end The input variables are the output from ode45, a vector T and a matrix M. The columns of M are the positions and velocities of the Sun and Jupiter, so X1 and Y1 get the coordinates of the Sun; X2 and Y2 get the coordinates of Jupiter. minmax is a vector of four elements which is used inside the loop to set the axes of the figure. This is necessary because otherwise MATLAB scales the figure each time through the loop, so the axes keep changing, which makes the animation hard to watch. Each time through the loop, animate_func uses clf to clear the figure and axis to reset the axes. hold on makes it possible to put more than one plot onto the same axes (otherwise MATLAB clears the figure each time you call plot). Each time through the loop, we have to call drawnow so that MATLAB actually displays each plot. Otherwise it waits until you finish drawing all the figures and then updates the display. draw_func is the function that actually makes the plot: function draw_func(x1, y1, x2, y2) plot(x1, y1, ’r.’, ’MarkerSize’, 50); plot(x2, y2, ’b.’, ’MarkerSize’, 20); end The input variables are the position of the Sun and Jupiter. draw_func uses plot to draw the Sun as a large red marker and Jupiter as a smaller blue one. Exercise 4
To make sure you understand how animate_func works,
try commenting out some of the lines to see what happens.
One limitation of this kind of animation is that the speed of the animation depends on how fast your computer can generate the plots. Since the results from ode45 are usually not equally spaced in time, your animation might slow down where ode45 takes small time steps and speed up where the time step is larger. There are two ways to fix this problem:
Exercise 5
Use animate_func and draw_func to vizualize your
simulation of Jupiter. Modify it so it shows one day of simulated
time in 0.001 seconds of real time—one revolution should take
about 4.3 seconds.
12.5 Conservation of EnergyA useful way to check the accuracy of an ODE solver is to see whether it conserves energy. For planetary motion, it turns out that ode45 does not. The kinetic energy of a moving body is m v^{2} / 2; the kinetic energy of a solar system is the total kinetic energy of the planets and sun. The potential energy of a sun with mass m_{1} and a planet with mass m_{2} and a distance r between them is
Exercise 6
Write a function called energy_func that takes the output of
your Jupiter simulation, T and M, and computes the total
energy (kinetic and potential) of the system for each estimated
position and velocity. Plot the result as a function of time and
confirm that it decreases over the course of the simulation. Your
function should also compute the relative change in energy, the
difference between the energy at the beginning and end, as a
percentage of the starting energy.
You can reduce the rate of energy loss by decreasing ode45’s tolerance option using odeset (see Section 11.1): options = odeset(’RelTol’, 1e5); [T, M] = ode45(@rate_func, [0:step:end_time], W, options); The name of the option is RelTol for “relative tolerance.” The default value is 1e3 or 0.001. Smaller values make ode45 less “tolerant,” so it does more work to make the errors smaller. Exercise 7
Run ode45 with a range of values for RelTol and confirm
that as the tolerance gets smaller, the rate of energy loss
decreases.
Exercise 8
Run your simulation with one of the other ODE solvers MATLAB provides
and see if any of them conserve energy.
12.6 What is a model for?In Section 7.2 I defined a “model” as a simplified description of a physical system, and said that a good model lends itself to analysis and simulation, and makes predictions that are good enough for the intended purpose. Since then, we have seen a number of examples; now we can say more about what models are for. The goals of a model tend to fall into three categories.
The exercises at the end of this chapter include one model of each type. 12.7 Glossary
12.8 ExercisesExercise 9
If you put two identical bowls of water into a freezer, one at
room temperature and one boiling, which one freezes first? Hint: you might want to do some research on the Mpemba effect. Exercise 10
You have been asked to design a new skateboard ramp; unlike a typical
skateboard ramp, this one is free to pivot about a support point.
Skateboarders approach the ramp on a flat surface and then coast up
the ramp; they are not allowed to put their feet down while on the
ramp. If they go fast enough, the ramp will rotate and they will
gracefully ride down the rotating ramp. Technical and artistic
display will be assessed by the usual panel of talented judges. Your job is to design a ramp that will allow a rider to accomplish this feat, and to create a physical model of the system, a simulation that computes the behavior of a rider on the ramp, and an animation of the result. Exercise 11
A binary star system contains two stars orbiting each other and sometimes planets that orbit one or both stars^{4}. In a binary system, some orbits are “stable” in the sense that a planet can stay in orbit without crashing into one of the stars or flying off into space. Simulation is a useful tool for investigating the nature of these orbits, as in Holman, M.J. and P.A. Wiegert, 1999, “LongTerm Stability of Planets in Binary Systems,” Astronomical Journal 117, available from http://citeseer.ist.psu.edu/358720.html. Read this paper and then modify your planetary simulation to replicate or extend the results.

Are you using one of our books in a class?We'd like to know about it. Please consider filling out this short survey.
