[SciPy-User] Laplacian of a matrix (MATLAB del2 function eqv. in NumPy)

[Jonathan Guyer]

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?