Get 1 week of unlimited access
Class Notes (1,053,408)
CA (602,327)
UTSC (35,382)
Lecture 8

# PSCB57H3 Lecture Notes - Lecture 8: Celestial Mechanics, Astronomical Object, Fiji

10 pages93 viewsSpring 2017

Department
Physical Sciences
Course Code
PSCB57H3
Professor
Hanno Rein
Lecture
8

This preview shows pages 1-3. to view the full 10 pages of the document.
Introduction to Scientiﬁc Computing
Lecture 8
Professor Hanno Rein
Last updated: November 7, 2017
1 N-body integrations
1.1 Newton’s law
We’ll discuss today an important real world application of numerical ODE integration: celestial me-
chanics. Celestial mechanics is the branch of astronomy that deals with the motions of any kind of
celestial object. These motions are governed by Newton’s law of gravitation. In some cases there are
small correction due to general relativity, tides or radiation forces, which we will not include in our
model.
Newtonian gravity says that a bodies of mass mifeels another body of mass mj, separated by the
vector ~r, with the following force
~
Fij =G mimj
|~rij |2·~rij
|~rij |
where Gis the gravitational constant. The last term is simply telling us in which direction the force
is directed (towards the other body). If there are Nbodies in total, then the force is a sum over all
bodies (except the body on which the force is exerted):
~
Fi=
N1
X
j=0
i6=j
G mimj
|~rij |2·~rij
|~rij |
To get an acceleration from a force, we use the well known law
~
F=m~a
so that we end up with the acceleration on the i-particle being
~ai=
N1
X
j=0
i6=j
G mj
|~rij |2·~rij
|~rij |
The acceleration is the second time derivate of the position, so we can also write the last equation as
¨
~ri=
N1
X
j=0
i6=j
G mj
|~rij |2·~rij
|~rij |
We have a total of Nparticles, thus we end up with a set of diﬀerential equation of second order
¨
~ri=F(~r)
1

Unlock to view full version

Only half of the first page are available for preview. Some parts have been intentionally blurred.

one for each i= 0...N 1 where F(~r) means that Fdepends on the rof all particles. All those
equations are coupled. We already discussed how to solve this equation using a trick to convert it so
ﬁrst order.
r
˙r0
=˙r
F(r)
Thus, we can simply use the Euler or midpoint method which we already know. Let’s try that.
2

Unlock to view full version

Only half of the first page are available for preview. Some parts have been intentionally blurred.

1.2 Initial conditions
As a test case, we use data from NASA’s Horizon system to simulate the Solar system. Besides
the major planets, we also include the comet 67P/Churyumov-Gerasimenko which is currently being
visited by the Rosetta spacecraft. The data is given in the ﬁle init.txt and has one line per body
and the format of name, followed by the x, y and z position and the x, y and z velocity. The units are
km for the position, km/s for the velocity and the mass is given in kg.
sun 1.988544 E+30 .4914033 47836458 E+05
3.632715592552171 E+05 1.049148558556447E +04
6.242516505718979E03 1.163260794114081E02
2.475374674040771E04
mercury 3.302 E+23 5.419423 38 5170247 E +07
1.625487416441774 E+07 6.231968220825911E +06
4.381057475129131E+00 4.891120140478765E+01
3.593097527852676E+00
venus 4.8685 E+24 2.256798710024889E +07
1.046561273591759E+08 2.760223024328955E+06
3.431401385797417 E+01 7.710006447409605E +00
1.875067390915588E+00
earth 5.97219 E+24 1.418287679667581E +08
4.126884716814923E+07 1.125285256157373E+04
8.830883522656004 E+00 2.868395352996171E +01
1.085239827735510E05
mars 6.4185E +23 2.40827 9029085545 E+07
2.311290271817935E+08 4.261631465386531E+06
2.317767379354534 E +01 4.5213 1275238314 1 E+00
6.638119814842468E01
jupiter 1.89813 E+27 7.778414202838075 E+08
2.244518699271340E+08 1.647441732604697E+07
3.784015648894145 E+00 1.193534649902395E +01
1.342303598962500E01
saturn 5.68319 E+26 2.818805400917871 E+08
1.321479657136385E+09 1.177469869993168E+07
9.963111372343775 E+00 2.038111913779875E +00
4.325059188729626E01
uranus 8.68103 E +25 2.668275 057715974 E+09
1.368207197729832 E+09 3.965456607834578E +07
3.057369253435637E+00 5.742539618944625E+00
1.841246675277081E02
neptune 1.0241 E+26 3.068964 88 3542690 E +09
3.290508769241064 E+09 2.965492893564129E +06
3.938602703572163E+00 3.739029327757723E+00
1.671116647134993E01
67P 0. 5.829716977332731 E+08 3.371344262751943E
+08 3.008099774720960E+07 1.302925387423112E
+00 1.173853622210837 E+01 7.989778520791224E
01
Code 1: The ﬁle init.txt
The plot below show the position of the planet on March 4th 2004.
3