Let's say I have a system of 1000 particles in a 2D box. When they hit a wall, they bounce back. The particles also affect each other with this force:

$$F={1 \over r^2}$$


$$F_x={\Delta x \over r^3}$$ $$F_y={\Delta y \over r^3}$$


$$\Delta x_i=x_i-x_j$$ $$\Delta y_i=y_i-y_j$$

And the total energy of the system is

$$E_k=\Sigma_i v_i^2$$

I expect $E_k$ remains constant during time. Am I making any mistakes?

I have written a C++ program and found that kinetic energy keeps increasing in time:


  • $\begingroup$ Please explain your setup with more detail. What are Delta_x and Delta_y? How do you simulate the interactions? $\endgroup$ – fffred Dec 23 '16 at 11:58
  • $\begingroup$ @fffred, Thanks a lot. Please see the update. This is the pure physics question. However if you are looking for the code, you can find it here. $\endgroup$ – ar2015 Dec 23 '16 at 12:02

A few things:

  • Make sure you have a notion of "timestep", i.e. a physical time between each increment of velocities and positions. A small timestep is necessary to avoid numerical heating. However, no matter how small you make your timestep, the simulation will overheat at some point. You can just delay that point.

  • Do not expect kinetic energy to be conserved. If you start your distribution as Maxwell-Boltzmann, then it should be conserved in average though. Otherwise, the distribution should tend towards Maxwell-Boltzmann and stabilize the kinetic energy.

  • $\begingroup$ The problem has been the step size. Do you mean that with any numerical method particles will finally accelerate in a crazy way? $\endgroup$ – ar2015 Dec 23 '16 at 12:33
  • $\begingroup$ Not any method. You are using the n-body technique. Some techniques conserve the energy well, but might be bad for other things. Some also get numerical cooling $\endgroup$ – fffred Dec 23 '16 at 12:37

You are neglecting the potential energy associated with the forces between the particles.

Consider, e.g., a positively and a negatively charged particle, initially at rest. They accelerate towards each other. Their kinetic energy increases - kinetic energy isn't conserved - but their kinetic plus potential energy is conserved, i.e. constant in time.

Potential energy, $V$, and force are related by $$ \vec F = - \vec\nabla V $$ where I've used the nabla operator. If you consider kinetic and potential energy, you should find that energy is conserved in your simulation.

  • $\begingroup$ But here the kinetic energy is increasing up to infinity. Something must be wrong. The gas inside a box finally reach to an equilibrium and it is impossible that all particles increase in speed up to infinity. $\endgroup$ – ar2015 Dec 23 '16 at 12:20
  • $\begingroup$ I also don't find this answer helpful. While you are right that the potential energy should be taken into account when it comes to energy conservation, still answer still doesn't adress the problem in the question. Let's say that there is a state where the potential energy of the particles is minimized, $U_{min}$, and the starting potential energy is $U_0$. The increase in kinetic enery can't exceed the difference between $U_0$ and $U_{min}$. $\endgroup$ – Quantumwhisp Dec 23 '16 at 12:37
  • $\begingroup$ The question states that e.g., total energy of the system is equal to KE. Yes, there may be further issues with the simulations, but surely pointing out the biggest, most glaring mistake is helpful. $\endgroup$ – innisfree Dec 23 '16 at 12:41

What's this bit

                    // F=1/r^2
                    double r=sqrt(r2);

looks like you have an inverse cube law not an inverse square law as stated in the comment?

Once this is resolved/understood there may still be errors, 1)you are using a finite time step (of one) and differentials to approximate the motion of the particles. Approximate being the operative word.

  • so I think you will find if you crank your time step down it will get better.
  • normally it's necessary to add a compensation factor to prevent the run-away.

2) your boundary conditions are a hard boundary, which is not real, flipping the velocity round introduces 2v worth of acceleration and therefore energy. There is no "equal and opposite re-action". - normally it's necessary to add a compensation factor to prevent the run-away.

Hope this helps.

  • $\begingroup$ Flipping velocities at the boundary doesn't add energy, does it? $\endgroup$ – innisfree Dec 23 '16 at 12:45
  • $\begingroup$ beside that, $\sqrt{f_x^2+f_y^2}=\sqrt{r^2/r^6}=1/r^2$ which is correct. $\endgroup$ – ar2015 Dec 23 '16 at 12:47
  • $\begingroup$ it should be fx+=dx/(r*r) $\endgroup$ – JMLCarter Dec 23 '16 at 12:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.