# [Numpy-discussion] Fwd: [Matplotlib-users] Vectorization

George Nurser gnurser@gmail....
Fri Jul 30 16:59:52 CDT 2010

```---------- 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.
```