37049 # Numpy divide by zero. Why?

```import pygame import random import math import numpy as np import matplotlib.pyplot as plt fx = np.zeros() fy = np.zeros() x = np.zeros() y = np.zeros() x = 11 y = 1 sigma = 1.01 e = 1.1 dt = 0.1 def LJ(x,y): for i in range(1,10): for j in range(1,10): rx = (x[i]-x[j]) ry = (y[i]-y[j]) fx[i] = 24*e*(((2/rx)*sigma/rx**12)-((1/rx)*sigma/rx**6)) fy[i] = 24*e*(((2/ry)*sigma/ry**12)-((1/ry)*sigma/ry**6)) print fx, fy ```

Why am I still getting errors `RuntimeWarning: divide by zero encountered in double_scalars` and

```RuntimeWarning: invalid value encountered in double_scalars ```

The result I am getting is

```[ 0. nan nan nan nan nan nan nan nan nan 0.] [ 0. nan nan nan nan nan nan nan nan nan 0.] ```

I tried to modify starting x and y but that gave no effect.

In this code:

```def LJ(x,y): for i in range(1,10): for j in range(1,10): ... ```

If `i == j`, you are comparing a particle to itself. Try skipping that iteration of the for-loops like so:

```def LJ(x,y): for i in range(1,10): for j in range(1,10): if i == j: continue rx = (x[i]-x[j]) ry = (y[i]-y[j]) fx[i] = 24*e*(((2/rx)*sigma/rx**12)-((1/rx)*sigma/rx**6)) fy[i] = 24*e*(((2/ry)*sigma/ry**12)-((1/ry)*sigma/ry**6)) ```

Also, you will need to put in actual values for the x and y lists, as they are all currently 0's. Two particles at exactly the same location exert an infinite amount of force according to that equation, so the divide by 0 is accurate in that scenario.

You should stay with the physics. The potential function is, reconstructing from your force formula,

```LJ(r) = 4*e*( (sigma/r)**12 - (sigma/r)**6 ) ```

The force corresponding to the potential is, for any potential,

```F = - LJ'(r)*(rx/r, ry/r) ```

Thus the procedure to compute the combined interaction forces should look like

```def LJ(x,y): for i in range(x.size): for j in range(x.size): rx = (x[i]-x[j]) ry = (y[i]-y[j]) r = sqrt(rx*rx+ry*ry+1e-12) f = 24*e/r*( 2*(sigma/r)**12 - (sigma/r)**6 ) fx[i] += f*(rx/r) fy[i] += f*(ry/r) print fx, fy ```

The condition is that at the call of the procedure the force array is initialized, either to zero or the gravitational forces only.

<hr>

The addition of the term `eps**2` in the computation of the distance serves to avoid zero distances and thus infinite or very large values in the force computation. This epsilon `eps` should of course be small against typical distances of particles resp. the dimension of the scene, here `sigma=1.01` is a typical distance, thus `eps=1e-6` is small against it.