IndexError: index out of bounds

Question

I am new to Python and appreciate any constructive feedback. My task is to solve a PDE using ODEINT where my spatial grid is defined as follows:

L = [8.0, 4.0, 8.0]   # Lenght of spatial zones 
dN = 1.0e2            # Number of grid points
N = [int(l*dN) for l in L] # Gives number of grid points over each spatial zone

By using the ODEINT function believe I must define an array of zeros for my solution, and an array of time points. I am struggling to find out how I write the correct size for the array's. I have tried:

dydt = np.zeros(shape=N) # For solution
t = np.linspace(0, 2, 2000) # for time

But that gives me the error message:

IndexError Traceback (most recent call last) in () 70 t = np.linspace(0, 2, lN) 71 ---> 72 U = odeint(numeric, y0, t, args=(h, D1, D2, N)) 73 74 Fexit = U/Np

/Users/severinlangeberg/anaconda/lib/python3.5/site-packages/scipy/integrate/odepack.py in odeint(func, y0, t, args, Dfun, col_deriv, full_output, ml, mu, rtol, atol, tcrit, h0, hmax, hmin, ixpr, mxstep, mxhnil, mxordn, mxords, printmessg) 213 output = _odepack.odeint(func, y0, t, args, Dfun, col_deriv, ml, mu, 214 full_output, rtol, atol, tcrit, h0, hmax, hmin, --> 215 ixpr, mxstep, mxhnil, mxordn, mxords) 216 if output[-1] < 0: 217 warning_msg = _msgs[output[-1]] + " Run with full_output = 1 to get quantitative information."

in numeric(y, t, h, D1, D2, N) 39 40 for i in [N[0] + 1]: # 800 = special 2, must be from begining of N[1] ---> 41 dydt[i] = D2 * (-3*y[i] + 4*y[i + 1] - y[i + 2] / h**2) # pluss 42 43 for i in range (N[0] + 2, N[0] + N[1] - 2): # from 801 to 1197, or 2 to 397

IndexError: index 801 is out of bounds for axis 0 with size 800

Which seems to come from my loops:
# ODE setup
def numeric(y, t, h, D1, D2, N):

    dydt = np.zeros((N))

    # End grid points
    dydt[0] = 0
    dydt[-1] = 0


    # Inner grid points
    for i in range (1, N[0] - 2): # 1 to 797
        dydt[i] = D1 * (y[i - 1] - 2*y[i]+ y[i + 1]) / h**2 # range i = j - 1

    for i in [N[0] - 2]: # 798 = special 1
        dydt[i] =  D1 * (-3*y[i] - 4*y[i - 1] + y[i - 2]) / h**2 # minus

    for i in [N[0] - 1]: # 799 = unknown, and last index in N[0]
        dydt[i] = ((D1*(-3*y[i] + 4*y[i + 1] - y[i + 2])) - (D2*(-3*y[i] - 4*y[i - 1] + y[i - 2])))

    for i in [N[0] + 1]: # 800 = special 2, must be from begining of N[1]
        dydt[i] = D2 * (-3*y[i] + 4*y[i + 1] - y[i + 2] / h**2) # pluss

    for i in range (N[0] + 2, N[0] + N[1] - 2): # from 801 to 1197, or 2 to 397
        dydt[i] = D2 * (y[i - 1] - 2*y[i] + y[i + 1]) / h**2 # range

    for i in [N[0] + N[1] - 1]: # at 1198, or 398 = special 3
        dydt[i] = D2 * (-3*y[i] - 4*y[i - 1] + y[i - 2]) / h**2 # range

    for i in [N[0] + N[1]]: # 1199, or 399 = unknown
        dydt[i] = ((D1*(-3*y[i] + 4*y[i + 1] - y[i + 2])) - (D2*(-3*y[i] - 4*y[i - 1] + y[i - 2])))

    for i in [N[0] + N[1] + 1]: # at 1200, or 1 = special 4
        dydt[i] =  D1 * (-3*y[i + 1] + 4*y[i + 2] - y[i + 3] / h**2)  # pluss

    for i in range (N[0] + N[1] + 2, N[0] + N[1] + N[2] - 1): # 1202 to 1999
        dydt[i] = D1 * (y[i - 1] - 2*y[i]+ y[i + 1]) / h**2 # range

    return dydt

Answer 1

L = [8.0, 4.0, 8.0] 
dN = 1.0e2 
N = [int(l*dN) for l in L] # create array [800, 400, 800] 

dydt = np.zeros(shape=N) # create 3 dimensional matrix 
                         # 1d has size 800 (from 0 to 799 including)
                         # 2d has size 400 (from 0 to 399 including)
                         # 3d same as 1st

In this code:

for i in[N[0]]: 
    dydt[i] = 0 # Unknown point

I is always 800 so it actually same as dydt[ 800 ] = 0. which is bigger then 799.

It's looks like that you use numpy arrays in wrong way.

Yes in low level they represents as contiguous array of ints or floats of whatever. But on python level np.array represents as array of arrays. (ie first element as 3d matrix mtx = np.zeros(shape=[800,400,800] in python will be mtx[0][0][0])).

So all code after for i in[N[0]]: will not work. It tries to access subarray that out of matrix bounds.

UPDATE:

for example you have 3d matrix (3x3x3)

 mtx=np.zero(shape=[3,3,3])
 print(mtx) 
 #[ [ [ 0.  0.  0.]
 #    [ 0.  0.  0.]
 #    [ 0.  0.  0.] ]

 #  [ [ 0.  0.  0.]
 #    [ 0.  0.  0.]
 #    [ 0.  0.  0.] ]

 #  [ [ 0.  0.  0.]
 #    [ 0.  0.  0.]
 #    [ 0.  0.  0.] ] ]

 mtx[0] = 1 
 print( mtx )
 #[ [ [ 1.  1.  1.]
 #    [ 1.  1.  1.]
 #    [ 1.  1.  1.] ]

 #  [ [ 0.  0.  0.]
 #    [ 0.  0.  0.]
 #    [ 0.  0.  0.] ]

 #  [ [ 0.  0.  0.]
 #    [ 0.  0.  0.]
 #    [ 0.  0.  0.] ] ]

 mtx[0][1] = 2
 #[ [ [ 1.  1.  1.]
 #    [ 2.  2.  2.]
 #    [ 1.  1.  1.] ]

 # others are zeros

 mtx[0][1][2] = 3
 #[ [ [ 1.  1.  1.]
 #    [ 2.  2.  3.]
 #    [ 1.  1.  1.] ]

 # others are zeros

Answer 2

From a quick look, it seems the problem is the for loop, instead of using for i in [N[0]], are you sure it's not for i in range(N[0])?

for i in [N[0]] is rather redundant because the list only contains one element, ie. N[0], so the list of size 800 becomes out of index when you try to seek dydt[800] <-- as the max index is 799.

Like so:

In [9]: N
Out[9]: [800, 400, 800]

In [10]: for i in [N[0]]:  # ie. i == N[0] == 800
   ....:     dydt[i] = 0
   ....:     
---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
<ipython-input-10-04daa45ec403> in <module>()
      1 for i in [N[0]]:  # ie. i == N[0] == 800
----> 2     dydt[i] = 0
      3 

IndexError: index 800 is out of bounds for axis 0 with size 800

In [11]: for i in range(N[0]):
   ....:     dydt[i] = 0
   ....: