Encountering divide by zero after a conditional preventing this

Question

I am writing a simple piece of code for a physics assignment where we are supposed to model a water molecule as a collection of charged particles. I am supposed to plot the electrostatic potential of the molecule, calculated by the formula: V = q/(k*r) where r is the distance from the center of the particle to a point.

Inside the particle however, the following formula applies: V = (r/R)^2 * q/(k*R) where R is the radius of the particle (which is defined as a non- zero number early in the code). For simplicity, k and q are set equal to 1.

The first formula has the possibility to encounter a divide by zero, but the second one does not. My code works fine, and the plot looks correct, but the code tells me that at some point it is running into a divide by zero in the following line of code:

V = np.where(r<R,
        V + rSquared/RSquared * q/(k*R),
        V + q/(k*r))

The conditional should ensure that if r is a very small number, it uses the formula that doesn't divide by r, theoretically making divide by zeros impossible.

To be more specific, I am using a for-loop that calculates V for each particle in a list of nested lists containing information about each particle such as their position, q and R. rSquared is a two- dimensional matrix generated by performing a calculation of the magnitude of the vector between the position of the particle and points on a meshgrid (but without taking the square root).

Does anyone know how this divide by zero could happen? My code runs and plots the field just fine, but this error message annoys me, since it shouldn't happen as far as I know.

Edit: Thanks for the answers! I was assuming that the numpy.where() function doesn't calculate q/(k* r) because of the conditional, but this is not true. Turns out it calculates both q/(k* r) and rSquared/RSquared * q/(k*R), and later assesses which result to use, which is why the divide by zero happens.

Answer 1

As commented by @WarrenWeckesser, the three arguments are evaluated, causing warnings to appear.

But Numpy does a good job of returning not-a-number or infinity where appropriate, which are then properly used in further calculations. For example:

>>> a = np.array([0, -1, 1])
>>> result = np.log(1 / a)
RuntimeWarning: divide by zero encountered in true_divide
RuntimeWarning: invalid value encountered in log
>>> result
array([inf, nan,  0.])

Results are valid, but warnings can be annoying. If you expect a section of your code to produce such warnings, you can silence them using errstate:

>>> with np.errstate(divide='ignore', invalid='ignore'):
...     result = np.log(1 / a)
>>> result
array([inf, nan,  0.])

You can use more convoluted methods to avoid warnings instead of silencing them, but this normally involves additional computations just to get the same effect.