Exclude the endpoints of an interval

Question

I have a 2D mesh with an interval of Xleft = -1.5 Xright = 1.5, Ydown = -1.5 Yup = 1.5. Nx and Ny are the variables that indicate how many intervals I want. The problem originates in 1D, so let's go to 1D to simplify.

My interval is [-1.5,1.5] and let's take N=12. My grid-step will be 3/12 = 0.25. I want to output an array with the inner points: (-1.5, 1.5). I use the following function.

for 2D

np.mgrid[intervalX,intervalY]

for 1D

np.mgrid[intervalX] = np.mgrid[start:stop:stepsize]

Because the start will be the literal start and the stop will be the stop - dx. I define the array of inner points as:

LeftX  = -1.5
RightX =  1.5

Nx = 12 

x_int = RightX - LeftX #interval in the x-direction

dx = (RightX - LeftX)  / Nx #grid step in x-direction

matrix0 = np.mgrid[LeftX+dx:RightX:(x_int/Nx)]

the output as expected is:

[-1.25 -1.   -0.75 -0.5  -0.25  0.    0.25  0.5   0.75  1.    1.25]

Most of the time it goes right.

But in some cases with a specific N, the program does something unexpected:

The output for N=9:

[-1.16666667 -0.83333333 -0.5        -0.16666667  0.16666667  0.5
  0.83333333  1.16666667  1.5]

You see that the right boundary point is included in the array. I would expect the same array expect for the most right 1.5. But this is unwanted. The same happens with N=10 or N=22.

• Question: Why does this happen?
• Question: How can I change my code so this doesn't happen anymore?

Answer 1

This is basically a numerical rounding error. Some more information.

What happens inside mgrid[start:end:step] (or arange) can be simplified as something like this:

def mgrid(start, end, step)
    a = [start]
    i = 1
    while start+i*step < end:
        a.append(start+i*step)
        i += 1
    return a

To use your numbers as an example and show where it fails:

LeftX = -1.5
RightX = 1.5
Nx = 9
dx = (RightX - LeftX) / Nx

#    start        + i * step
a8 = (LeftX + dx) + 8 * dx 
# 1.4999999999999998
a8 < RightX
# True

so a8 is numerically slightly smaller than your endpoint 1.5 and is therefore added to the interval. Nevertheless it is displayed as 1.5 because it is pretty close.

As a solution you could use the np.linspace() function. The result includes by default the start and endpoint and you could use just the elements x[1:-1], or shift the start and end point:

import numpy as np
x1 = np.linspace(LeftX, RigthX, Nx+1)x[1, :-1]
x2 = np.linspace(LeftX+dx, RightX-dx, Nx-1)

By including the endpoints explicitly you don't have to worry if the check start + i*step < end gives a wrong result for your last index due to a rounding error.

Alternatively you could use a small epsilon for which you know it is definitely smaller than your step size and do something like this:

eps = 1e-10
np.mgrid[LeftX+dx:RightX-eps:dx]

By subtracting this epsilon you make sure the check start + i*step < end always gives the expected result.