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