Skip to content Skip to sidebar Skip to footer

Numpy Divide By Zero. Why?

import pygame import random import math import numpy as np import matplotlib.pyplot as plt fx = np.zeros([11]) fy = np.zeros([11]) x = np.zeros([11]) y = np.zeros([11]) x[0] = 11 y

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?"