Midterm Project: TheThree-Body Problem
ENGR 3
Due: 11:59 PM, Sunday, May 3
1 Introduction
You have probably learned, perhaps repeatedly, the canonical equation of Newtonian physics: F = ma. This equa-
tion underlies essentially every dynamical system, from the interaction between atoms to the movement of entire
galaxies. Sometimes, when we’re talking about a particularly simple solution, we can use this simple rule to write
out the trajectory of an object. For example, if we throw a ball through the air and ignore air resistance, we know
that its height as a function of time can be represented as:
y (t) = y0 + v0t+
1
2
gt2.
This is known as a closed-form solution—a solution where we can analytically solve for the position at any time just
by looking at the equations.
Let’s stick with gravity but up the scale: think about the Earth orbiting the Sun. The acceleration of an object
of massm under the influence of gravity from an object of massM is given by:
a =
GM
r2
, (1)
where G is a constant and r is the distance between the two objects’ centers. If we consider just the Earth and the
Sun, we can again write out a closed-form solution for the Earth’s position at any time in the future!
Our solar system is, of course, not just the Sun and the Earth, so let’s also add in the Moon as well. The
Moon is quite small compared to these other two bodies, so we expect that it shouldn’t have much impact on the
trajectories. And this is true, but it turns out the effect it has is, while small, monumental in what it does to our
ability to solve for the solution. That’s because as soon as we add the Moon—our third body—there no longer exists
a closed-form solution for the trajectories of any of the three objects.
Since we can’t write out an analytic expression for the trajectory anymore, we instead turn to numerical
methods, simulating the trajectory one time step at a time to construct the resulting trajectory. We’re going to learn
more in the second half of the course about how to utilize some built-in tools in MATLAB to numerically simulate
such a system, but for the midterm project, you’re going to create a rudimentary solver and use it to simulate a model
of the motion of the Sun, Earth, and Moon. It may sound complicated, but we will tell you step by step what you
need to do.
2 Problem Statement
You can find on Gauchospace a file called “3body.mat.” Loading the file in MATLAB will provide you with five 1× 3
vectors, corresponding to values for three different objects (so the first value in each vector is for the first object,
etc). They are:
• masses: the masses of the objects;
• x: the initial x-coordinates of the objects;
• y: the initial y-coordinates of the objects;
• vx: the initial x-velocities of the objects; and
• vy: the initial y-velocities of the objects.
Using these provided values, you will be responsible for simulating the system. You will do this by writing two
functions and a script, the details of which can be found below. The files you will be expected to submit are:
1
• compute_acceleration.m: a function to calculate the gravitational accelerations on each object;
• euler_step.m: a function to take a single time step in the simulation; and
• threebody.m: a script to load in the initial values and simulate the system, plotting the positions of the three
objects at fixed intervals.
Finally, Section 4 of this homework asks you to answer some questions. Please submit this as a PDF file generated
by Word in addition to your Matlab files. If you don’t have Microsoft Word on your computer, not to worry. It is
easy to create a Word document on Box. Instructions on how to do this are included for download in this midterm
project on GauchoSpace.
3 Program Specifications
Detailed instructions for each component of the midterm project are provided. If you have any questions, don’t
hesitate to ask!
3.1 compute_acceleration.m
This function should take as input the masses, x-coordinates, and y-coordinates of the three bodies and return two
vectors, representing accelerations in the x direction and y direction for each of the 3 masses. The acceleration on an
object is equal to the sum of the contributions from the other two objects, according to the gravity equation. Since
we need the components of the acceleration in each direction, it is helpful to consider the distances in each direction.
Let’s consider the acceleration applied to object 1 with massm1 by object 2, with massm2. If the distance between
object 1 (at coordinates (x1, y1)) and object 2 (at coordinates (x2, y2)) in the x-direction is ∆x = x2 − x1, and the
distance between them in the y-direction is∆y = y2 − y1, then the total distance between them is:
r =
√
∆x2 +∆y2, (2)
and the acceleration in each direction is given by:
ax =
∆x m2
r3
, (3)
ay =
∆y m2
r3
. (4)
The total acceleration on object 1 would be this plus the acceleration induced by object 3:
ax,1 =
(x2 − x1)m2√
(x2 − x1)2 + (y2 − y1)2
3 +
(x3 − x1)m3√
(x3 − x1)2 + (y3 − y1)2
3 . (5)
Now, using the same reasoning, write down the total acceleration on object 2 and on object 3 and then use all of the
accelerations to write your compute_acceleration.m function.
To verify your code, we will provide the results for the first acceleration computation:
Object 1 2 3
ax 0.000100907029478458 -0.0296000000000000 -0.0672108843537415
ay 0 0 0
3.2 euler_step.m
This function takes as input the x and y coordinates, velocities, and accelerations as well as a time step dt, and returns
(outputs) the new x and y coordinates and velocities. The new position is calculated by approximating the velocity
as a constant over the small time interval dt, and the same is done with the acceleration; this method of numerical
2
simulation is called the forward Euler method. (By the way, Euler is a German name that is pronounced like ‘Oiler’.)
The new position is given by:
xnew = x+ dt ∗ vx, (6)
and the new velocity is given by:
vxnew = vx + dt ∗ ax. (7)
Using the same reasoning, derive similar equations to update y and vy . Do this for all 3 objects. After the first update
using the provided value for dt (see below), the new values should be:
Object 1 2 3
x 0 10 10.5000000000000
y 0 0.000273861278752583 0.000344571956871238
vx 5.04535147392290e-06 -0.00148000000000000 -0.00336054421768708
vy 0 0.547722557505166 0.863950323522004
3.3 threebody.m
This script will first load in the initial values. This script should also set the value of dt, which you will initially set
at dt = 0.0005, and should set t = 0. It should use your function euler_step.m to simulate the 3-body problem from
one t to the next, t+ dt. You should continue the simulation until one of two conditions is met:
• t reaches 1000, or
• the distance between objects 2 and 3 exceeds a value of 3 units.
In addition, every 0.5 time units, you should call the provided function “draw_3body” to animate the system. This
function takes as inputs the 3× 1 vectors representing the x- and y-coordinates.
4 Analysis
Once you have gotten your code working for the provided value of dt, try increasing and decreasing the value of
dt. How does changing this affect the simulation? What does this suggest to you about the forward Euler method?
Write a brief response to these questions. Later in the course, we will utilize numerical solvers that avoid some of
the problems that you have (hopefully!) identified with the forward Euler method.