3
$\begingroup$

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?

$\endgroup$
7
  • 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 '18 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 '18 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 '18 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 '18 at 4:07
  • 1
    $\begingroup$ Runge-Kutte 4 order, shorten the time step didn't change anything. $\endgroup$ – אבנר יעקב Aug 16 '18 at 6:20