Harmonic function: definitions and properties


This is the second part of my tutorial series on bounded harmonic functions. For a quick introduction and examples of use of harmonic functions read the first part. In this part I define harmonic functions and their properties. This is the hard part with a lot of mathematics. But it's a mandatory step to understand how harmonic functions work. This will allow you to apply them in a broader context and understand many scientific papers relying on these. In addition this will be my entry point to introduce Finite Element Method in future posts. So hang on it's worth it!

Harmonic function definition and notations

In this section I will introduce the vocabulary and notions needed for the rest of the post.

So what is an harmonic function? It's a function \( f: \mathbb R^n \rightarrow \mathbb R\) which second derivatives sum to zero. So it must satisfy the equation:

\begin{equation} \frac{\partial^2f}{\partial x_1^2} + \frac{\partial^2f}{\partial x_2^2} + \cdots + \frac{\partial^2f}{\partial x_n^2} = 0 \label{eq:laplace} \end{equation}

Which for a 2D function \( f(x,y) \) would be \( \frac{\partial^2 f}{\partial x^2} + \frac{\partial^2 f}{\partial y^2} = 0 \) or \(f_{xx} + f_{yy} = 0\) depending on how you like to write your derivatives. This equation is a differential equation.  However it is a special type of differential equation. Because \( f \) has multiples variables it is called a partial differential equation (PDE). Among the numerous PDEs the above equation is also named the Laplace equation.

I showed earlier how to compute with a diffusion process discrete harmonic functions. In these examples if you pick any point \( (x, y) \) then \( f(x,y) \) has to verify the Laplace equation \(f_{xx} + f_{yy} = 0\). We'll see every points do, even if for some it is not so obvious, but we'll get back on that.

What we have done using the diffusion process is: solving numerically the Laplace equation \(f_{xx} + f_{yy} = 0\) (a PDE) for every points of the grid. It is a numerical solution since the solution \( f \) is represented with a discrete grid of finite numerical values. 

Values at the boundary are sometime called constraints,  they're fixed. By opposition values inside the boundary (interior domain) are free or unknown.

Laplace operator

Usually the Laplace equation (\ref{eq:laplace}) is written using the Laplace operator denoted \( \Delta \). So the Laplace equation is equivalent to:

$$ \Delta f = 0 $$

\( \Delta \) is called a differential operator. Much like a minus sign '\( - \)' is applied to a value to negate it, differential operators are applied to functions and differentiate them. And so, \( \Delta = \frac{\partial^2}{\partial x_1^2} + \frac{\partial^2}{\partial x_2^2} + \cdots + \frac{\partial^2}{\partial x_n^2}\) and when applied \( \Delta f = \frac{\partial^2f}{\partial x_1^2} + \frac{\partial^2f}{\partial x_2^2} + \cdots + \frac{\partial^2f}{\partial x_n^2}\)

N.B: Laplace operator \( \Delta \) can be defined with the gradient \( \nabla \) and divergence \( \nabla \mathbf . \) (please notice the dot) differential operators:

$$ \nabla \mathbf . \nabla f = \Delta f $$

The gradient \( \nabla = \begin{bmatrix} \frac{\partial }{\partial x_1}, & \frac{\partial}{\partial x_2}, & \cdots &, \frac{\partial}{\partial x_n} \end{bmatrix} \) is a vector of first derivatives. for \( f: \mathbb R^n \rightarrow \mathbb R\) its dimension equal the number of arguments. So when applied we have: \( \nabla f: \mathbb R^n \rightarrow \mathbb R^n \).

The divergence \( \nabla \mathbf . \) (don't forget the dot) return a scalar value. It is meant to be applied to a vector field \( f: \mathbb R^n \rightarrow \mathbb R^n\):

$$f(x_1, \cdots, x_n) =
    f_1(x_1, \cdots, x_n) \\
    \vdots \\
    f_n(x_1, \cdots, x_n)  
\end{bmatrix} $$

It computes \( \nabla \mathbf . f = ( \frac{\partial f_1 }{\partial x_1} + \frac{\partial f_2 }{\partial x_2} + \cdots + \frac{\partial f_n }{\partial x_n} ) \).

Some properties

The following will present some properties of harmonic functions. This will help you foresee and understand the shape of the function given a specific boundary condition.

Minimal surfaces

My preferred way to understand and predict the shape of an harmonic function is the analogy with soap bubbles or more precisely soap films:

helicoid soap
catenoid soap
Miscelaneous soap film
Helicoid Catenoid Miscellaneous
Courtesy of © The Exploratorium.
Courtesy of soapbubble.dk.
(license CC BY-NC 3.0)
Courtesy of math.hmc.edu

 The soap film have the particularity to be minimal surfaces. Depending on the structure the soap film is spread onto, the film shape will try to have the smallest possible surface in order to minimize the overall tension.

So what does it has to do with harmonic functions? Well minimal surfaces have their mean curvature null everywhere i.e they respect the Laplace equation (\ref{eq:laplace}). When you visualize in three dimension a 2D harmonic function (Figure 1 or 2) you will see that indeed harmonic functions are minimal surfaces and look exactly like soap films.

You like the analogy see this great video about minimal surfaces.

1D functions

Lets take a look at the solution of the Laplace equation in 1D to better understand how it behaves. In 1D we seek \( f : \mathbb R \rightarrow \mathbb R \) which satisfy \( f''(x) = 0 \). It isn't that hard to find out:

$$ \begin{eqnarray*} \iint f''(x) & = & \iint 0 \\ \int f'(x) & = & \int a \\ f(x) & = & ax + b \\  \end{eqnarray*} $$

This analytical solution tells us something interesting: 1D harmonic functions are always linear. This explains our earlier experiment (Figure 1) when diffusing values on a 1D array it converges to a straight line. Another way to see it, is that a 1D function is a curve, and having \( f''(x) = 0 \) everywhere means it has no curvature, hence the straight line.

2D functions

The 2D version is a bit more tricky. To satisfy \(f_{xx} + f_{yy} = 0\) the obvious solution is for \( f \) to be linear in its principal directions \( x \) and \( y \). Indeed if \( f \) is linear in \(x\) and \(y\) directions then \( f_{xx} = 0 \) and \( f_{yy} = 0 \). Figure 2 is the illustration of this case. We can see that the variations of \( f \) from bottom-up is linear (visualized in three dimensions you would see a tilted plane). You can deduce from this observation that the functions of the form:

$$ f(x,y) = ax +by $$

Are all harmonic functions. Indeed every plane have a null mean curvature.

But this doesn't mean this is the only type of harmonic functions in 2D. The definition is \(f_{xx} + f_{yy} = 0\) not \( f_{xx} = 0 \) and \( f_{yy} = 0 \). So the other possibility is that the second derivatives \( f_{xx} \) and \( f_{yy} \) cancel each other out without individually being equal to zero. Take an analytical example such as \( f(x,y) = x^2 - y^2 \). Second derivatives are \( f_{xx} = 2 \) and \( f_{yy} = 2 \) and cancel out. This is a saddle function of the form:

$$ f(x,y) = ax^2 - by^2 $$ IMAGE Saddle

The saddle has two concavity with opposite orientation or opposite curvatures which cancels out. This is not the most practical example but it is good illustration for what follows. When \( f \) is harmonic and not linear then it has to have this saddle shape locally everywhere. This is the only way second derivatives can cancel each other out.

This bring us to the last example Figure 3 where it isn't so clear that \( \Delta f = 0 \) everywhere (at least to me the first time I saw it). Again here an analytical example can help. The function \( f(x,y) = \ln{ \sqrt{x^2 + y^2} } \) look a lot like our example:


If you trust me on the computation of the second derivatives they are as follows:

$$ f_{xx}(x, y) = \frac{-2x^2}{ (x^2 + y^2)^2 } + \frac{1}{x^2 + y^2} \text{ and } f_{yy}(x, y) = \frac{-2y^2}{ (x^2 + y^2)^2 } + \frac{1}{x^2 + y^2}$$

and do cancel each other out:

$$ \begin{eqnarray*} f_{xx} + f_{yy} & = & \frac{-2x^2}{ (x^2 + y^2)^2} + \frac{1}{x^2 + y^2} + \frac{-2y^2}{ (x^2 + y^2)^2} + \frac{1}{x^2 + y^2} \\ & = & \frac{-2(x^2 + y^2)}{ (x^2 + y^2)^2 } + \frac{2}{x^2 + y^2} \\ & = & \frac{-2}{ x^2 + y^2} + \frac{2}{x^2 + y^2} \\ & = & 0 \end{eqnarray*} $$


TODO: unfortunately this tutorial is still work in progress. I'll leave the outline here as a reminder for me and for your reference, who knows it my help you.

The maximum principle

The maximum principle states that a non-constant harmonic function cannot attain a maximum (or minimum) at an interior point of its domain. This result implies that the values of a harmonic function in a bounded domain are bounded by its maximum and minimum values on the boundary. In other words: "no interior local extrema in f"

Uniqueness Principle:

If we know the values of the function along the boundary of a region, there is only one harmonic function defined in that region with those boundary values.

mean-value property

"it streches"

The value at any point \( x_k \) is equal to the average value over any sphere around that point.

On a 2D domain (hot plate examples the average temp on the border of a circle is the same at its center):

\[ u(x_k) = \frac{1}{2} \int_0^{2\pi} { u(x_k + r e^{i \theta}) d \theta } \]

Hot plate illustration:

Boundary continuity


(todo: image C0 function, mesh deformation)

Dirichlet energy

The Dirichlet energy is a measure of how variable a function is. It's defined as: 
\[ E[u] = \frac 1 2 \int_\Omega \| \nabla u(x) \|^2 \, dx \]

- \( \Omega \in \mathbb R^n \) an open set
- \(u : \Omega \rightarrow \mathbb R\) a scalar function

An harmonic function (i.e. that respect the Laplace equation \( \Delta u = 0 \)) minimizes the Dirichlet energy!


Solving the Laplace equation


Setup the linear system

Solving the system

Diffuse values with Gauss-Seidel

Multi grid approach

LU decomposition

Code to solve

LLT decompostion

Code to solve


No comments

(optional field, I won't disclose or spam but it's necessary to notify you if I respond to your comment)
All html tags except <b> and <i> will be removed from your comment. You can make links by just typing the url or mail-address.
Anti-spam question: