- #1

- 3

- 0

## Main Question or Discussion Point

Hey guys!

Working on a side project using machine learning and n-body simulations. I could just use the library Rebound (http://rebound.readthedocs.io/en/latest/) but I decided to first try building my own simulator before taking shortcuts.

Simulation uses the masses of the bodies (from Sun to Pluto) in kgs, gets the positions and velocities for a certain date from NASA's Horizons (positions in AU and velocities in AU/hr)

- I convert G to the appropriate units.

- I use a timestep of 0.001

- I run till time = 2000.

- Integrating using leapfrog method.

- Using numpy's longdouble for better precision.

- Calculating acceleration using following code (brute particle to particle method).

t = 0

t = 2000

Questions.

- If I plot it for a longer time duration, the bodies are moving further and further away from the Sun. Am I doing something wrong? Is it because Python cannot handle the right precision, my timestep is not small enough, leapfrog not good enough?

If you could give me some direction or hints, that would be great!

Working on a side project using machine learning and n-body simulations. I could just use the library Rebound (http://rebound.readthedocs.io/en/latest/) but I decided to first try building my own simulator before taking shortcuts.

Simulation uses the masses of the bodies (from Sun to Pluto) in kgs, gets the positions and velocities for a certain date from NASA's Horizons (positions in AU and velocities in AU/hr)

- I convert G to the appropriate units.

- I use a timestep of 0.001

- I run till time = 2000.

- Integrating using leapfrog method.

- Using numpy's longdouble for better precision.

- Calculating acceleration using following code (brute particle to particle method).

Python:

```
# Distances between the two bodies.
dx = other.x - body.x
dy = other.y - body.y
dz = other.z - body.z
r = math.sqrt(dx * dx + dy * dy + dz * dz)
# Gravitational force between particles.
f = -G / (r *r * r) * other.m
body.ax += f * dx
body.ay += f * dy
body.az += f * dz
```

Code:

```
SUN, m=1.00002250368 x=-0.00379862800765 y=0.00681683726634 z=0.0029805622517
MERCURY, m=1.66054877697e-07 x=-0.0671159633132 y=-0.40333325859 z=-0.209555180056
VENUS, m=2.44832880699e-06 x=0.719400450883 y=0.0713573963636 z=-0.0137387980893
EARTH, m=3.04135120774e-06 x=-0.170126299491 y=0.895975454466 z=0.388429986171
MARS, m=3.22781112204e-07 x=-1.32389966792 y=-0.811358049279 z=-0.336669172693
```

Code:

```
SUN, m=1.00002250368 x=-0.00449593006396 y=0.00665448086074 z=0.00293125410067
MERCURY, m=1.66054877697e-07 x=2.20774164708 y=-1.9084499026 z=-1.24941970561
VENUS, m=2.44832880699e-06 x=1.52493838842 y=2.00993141858 z=0.80753519142
EARTH, m=3.04135120774e-06 x=-1.93442316424 y=1.25283838345 z=0.543155283391
MARS, m=3.22781112204e-07 x=-0.886516277866 y=-1.83209935391 z=-0.81665499933
```

- If I plot it for a longer time duration, the bodies are moving further and further away from the Sun. Am I doing something wrong? Is it because Python cannot handle the right precision, my timestep is not small enough, leapfrog not good enough?

If you could give me some direction or hints, that would be great!