C++
Exercise 7
TASK DESCRIPTION
You will compute the two-dimensional orbits of multiple bodies under the force of gravity, using the Euler method of numerical integration. From the estimated trajectories, you will determine when the bodies are aligned.
INPUT FORMAT
Your program should read the descriptions of the bodies to simulate from standard input in the following format:
n x1 y1 dx1 dy1 m1 x2 y2 dx2 dy2 m2 ... xn yn dxn dyn mn
The integer n dictates how many objects will follow. For each of the n objects, five real numbers denote, respectively, the starting x position, y position, x velocity, y velocity, and mass of the object.
GRAVITY
You will apply the force of gravity to these objects over time, using the following (simplified) equation of gravitational attraction, where d is the displacement vector between the objects:
Fij = mi·mj ⋅ d ||d||2 ||d||
That is, the force vector drawing objects i and j towards each other is proportional to the product of their masses and inversely proportional to their squared distance from each other. The vector sum of all the forces exerted on an object at a given instant, divided by the object's mass, yields its instantaneous acceleration.
In two dimensions, for bodies A and B with positions (xA, yA) and (xB, yB), and with masses mA and mB b, we have:
| d | = | (xA - xB, yA - yB) | |
| FB | = | mA⋅mB⋅ | d |
| ||d||3 | |||
| FA | = | −FB | |
THE EULER METHOD
Given a differential equation y' = f(y), the Euler method approximates the function locally by a line, and takes a step along that line. For a fixed time step size of h, the Euler method shows how to approximate the next value of the function given the current one:
yn+1 = yn + h ⋅ f(yn)
In this exercise, the state y consists of the position and velocity of all the bodies; y0 is the starting state your program reads from input. The gravitational force equation determines the acceleration of the bodies -- the derivative of the velocities over time. To calculate the acceleration at each time step, you need to evaluate the gravitational attraction between all pairs of objects, adding up the total acceleration per object as you go. Keep in mind that the equation for gravity given above is directional -- it will produce forces of equal magnitude on bodies A and B, pulling A towards B and B towards A.
Once you have computed the accelerations for the current time step, you can update the positions and velocities. If we call an object's x and y accelerations ax and ay, the Euler update is given by:
xt+h = xt + h ⋅ dxt yt+h = yt + h ⋅ dyt dxt+h = dxt + h ⋅ ax dyt+h = dyt + h ⋅ ay
ALIGNMENT
For any three bodies A, B, and C, you can determine the angular deviation of C from the line between A and B using the dot product. If the position vector for body X is denoted pX, we have:
cos θ = (pC - pA)·(pB - pA) ||pC - pA||⋅||pB - pA||
In the case of two dimensional positions, the formula becomes:
cos θ =
(xC - xA)(xB - xA) + (yC - yA)(yB - yA) [(xC - xA)2 + (yC - yA)2]½ ⋅ [(xB - xA)2 + (yB - yA)2]½
Note that many of the terms above are used more than once, so can be stored in variables to make everything clearer.
To test that three objects are aligned within ε radians, you need only test that (1 - |cos &theta|) < (1 - cos ε).
In C++, abs(x) yields the absolute value of the real number x. For the purposes of this exercise, you need only test alignment of the first three bodies in the order specified.
TEST CASE
You should write a program that uses the Euler method to compute the trajectories of multiple objects over discrete time steps. To evaluate your program you should run it with the following three body input:
3 0 0 0 0 1 0.9 0 0 -1.0 0.1 0 1.2 1.0 0 0.01
Use a time step of h = 0.01, and run from t=0 to t=10. At each time step, you should test whether the three bodies are aligned to within 0.01 radians (cos(0.01) ≈ 0.99995), and output the value of 't' when they are. Record the outputted time values. Keep in mind that if you update the bodies at time t and then check for alignment, then the time of alignment is (t+h).
SUGGESTIONS
Your program should roughly follow these steps:
- Read the body descriptions and store them.
- Set the total acceleration on each object to zero.
- For each pair of bodies, compute the gravitational force between them according to their positions and the equation of gravity given above. Add to the total accelerations for both bodies.
- Update each body's position and velocity according to its current velocity and total acceleration due to gravity, and the time step h.
- Check whether the bodies are aligned, and output t if they are.
- Increase t by h, and repeat from step 2 until t reaches the desired limit.
Hence your program should look something like:
int main() {
// read body descriptions and store them
// for each time step {
// set the acceleration for each body to zero
// For each pair of bodies {
// compute the gravitational force between them according to their positions and the equation of gravity.
// Add to the running total accelearation for each body.
// }
// For each body {
// update its position acording to its velocity and the time step h
// update its velocity according to its acceleration and the time step h
// }
// Check whether the bodies are aligned
// if they are, output t
// increase t by h and repeat for next time step
// }
}
For storing the state of the simulated bodies, you should use a class. This allows you to group multiple variables together into a compound data type. You'll want something like this:
// Declare a data type
class Body
{
public:
double x, y; // position
double dx, dy; // velocity
double mass; // mass
double ax, ay; // computed acceleration for Euler
};
You can access the class members using the '.' operator:
Body a, b; a.y = 0; b.mass = 47.0;
You can keep all of the bodies in a vector (remember to #include <vector>):
int n;
cin >> n;
// Declare a vector with n elements
vector<Body> bodies(n);
for (int i=0; i<n; ++i)
{
cin >> bodies[i].x
>> bodies[i].y
>> bodies[i].dx
>> bodies[i].dy
>> bodies[i].mass;
}
You can save the input text to a file (bodies.in, for instance) and pipe it into your program. Here's what the contents of the file might look like:
3 0 0 0 0 1 0.9 0 0 -1.0 0.1 0 1.2 1.0 0 0.01
Here's how to tell your program to read input from the file:
user@tw:~> ./ex4 < bodies.in
In order to compute gravity between all pairs of bodies, you'll need a nested loop over the bodies:
for (int i=0; i<n; ++i)
{
bodies[i].ax = bodies[i].ay = 0; // start with zero acceleration for bodies[i]
for (int j=0; j<i; ++j)
{
// compute gravity between bodies[i] and bodies[j],
// and add to bodies[i].ax, bodies[i].ay, bodies[j].ax, and bodies[j].ay
}
}
Note how the bounds on i and j ensure that each pair is considered exactly once. Don't forget to divide by an object's mass to convert force to acceleration.
You might want to use the standard sqrt (square root) and acos (arc cosine) functions for computing gravity and testing for alignment of bodies. They are both available with #include <cmath>. The acos function returns a value between 0 and π. Note that you don't really need to use acos; you can compare the dot product to cos θ directly.
VISUALISATION
To see the motion of the bodies, you can output the positions and masses of the bodies at each time step and feed them to a program called show_sim, available on the teaching system.
In order to use show_sim, your program should output, at each time step, the number of bodies followed by their positions and masses. For instance, the output corresponding to t=0 for the test case given above would be:
3 0 0 1 0.9 0 0.1 0 1.2 0.01
The following code will print out this information for all of the bodies, assuming the same variable names and class as the example code shown above:
cout << n << endl;
for (int i=0; i<n; ++i) {
cout << bodies[i].x << " " << bodies[i].y << " " << bodies[i].mass << endl;
}
At the command line, you can pipe your program's output directly into the visualiser:
user@tw:~> ./ex4 < bodies.in | /export/teach/1BComputing/show_sim
Alternatively, you can save your output to a file and then run the visualiser:
user@tw:~> ./ex4 < bodies.in > my_output user@tw:~> /export/teach/1BComputing/show_sim < my_output
Be careful not to output anything else to standard output when using the visualiser. If you want information to be printed to the console instead of piped to the visualiser, use cerr instead of cout.
- Back to the main 1B C++ page
- © 2010 University of Cambridge Department of Engineering
Information provided by Roberto Cipolla and Ethan Eade (Last updated: December 2010) - Privacy
- Accessibility