TITLE: Energy conservativation of a system of particles
QUESTION [0 upvotes]: 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:

REPLY [0 votes]: 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.

REPLY [0 votes]: 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.

REPLY [0 votes]: 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.
{"set_name": "stack_exchange", "score": 0, "question_id": 300450}