Help understanding how to compute Laplace operator? I'm working on a diffusion simulation, and I am struggling with how to calculate the Laplace operator for my density function as it applies to Flick's second law.
Conceptually I believe it is clear: given a density field, the rate of change in density at each point in the field is proportional to the current density times the Laplacian value at that point.
This makes intuitive sense because the Laplacian is a function of the gradient: i.e. if the surrounding points on the gradient point toward a given point, that point will have a decreasing density.  In other words, the Laplacian is positive when the neighborhood of a point is more dense, and negative when the neighborhood is less dense, which is exactly what we want for modeling diffusion.
So where I'm blocked is in how I can implement this.  
My density function is modeled by a 2D array of values.  At each time point I want to calculate the rate of change at each point in the array, and use that to compute the value at the next time point.
Computing the gradient seems easy enough: I can get a directional difference in density for each of the neighbors and add them together.  Where I'm lost is how to compute the laplacian from the gradient values.
 A: Assuming you are using a Cartesian coordinate system, the 2 dimensional Laplacian operator is
$$\nabla^2 f = \left( \frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2} \right)f$$
Since you are doing a simulation, you will need to numerically estimate this.
We can estimate a derivative as a finite difference.
$$\frac{\partial}{\partial x}f(x_0,y_0) \approx \frac{1}{\delta}\left(f(x_0+\delta,y_0) - f(x_0,y_0)\right)$$
or centered on the point
$$\frac{\partial}{\partial x}f(x_0,y_0) \approx \frac{1}{\delta}\left(f(x_0+\delta/2,y_0) - f(x_0-\delta/2,y_0)\right)$$
For the Laplacian you need the second derivative, so the "difference of the differences".
$$\frac{\partial^2}{\partial x^2}f(x_0,y_0) \approx \frac{1}{\delta}\left(\frac{\partial f}{\partial x}(x_0+\delta/2,y_0) - \frac{\partial f}{\partial x}(x_0-\delta/2,y_0)\right)$$
$$ \approx \frac{1}{\delta}\left(\frac{1}{\delta}\left(f(x_0+\delta,y_0) - f(x_0,y_0)\right) - \frac{1}{\delta}\left(f(x_0,y_0) - f(x_0-\delta,y_0)\right)\right)$$
$$ = \frac{1}{\delta^2}\left(f(x_0+\delta,y_0) - 2 f(x_0,y_0)+ f(x_0-\delta,y_0)\right) $$
Do similarly for the second derivative in the $y$ direction.
