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 vectors
12.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 force1.
These quantities cannot be described with a single number because they contain multiple components. For example, in a 3-dimensional 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 2-dimensional 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, vx and vy, which are the components of velocity in the x and y directions.
>> V = [30, 40] % velocity in m/s
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 Cd, 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 name2 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 Fd 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 Fd 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/s2.
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 products
Multiplying 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 3-dimensional vectors; the result is a 3-dimensional vector.
A common use of cross is to compute torques. If you represent a moment arm R and a force F as 3-dimensional 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 Newton-meters.
12.3 Celestial mechanics
Modeling celestial mechanics is a good opportunity to compute with spatial vectors. Imagine a star with mass m1 at a point in space described by the vector P1, and a planet with mass m2 at point P2. The magnitude of the gravitational force3 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 m2 / kg2. 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 P1 is in the direction toward P2. 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 × 1030 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.
Animation 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) % animate the positions of the planets, assuming that the % columns of M are x1, y1, x2, y2. 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 Energy
A 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 v2 / 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 m1 and a planet with mass m2 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', 1e-5); [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 1e-3 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.
Exercise 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.
A binary star system contains two stars orbiting each other and sometimes planets that orbit one or both stars4. 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, “Long-Term 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.