Due: Tuesday, March 19 at 11:59pm
Submission: On Courseworks
Skeleton Code: hw4.zip
Please read the instructions carefully. An incorrectly formatted submission may not be counted.
In this homework, we will start from some skeleton code. You are given a module called nbody with python files inside. To submit, place this module inside a folder called uni-hw4 (for for me this would be tkp2108-hw4). Then compress this folder into a .zip file. For most operating systems, you should be able to right click and have a “compress” option in the context menu. This should create a file called tkp2108-hw4.zip (with your uni). Submit this on courseworks.
In this homework, we will start working towards the implementation of our N-Body simulator. In this simulation, we will attempt to model the movement of a number of bodies interacting with and on each other. These simulations are applicable in fields like astrophysics, particle physics, and molecular dynamics, where a number of bodies in motion influence each other with gravitational, electromagnetic, and/or nuclear forces.
The motion of these bodies is dictated by continuous equations - that is, as each body moves, the forces induced by every other body are in flux, altering those forces and influencing the movement. Such equations are extremely difficult if not impossible to solve analytically, so exact solutions for the position of each body over the course of the interaction is impossible to solve as a pure function of time.
In the lectures we discussed simulations and how we can reduce the complexity of our math by freezing time and relying on the values of much simpler instanteneous equations over short time steps. In layman’s terms, rather than attempt to solve the above problem, we:
∆t, assuming the force is constant over that small interval∆t and repeatThus, instead of requiring very complicated analytic solutions, we can rely on the much simpler instantaneous equations (and accumulating some error in the process).
For the next few homeworks, we will be building a simulator to model an arbitrary number of objects in two dimensions. Each object will be randomly assigned a position, mass, and initial velocity. Using the equation for gravitational force we covered in the lectures on simulations, we will then simulate each objects motion over a certain simulation duration. Over the course of the simulation, we will collect statistics about what is going on during the simulation, and we will animate (in the form of a sequence of plots) the motion of our objects. The end result should look like a much more interesting and complicated version of the gif above.
G = 6.6742 x 10^-11 is the gravitational constant
m1 = Mass of Object 1
m2 = Mass of Object 2
r = Distance between them
F = G·m1·m2 / r^2
Before we write the simulation, we will familiarize ourselves with some of the libraries and tools we will use. In this first homework on N-Body simulations, we will get comfortable using numpy arrays and matplotlib plotting functionality, as well as put our knowledge of proper python module format into practice by building of first proper library.
We want to be able to generate randomly positioned objects. In this part, we will write three functions to generate a collection of random mass objects in 2d space, plot them, and be able to calculate the force between any two of them:
def generateObjects(count, center=0, width=10, massmult=100):
'''
Generate a `count` x 3 matrix of (x, y) coordinates corresponding
to `count` number of arbitrarly placed elements
Args:
count (int): the number of (x, y) pairs to generate
center (int): the center of our distribution (average if `method`=normal)
width (int): the width of our distribution (covariance if `method`=normal)
Returns:
a (`count`, 3) array of positions
The return value should look like:
[[x1, y1, mass1],
[x2, y2, mass2],
...
]
'''
def plotObjects(objects):
'''Draw a matplotlib chart of our objects. `objects` argument is an Nx3 matrix of (x,y,mass) coordinates'''
def forceVec(vec1, vec2):
'''Calculate the force vector between vec1 and vec2, where vec1 and vec2 are like [x, y, mass]
The force vector should be of the form [Fx, Fy]
'''
This function is is pretty straightforward. We should generate a list of lists (or 2 dimensional numpy array if you want) where the inner lists are the form: (x coordinate, y coordinate, mass). The only trickiness is making our random numbers conform to the center, width, and massmult arguments provided to the function. center and width refer to where we center our coordinates and how wide they can span, and massmult refers to how we randomly generate the masses.
First, lets think about the behavior of random.random() from python. This function picks a random number unifomly between [0, 1).
If we want to scale this range (e.g. if we want our masses to be from 0 to massmult), we can multiple the random number by a scalar (massmult in our case).
If we want to shift this range (e.g. to center it around center with width width as for our positions), we can scale first (as before) then subtract (width/2 - center).
With the default arguments (center=0, width=10, massmult=100), the output should look like:
[[ 4.298671314866175, -3.5964251775075926, 10.108602806255728 ],
[ -1.2199514255838748, -3.075078260746047, 71.06973748368017 ],
[ 0.43627207943535673, -0.3028974104491997, 34.09480628094399 ],
[ 2.298497565434139, 4.638535489191206, 92.87857278297581 ],
[ -3.7179665494764325, 3.794670456131353, 64.25410023288761 ],
[ -0.24351285066161843, -2.056854085174323, 22.845317256943055 ],
[ -3.67456000843141, 2.0043744330907822, 39.96689798311353 ],
[ -2.1476591013188386, -1.6020214751604103, 31.807598762388658 ],
[ 4.927651653767887, 2.6000126712950706, 37.30082584067303 ],
[ 4.497739723733117, 0.9902447800327749, 9.366348765977428 ]]
The first column is our data points' x coordinates, the second column is the y coordinates, and the third column is the mass.
If we plot our uniformly distributed and random-mass points with the default arguments, our result should look like:
Notice how the center of our data points is around 0, the width in x and y coordinates is 5 (width / 2) and the size corresponds to our mass multiplier. You are also given a function that generates the same type of dataset, but instead of uniformly distributed, it samples from a 2 dimensional normal distribution (e.g. a 2-d bell curve).
def generateObjectsNormal(count, center=0, width=10, massmult=100):
import numpy as np
return np.random.multivariate_normal([center, center, 100], [[width/2, 0, 0], [0, width/2, 0], [0, 0, massmult*10]], count)
If we were to plot the same distribution, rather than randomly distributed it would look like a 2d bell curve.
We will spend more time over the course of the semester getting comfortable with matplotlib. For this assignment, it is simple enough for us to rely on just 3 matplotlib commands.
import matplotlib.pyplot as plt
# Creates a new "figure", or plot
plt.figure(figsize=(14, 8))
# Scatter X_data vs Y_data with size S_data and transparency 0.5
plt.scatter( X_data, Y_data, s=Size_data, alpha=0.5)
# Show the resulting plot
plt.show()
So looking at the function we need to implement:
def plotObjects(objects):
'''Draw a matplotlib chart of our objects. `objects` argument is an Nx2 matrix of (x,y) coordinates'''
All we need to do is use these three commands with the appropriate slices of our objects list-of-lists corresponding to the first column (x), second column (y), and third column (mass).
The last part of our problem is to calculate the force between two rows of our dataset. This is straightforward, but requires us to remember a bit of trigonometry from high school. Each row of our dataset is of the form (x, y, mass), and lets assume we have two of these objects vec1 = (x1, y1, mass1), vec2 = (x2, y2, mass2). Visualizing our two data points would look like this:
To calculate the force, we need to know:
vec1[2])vec2[2])The distance is the only one we don’t know, but we can use the Pythagorean theorem: distance = sqrt((x1 - x2)^2 + (y1 - y2)^2). Of course, let’s write a little helper for this (this is just a suggestion, you may implement however you like):
def distance(s1, s2):
'''s1 and s2 are coordinate pairs (x1, y1) and (x2, y2)'''
With the two masses and the distance, we can now calculate force. To make our lives easier with simpler numbers, let’s assume all our masses are in the millions (e.g. a weight of 5 means 5,000,000 kg). Our formula for force is then simply:
# If our units are in the millions
G = 6.6742*10**-5
def force(m1, m2, s):
'''m1 is mass of first object, m2 is mass of second object, s is distance between'''
return (G * m1 * m2) / s**2
The last step requires our trigonometry…ugh. Since we eventually want to calculate the net force on our objects due to every other object, and all the forces will be in different directions, it makes sense for us to break each force into x and y components and sum up these independently, since they will all align along the x and y axis. If we remember from trigonometery, if our two points are at angle θ:
Thankfully, it is pretty easy to calculate θ:
def angle(vec1, vec2):
'''Given vec1 and vec2 vectors like: (x1, y1) and (x2, y2)'''
return math.atan2(vec2[1] - vec1[1], vec2[0] - vec1[0])
Our final function, and the one we have to implement for the assignment ties these concepts together. Here is what my pseudocode looks like:
def forceVec(vec1, vec2):
'''Calculate the force vector between vec1 and vec2, where vec1 and vec2 are like [x, y, mass]
The force vector should be of the form [Fx, Fy]
'''
# Get the mass of the first vector
# Get the mass of the second vector
# Calculate the force between the two objects
# Calculate ethe angle between the two objects
# Calculate the X component of the force vector
# Calculate the Y component of the force vector
# return X component, Y component
Skeleton code is provided with some function signatures provided, but not implemented. Your job is to implement these functions however you like, but such that the results can be used according to the instructions above. Here is a sample python session from a Jupyter notebook using the library.
Skeleton Code: hw4.zip
from nbody import generateObjects, generateObjectsNormal, plotObjects, forceVec
uniform_data = generateObjects(count=1000)
normal_data = generateObjectsNormal(count=1000)
plotObjects(uniform_data)
plotObjects(normal_data)
forceVec(uniform_data[0], uniform_data[1])
[ 0.00036356, -0.00046474]
# Some nice numbers to check
# data point at 0,0 with mass 1935 and data point at 3,4 with mass 1936
forceVec([0, 0, 1935], [3, 4, 1936])
[6.00062906, 8.00083874]