How to compute laplacian of a field?

Question

I'm trying to compute the laplacian of a 2d field A using scipy.ndimage.convolve.

stencil = numpy.array([[0, 1, 0],[1, -4, 1], [0, 1, 0]])
scipy.ndimage.convolve(A, stencil, mode='wrap')

This doesn't seem to give me the right answer though. Any ideas where I'm going wrong, or are there better ways of computing the laplacian in numpy?

Answer 1

did you try another laplace convolution kernel, like [[1,1,1][1,-8,1][1,1,1]] ?

Answer 2

I tried your example, and it does work with the latest SciPy on my machine.

I would suggest that you plot image-ndimage.convolve(…) in order to see how the convolution changes your image.

Example:

image = vstack(arange(-10, 10)**2 for _ in range(20))
ndimage.convolve(image, stencil, mode='wrap')

yields

array([[-38,   2,   2,   2,   2,   2,   2,   2,   2,   2,   2,   2,   2,...)

which is quite correct (the second order derivative of x**2 being 2)—except for the left border.

Answer 3

I got another idea: did you take into account that your stencil, in order to approximate the Laplacian, should be divided by step**2, where step is the step size of your grid? Only then can you compare the ndimage.convolve result with the analytical result.

In fact, with a Gaussian, I obtain results that indicate that ndimage.convolve works quite well:

from scipy import ndimage

stencil = numpy.array([[0, 1, 0],[1, -4, 1], [0, 1, 0]])
x = linspace(-10, 10, 100)
y = linspace(-10, 10, 100)
xx, yy = meshgrid(x, y)
image = exp(-xx**2-yy**2)  # Standard deviation in x or y: 1/sqrt(2)

laplaced = ndimage.convolve(image, stencil)/(x[1]-x[0])**2  # stencil from original post
expected_result = -4*image + 8*(xx**2+yy**2)*image  # Very close to laplaced, in most points!