# Precision error in Nbody simulation

**URL:** <https://discourse.julialang.org/t/precision-error-in-nbody-simulation/79061>\
**Category:** General Usage\
**Created:** [April 5, 2022, 4:30pm UTC](https://discourse.julialang.org/t/precision-error-in-nbody-simulation/79061 "2022-04-05T16:30:40Z")\
**Posts on this page:** 12\
**Page:** 1

<div class="post-metadata">

**Author:** ![staufenbach](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/staufenbach/32/35221_2.png) [@staufenbach](https://discourse.julialang.org/u/staufenbach)\
**Post date:** [April 5, 2022, 4:30pm UTC](https://discourse.julialang.org/t/precision-error-in-nbody-simulation/79061/1 "2022-04-05T16:30:40Z")

</div>

Hi,

I’m currently writing a Nbody simulation using the Barnes-Hut algorithm. In order to validate what i wrote, i want to compare it to the brute force calculation. (setting the theta parameter in BH makes it degenerate in a brute force calculation)

However,for reason i dont fully understand there a discrepancy between the two result (at machine precision level) which end up accumulating and diverging.

Except that the BH algorithm is written in a recursive form in a separate module, i see no reason why i shouldn’t get the same result.

does anyone has any insight about what is happening ? and any workaround if possible ?

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [April 5, 2022, 4:32pm UTC](https://discourse.julialang.org/t/precision-error-in-nbody-simulation/79061/2 "2022-04-05T16:32:40Z")

</div>

> [@staufenbach](#):
>
> However,for reason i dont fully understand there a discrepancy between the two result (at machine precision level) which end up accumulating and diverging.

> [@PSA: floating-point arithmetic](https://discourse.julialang.org/t/psa-floating-point-arithmetic/8678):
>
> Sometimes people are surprised by the results of floating-point calculations such as julia\> 5/6 0.8333333333333334 # shouldn't the last digit be 3? julia\> 2.6 - 0.7 - 1.9 2.220446049250313e-16 # shouldn't the answer be 0? These are not bugs in Julia. They’re consequences of the IEEE-standard 64-bit binary representation of floating-point numbers that is burned into computer hardware, which Julia and many other languages use by default. Brief explanation You can t…

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [April 5, 2022, 4:33pm UTC](https://discourse.julialang.org/t/precision-error-in-nbody-simulation/79061/3 "2022-04-05T16:33:36Z")

</div>

> [@staufenbach](#):
>
> Except that the BH algorithm is written in a recursive form in a separate module, i see no reason why i shouldn’t get the same result.

In particular, floating-point arithmetic is not associative, so even adding up the same numbers in a different order (e.g. [recursively](https://en.wikipedia.org/wiki/Pairwise_summation)) will generally give a slightly different result.

---

<div class="post-metadata">

**Author:** ![Sukera](https://avatars.discourse-cdn.com/v4/letter/s/ce7236/32.png) [@Sukera](https://discourse.julialang.org/u/Sukera)\
**Post date:** [April 5, 2022, 4:33pm UTC](https://discourse.julialang.org/t/precision-error-in-nbody-simulation/79061/4 "2022-04-05T16:33:39Z")

</div>

> [@staufenbach](#):
>
> Except that the BH algorithm is written in a recursive form in a separate module, i see no reason why i shouldn’t get the same result.

In addition to the links posted by @stevengj, the error is probably due to that - changing the order of operations when dealing with floating point numbers can change the result.

* * *

I see you were faster 😂

---

<div class="post-metadata">

**Author:** ![goerch](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/goerch/32/29122_2.png) [@goerch](https://discourse.julialang.org/u/goerch)\
**Post date:** [April 5, 2022, 4:34pm UTC](https://discourse.julialang.org/t/precision-error-in-nbody-simulation/79061/5 "2022-04-05T16:34:07Z")

</div>

> [@staufenbach](#):
>
> Except that the BH algorithm is written in a recursive form in a separate module, i see no reason why i shouldn’t get the same result.

Same as in `==` or in `isapprox`? And N-body is hard AFAIU.

Edit: the BH algorithm gives an approximate solution? And you are comparing it to an better approximation via brute force?

---

<div class="post-metadata">

**Author:** ![PeterSimon](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petersimon/32/25193_2.png) [@PeterSimon](https://discourse.julialang.org/u/PeterSimon)\
**Post date:** [April 5, 2022, 5:24pm UTC](https://discourse.julialang.org/t/precision-error-in-nbody-simulation/79061/6 "2022-04-05T17:24:59Z")

</div>

The problem may be ill-conditioned. I would try gradually increasing the precision used in each algorithm and see if the answers begin to converge.

---

<div class="post-metadata">

**Author:** ![staufenbach](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/staufenbach/32/35221_2.png) [@staufenbach](https://discourse.julialang.org/u/staufenbach)\
**Post date:** [April 5, 2022, 5:30pm UTC](https://discourse.julialang.org/t/precision-error-in-nbody-simulation/79061/7 "2022-04-05T17:30:33Z")

</div>

oh, i forgot about the impact of the order of operations. its effectively not the same.

@goerch yes BH is approximate but the level of approximation can be control. and you can make it so its exactly equivalent to a brute force approach. also isapprox is true for the ~20 first iterations.

thanks for the answer. i probably round of the acceleration after each iteration. its will be slower but will be easier to see if everythin works fine

---

<div class="post-metadata">

**Author:** ![goerch](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/goerch/32/29122_2.png) [@goerch](https://discourse.julialang.org/u/goerch)\
**Post date:** [April 5, 2022, 5:33pm UTC](https://discourse.julialang.org/t/precision-error-in-nbody-simulation/79061/8 "2022-04-05T17:33:16Z")

</div>

> [@staufenbach](#):
>
> i probably round of the acceleration after each iteration.

You could also parameterize `isapprox`

```julia
help?> isapprox
search: isapprox

  isapprox(x, y; atol::Real=0, rtol::Real=atol>0 ? 0 : √eps, nans::Bool=false[, norm::Function])
[...]

```

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [April 5, 2022, 5:35pm UTC](https://discourse.julialang.org/t/precision-error-in-nbody-simulation/79061/9 "2022-04-05T17:35:21Z")

</div>

> [@staufenbach](#):
>
> which end up accumulating and diverging.

Note that N-body dynamics are typical chaotic, so even tiny differences in roundoff errors will eventually lead to [exponentially differing results](https://en.wikipedia.org/wiki/Lyapunov_exponent).

The usual rejoinder is that any quantity that is exponentially sensitive to initial conditions or other perturbations is not physically meaningful to begin with, and you want to instead look at statistical properties of the solution that are robust to tiny perturbations.

---

<div class="post-metadata">

**Author:** ![goerch](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/goerch/32/29122_2.png) [@goerch](https://discourse.julialang.org/u/goerch)\
**Post date:** [April 5, 2022, 5:40pm UTC](https://discourse.julialang.org/t/precision-error-in-nbody-simulation/79061/10 "2022-04-05T17:40:22Z")

</div>

> [@stevengj](#):
>
> Note that N-body dynamics are typical chaotic.

Not sure about this, isn’t the solar system N-body?

Edit: dangerously off topic…

---

<div class="post-metadata">

**Author:** ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)\
**Post date:** [April 5, 2022, 5:41pm UTC](https://discourse.julialang.org/t/precision-error-in-nbody-simulation/79061/11 "2022-04-05T17:41:43Z")

</div>

> **[Is the Solar System Stable?](https://www.ias.edu/ideas/2011/tremaine-solar-system)**
>
> The stability of the solar system is one of the oldest problems in theoretical physics, dating back to Isaac Newton. After Newton discovered his famous laws of motion and gravity, he used these to determine the motion of a single planet around the...

---

<div class="post-metadata">

**Author:** ![xzackli](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/xzackli/32/38301_2.png) [@xzackli](https://discourse.julialang.org/u/xzackli)\
**Post date:** [April 5, 2022, 6:13pm UTC](https://discourse.julialang.org/t/precision-error-in-nbody-simulation/79061/12 "2022-04-05T18:13:02Z")

</div>

For some arguments on the origin of the instabilities, see i.e.

> **[On the Exponential Instability of N-Body Systems](https://ui.adsabs.harvard.edu/abs/1993ApJ...415..715G/abstract)**
>
> We reconsider the old problem of the growth of numerical errors in N-body integrations. We analyze the effects of successive encounters and show that these tend to magnify errors on a time scale which is comparable with the crossing time. This...