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 Mechanics

13.1 Celestial mechanics

Modeling celestial mechanics is a good opportunity to compute with spatial vectors. Imagine a star with mass m₁ at a point in space described by the vector P₁, and a planet with mass m₂ at point P₂. The magnitude of the gravitational force between them is

F_g = G

m₁ m₂

r²

(1)

where r is the distance between them and G is the universal gravitational constant, which is about 6.67 × 10⁻¹¹ N m² / kg².¹ 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₁ is in the direction toward P₂. We can compute relative direction r by subtracting vectors; if we compute r = P₂ − P₁, then the direction of r is from P₁ to P₂.

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₁₂, a vector that represents the force on the star due to the planet, and F₂₁, 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³⁰ 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 Animation

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

When you call ode45 you can give it a vector of points in time where it should generate estimates. Here is an example:
```
end_time = 1000;
step = end_time/200;
[T, M] = ode45(@rate_func, [0:step:end_time], W);
```
The second argument is a range vector that goes from 0 to 1000 with a step size determined by step. Since step is end_time/200, there will be about 200 rows in T and M (201 to be precise).
This option does not affect the accuracy of the results; ode45 still uses variable time steps to generate the estimates, but then it interpolates them before returning the results.
You can use pause to play the animation in real time. After drawing each frame and calling drawnow, you can compute the time until the next frame and use pause to wait:
```
dt = T(i+1) - T(i);
pause(dt);
```
A limitation of this method is that it ignores the time required to draw the figure, so it tends to run slow, especially if the figure is complex or the time step is small.

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 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 v² / 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₁ and a planet with mass m₂ and a distance r between them is

U = −G

m₁ m₂

(2)

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', 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 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.

1: See https://en.wikipedia.org/wiki/Gravity.

Are you using one of our books in a class?

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