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 13 Celestial Mechanics13.1 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 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}.^{1} Remember that when using this value of G it is easiest to directly specify masses in kilograms and distances in meters, and forces are calculated 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 r by subtracting vectors; if we compute r = P_{2} − P_{1}, then the direction of r is from P_{1} to P_{2}. 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 F_{12}, a vector
that represents the force on the star due to the planet, and
F_{21}, the force on the planet due to the star. Encapsulate
these statements in a function named gravity_force_func that
takes P1, m1, P2, and m2 as input variables
and returns F12 and F21.
Exercise 2
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
https://en.wikipedia.org/wiki/Jupiter. Confirm that it takes
about 4332 days for Jupiter to orbit the Sun.
13.2 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 2: function animate_func(T,M) % 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 3
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 4
Use animate_func and draw_func to visualize 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.
13.3 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 5
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.9): 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 6
Run ode45 with a range of values for RelTol and confirm
that as the tolerance gets smaller, the rate of energy loss
decreases.
Exercise 7
Run your simulation with one of the other ODE solvers MATLAB provides
and see if any of them conserve energy.

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