Many problems in statistical physics and astrophysics requiring solving problems consisting of many particles at once (sometimes on the order of thousands or more!). This can't be done by the traditional pen and paper techniques you are all learning in your physics classes. Instead, we must impliment numerical solutions to these problems.
Today, you will create your first of many numerical simulation for a simple problem is that solvable by pen and paper already, the 2 body problem in 2D. In this problem, we will describe the motion between two particles that share a force between them (such as Gravity). We'll design the simulation from an astronomer's mindset with their astronomical units in mind. This simulation will be used to confirm the general motion of the earth around the Sun, and later will be used to predict the motion between two stars within relatively close range.
We will guide you through the physics and math required to create this simulation. The problem here is designed to use the knowledge of scientific python you have been developing this week.
Like any code in python, The first thing we need to do is import the libraries we need. Go ahead and import Numpy and Pyplot below as np and plt respectfully. Don't forget to put matplotlib inline to get everything within the notebook.
In [ ]:
Now we will define the physical constants of our system, which will also establish the unit system we have chosen. We'll use SI units here. Below, I've already created the constants. Make sure you understannd what they are before moving on.
In [3]:
#Physical Constants (SI units)
G=6.67e-11
AU=1.5e11 #meters. Distance between sun and earth.
daysec=24.0*60*60 #seconds in a day
Next, we will need parameters for the simulation. These are known as intial condititons. For a 2 body gravitation problem, we'll need to know the masses of the two objects, the starting posistions of the two objects, and the starting velocities of the two objects.
Below, I've included the intial conditions for the earth (a) and the Sun (b) at the average distance from the sun and the average velocity around the sun. We also need a starting time, and ending time for the simulation, and a "time-step" for the system. Feel free to adjust all of these as you see fit once you have built the system!
a note on dt:
As already stated, numeric simulations are approximations. In our case, we are approximating how time flows. We know it flows continious, but the computer cannot work with this. So instead, we break up our time into equal chunks called "dt". The smaller the chunks, the mroe accurate you will become, but at the cost of computer time.
In [4]:
#####run specfic constants. Change as needed#####
#Masses in kg
Ma=6.0e24 #always set as smaller mass
Mb=2.0e30 #always set as larger mass
#Time settings
t=0.0 #Starting time
dt=.01*daysec #Time set for simulation
tend=300*daysec #Time where simulation ends
#Intial conditions (posistion [m] and velocities [m/s] in x,y,z coorindates)
#For Ma
xa=1.0*AU
ya=0.0
vxa=0.0
vya=30000.0
#For Mb
xb=0.0
yb=0.0
vxb=0.0
vyb=0.0
It will be nice to create a function for the force between Ma and Mb. Below is the physics for the force of Ma on Mb. How the physics works here is not important for the moment. Right now, I want to make sure you can transfer the math shown into a python function. I'll show a picture on the board the physics behind this math for those interested.
and
If we break Fg into the x and y componets we get:
$$Fx=\frac{-GM_aM_b}{r^3}x$$
$$Fy=\frac{-GM_aM_b}{r^3}x$$
So, $Fg$ will only need to be a function of xa, xb, ya, and yb. The velocities of the bodies will not be needed. Create a function that calculates the force between the bodies given the posistions of the bodies. My recommendation here will be feed the inputs as seperate componets and also return the force in terms of componets (say, fx and fy). This will make your code easier to make and easier to read.
In [1]:
#Function to compute the force between the two objects
def FG(xa,xb,ya,yb):
#Computer rx and ry between Ma and Mb
rx=xb-xa
ry=#Write it in
#compute r^3
r3=#Write in r^3 using the equation above. Make use of np.sqrt()
#Compute the force in Newtons. Use the equations above as a Guide!
fx=#Write it in
fy=-#Write it in
return #What do we return?
Now that we have our function, we need to prepare a loop. Before we do, we need to intialize the loop and choose a loop type, for or while. Below is the general outline for how each type of loop can gp.
For loop:
np.zeros or np.linspace for the amount of steps needed to go through the simulation (which is numSteps=(tend-t)/dt the way we have set up the problem). The for loop condition is based off time and should read rough like: for i in range(numSteps)
np.array([]) and use np.append() to tact on new values at each step like so, xaArray=np.append(xaArray,NEWVALUE). The while condition should read, while t<tendMy preference here is While since it keeps my calculations and appending seperate. But, feel free to use which ever feels best for you!
In [6]:
#Run a loop for the simulation. Keep track of Ma and Mb posistions
#Intialize vectors
xaAr=np.array([])
yaAr=np.array([])
xbAr=#Write it in for Particle B
ybAr=#Write it in for Particle B
Now for the actual simulation. This is the hardest part to code in. The general idea behind our loop is that as we step through time, we calculate the force, then calculate the new velocity, then the new posistion for each particle. At the end, we must update our arrays to reflect the new changes and update the time of the system. The time is super important! If we don't (say in a while loop), the simulation would never end and we would never get our result.
Calculate the force with the last known posistions (use your function!)
Calculate the new velocities using the approximation: vb = vb + dt*fg/Mb and va= va - dt*fg/Ma Note the minus sign here, and the need to do this for the x and y directions!
Calculate the new posistions using the approximation: xb = xb + dt*Vb (same for a and for y's. No minus problem here)
Update the arrays to reflect our new values
Update the time using t=t+dt
Now when the loop closes back in, the cycle repeats in a logical way. Go one step at a time when creating this loop and use comments to help guide yourself. Ask for help if it gets tricky!
In [7]:
#Your loop here
#using while loop method with appending. Can also be done with for loops
while #What is our condition for ending?:
#Compute current force on Ma and Mb. Ma recieves the opposite force of Mb
fx,fy=Fg(xa,ya,xb,yb)
#Update the velocities and posistions of the particles
vxa=vxa-fx*dt/Ma
vya=#Write it in for y
vxb=#Write it in
vyb=#Write it in
xa=xa+vxa*dt
ya=#Wite it in
xb=#Write it in
yb=#Write it in
#Save data to lists
xaAr=np.append(xaAr,xa)
yaAr=#How will I save it for yaAr?
xbAr=np.append(xbAr,xb)
ybAr=np.append(ybAr,yb)
#update the time by one time step, dt
#How do I update the time?
Now for the fun part (or not so fun part if your simulation had an issue), plot your results! This is something well covered in previous lectures. Show me a plot of (xa,ya) and (xb,yb). Does it look sort of familiar? Hopfully you get something like the below image (in units of AU).
In [9]:
from IPython.display import Image
Image("Earth-Sun-averageResult.jpg")
Out[9]:
In [10]:
#Your plot here
plt.plot(#Particle A plot
plt.plot(#Partcile B plot
Now that the code is made, the fun can begin! This is where the code stops and the science begins. Start experimenting with different sinarios. Here are some ideas to get you started:
-We just ran for the average distance and average speed around the sun for more then a year but we didn't complete the oribit. Why? Let's try and look up the maximum speed and smallest distance (or the opposite) and run that for the same time? Do we complete the orbit now? Why or why not?
-Let's do a binary star system. Let's do one Sun and one star that is 2 times larger. Set one at the origin and one at 1AU away. Try different starting speeds of one of the stars. When do you get a repeating pattern? How does that pattern change with speed?
In [ ]: