---------- Forwarded message ---------- From: George Nurser <gnurser@gmail.com> Date: 30 July 2010 22:37 Subject: Re: [Matplotlib-users] Vectorization To: Nicolas Bigaouette <nbigaouette@gmail.com>, Discussion of Numerical Python <numpy-discussion@scipy.org>

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             )

If dx ,dy, dz vary, I don't think this is quite correct. To get dphi/dz at the interface between k & k+1 you need the distance between the psi-points inside the kth and k+1th gridbox. If you assume the psi-points lie at the centres of the gridboxes (other assumptions are possible), this distance is  .5*(dz[k]+dz[k+1]) So, writing rdzw[k] = 2./(dz[k]+dz[k+1]) d/dz(dphi/dz))[k] =  {(phi[k+1]-phi[k])*rdzw[k] - (phi[k]-phi[k-1])*rdz[k-1]} / dz[k]                         =  {phi[k+1]*rdz[k]- phi[k]*(rdz[k] + rdz[k-1]) + phi[k-1]*rdz[k-1]} / dz[k] & similarly for the d2phi/dx2 & d2phi/dy2 terms There are other metric terms which appear in the discretization of the Laplacian if dz depends on i or j, or dy depends on k or i, or dx depends on j or k, but these do not appear in your problem where dz depends only on k, dy only on j and dx only on i. HTH. George Nurser.