---------- Forwarded message ----------
From: George Nurser <gnur...@gmail.com>
Date: 30 July 2010 22:37
Subject: Re: [Matplotlib-users] Vectorization
To: Nicolas Bigaouette <nbigaoue...@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.
_______________________________________________
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion
