0

(python 2) I have this code that takes my arrays, x1,y1,z1, vx1, vy1,vz1, and operates on them (this is the bulk of the code), and at the end I'm left with new arrays x2,y2,z2, vx2,vy2,vz2. What I want to do is figure out how to loop the whole code, but update x1,y1,...etc with x2,y2,... giving me x3,y3.... etc until I have xn,yn,... etc.

I tried to do this by using, 

 for timestep in xrange(0,1000):

but this just runs the whole program 1000 times giving me x2,y2,... 1000 times. In addition, it's really slow. However, what I want to do is get x1000,y1000,.... preferably using Numpy so my code also runs more quickly.

#x1,y1,z1,vx1,vy1,vz1,mass1,and n1 are all arrays with npoints elements

for timestep in xrange(0,1000):     
     M = np.zeros((npoints,npoints))

     def do_work(xn, step):

         #Not my actual function, but to give you an idea
         M[xn,step] = x1[step]+y1[step]+z1[step]+x1[xn]+y[xn]+z[xn]

     #the rest of this is all just more code to create my new arrays.
     [do_work(xn, step) for (xn,step) in itertools.product(xrange(npoints), xrange(npoints))]



     a=[np.sum(arr) for arr in M]


     a = np.array(a)
     vxx = np.array(vx1)
     vyy=np.array(vy1)
     vzz=np.array(vz1)
     vx=vxx[0:npoints]
     vy=vyy[0:npoints]
     vz=vzz[0:npoints]

     vx2 = vx + (a +a)/2  * dt
     vy2 = vy + (a +a)/2  * dt
     vz2 = vz + (a+a)/2   * dt

     xx = np.array(x1)
     yy = np.array(y1)
     zz = np.array(z1)

     x=xx[0:npoints]
     y=yy[0:npoints]
     z=zz[0:npoints]


     x2= (x+vx2*dt) + (a*dt**2)/2
     y2= (y+vy2*dt) + (a*dt**2)/2
     z2= (z+vz2*dt) + (a*dt**2)/2


#plotting and printing 


#print x1
#plt.scatter(x2,y2)

plt.show()
Hoog
  • 2,280
  • 1
  • 14
  • 20
Mckenzie
  • 25
  • 1
  • 1
  • 8

3 Answers3

2

Your problem is that you are overwriting your x2, etc in every iteration. I suggest to only calculate your deltas in the loop and add them either to your x1, etc or create a copy of x1, etc before the loop e.g. as x_final, etc and add up your deltas cumulatively in the loop.

#x1,y1,z1,vx1,vy1,vz1,mass1,and n1 are all arrays with npoints elements

#x_final, y_final as copy

for timestep in xrange(0,1000): 
    # do your calculations

    x_final = x_final + dx
M. Schlenker
  • 81
  • 3
1

for timestep in xrange(0,1000):
M = np.zeros((npoints,npoints))

The start of your code above is initializing M every for loop iteration. Put M outside of the for loop if you are intending to write to it as your answer to avoid this.

Consider using the enumerate to drill down to where you want to calculate. I find it to be useful because it provides access to the index you are on, which is helpful for calculating prior and current timesteps.

See this stack overflow post: How to iterate 1d NumPy array with index and value

An example of this construct you would use:

npoints = 1000
M = np.zeros((npoints,npoints))

def add_to_number(prev_num):
    new_num = prev_num + 1
    return new_num

for i, row in enumerate(M[0:-1]):
    # this is a row slice of M
    for j, value in enumerate(row[0:-1]):
        # this is a value slice of the above row
        # having position i, j
        M[i+1, j+1] = add_to_number(value)

you will see that it will go through M row by row, value by value, and add numbers along M. Hope that helps!

Kenan
  • 41
  • 1
  • 2
0

If you don't need all of the x2, y2, z2, x3, ... in memory at the end you could add in a few lines at the end

x1 = x2
y1 = y2
z1 = z2

This way on the second pass through the loop you will have updated X1, Y1, Z1 values as input, and you will just write over x2, y2, z2 with (what one would consider) x3, y3, z3. Those guys are loaded into x1, y1, z1 and the loop continues.

If you do need all the intermediate x, y and z values it may be worth printing them to file at each step so you can access them later.

Hoog
  • 2,280
  • 1
  • 14
  • 20
  • For some reason when I do this, it changes the size of some arrays. I get the error: index 999 is out of bounds for axis 0 with size 999 – Mckenzie Feb 08 '19 at 19:17