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

[Juan Luis Cano]

On Thu, Apr 11, 2013 at 4:48 PM, Jonathan Guyer <guyer at nist.gov> wrote:

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

Thank you for all of this insight, actually I am asking this question for
another person so I don't know all the details but this is very useful.
Anyway, even though 2nd order is more accurate I was not able to equal the
result of this last example (MATLAB):

http://www.mathworks.es/es/help/matlab/ref/del2.html#f81-998436


> As asked on the StackOverflow thread, why is it important to get the same
> answer that Octave gives?
>

Point is most of the times porting a code from MATLAB, as usual. Anyway
being able to specify the spacing is more important than the boundary
itself.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20130413/4ff700a7/attachment.html>