On Jul 2, 2010, at 2:15 PM, Nicolas Bigaouette wrote: > Hi all, > > I don't really know where to ask, so here it is. > > I was able to vectorize the normalization calculation in quantum > mechanics: <phi|phi>. Basically it's a volume integral of a scalar > field. Using: > norm = 0.0 > for i in numpy.arange(len(dx)-1): > for j in numpy.arange(len(dy)-1): > for k in numpy.arange(len(dz)-1): > norm += psi[k,j,i]**2 * dx[i] * dy[j] * dz[k] > if dead slow. I replaced that with: > norm = (psi**2 * dx*dy[:,numpy.newaxis]*dz > [:,numpy.newaxis,numpy.newaxis]).sum() > which is almost instantanious. > > I want to do the same for the calculation of the kinetic energy: > <phi|p^2|phi>/2m. There is a laplacian in the volume integral which > complicates things: > K = 0.0 > for i in numpy.arange(len(dx)-1): > for j in numpy.arange(len(dy)-1): > for k in numpy.arange(len(dz)-1): > K += -0.5 * m * phi[k,j,i] * ( > (phi[k,j,i-1] - 2.0*phi[k,j,i] + phi[k,j,i+1]) / > dx[i]**2 > + (phi[k,j-1,i] - 2.0*phi[k,j,i] + phi[k,j+1,i]) / > dy[j]**2 > + (phi[k-1,j,i] - 2.0*phi[k,j,i] + phi[k+1,j,i]) / > dz[k]**2 > ) > > My question is, how would I vectorize such loops? I don't know how > I would manage the "numpy.newaxis" code-foo with neighbours > dependency... Any idea? >
I would first create a 3d array of the integrand, probably using scipy.signal.convolve to convolve phi with a kernel such as [[[0,0,0],[0,1,0],[0,0,0]], [[0,1,0],[1,-6,1],[0,1,0]], [[0,0,0],[0,1,0],[0,0,0]]] Then just multiply by whatever factors of dx, dy, dz, and m, and sum the 3d integrand. If dx,dy,dz are non-uniform, it is a harder problem... Hope that helps, -Jeff P.S. Careful, the code you wrote will multiply by the mass instead of dividing by it. ------------------------------------------------------------------------------ This SF.net email is sponsored by Sprint What will you do first with EVO, the first 4G phone? Visit sprint.com/first -- http://p.sf.net/sfu/sprint-com-first _______________________________________________ Matplotlib-users mailing list Matplotlib-users@lists.sourceforge.net https://lists.sourceforge.net/lists/listinfo/matplotlib-users
