n-Body simulation for Physics students
Author
Katja Della Libera
Title
n-Body simulation for Physics students
Description
an overview of use cases of the n-Body simulation function and neat examples
Category
Educational Materials
Keywords
physics, n-body problem, particles, planets
URL
http://www.notebookarchive.org/2019-06-7vbuqpm/
DOI
https://notebookarchive.org/2019-06-7vbuqpm
Date Added
2019-06-17
Date Last Modified
2019-06-17
File Size
4.5 megabytes
Supplements
Rights
Redistribution rights reserved



N-Body Simulations
N-Body Simulations
At different points in NS110U you have come across particles with either mass or charge that are subject to the forces of gravity and electromagnetism. You know how to calculate the force exerted by these charges on each other and for a given moment can calculate the acceleration they would experience.
But what would happen if you could observe them moving around? Since you can calculate the acceleration at every individual point, even after the charges or masses moved, you would have to do a complicated integral calculation.
Fortunately, you can use the function nBodySimulation to do this for you.
Coulomb’s Law
Coulomb’s Law
The following is a drawing from Chabay and Sherwood’s textbook used for the class. I will not solve the actual exercises (P47, page 542), but use the setup from the drawing as an example
I will assume the objects are starting from rest. For the positions, I will set Q2 as 0 in the coordinate system and assume it is a two-dimensional space
In[]:=
coulombData= NBodySimulation["Coulomb", <|"Mass"1,"Position"{-0.4,0},"Velocity"{0,0},"Charge"-4×|>,(*Q1*) <|"Mass"1,"Position"{0,0},"Velocity"{0,0},"Charge"3×|>,(*Q2*) <|"Mass"1,"Position"{0,-0.3},"Velocity"{0,0},"Charge"-2×|>,(*Q3*) 5]
-6
10
C
-6
10
C
-6
10
C
Out[]=
NBodySimulationData
| ||||
Data not saved. Save now |
So what did this actually calculate? Let’s look at all the parameters known about the state of the 3 particles at time 4:
In[]:=
coulombData[All,4]
Out[]=
Mass1,Position{-0.968502,0.491381},Velocity{-0.219485,0.132175},Charge,Momentum{-0.219485,0.132175},Mass1,Position{0.264546,-0.465026},Velocity{-0.082725,-0.00452082},Charge,Momentum{-0.082725,-0.00452082},Mass1,Position{0.303956,-0.326354},Velocity{0.30221,-0.127654},Charge,Momentum{0.30221,-0.127654}
-
1
250000
C
3
1000000
C
-
1
500000
C
Easier to see is the behavior in a plot.
Below is the position of the particles plotted from time 0 to 2s.
Below is the position of the particles plotted from time 0 to 2s.
In[]:=
ParametricPlot[Evaluate[coulombData[All,"Position",t]],{t,0,2}]
Out[]=
Try changing the time interval. What happens to each of the particles?
A fancier version of this could even add an interactive element, where you can observe the position of the celestial bodies at a given moment.
Try moving the slider below (you will have to enable dynamic content).
Try moving the slider below (you will have to enable dynamic content).
In[]:=
Manipulate[ Show[ ParametricPlot[Evaluate[coulombData[All,"Position",t]],{t,0,2}], ListPlot[ coulombData[All,"Position",time], PlotStyleDirective[PointSize[Medium],Red]]], {time,0,2}, SaveDefinitionsTrue]
Out[]=
| |||||||
|
We can also look at the other parameters, or at individual particles. Try understanding the plot of the velocity of particle one shown below.
Why do I need to use the Sqrt function?
How does this correspond to the path of the particle in the plot above? You can try playing around with the slider to find the peek you see below and explain it.
Why do I need to use the Sqrt function?
How does this correspond to the path of the particle in the plot above? You can try playing around with the slider to find the peek you see below and explain it.
In[]:=
Plot[Sqrt[+],{t,0,2},PlotRangeAll]
2
Evaluate[coulombData[1,"Velocity",t]][[1]]
2
Evaluate[coulombData[1,"Velocity",t]][[2]]
Out[]=
Run the line below to get the velocity plot for particle two. What do you expect?
In[]:=
Plot[Sqrt[+],{t,0,2},PlotRangeAll]
2
Evaluate[coulombData[2,"Velocity",t]][[1]]
2
Evaluate[coulombData[2,"Velocity",t]][[2]]
Try plotting at least 1 more parameter for any combination of bodies yourself.
Newtonian Gravity
Newtonian Gravity
Of course, the simulation doesn’t just work for particles, but also celestial bodes or any mass interacting through gravity.
In[]:=
celestialData = NBodySimulation["Newtonian", <|"Mass"Quantity[1000,"MetricTons"], "Position"Quantity[{0,0,0},"Meters"], "Velocity"Quantity[{0,0,0},"Meters"/"Seconds"]|>, <|"Mass"Quantity[100,"MetricTonsTons"], "Position"Quantity[{1,1,0},"Meters"], "Velocity"Quantity[{0,0,0},"Meters"/"Seconds"]|>, <|"Mass"Quantity[200,"MetricTons"], "Position"Quantity[{0,1,1},"Meters"], "Velocity"Quantity[{0,0,0},"Meters"/"Seconds"]|>, <|"Mass"Quantity[300,"MetricTons"], "Position"Quantity[{1,0,1},"Meters"], "Velocity"Quantity[{0,0,0},"Meters"/"Seconds"]|>, <|"Mass"Quantity[400,"MetricTons"], "Position"Quantity[{-1,-1,0},"Meters"], "Velocity"Quantity[{0,0,0},"Meters"/"Seconds"]|>, <|"Mass"Quantity[500,"MetricTons"], "Position"Quantity[{0,-1,-1},"Meters"], "Velocity"Quantity[{0,0,0},"Meters"/"Seconds"]|>, <|"Mass"Quantity[600,"MetricTons"], "Position"Quantity[{-1,0,-1},"Meters"], "Velocity"Quantity[{0,0,0},"Meters"/"Seconds"]|>, 20]
Out[]=
NBodySimulationData
| ||||
Data not saved. Save now |
In[]:=
Manipulate[ Show[ ParametricPlot3D[Evaluate[celestialData[All,"Position",t]],{t,0,20}], ListPointPlot3D[ celestialData[All,"Position",time], PlotStyleDirective[PointSize[Medium],Red]]], {time,0,20}, SaveDefinitionsTrue]
Out[]=
| |||||||
|
Oscillations
Oscillations
Another type of n-body problem we talked about is masses attached to springs that are oscillating. To find the exact position at any moment, you can also use the nBodySimulation function.
In[]:=
oscillationData=NBodySimulation["Harmonic",{<|"Mass"1,"Position"{-1/3},"Velocity"{0}|>,<|"Mass"4,"Position"{1/3},"Velocity"{0}|>},10]
Out[]=
NBodySimulationData
|
Can you tell from the magnitude of oscillations which body is heavier?
In[]:=
Plot[Evaluate[oscillationData[All,"Position",t]],{t,0,10}]
Out[]=
This last example is the first example in the official documentation for the nBodySimulation, which you can check out here: .


Cite this as: Katja Della Libera, "n-Body simulation for Physics students" from the Notebook Archive (2019), https://notebookarchive.org/2019-06-7vbuqpm

Download

