Numpy Divide By Zero. Why?
Solution 1:
In this code:
defLJ(x,y):
for i inrange(1,10):
for j inrange(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.
Solution 2:
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.
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.
Post a Comment for "Numpy Divide By Zero. Why?"