This is the second C++ assessment for 2017 / 2018. It is worth 40% of the total mark available for this
module. The submission date is in 30 April 2018 by 12:00.
Simulating 𝑵-body interaction
An 𝑁-body problem refers to modelling the interaction of 𝑁 astronomical bodies in three dimensions
through time. The assessment task is to write a general purpose 𝑁-body simulation programme and run it
under three scenarios. Assume a Newtonian universe in which relativistic effects can be ignored.
Two different simulation methods must be implemented and their effectiveness, in terms of run-
times and accuracy, compared. The methods are
1) A naive Euler discretization.
2) A Runge-Kutta discretization.
Details of these discretizations are given in the appendix. (Note: neither of these discretization methods
are accurate over periods even as brief as hundreds of thousands of years.)
You must produce simulation results for the following scenarios
i) A system containing only Sun-Earth (based on the solar system).
ii) A system containing only Sun-Venus-Earth-Mars-Jupiter (based on the solar system).
iii) A system containing only Sun-Venus-Earth-Mars perturbed by a passing neutron star.
Initial data and evolution equations 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.
For scenarios:
i) You must comment on the effectiveness of the discretizations,
and the stability of the system, if possible, over millions of years, in this simplest case.
ii) You must comment on the effectiveness of the discretizations,
and the stability of the system.
iii) You should alter the initial velocity (speed and direction) of the neutron star.
Draw conclusions about how destructive the passage of the neutron star is.
The application must be object oriented with a clear client interface and input validation. The code you
submit must to be able to run without addition or modification with the DevCpp installation on the
machines in the teaching area.
The report
Your report must include a clear description of your application and an assessment of its performance.
You will be assessed on
a) The quality of your code in terms of readability, maintainability, clarity,
sophistication, generalizability, and general presentation.
b) The degree to which your code has the required functionality.
c) Your assessment of the efficacy of your implementation.
d) The presentational quality of the assessment as a whole.
Group work is not permitted
The code you submit must be yours alone. Clear similarities in code, in the write-up or in results, will
be taken as evidence of group work (for instance if identical or very similar screen shots as given, 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
C++_ass_Mar_2018.docx 2
Appendix: an 𝑵-body system interacting under gravity
Let there be 𝑁 (point) bodies, {𝐵𝑖}𝑖=1,…,𝑁 interacting under gravity. Each body 𝐵𝑖 = {𝑚𝑖,𝑃𝑖,𝑉𝑖} has a
mass 𝑚𝑖, a position 𝑃𝑖 = (𝑥𝑖,𝑦𝑖,𝑧𝑖) ∈ ℝ3, and a velocity 𝑉𝑖 = (𝑢𝑖,𝑣𝑖,𝑤𝑖) ∈ ℝ3. Bodies interact under
the influence of Newtonian gravity. 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 𝐴𝑖(𝑡) is acceleration due to gravity,
𝐴𝑖(𝑡) = ∑ 𝐺 𝑚𝑘|𝑃
𝑘(𝑡) −𝑃𝑖(𝑡)|3
(𝑃𝑘(𝑡) −𝑃𝑖(𝑡))
𝑘=1,…,𝑁,
𝑘≠𝑖
,
and (in MKS units) 𝐺 = 6.674 ×10−11 m3kg-1s-2 is the constant of gravitational attraction.
Here 𝑃𝑘 −𝑃𝑖 ∈ ℝ3 is a 3-vector, and |𝑃𝑘 −𝑃𝑖| is the ordinary Euclidian distance,
|𝑃𝑘 −𝑃𝑖| = √(𝑥𝑘 −𝑥𝑖)2 +(𝑦𝑘 −𝑦𝑖)2 +(𝑧𝑘 −𝑧𝑖)2.
Table 1 gives initial data. This is based on approximate data for the solar system.
Table 1. Initial data
𝑚𝑖 (kg) 𝑃𝑖0 (m) 𝑉𝑖0 (m/s)
Sun 1.989×1030 (0,0,0) (0,0,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)
If using MKS units you may need to take Δ𝑡 = 3600×24×365.25100,000 , one hundred thousandth of a year (in
seconds), in the Euler discretization. In the Runge-Kutta discretization a larger time step might be
possible.
It may be convenient to use one year as the unit of time, instead of seconds, in which case Δ𝑡 = 1100,000,
and the astronomical units, 𝐴𝑈 = 1.495978707 ×1011 m, as the unit of length, instead of the metre.
(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.
C++_ass_Mar_2018.docx 3
Appendix: the Euler and Runge-Kutta discretizations
Consider a system of equations of the form.
𝑑𝑦𝑖
𝑑𝑡 = 𝑓𝑖(𝑡,𝑦1,…,𝑦𝑚), 𝑖 = 1,…,𝑚.
Set 𝒚 = (𝑦1,…,𝑦𝑚)′, 𝒇 = (𝑓1,…,𝑓𝑚)′, so in vector notation one can write
𝑑𝒚
𝑑𝑡 = 𝒇(𝑡,𝒚).
In our case 𝑚 = 2𝑁, 𝑦𝑖 ∈ ℝ3, and
𝑦𝑖 = {𝑃𝑖, 𝑖 = 1,…,𝑁,𝑉
𝑖, 𝑖 = 𝑁 +1,…,2𝑁,
𝑓𝑖(𝑡,𝑦1,…,𝑦𝑚) = 𝑓𝑖(𝑃1,…,𝑃𝑁,𝑉1,…,𝑉𝑁) = {
𝑉𝑖, 𝑖 = 1,…,𝑁,
∑ 𝐺 𝑚𝑘|𝑃
𝑘 −𝑃𝑖−𝑁|3
(𝑃𝑘 −𝑃𝑖−𝑁)
𝑘=1,…,𝑁,
𝑘≠𝑖−𝑁
, 𝑖 = 𝑁 +1,…,2𝑁.
Note that the fact that 𝑦𝑖 ∈ ℝ3 causes no complications; one looks at each component of 𝑦𝑖 individually.
Time evolves in ticks of length ℎ = Δ𝑡 starting from time 𝑡0 = 0, so that 𝑡𝑗 = 𝑡𝑗−1 +h = 𝑗ℎ. Write 𝒚𝑗
and 𝒇𝑗 for the values of 𝒚 and 𝒇 at time 𝑡𝑗.
A naive Euler discretization sets
𝒚𝑗+1 = 𝒚𝑗 +ℎ𝒇𝑗.
In our case, writing 𝑃𝑖𝑗 and 𝑉𝑖𝑗 for the position and velocity of body 𝐵𝑖 at time 𝑡𝑗, the Euler discretization
sets
𝑃𝑖,𝑗+1 = 𝑃𝑖𝑗 +ℎ𝑉𝑖𝑗, 𝑖 = 1,…,𝑁,
𝑉𝑖,𝑗+1 = 𝑉𝑖𝑗 +ℎ ∑ 𝐴𝑖𝑘
𝑘=1,…,𝑁,
𝑘≠𝑖
, 𝑖 = 1,…,𝑁,
where
𝐴𝑖𝑘 = 𝐺 𝑚𝑘
|𝑃𝑘𝑗 −𝑃𝑖𝑗|3
(𝑃𝑘𝑗 −𝑃𝑖𝑗).
Note that the time step must be very small for this approximation to be useful.
C++_ass_Mar_2018.docx 4
A Runge-Kutta discretization is more accurate than the Euler discretization. The standard fourth order
Runge-Kutta method sets
𝒅1 = ℎ𝒇(𝑡𝑗,𝒚𝑗),
𝒅2 = ℎ𝒇(𝑡𝑗 + 12ℎ,𝒚𝑗 + 12𝒅1),
𝒅3 = ℎ𝒇(𝑡𝑗 + 12ℎ,𝒚𝑗 + 12𝒅2),
𝒅4 = ℎ𝒇(𝑡𝑗 +ℎ,𝒚𝑗 +𝒅3),
and finally sets
𝒚𝑗+1 = 𝒚𝑗 + 16𝒅1 + 13𝒅2 + 13𝒅3 + 16𝒅4.
(See, for example, Numerical Recipes.) One expects a larger time step can be used with a Runge-Kutta
method to achieve the same accuracy as the Euler method.