GSEL GSEL - 5 months ago 31
Python Question

IndexError: index out of bounds

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

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
   ....: