Intuitive interpretation of the Laplacian Operator Just as the gradient is "the direction of steepest ascent", and the divergence is "amount of stuff created at a point", is there a nice interpretation of the Laplacian Operator (a.k.a. divergence of gradient)?
 A: Sure.  One way is to interpret the laplacian of a function at a point, $\Delta f(p)$, as a measurement of the "average bending" of the image at $f(p)$.
Note that in higher dimensions, this is not really very intuitive, due to the "average" part of the "average bending": a function may have all second derivatives negative except one, which is very large and positive, and nevertheless have positive Laplacian at that point.  In such situations one would be better off considering the full Hessian of the function, or the eigenvalues of it, in turn.
A: It behaves like a local averaging operator. This can be seen easily when considering a finite-difference approximation to the Laplacian:
$$\nabla^2 f(x,y) \approx \frac{f(x+h,y) + f(x-h,y) + f(x,y+h) + f(x,y - h) - 4f(x,y)}{h^2}$$
You can rewrite this in vector notation as:
$$\nabla^2 f(x,y) \approx \frac{1}{h^2} \sum_{\mathbf{h}} f(\mathbf{x+h}) - f(\mathbf{x})$$
where the sum is over vectors in the $x$, $-x$, $y$ and $-y$ directions, and $h=||\mathbf{h}||$. Higher order approximations to the Laplacian will involve averaging the rates of change in more directions.
So you can think of the Laplacian as behaving like an 'average rate of change'. As pointed out in Glen Wheeler's answer, the average rate of change can be zero even when there is significant curvature at a point, for example as in the function $f(x,y)=x^2-y^2$.
In image processing, a discrete Laplacian (where $h$ is one pixel in the definition I gave above) can be used as a crude edge-detection filter. It is close to zero in regions where the image is varying smoothly, and has large values in regions where the image has sharp transitions from low to high intensity.
In physics, the Laplacian is interpreted as a diffusion operator, as in the equation
$$\frac{\partial u}{\partial t} = \nabla^2 u$$
This says that the rate of change of $u$ in time is given by the average rate of change of $u$ in space. If we interpret $u$ as a temperature (and hence $\partial u/\partial t$ is the rate of change of temperature) then we see that there is more heat exchange in regions where the temperature is highly variable, and less heat exchange when the temperature varies smoothly.
A: Although this thread is quite old, I think a quick picture would add up to the excellent answers here. An explanation is also below. I've tried to be as non-technical as I could get, since this would also add up intuition to the other more technical answers. Anyway, let's see how come that the Laplacian is a local averaging operator.

