# Energy conservativation of a system of particles

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}$$

So:

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

where

$$\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:

• Please explain your setup with more detail. What are Delta_x and Delta_y? How do you simulate the interactions? Dec 23, 2016 at 11:58
• @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. Dec 23, 2016 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.

• The problem has been the step size. Do you mean that with any numerical method particles will finally accelerate in a crazy way? Dec 23, 2016 at 12:33
• 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 Dec 23, 2016 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.

• 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. Dec 23, 2016 at 12:20
• 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}$. Dec 23, 2016 at 12:37
• 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. Dec 23, 2016 at 12:41

What's this bit

                    // F=1/r^2
double r=sqrt(r2);
fx+=dx/(r*r*r);


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.

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