When I have those kind of errors I always try to print out the len of more or less everything. Check if the length of Z[1:,1:-1] is correct, and so on.
–
Daniel Thaagaard AndreasenOct 19 '12 at 7:27

Did you try to profile your code with cProfile or something? I think you should identify your real bottleneck first.
–
Spirit ZhangOct 19 '12 at 7:28

The bottleneck is in the call to heat, which takes 99.9% of the program time. So yes, this is where I need to be optimizing and I have read that vectorizing will speed up my code greatly. I'm just off by 1 on the dimensions :/
–
ZackOct 19 '12 at 7:33

@user1493692, sorry its not that trivial to optimize directly (afterall). The for loop relies on previously changed values (is this wanted?), this will not happen if vectorized like this, before doing bigger hacks its, easier to use cython if you don't mind the small compiling step...
–
sebergOct 19 '12 at 9:38

2 Answers
2

You cannot vectorize that diffusion calculation in the time dimension of the problem, that still requires a loop. The only obvious optimization here is to replace the Laplacian calculation with a call to the numpy.diff function (which is precompiled C), so your heat equation solver becomes:

I'm trying to use technicaldiscovery.blogspot.com/2011/06/… as a guide, where they vectorize a double for loop of what looks a lot similar to mine. What does D*q*np.diff(Z[:,i-1], 2) actually take place of in my original code.
–
ZackOct 19 '12 at 8:06

The second difference calculated with diff replaces the inner loop over Z[j-1,i-1]-2*Z[j,i-1]+Z[j+1,i-1]). The blog post you refer to isn't vectorizing a loop like yours, despite the superficial similarities.
–
talonmiesOct 19 '12 at 8:19

To get back to your code:
You are filling Z[1,1:-1] with (first term only) Z[1:-1,:-1]. The mismatch in shape is readily apparent here.

Disregarding the second index (as you will have to loop anyway), your vectorized solution uses a different assumption than the non-vectorized approach: in the first version you have one side u0 (Z[:,0]) and two sides 0 (Z[0,:] and Z[-1,:]) while in the vectorized solution you do try to set Z[-1,:] to something other than 0 by filling Z[1:,i]. Which situation are you trying to simulate?!