4

This is a bit of a strange issue for me and I wasn't sure how to properly title the question. I have the follow MWE which simply generates a list of coordinate points (x,t) and performs some check to see if they lie on a user-prescribed boundary. In particular, if x[i] == 1.0 and t[i] != 0.0 then the program should print a statement stating so. I can't seem to figure out why the if conditional is never entered here. I have printed out the pairs of values x[i] and t[i] to verify that there are indeed pairs satisfying the conditional...

#Load Modules
import numpy as np
import math, random
from pylab import meshgrid

# Create the arrays x and t on an evenly spaced cartesian grid
N = 10
xa = -1.0;
xb = 1.0;

ta = 0.0;
tb = 0.4;

xin = np.arange(xa, xb+0.00001, (xb-xa)/N).reshape((N+1,1))
tin = np.arange(ta, tb+0.00001, (tb-ta)/N).reshape((N+1,1))

X_tmp,T_tmp = meshgrid(xin,tin)
x = np.reshape(X_tmp,((N+1)**2,1))
t = np.reshape(T_tmp,((N+1)**2,1))

# create boundary flags
for i in range(0,(N+1)**2):
    if (x[i] == xb and t[i] != ta):
        print("We are on the right-side boundary")
John Kugelman
  • 349,597
  • 67
  • 533
  • 578
user1799323
  • 649
  • 8
  • 25

1 Answers1

5

I think you're running into floating point precision issues. So while x[i] is very near, it is not exactly equal to xb. Perfect equality tests cause trouble like this with floating-point numbers. What you want is to test that the difference between these values is small. Try this:

ep = 1e-5 # choose this value based on how close you decide is reasonable
for i in range(0,(N+1)**2):
    if (abs(x[i] - xb) < ep and abs(t[i] - ta) > ep):
       print("We are on the right-side boundary")

Also I just learned Python 3.5 added isclose functions useful for such cases! See this question/answer for more discussion. Also note that if you want to do this for an array, NumPy provides the allclose function.

aganders3
  • 5,838
  • 26
  • 30