So, let's imagine the points $A$, $B$, $C$, $D$ and $O$ are some points on a metal sheet and we have an associated temperature at each (call it $a$, $b$, $c$, $d$, $o$).
Cool.
Now we want to learn something about how this temperature is distributed on the sheet and specifically on these five points.
Specifically, how fast the temperature changes from $A$ to $O$? Well, this should be the total change in temperature from $A$ to $O$ (call it $o-a$) divided by the distance from $A$ to $O$ (call it $h$, and the distance from $A$ to $C$ - $2h$). Thus, as we are going from $A$ to $O$, the temperature changes with a rate of $(o-a)/h$ per unit distance. Similarly, we can calculate $(b-o)/h$, $(c-o)/h$ and $(o-d)/h$. For brevity, let these be $a', b', c'$ and $d'$.
What these tell us is whether the temperature is increasing in a given direction (and how fast). For example, if $a'>0$ then the temperature at $O$ is higher than the temperature at $A$, and it is increasing as we go from $A$ to $O$ at a rate of $(o-a)/h$. Obviously, as $h$ tends to zero, this becomes the first derivative.
...but we can do more.
What if we wanted to know how the change of temperature changes as we move on the metal sheet? In other words, the "second derivatives"? (Formally, we are still "discrete", so these are not exactly derivatives, but let's call them so for convenience.)
Now let's find this for the temperature between $A$ and $C$ by taking the average change of $a'$ and $c'$: $$\frac{(a'+b')}{2h} = \frac{1}{2}\frac{\frac{c-o}{h}-\frac{o-a}{h}}{h} = \frac{1}{2}\frac{a-2o+c}{h^2}$$
We can calculate the second derivative between $B$ and $D$ analogously.
This doesn't seem too big of a step toward the Laplacian, but let's continue by adding what happens in the $y$ direction. Which is, let's add the second derivative in the y-direction to the second derivative in the x-direction, giving: $$\frac{1}{2}\frac{a-2o+c}{h^2}+\frac{1}{2}\frac{b-2o+d}{h^2}=\frac{1}{2h^2}(a+b+c+d-4o)$$
The last thing is the sum of the two second partial derivatives, which is - the Laplacian. Let's see what happens when it is equal to zero... $$\frac{1}{2h^2}(a+b+c+d-4o)=0$$discarding the $1/2h^2$ term: $$a+b+c+d-4o=0$$
...which gives:$$o=\frac{a+b+c+d}{4}$$
Wait... but this means that saying that the Laplacian being equal to zero is the same as saying that the temperature at $O$ is equal to the average of the temperatures at $A,B,C$ and $D$! Exactly as we wanted.
A: The laplacian is also the trace of the hessian matrix (the matrix of second-order partial derivatives). Since the trace of a matrix is invariant under a change of basis, then the laplacian does not change if you do a change of basis. For instance if you work in $\mathbb{R}^2$, then $$\Delta f=\frac{\partial^2 f}{\partial x^2}+\frac{\partial^2 f}{\partial y^2}$$ in cartesian coordinates and $$\Delta f=\frac{1}{r}\frac{\partial f}{\partial r}+\frac{\partial^2 f}{\partial r^2}+\frac{1}{r^2}\frac{\partial^2 f}{\partial\theta^2}$$ in polar coordinates but at each point $(x,y)$ or $(r,\theta)$ this is the same value.
[The hessian matrix contains information about the way your function bends around some point. Since it is symmetric (by the Schwarz theorem) and real, you can diagonalize it and so the laplacian is also the sum of the eigenvalues of the hessian, and these eigenvalues are important to detect which kind of bending you have locally.]
A: Assume the function $f$ is $C^2$ in a neighbourhood of ${\bf 0}\in{\mathbb R}^n$. Using the Taylor expansion of $f$ at ${\bf 0}$ one can prove the following:
$$\Delta f({\bf 0}) =\lim_{r\to 0+}\ {2n\over r^2}{1\over \omega(S_r)} \int\nolimits_{S_r}\bigl(f({\bf x})-f({\bf 0})\bigr)\ {\rm d}\omega({\bf x})\ .\qquad (*)$$
This formula says that $\Delta f({\bf 0})$ is "essentially" (i.e., up to the scaling factor ${2n \over r^2}$) equal to the average difference $f({\bf x})-f({\bf 0})$ over small spheres around ${\bf 0}$.
Using this interpretation one gets, e.g., an intuitive understanding of the heat equation
$${\partial u\over\partial t}=a^2\ \Delta u\ ,$$
namely: If averaged over small spheres around a point ${\bf p}$ it is hotter than at ${\bf p}$ itself, then in the next second the temperature at ${\bf p}$ will rise.
Given the interest in the above formula $(*)$, here are some hints for the proof: By Taylor's theorem one has
$$f({\bf x})-f({\bf 0})= \sum_{i=1}^n f_{.i} x_i +{1\over2}\sum_{i,k} f_{.ik} x_i x_k + o(|{\bf x}|^2)\qquad({\bf x}\to{\bf 0}) .$$
Here the $f_{.i}$ and the $f_{.ik}$ are the partial derivatives of $f$ evaluated at ${\bf 0}$, whence constants. Now we integrate this over $S_r$ with respect to the surface measure ${\rm d}\omega$  and  obtain
$$\int\nolimits_{S_r}\bigl(f({\bf x})-f({\bf 0})\bigr)\ {\rm d}\omega({\bf x})={1\over2}\sum_{i}f_{.ii}\int\nolimits_{S_r}x_i^2\ {\rm d}\omega({\bf x}) +o\bigl(r^{2+(n-1)}\bigr) \qquad(r\to0+)\ ,$$
because all other terms are odd in at least one variable. The integrals $\int\nolimits_{S_r}x_i^2\ {\rm d}\omega({\bf x})$ are all equal; therefore we have 
$$\int\nolimits_{S_r}x_i^2\ {\rm d}\omega({\bf x})={1\over n}\int\nolimits_{S_r}\sum_k x_k^2\ {\rm d}\omega({\bf x})={r^2\over n}\omega(S_r)\qquad(1\leq i\leq n)\ .$$
Putting it all together we obtain
$$\int\nolimits_{S_r}\bigl(f({\bf x})-f({\bf 0})\bigr)\ {\rm d}\omega({\bf x})={r^2\over 2n}\omega(S_r)\Delta f({\bf 0}) +o(r^{n+1})\qquad(r\to 0+)\ ,$$
and solving for $\Delta f({\bf 0})$ we get the stated formula.
For a proof using Gauss' theorem see here:
Nice way of thinking about the Laplace operator... but what's the proof?
A: From a probabilistic point of view the Laplacian is of great importance as it arises a the (infinitesimal) generator of the Brownian motion (even in manifolds).
A $d$-dimensional Brownian motion (or Wiener process) is a(n almost surely) continuous stochastic process with stationary and independent increments such that $B_t - B_s \sim \mathcal N(0,(t-s))^{\otimes d}$ are normally distributed, i.e., have a Gaussian density
$$p_{t-s}(x,\mathrm dy) = \frac 1{\sqrt{2\pi(t-s)}^d} \exp\left( -\frac12 \frac{|x-y|^2}{t-s} \right) \mathrm dy.$$
The (infinitesimal) generator of a strongly continuous contraction semigroup $(P_t)_{t\geq 0}$ on on a Banach space is the operator
\begin{align}
        \mathbf A u := \lim_{t\to 0} \frac{P_t u - u}{t},
\end{align}
where the domain $\mathcal D(\mathbf A)$ of $\mathbf A$ is the set of functions $u$ where the limit exists. By Taylor's formula this can be alternatively written as
\begin{align*}
    P_t u = u + t \mathbf A u + o(t),
\end{align*}
so that the generator $\mathbf A$ gives the first-order approximation to $P_t$ for small $t$. One can show that
\begin{align}
    \frac{\mathrm d}{\mathrm d t} P_t u(x) = \mathbf A P_t u(x) = P_t \mathbf A u(x),
\end{align}
i.e. the generator is the time derivative of the mapping $t\mapsto P_t u(x)$. Reading this as a partial differential equation we see that $u(t,x) := P_t f(x)$ is a solution of
\begin{align}
    \partial_t u(t,x) = \mathbf A u(t,x), \quad u(0,x) = f(x).
\end{align}
From a probabilistic point of view one defines the transition semigroup 
$$P_t u(x) := \mathbb E^x u(B_t) = \frac 1{\sqrt{2\pi t}^d} \int_{\mathbb R^d} \exp\left( - \frac{|x-y|^2}{2t} \right) \mathrm dy.$$
Thus $P_t u(x) = u * p_t(x)$ is a convolution of the heat kernel, thus a solution of the heat equation. In other words
\begin{align}
    \mathbb E^x u(X_t) \approx u(x) + \mathbf A P_t u(x).
\end{align}
So, essentially, the generator describes the movement of the process in an infinitesimal time interval.
One can show that for every Feller semigroup the limit exists in $(C_\infty, \Vert \cdot \Vert_\infty)$ and hence for every Markov process. The more important probabilistic reason is that for every adapted Feller process $(X_t, \mathcal F_t)_{t\geq 0}$ on $\mathbb R^d$ the process 
\begin{align}
    M_t^u := u(X_t) - u(x) - \int_0^t \mathbf A u(X_r) \mathrm d r \quad ( u \in \mathcal D(\mathbf A), {t\geq 0} )
\end{align}
is an $\mathcal F_t$-martingale. Since martingales have constant expectation it follows that
\begin{align*}
    \mathbb E^x u(X_t) - u(x) &= \mathbb E^x \int_0^t \mathbf A u(X_r) \mathrm d r\\
                        &= \int_0^t \mathbb E^x(\mathbf A u)(X_r) \mathrm d r\\
    \Longleftrightarrow P_t u(x) - u(x) &= \int_0^t P_r \mathbf A u(x) \mathrm d r,
\end{align*}
from which we easily attain Dynkin's formula by Doob's optional stopping theorem: Let $\sigma$ be an $\mathcal F_t$ stopping time with $\mathbb E^x \sigma < \infty$ then
\begin{align}
    \mathbb E^x u(X_\sigma) - u(x) = \mathbb E^x \int_0^\sigma \mathbf A u(X_r) \mathrm d r    \quad (u \in \mathcal D(\mathbf A)).
\end{align}
This connection allows us to solve classical partial differential equations, like the Dirichlet problem, with martingale methods in an elegant probabilistic way. 
Finally, back to the Brownian motion. For simplicity $d=1$. Then $\mathbb E^x(B_t - x) = 0$ and $\mathbb E^x(B_t -x)^2 = t$ and by Taylor's formula
\begin{align*}
    \mathbb E^x u(B_t) \approx \mathbb E^x\left( u(x) + u'(x)(B_t-x) + \frac 12 u''(x)(B_t-x)^2 \right) = u(x) +  0 + \frac 12 t u''(x).
\end{align*}
So we might assume that it holds $(\mathbf A, \mathcal D(\mathbf A)) = \left( \frac 12 \Delta, \mathcal C_\infty\right)$, where $\mathcal C_\infty :=: \mathcal C_\infty(\mathbb R^d)$ denotes the family of all continuous functions $f : \mathbb R^d \to \mathbb R$ vanishing at infinity. For $d=1$ that is true.
    For Brownian motion in higher dimensions, it can still be shown that $\mathbf A u = \frac 12 \Delta u$ where $u \in \mathcal D\left(\frac 12 \Delta \right)$. However, we only get $C_\infty \subsetneq \mathcal D\left( \frac 12 \Delta \right)$.
Moveover, we see that
    \begin{align*}
        M_t^u := u(X_t) - u(x) - \int_0^t u''(B_r) \mathrm d r \quad ( u \in \mathcal D(\mathbf A), {t\geq 0} )    
    \end{align*}
is a martingale. 
Thus, we see that $u(t,x) := \mathbb E^x f(B_t)$ is the unique solution of the heat equation
    \begin{align*}
        \partial_t u(t,x) = \frac{1}{2}\partial_x^2 u(t,x), \quad u(0,x)=f(x).
    \end{align*}
Let $D$ be an open, bounded and connected domain, $\tau_{D^c} := \inf\{ t>0 : B_t \not\in D\}$ the first exit time of the set $D$ and $f : \partial D \to \mathbb R$ continuous. If $\partial D$ suffices some regularity condition, then the Dirichlet problem 
    \begin{align*}
     \Delta u(x) &= 0  &\forall x \in D\\
     u(x) &= f(x)   &\forall x \in \partial D\\
     u(x) &\text{ continuous} &\forall x \in \overline D
    \end{align*}
    has the unique solution $u(t,x) := \mathbb E^x f\big(B_{\tau_{D^c}}\big)$.
The Laplacian as generator also plays a key role when defining Brownian motion on manifolds.
A: In addition to other answers, you also should pay attention to the fact that there is a connection between the Laplacian and the Gaussian kernel of the form 
$$
\frac{1}{\sqrt{2\pi t}}e^{-\|x\|^2/2t}
$$
which in fact one of solutions of the equation $u_t = \Delta u$ mentioned by Chris. This kernel has a lot of applications to the averaging since it has a well-known "Bell shape" and is a density function for a normal (Gaussian) distribution.
