[SciPy-User] Laplacian of a matrix (MATLAB del2 function eqv. in NumPy)
Jonathan Guyer
guyer@nist....
Thu Apr 11 09:48:34 CDT 2013
On Apr 11, 2013, at 8:48 AM, Juan Luis Cano wrote:
> Hello everybody, I would like to calcute de Laplacian Operator of a matrix with spacing between points and if it were possible the same boundary conditions as the function del2 does in Matlab. I am already aware of this SO question
>
> http://stackoverflow.com/q/4692196/554319
>
> but scipy.ndimage.filters.laplace() is not suitable because a) none of its modes contemplates extrapolation on the boundary and b) doesn't allow change spacing between the points. Another person suggested a direct translation of the code but doesn't work correctly. This translation option is a bit bothering because the original code is somewhat difficult to follow. Is there any other possibility?
>
> I wish you could help me. Thanks in advance.
It's overkill for this, but FiPy can calculate a Laplacian and will let you adjust the spacing between points, but it calculates the same laplacian that scipy.ndimage.filters.laplace() does, i.e.,
>>> import fipy as fp
>>> from fipy import numerix as nx
>>> mesh = fp.Grid2D(nx=4, ny=4)
>>> var = fp.CellVariable(mesh=mesh, value=[3, 4, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19])
>>> print var.faceGrad.divergence.reshape((4,4))
[[ 6. 6. 3. 3.]
[ 0. -1. 0. -1.]
[ 1. 0. 0. -1.]
[-3. -4. -4. -5.]]
The boundary cell values result from a zero gradient condition on the bounding faces.
The documentation for del2() at http://octave.sourceforge.net/octave/function/del2.html says "Boundary points are calculated from the linear extrapolation of interior points", which is a zero curvature condition.
The del2() documentation also says "For a 2-dimensional matrix M this is defined as
1 / d^2 d^2 \
D = --- * | --- M(x,y) + --- M(x,y) |
4 \ dx^2 dy^2 /"
The division by 4 indicates they're calculating a 1st order laplacian, whereas ndimage and FiPy are calculating a 2nd order laplacian. 2nd order is more accurate.
As asked on the StackOverflow thread, why is it important to get the same answer that Octave gives?
More information about the SciPy-User
mailing list