HW6: N-Body Simulation

1. Introduction
This project is a discrete simulation solution of an N-body gravitational force problem.
Given some number of bodies in empty space, the bodies are allowed to move freely and
interact with each other under gravitational force, although in this simulation you will not
simulate collisions between the bodies. The result yields an interesting animation.
2. Value
This program is worth a maximum of 20 points. See the grading rubric at the end for a
breakdown of the values of different parts of this project.
3. Due Date
This project is due by 11:59PM on Sunday, October 28,
4. Objectives
The purpose of this assignment is to give you more practice with the things you have
already learned.
Keywords: parallel vectors, loops, functions, parameters, return values, plotting, global
variables, pythagorean theorem
2
5. Background
Be sure you understand how to write and use for loops and while loops.
Be sure you understand how to write and use functions, including parameters, return
values, and arguments.
Read all you want about the N-body problem, but you won't need to know any of that
information for this project.
6. Assignment
In your CSE1010 folder in Matlab create a new folder for HW6.
6.1. Script nbody.m
Start with the main script file this time. Add your personal information to the comments.
% N-Body Simulation
% CSE1010 Project 6, Fall 2012
% (your name goes here)
% (the current date goes here)
% TA: (your TA's name goes here)
% Section: (your section number goes here)
% Instructor: Jeffrey A. Meunier
clc % clear command window
clear % clear all variables
clf % clear plot area
axis equal
hold on
Enter these statements into the main script file. These values will determine the
characteristics of the simulation.
numBodies = 4;
global G
G = 100; % not -9.8 in this project!
massFactor = 40;
minMass = 5;
distanceFactor = 1000;
velocityFactor = 3;
Now create the following vectors:
• Masses: A row vector of numBodies random values in the range minMass to
minMass + massFactor.
3
• Xs: A row vector of numBodies random values in the range -distanceFactor/2 to
+distanceFactor/2.
• Ys: A row vector of numBodies random values in the range -distanceFactor/2 to
+distanceFactor/2.
• Dxs: A row vector of numBodies random values in the range -velocityFactor/2 to
+velocityFactor/2.
• Dys: A row vector of numBodies random values in the range -velocityFactor/2 to
+velocityFactor/2.
The Masses vector lists the masses of all the bodies.
The Xs and Ys vectors list the locations on the screen where the bodies are.
The Dxs and Dys vectors list the velocities of the bodies.
There's a lot more to put in the main script, but you should start by writing a few functions
first.
6.2. Function circle
Believe it or not, Matlab does not have a way to draw a circle of an arbitrary size or color.
You will write this function.
Create a new function called circle. This function has 4 parameters:
1. x: this is the x coordinate of the center of the circle
2. y: this is the y coordinate of the center of the circle
3. radius: this is the radius of the circle
4. color: this is the color of the circle
One way to draw a circle in Matlab is to plot a bunch of points in a circle.
1. First you must create a vector of numbers in the range 0 to 2pi in 0.01 increments.
Call this vector Theta. Use the colon notation.
2. Generate a vector of X values for the circle by calling the cos function on the Theta
vector and multiplying the result by radius, and adding x to the vector.
3. Generate a vector of Y values for the circle by calling the sin function on the Theta
vector and multiplying the result by radius, and adding y to the vector.
4. Plot all the points at once using the specified color: plot(X,Y,color)
Write comments in the function. Mine look like this. You can copy & paste for this one,
but for the rest of the functions you must write your own comments.
% Draws a circle at x,y having given radius and color.
% Use: circle(x, y, radius, color)
4
Make sure the help feature works with this function:
help circle
Draws a circle at x,y having given radius and color.
Use: circle(x, y, radius, color)
Test the function
Back in the main script file, you will test the circle function by plotting the first two bodies
in the vectors that you generated randomly. Add these lines at the bottom and then run the
program:
circle(Xs(1), Ys(1), Masses(1), 'k');
That statement draws a circle centered at Xs(1), Ys(1) with mass (or radius) Masses(1)
using black. Draw the next one:
circle(Xs(2), Ys(2), Masses(2), 'k');
Those two statements plot the first two random bodies, showing two circles of different
sizes in different locations. Run the program several times to be sure it works correctly.
Here's what mine looked like one time:
6.3. Function plotBodies
Instead of calling the circle function once for each body, it would be nice to write a
function that calls the circle function for all the bodies.
Create a new function called plotBodies that has these parameters:
1. Xs: this is the vector of all the bodies' x coordinates
5
2. Ys: this is the vector of all the bodies' y coordinates
3. Masses: the vector of masses
4. color: the color to plot all the bodies
Determine the number of bodies that are to be plotted (it's the length of either the Xs, Ys,
or Masses vectors) and store it in a variable.
Start a for loop that goes from 1 to the number of bodies. Use a loop control variable
called something like b or body, because you're iterating over the number of bodies. I
used b in my program.
Inside the loop, call the circle function. You need to supply 4 arguments when you call
that function: the x coordinate, the y coordinate, the radius, and the color.
Recall how you tested the circle function:
circle(Xs(1), Ys(1), Masses(1), 'k');
This function call draws the circle for mass 1. Now you want to drawn the circle for mass
1, you want to draw a circle for mass b.
End the for loop.
Write comments in the function.
Test the function
In the main script, replace the two calls to the circle function with a single call to the
plotBodies function using the correct arguments. Run the program. This is what my plot
looks like. Yours will be slightly different, but there should be 4 circles of different sizes
placed randomly in the range +/-500 on the x and y axes:
6
6.5. Function moveBodies
Now that you can draw the n bodies on the screen, you need to be able to move them.
The velocities are specified in the Dxs and Dys vectors. In order to move the bodies, the
position of each body will be changed by its dx & dy values found in the vectors, then the
bodies will be re-drawn, then moved, then re-drawn, and so on.
Create a new function called moveBodies. It has 4 parameters:
1. Xs: these are the x positions of all the bodies
2. Ys: these are the y positions of all the bodies
3. Dxs: these are the x velocities of all the bodies
4. Dys: these are the y velocities of all the bodies
This function has two separate return variables, Xs and Ys. This function is responsible for
changing each element in Xs by the corresponding amount in Dxs, and the amount in Ys
by the corresponding amount in Dys. This function returns the two vectors Xs and Ys after
the vectors have been modified.
Recall that it is allowed to use the same variable name (or names) for a parameter and a
return variable.
Notice that the Dxs vector contains the amount that the Xs vector needs to change. For
example, if Dxs(1) is 2.3, then Xs(1) needs to have 2.3 added to it. If Dxs(1) is -0.02, then
Xs(1) needs to have -0.02 added to it. The same is true for location 2, location 3, and the
rest of the vectors. The same is also true for the Ys & Dys vectors. This is actually much
simpler than it sounds, and it does not require the use of a loop.
Notice that you are changing the values of the Xs and Ys vectors. Recall that when you
7
change Xs and Ys inside this function, you're changing copies of those vectors. Thus, these
two vectors must be returned from this function.
Write comments in the function.
Test the function
In the main script, after you call plotBodies, call the moveBodies function and store the
result into Xs and Ys, like this:
[Xs Ys] = moveBodies(Xs, Ys, Dxs, Dys);
and then call plotBodies again right after that to show the result of moving all the bodies.
If you run this, it's difficult to see that the bodies have moved at all, so put those two move
& plot function calls inside a for loop that goes from 1 to 100. You'll end up writing
something like this in the main script:
for i = 1:100
moveBodies
plotBodies
end
When I run my program this is what I see:
Each body turns into a black streak. This can be fixed by doing this:
• Just before each body is moved, plot the bodies using a color other than black (this
draws over the black circles with green ones).
• Move the bodies.
• Plot the bodies using black.
8
for i = 1:100
plotBodies using green
moveBodies
plotBodies using black
end
Now it actually looks kind of nice.
6.6. Function calculateForce
This function calculates the gravitational force on a body by another body.
This function needs the following parameters:
1. The vector index number of a body. Let's call it body1.
2. The vector index number of another body. Let's call it body2.
3. Xs: these are the x positions of all the bodies.
4. Ys: these are the y positions of all the bodies.
5. Masses
This function returns two values: fx and fy, which are the x and y components of the
resultant force.
The general formula for gravitational force is this:
9
where G is the universal gravitational constant, M1 and M2 are the masses of the two
bodies, and r is the distance between the two bodies.
In this program, G is a global variable that you defined in the main script file. I have found
a value for G that, while not realistic, does create a nice simulation. This function can have
access to that variable by using this statement:
global G
But do not re-assign a value to this variable! You mean to use the variable G that already
exists, not create a new one.
M
1 and M2 can be determined by reading the values out of the Masses vector at locations
body1 and body2.
r is the distance between the two bodies. You know the locations of the two bodies: Xs
(body1), Ys(body1) and Xs(body2), Ys(body2). Use the Pythagorean theorem to determine
the distance between the bodies. It would be useful to create two other variables first: dx
and dy, where dx is the difference between the x values and dy is the difference between
the y values. You will need to use these variables again below.
Since r will be used as the denominator in the formula, you should prevent r from
becoming too small or else the resultant force f will become artificially inflated. In reality,
two bodies will never become closer than the sum of their radii (because that's when they
crash into each other), but this program does nothing to prevent two bodies from
overlapping and even occupying the same point. If this were to happen the distance
between them approaches 0 and the force approaches infinity.
This is how you will solve the problem:
Consider the radius of a body to be the same as the body's mass. For example, a mass of
50 has a radius of 50.
Create an if statement that checks if r (in this case it's the distance between the centers of
the two masses) is less than the sum of the two masses' radii. If it is, set r equal to the sum
of the radii.
Now calculate the force f using G, the mass of body 1, the mass of body 2, and the
distance squared.
Now that you have the resultant force, you need to decompose it into its x and y
10
components.
Calculate the angle of the force by calling atan2(dy, dx) where dy and dx are the variables
you created previously.
Assign to fx the value of calling the cos function on the angle.
Assign to fy the value of calling the sin function on the angle.
Write comments in the function.
Test the function
Put this function call in your main script file and then run the program.
calculateForce(1, 2, [0 400], [0 400], [100 100])
My program displayed this answer:
ans =
2.2097
After you verify your function, remove the function call from the main script.
6.7. Function accelerateBody
This function will apply the gravitational attraction of all other bodies to a single body. This
will change the velocity of the body. This is what accelerate means in this context: not
necessarily increase the velocity of a body, but change its velocity and direction.
This function has these parameters:
1. A body number. I called mine b.
2. Xs: these are the x positions of all the bodies
3. Ys: these are the y positions of all the bodies
4. Dxs: these are the x velocities of all the bodies
5. Dys: these are the y velocities of all the bodies
6. Masses
This function returns two values: dx and dy, which is the new velocity of the body.
To calculate the total force on a body, you must add up all the separate forces on a body
by all bodies other than itself.
11
Start a for loop that does this:
For each body number from 1 to the number of bodies:
• If the body number is not equal to b:
⁃ Calculate the force on the body by calling calculateForce.
⁃ Calculate the x & y accelerations. The acceleration is the force divided by
the mass. The mass is the mass of body b.
⁃ Add the x & y accelerations to the x & y velocities of body b.
You must check to see if the body number is equal to b so that you don't calculate the
force between body b and itself.
Make sure that this function returns the x & y velocities of body b.
Write comments in the function.
Test the function
In the main script, call the accelerateBody function on body number 1 and store the return
values back into the Dxs and Dys vectors:
[Dxs(1) Dys(1)] = accelerateBody(1, ___, ___, ___, ___, ___);
Put that statement inside the loop somewhere. I made it the last statement in the loop.
Now run the program to observe the effect. Body 1 should no longer move in a straight
line because now it is under the gravitational effects of all the other bodies.
You can tell which body is body 1 because it's the one with the curved trace. Here you
can see body 1 accelerating (actually moving faster) and being pulled into an orbit around
another more massive body:
12
I ran the program several times and most plots were usually not this dramatic. In many of
them it was difficult to determine if the path of body 1 was changing at all.
Wow, here's an even better one. Body 1 is a very small body that is being thrown around
by two very massive bodies:
6.8. Function accelerateBodies
You're almost done! Now that you can accelerate a single body, you can apply the
accelerateBody function to all bodies.
This function will apply the gravitational attraction to all bodies, changing the velocities of
each body. The function has these parameters:
1. Xs
2. Ys
3. Dxs
4. Dys
5. Masses
The function returns two values: Dxs and Dys.
Write a for loop that iterates over all the bodies. It calls the accelerateBody function on
each body, and stores the results back into the Dxs and Dys vectors for that body.
Write comments in the function.
Change the main script to call this function instead of the accelerateBody function. Make
sure you store the results back into the Dxs and Dys vectors:
13
[Dxs Dys] = accelerateBodies(___, ___, ___, ___, ___);
Here's what mine looks like:
6.9. Animate the program
This static display is nice, but wouldn't it be nice to watch these bodies move around in
real time?
Do this in the main script:
• Change the for loop to a while true loop. Do not run the program yet. You will not
see any output.
• Add the statement pause(eps) to the bottom of the loop, inside the loop. The eps
value is a built-in constant that is the smallest number that Matlab can represent
accurately. This causes Matlab to update the display and pause for about 2 x 10-16
second. While your program pauses, Matlab will update the display.
How to stop the program
Since this program runs in an unbounded while loop, you must interrupt the program
manually by typing Control-C (hold the Control key and press c key) in either the plot
window or the command prompt window.
Here's what mine looks like after running for a minute or so.
14
You'll notice that the program gets slower and slower the longer it runs. This has
something to do with the way Matlab stores and draws graphics on the screen.
Furthermore, not only is Matlab getting slower, but it's consuming significantly more
memory the longer it runs.
You don't have to let the program run for very long, but test it several times to be sure it's
working properly.
6.10. Comment all the files
Make sure you place "help" comments inside each function, like this:
function circle(___)
% Draws a circle...
% Use: circle(___)
...
end
7. Report
Create a Microsoft Word or OpenOffice Word file (or some other format that your TA
agrees on -- ask him or her if you are not sure). Save the file with a .doc or .docx format.
At the beginning of the document include this information:
N-Body Simulation
CSE1010 Project 6, Fall 2012
Your name goes hereT
he current date goes here
TA: Your TA's name goes here
Section: Your section number goes here
15
Instructor: Jeffrey A. Meunier
Be sure to replace the parts that are underlined above.
Now create the following sections in your document.
1. Introduction
In this section copy & paste the text from the introduction section of this
assignment. You do not need to copy the plot images, but you can if you want to.
Rewrite it a bit to make it flow better. Keep in mind that the introduction should
explain the project to someone who does not have the actual assignment document
handy: this could be you next semester or next year. (Do not discard any of your
work!)
2. Output
Run the program 4 times. Each time let the program run for a minute or so, then
save an image of the plot window. Paste the 4 plot images into the document.
3. Source code
Copy & paste the contents of your .m file(s) here. You do not need to write
anything.
8. Submission
Submit the following things things on HuskyCT:
1. All the .m files for this project stored in a single zip file.
2. The MS Word document.
If for some reason you are not able to submit your files on HuskyCT, email your TA before
the deadline. Attach your files to the email.
9. Notes
Here are some notes about working on this project:
• I can't think of anything right now.
10. Grading Rubric
Your TA will grade your assignment based on these criteria:
• (2 points) The comment block at the beginning of the script file is correct, and each
function has a brief comment section in it.
• (12 points) The program works correctly and contains all the necessary functions.
• (3 points) The program is formatted neatly. Follow any example in the book, or see
the web site here: https://sites.google.com/site/jeffscourses/cse1010/matlabprogram-formatting
• (3 points) The document contains all the correct information and is formatted
16
neatly.