Numerical solutions assessment 4, May 2020
This is the fourth Numerical Solutions assessment. It is worth 40% of the total mark available for this module. The deadline for submission is Sunday 7th June 2020 by 23:59. Please submit a .pdf and a .zip to BlackBoard. The assessment is conducted in pairs.
Simulating N-body interactions: the stability of a simplified solar system
An 𝑁-body problem refers to modelling the interaction of 𝑁 astronomical bodies in three dimensions through time. This is modelled by a system of ODEs. The task is to write a general purpose 𝑁-body simulation programme, investigate the behaviour of a set of 𝑁-body systems based on the solar system, and write up the results. There are three scenarios:
i) A system containing only Sun-Earth.
ii) A system containing only Sun-Venus-Earth-Mars-Jupiter.
iii) A system containing only Sun-Venus-Earth-Mars perturbed by a passing neutron star.
The evolution equations and initial values are given in the appendix. You need to be able to output a representation of the positions of the bodies in the system at intervals you can specify.
You are free to choose the solution method. However for purposes of computational efficiency it is recommended that an appropriate higher order method be used. Whatever the method you choose you will need to take a reasonably large number of time steps per year.
For the scenarios you must:
Scenario (i):
Scenario (ii): Scenario (iii):
Comment on the behaviour and stability of the system over many years (see notes) in this simplest case.
Comment on the behaviour and stability of the system over a large number of years. Draw conclusions about how destructive the passage of the neutron star is.
You should assess the effects of the neutron star starting from various initial positions and velocities (speed and direction). You need evolve the system only until the effects on it of the passage of the neutron star are clear.
In each case comment on the accuracy of your numerical solution.
The submission
You must submit two files: (a) a report in .pdf format and (b) a .zip file containing code. The report, as a single pdf file, must contain:
i) A description of the behaviour and stability of the system in each scenario.
ii) An assessment of the accuracy of your numerical solution.
iii) A code listing.
For parts (i) and (ii) you may wish to include graphical representations. Part (ii) should include results of checks you have made on the accuracy of your simulations, perhaps by using different numbers of time steps per year, and referring to plots or tables. Between them, parts (i) and (ii) should be less than five sides excluding plots and tables.
You do not need to ‘re-invent the wheel.’ If you want you may base your code on one of my applications, re-using my code freely as appropriate. Should you do so you must state which version of my code you base yours on. You must list (a) code that you have ‘touched’ eg the new objects you have written, and objects of mine that you have altered and (b) sufficient other application code so that I can understand the programme flow through your application. If in doubt as to what you should or should not include in your listing, please contact me.
The zip file must contain your C++ code according to guidelines issued elsewhere. Your application code must be object oriented with clear client interface. The code you submit must be the complete application, able to run without addition or modification with DevCpp on my standard Windows laptop. the .zip must contain only the .cpp, .h, and .dev files (plus, if used, any I/O files).
You will be assessed on
a) Your assessment and discussion of the behaviour of the system in the three scenarios.
b) Your assessment and discussion of the accuracy of your numerical solutions.
c) The quality of your code in terms of functionality, readability, maintainability, clarity,
efficiency, sophistication, generalizability, and general presentation.
d) The presentational quality of the report as a whole.
NS_ass_4_May_2020.docx 1
Notes
The period of evolution
‘Many’ or a ‘large number’ of years in scenarios (i) and (ii) means a sufficient length of time to establish the long term behaviour of the systems. This may be hundreds of thousands, or millions, of years. The period you choose will also be determined by the run time of your code.
Code efficiency
It is important that your code be reasonably efficient. Simulating scenarios (i) and (ii) over a reasonable period of time (perhaps hundreds of thousands of years), with a time step small enough for numerical solutions to be reasonably accurate over the entire period, may take an excessive amount of time if your code does not take efficiency into account. For instance you should avoid excessive construction and copying of vectors; these operations are relatively expensive.
With reasonably efficient code you should expect a run time for a single simulation of the order of minutes, not seconds, but significantly under an hour.
Plots
You are strongly advised to include plots representing planetary or solar movement in your report. It may be convenient to import your results into Excel or Matlab to plot them there.
Stability
When investigating the stability of the system one should look at eg the Sun-Earth distance, not the distance of the earth from (0, 0, 0). This is because the whole system may drift. You may wish to compare true orbital periods, and other real life properties, with those obtained in your simulation (although you should to expect these to be identical).
Collaboration between pairs is not permitted
The code you submit must be that of your pair alone. Clear similarities in code, in the write-up, or in results, will be taken as evidence of collaboration (for instance if identical or very similar results are given, if there are strong similarities in the write-up, or if code is clearly shared.)
Code closely modelled on that from other sources is not permitted
Code to perform components of this exercise could be found in various places and might, for instance, be downloadable from the web. Use of code taken from, or altered from, such sources is not permitted. You must devise your code yourself. Specifically excluded from this prohibition is code I have declared that you may use. In particular you may use the library code I have made available, without attribution.
Nick Webber May 2020
NS_ass_4_May_2020.docx 2
Appendix: an N-body system interacting under Newtonian gravity
Let there be 𝑁 (point) bodies, {𝐵 } interacting under gravity. Each body 𝐵 = {𝑚 , 𝑃 , 𝑉 } has a 𝑖𝑖=1,…,𝑁 𝑖 𝑖𝑖𝑖
mass 𝑚 , a position 𝑃 = (𝑥 , 𝑥 , 𝑥 ) ∈ R3, and a velocity 𝑉 = (𝑣 , 𝑣 , 𝑣 ) ∈ R3. Bodies interact 𝑖 𝑖 𝑖1 𝑖2 𝑖3 𝑖 𝑖1 𝑖2 𝑖3
under the influence of Newtonian gravity – any relativistic effects can be ignored. They do not spin, and experience no tidal forces. They behave like point masses.
The equations of motion under Newtonian gravity are
𝑑𝑃 (𝑡) 𝑖
𝑑𝑡
𝑑𝑉(𝑡) 𝑖
𝑑𝑡
𝑘=1,…,𝑁 𝑘≠𝑖
=𝑉(𝑡), 𝑖
𝑖 = 1,…,𝑁,
=𝐴𝑖(𝑡),
where 𝐴𝑖(𝑡) ∈ R3 is acceleration due to gravity of the 𝑖th body,
𝑖 = 1,…,𝑁, 𝐴(𝑡)= ∑ 𝐺 𝑚𝑘 (𝑃(𝑡)−𝑃(𝑡)),
𝑖
|𝑃 (𝑡) − 𝑃 (𝑡)|3 𝑘 𝑖 𝑘𝑖
and (in MKS units) 𝐺 = 6.674 × 10−11 m3kg−1𝑠−2 is the constant of gravitational attraction. Here 𝑃 (𝑡) − 𝑃 (𝑡) ∈ R3 is a 3-vector, and |𝑃 (𝑡) − 𝑃 (𝑡)| is the ordinary Euclidian distance,
𝑘𝑖 𝑘𝑖
|𝑃 (𝑡)−𝑃(𝑡)|=√(𝑥 −𝑥 )2 +(𝑥 −𝑥 )2 +(𝑥 −𝑥 )2.
𝑘 𝑖
𝑘1 𝑖1 𝑘2 𝑖2
𝑘3 𝑖3
Note that this system is of the form
where𝑢 ≡{𝑃}, 𝑢 ≡{𝑉}, and𝑓(𝑢)≡{𝐴}with𝐴 =∑𝑘=1,…,𝑁𝐺
𝑢′ = 𝑢 , 12
𝑢′ = 𝑓(𝑢 ), 21
𝑚𝑘 (𝑃(𝑡)−𝑃(𝑡)) Here 1 𝑖 2 𝑖 1 𝑖 𝑖 𝑘≠𝑖 |𝑃𝑘(𝑡)−𝑃𝑖(𝑡)|3 𝑘 𝑖
𝑢1 and 𝑢2 each have 3𝑁 components.
Table 1 gives initial data (in MKS units). This is based on approximate data for the solar system.
Table 1. Initial data
The sun is given an initial velocity. With an initial zero velocity for the Sun, simulating over 106 years, it drifts around 18 AU away from its initial position. However, giving the Sun an initial velocity of (0, -8.94E- 02, 0), chosen so that the momentum of the Sun-Earth system is roughly zero, it drifts only 0.003 AU in this time. In both cases the Sun-Earth distance is stable between 1.003 and 0.996 AU.
You are strongly advised to use the year as the unit of time, not seconds, and astronomical units, 𝐴𝑈 = 1.495978707 × 1011 m, as the units of length, not metres. One AU is roughly the average distance between the Earth and the Sun. In AU-kg-year units the gravitational constant is 𝐺 = 1.9853 × 10−29 AU3 kg-1 yr-2.
𝑚𝑖 (kg)
𝑃 (m) 𝑖0
𝑉 (m/s) 𝑖0
Sun
1.989 × 1030
(0,0,0)
(0, −8.94 × 10−2, 0)
Venus
4.867 × 1024
(−1.082 × 1011, 0, 6.411 × 109)
(0, −3.502 × 104, 0)
Earth
5.972 × 1024
(1.496 × 1011 , 0,0)
(0, 2.978 × 104, 0)
Mars
6.417 × 1023
(0, 2.279 × 1011, 7.359 × 109)
(−2.401 × 104 , 0, 0)
Jupiter
1.898 × 1027
(0, −7.785 × 1011, −1.77 × 1010)
(1.307 × 104, 0, 0)
Neutron star
4.0 × 1030
(0,0,1.0×1015)
(1.2 × 102, 0, −1.0 × 105)
NS_ass_4_May_2020.docx 3