I wrote a code that follows the Barnes-Hut algorithm for gravitational dynamics. Everything looks pretty good, except that when I plot the total energy $$E_{tot}\equiv E_k+E_u=const$$ I get this result:

enter image description here

and you can see that although it seems constant for most of the time, it has some increasing transient in the beginning. I couldn't find a bug in the code that could be responsible for this, is there may be other explanation that can make sense of this?

  • 2
    $\begingroup$ We cannot see your code, so its impossible to answer. And actual coding questions should be placed on Stack Overflow $\endgroup$
    – user1569
    Aug 15, 2018 at 7:39
  • $\begingroup$ Constants of motion not staying constant is a common problem with numerical integration codes. It needn't come from a very obvious bug, but merely from an inappropriate choice of algorithm or implementation. I'd suggest you search the literature, as it's very unlikely you're the first to run into this. $\endgroup$
    – Mark Olson
    Aug 15, 2018 at 13:38
  • $\begingroup$ The Barnes-Hut has major problems with discontinuity. If a star moves from an area into another, what will happen? $\endgroup$
    – peterh
    Aug 15, 2018 at 21:15
  • $\begingroup$ This is a simulation question, not a programming question; don't migrate to Stack Overflow but consider Computational Science SE. Anyway, what happens if you shorten the timestep? $\endgroup$
    – Mike G
    Aug 16, 2018 at 4:07
  • 1
    $\begingroup$ Runge-Kutte 4 order, shorten the time step didn't change anything. $\endgroup$ Aug 16, 2018 at 6:20


Browse other questions tagged .