Home - Lesson 1 - Lesson 2 - Lesson 3 - Lesson 4 - References.

Lesson 2 - Numeric Algorithms

Brian Bockelman

Abstract:

This lesson starts off with two simple numeric application - explicitly constructing all the harmonic functions on the $ SG$ , and doing numeric integrations. We then discuss the theory behind how a numeric solver may be implemented, and present one that actually has been made. The reader is encouraged to download the code, and produce the results in this lesson on their own. Finally, in order to demonstrate how the solver works and the general correctedness of the solver, we discuss the eigenvalues and eigenfunctions of the Laplacian on the $ SG$ .

Harmonic Functions, Revisited

I ended the section on harmonic functions from the last lesson rather cryptically, saying there exists a explicit construction of any harmonic. This served partly as a teaser, and partly because the algorithm is a great way to introduce the topic of this lesson: computations on the Sierpinski Gasket. Specifically, by the end of this lesson, you should be able to solve dynamic equations on the Sierpinski Gasket (or at least understand the steps), one of the main goals of this tutorial.

One obvious way to create harmonic function on the $ SG$ is to create a function whose Laplacian is identically zero on each approximation graph, $ \Gamma_k$ . That way, when we take the limit $ k \rightarrow \infty$ , we would get the limit of the Laplacian of the function is zero.

Start off by picking the value of the function at the three boundary points, $ v_0$ , $ v_1$ , and $ v_2$ . Then, there are three points in $ V_1$ whose value is unknown; $ u_0$ , $ u_1$ , and $ u_2$ as labeled below:

We want the Laplacian at $ u_i$ to be zero - this leads to the three equations:

$\displaystyle \triangle_1 f(u_i) =\frac32 5^1 \sum_{y \in N(y)} (f(y) - f(u_i)) = 0, i=0,1,2 $

Expanding these three equations (with three unknowns), we get:

$\displaystyle f(v_0) + f(v_1) + f(u_0) + f(u_1) - 4 f(u_2) = 0 $

$\displaystyle f(v_1) + f(v_2) + f(u_1) + f(u_2) - 4 f(u_0) = 0 $

$\displaystyle f(v_2) + f(v_0) + f(u_2) + f(u_0) - 4 f(u_1) = 0 $

I leave it up to you to solve for $ f(u_i)$ ; the answers are:

$\displaystyle f(u_0) = \frac15 f(v_0) + \frac25 f(v_1) + \frac25 f(v_2) $

$\displaystyle f(u_1) = \frac15 f(v_1) + \frac25 f(v_2) + \frac25 f(v_0) $

$\displaystyle f(u_2) = \frac15 f(v_2) + \frac25 f(v_0) + \frac25 f(v_1) $

So, the three inner vertices of level $ 1$ are totally determined by the vertices of level 0. If we would now split level $ 1$ up into the three cells, and do the calculations for level 2, we would get similar formulas.
So, the values of a harmonic function on level $ m+1$ is totally determined by level $ m$ by what we call the $ \frac15-\frac25-\frac25$ rule, illustrated in the following diagram:

\includegraphics{152525rule.eps}

Figure illustrating the 1/5-2/5-2/5 rule. Notice how u_1 is 2/5 of its closest neighbors, and 1/5 of the far one.

Because the values in level $ m+1$ are completely determined by those of $ m$ for harmonic functions, and $ V_0$ contains three points, the set of all harmonic functions forms a vectors space of dimension 3. One convenient basis is the functions $ h_i$ , for $ i=0,1,2$ , defined by

$\displaystyle h_i(v_j) = \delta_{ij}, $

where $ \delta_{ij}$ is the delta function; it is $ 1$ when $ i = j$ and 0 otherwise. That is, $ h_i$ is $ 1$ on the boundary point $ v_i$ and 0 otherwise. So, if $ f$ is a harmonic function, $ f(v_0) = a$ , $ f(v_1) = b$ , and $ f(v_2) = c$ , then we can write $ f$ as

$\displaystyle f(x) = a h_0(x) + b h_1(x) + c h_2(x) $

A picture of one of the basis elements follows:

An element of the level 2 spline basis.

Notice that the other two basis functions are just rotations of this one.

Generalized Harmonic Functions and Splines

Harmonic functions were our first step to the $ SG$ 's equivalent for piecewise linear functions. A piecewise harmonic function is a continuous function which is harmonic on each of its cells. Suppose that $ w_1, \ldots, w_n$ are $ n$ words such that

$\displaystyle SG = \bigcup_{i=1}^n F_{w_i}(SG). $

That is, $ \delta_i := F_{w_i}(SG)$ are all cells, and Sierpinski Gasket is decomposed into these cells. Let $ \delta_i^0 := F_w(V_0)$ . Then, for a function $ f$ to be piecewise harmonic, we need for $ \triangle f\vert _{\delta_i \setminus \delta_i^0}$ to be zero. We also impose a matching condition:

    If $ F_{w_i}(v_j) = F_{w_k}(v_\ell)$ , then $ f(F_{w_i}(v_j)) = f(F_{w_k}(v_\ell))$ $\displaystyle $

This condition simply requires that the different 'pieces' match up in a continuous fashion.

The space of all piecewise harmonic functions which are zero on the boundary points is called the spline space. The function must be zero on the boundary because the spline spaces will be used to solve equations with Dirichlet boundary values, $ f(V_0) \equiv 0$ . If we further require the cells to be all in level $ k$ , then the spline space is labeled $ \mathcal{S}_k$ . Each $ \mathcal{S}_k$ has a very natural basis. Label the points in $ V_k \setminus V_0$ as $ x_1, \ldots, x_n$ . Then, $ \phi_i$ is the piecewise harmonic function that is zero at all of the $ x_j$ for $ i \neq j$ and $ 1$ at $ x_i$ . Once we have defined $ \phi$ at these points, we have defined $ \phi$ precisely at each of the points in the boundaries of each cell of level $ k$ . Thus, extend the function into each cell using the harmonic algorithm discussed above, and $ \phi_i \in \mathcal{S}_k$ . I leave it up to you to see that $ \{ \phi_i \}_{i=1}^n$ is a basis for the vector space $ \mathcal{S}_k$ . Here's a picture of one of these spline basis functions for $ k=2$ :

We know that, in the case of the real interval, we can gain a lot of accuracy by approximating with quadratics instead of approximating with lines.

There are similar functions in the case of $ SG$ , the biharmonic functions. In general, a function is n-harmonic if

$\displaystyle \triangle^n f \equiv 0 $

That is, the Laplacian of $ f$ , applied $ n$ times, is identically zero:

$\displaystyle \underbrace{\triangle ( \triangle ( \cdots \triangle}_{\text{n times}} f ) \cdots ) = 0.$

Here's an example of a tri-harmonic function.

So, a biharmonic function is a function such that $ \triangle^2 f \equiv 0$ ; these are the $ SG$ 's equivalent of a quadratic function. The spline space of functions that are biharmonic on all the cells of level $ k$ and 0 on the boundary is called $ \mathcal{S}_k^2$ . We must add a new matching condition; the normal derivatives must match up at the boundaries of the cells in order to make the spline function smooth. This spline space has a ``nice" basis also. Let $ x_1, \ldots, x_n$ be the points of $ V_k \setminus V_0$ . Then, define

$\displaystyle f_i(x_j) = \delta_{ij} $

$\displaystyle \partial f_i(x_j) \equiv 0 $

$\displaystyle g_i(x_j) \equiv 0 $

$\displaystyle \partial g_i(x_j) = \delta_{ij}. $

Remember, $ \delta_{ij}$ is defined to be

$\displaystyle \delta_{ij} = \begin{cases}1, & i = j \\ 0, & i \neq j \end{cases}, $

and $ \partial$ is the normal derivative.
If we extend each of the $ f_i$ and $ g_i$ in a biharmonic manner, we end up with a basis for $ \mathcal{S}_k^2$ .

Integration

A second important, but simple example of computation on the Sierpinski Gasket is that of numerical integration. For those of you who are familiar with measure theory, you already know that there is a clear and ``natural" way to define integrations on the $ SG$ , but I conveniently left out any algorithms for doing so.

It is clear that, like last section, we cannot calculate the exact integral of an arbitrary function on the $ SG$ , just like we cannot calculate the precise integral of every algebraic function. To continue our analogy with the real line, closed forms are known for certain algebraic functions (such as all polynomials), but a numeric approximation must be made for the more complex ones.

Remember our approximation of the real Laplacian on the unit interval by using the didactic points? We will construct our integral on the $ SG$ in a similar way. Let's say we would like to do approximate a Riemann integral on the level 1 of the unit interval. So, we break up the interval into $ [0,\frac12]$ and $ [\frac12,1]$ , and do the Riemann sum of this partition:

$\displaystyle \int f(x) dx \approx f(0) (\frac12 - 0) + f(\frac12)(1 - \frac12)... ...rac12 f(0) + \frac12 f(\frac12) = \frac12 \sum_{x \in D_1 \setminus\{1\}} f(x) $

What we are really doing is taking each interval, multiplying its length / size by the function value, and summing over each interval. By looking at more intervals, we can get a better approximation of the integral. So, for $ k = 3$ , $ D_k$ has 9 points, and 8 intervals. Each interval has length $ \frac18$ , so

$\displaystyle \int f(x) dx \approx \frac18 \sum_{x \in D_1 \setminus \{1\}} $

If we want a really good approximation of the integral, we might decide to pick $ k = 6,7$ , so we sample the function at $ 2^k-1$ points.

We want to perform numerical integrations on the $ SG$ in the exact same way. One stumbling block remains - we naturally knew the ``size" of each piece of the unit interval. What should be the ``size" of one of the cells of $ SG$ ?

To start out with, we want to think of the normal $ SG$ as a single unit in size - so, it would make sense that the integral of the function $ f: SG \rightarrow \mathbb{R}$ defined by $ f(x) = 1$ to be one. That is, we would like

$\displaystyle \int_{SG} 1 dx = 1 $

In our approximation, when we compare the area of the triangle in $ S_0$ to the area of one of the three triangles in $ S_1$ , we find that the area decreases by $ \frac13$ . So, we want the area of each triangle to decrease by $ \frac13$ as we move down each level. For a cell in level $ m$ , we set the size of the cell in the approximation to be $ 3^{-m}$ .

An example of a numeric integration of a function

For an example of an exact integration, consider the three harmonic basis functions, $ h_i$ , $ i=0,1,2$ . Notice that the sum $ h_0 + h_1 + h_2$ is identically one. Further, the integral of each $ h_i$ should be exactly the same, since the functions are identical except for ordering of the vertices. Thus,

$\displaystyle \int_{SG} h_0 dx = \int_{SG} h_1 dx = \int_{SG} h_2 dx $

$\displaystyle \int_{SG} h_0 dx + \int_{SG} h_1 dx + \int_{SG} h_2 dx = \int_{SG} (h_0 + h_1 + h_2) dx = \int_{SG} 1 dx = 1 $

Hence,

$\displaystyle \int_{SG} h_i dx = \frac13, \, i = 0,1,2 $

Now, for any harmonic function, $ f$ ,

$\displaystyle f(x) = f(v_0) h_0(x) + f(v_1) h_1(x) + f(v_2) h_2(x). $

Thus,

$\displaystyle \int_{SG} f(x) = \int_{SG} ( f(v_0) h_0(x) + f(v_1) h_1(x) + f(v_2) h_2(x) ) dx $

$\displaystyle = f(v_0) \int_{SG} h_0(x) dx + f(v_1) \int_{SG} h_1(x) dx + f(v_2) \int_{SG} h_2(x) dx $

$\displaystyle \int_{SG} f(x) = \frac{f(v_0)}{3} + \frac{f(v_1)}{3} + \frac{f(v_2)}{3}$ (1)

Trapezoidal Rule and Simpson's Method

Just like with the real line, the naive implementation of numerical integration isn't the best. There are two other methods that, at very little computation cost, can increase the accuracy of your integrations.

First, there is the so-called trapezoidal rule. Instead of evaluating at one point per cell, we can evaluate at all three boundary points of the cell, and average them out. This method is rather straightforward, and illustrated below:

An example of how the trapezoidal rule is used on a function

Secondly, there is Simpson's Rule. Here, we use the next level, $ k+1$ to help us estimate the value to use for a cell in level $ k$ . Simpson's Rule is illustrated below:

An example of how the Simpson's rule is used on the same function as above. Compare the two results.

Energy

We have one last theoretical hurdle before we can describe our methods of solution. We must construct a way to measure the energy of a function on $ SG$ . The energy of a function on the real line is

$\displaystyle \mathcal{E}(f,f) = \int (f')^2 dx $

It will be a bit more difficult to derive an exact analogy for the $ SG$ just from this equation, as we have no analogue for the derivative. Instead, think of the energy form as somehow measuring the amount of change in the function $ f$ . Hence, in our construction, we want the energy of a constant function to be zero, just like in the case of energy on the real line.
To mimic the term inside the integral, $ (f')^2$ , we use $ (f(x) - f(y))^2$ where $ x,y$ are neighbors in level $ k$ . In order to simplify notation, we say that $ x \sim_k y$ if $ x$ and $ y$ are neighbors on level $ k$ ; another way to look at this is to say $ x \sim_k y$ if and only if $ x$ and $ y$ are two endpoints of a edge in the graph $ \Gamma_k$ . Hence, we want our energy to have the form

$\displaystyle E_k(f,f) = \sum_{x \sim_k y} c_k (f(x) - f(y))^2$ (2)

where the $ c_k$ is some constant. By taking the sum over $ x \sim_k y$ , we mean sum over each neighbor-pair once, or we can think of the sum as taken over all the edges of $ \Gamma_k$ .
We now impose conditions so the $ c_k$ can be determined. We want a renormalization identity to hold for all $ k$ :

$\displaystyle E_k(f,f) = \sum_{i=0}^2 r_k E_{k-1}(f \circ F_i, f \circ F_i)$ (3)

We also want the following:

$\displaystyle E_1(\bar f,\bar f) = E_0(f,f),$ (4)

where $ \bar f$ is the harmonic extension of $ f$ from $ \Gamma_0$ to $ \Gamma_1$ . Equations 2 - 4 set up a rather complicated system to solve for the constants. However, if one draws a picture and works through all the conditions, we find that these conditions are consistent, and lead to the following equation for the energy:

$\displaystyle \mathcal{E}_k(f,f) = \left(\frac53\right)^k \sum_{x ~_k y} (f(x) - f(y))^2 $

and

$\displaystyle \mathcal{E}(f,f) = \lim_{k \rightarrow \infty} \mathcal{E}_k(f,f) $

In fact, we can extend the definition a little more to be able to examine the energy of two functions:

$\displaystyle \mathcal{E}_k(f,g) = \left( \frac53 \right)^k \sum_{x ~_k y} (f(x) - f(y))(g(x) - g(y)) $

$\displaystyle \mathcal{E}(f,g) = \lim_{k \rightarrow \infty} \mathcal{E}_k(f,g)$

Here's an example of how to approximate the energy of a function; we use a level 1 approximation.

We have the following identity, called the weak formulation of the Laplacian, which relates the energy, integral, and Laplacian of a function, $ u$ :

$\displaystyle \mathcal{E}(u,v) = -\int_{SG} v(x) \triangle u(x) dx$ (5)

for all $ v$ which vanish on the boundary ( $ v(V_0) \equiv 0$ ) and have finite energy. In fact, in more rigorous discussions of the $ SG$ , the Laplacian is defined to be the function $ \triangle u(x)$ such the above formula holds true for all $ v$ . It turns out that this formula is the key for the finite element method, as it allows us to relate the energy to the Laplacian.

The Finite Element Method

Overall, methods for solving differential equations fall into one of two categories: finite difference methods and finite element methods. It turns out that we have all the background needed to do the second, a finite element method.

We want to solve the equation

$\displaystyle - \triangle u + qu = f$ (6)

with what is called Dirichlet initial conditions, that is,

$\displaystyle u\vert _{V_0} = 0. $

It may not be easy to solve this problem analytically, but we can do an approximation. Suppose $ u$ solves the above problem. Then,

$\displaystyle \tilde u = \sum_{i=1}^{j_k} c_i \phi_i $

is the approximation of $ u$ in the spline space $ \mathcal{S}_k$ , where $ \phi_i, i = 1, \ldots, j_k$ are the basis elements of $ \mathcal{S}_k$ . We multiply equation 6 by an arbitrary element, $ v$ of $ \mathcal{S}_k$ , and integrate:

$\displaystyle - \int_{SG} \triangle \tilde uv + \int_{SG} quv = \int_{SG} fv $

$\displaystyle - \int_{SG} \sum_{k=1}^{j_k} c_k \phi_i v + \int_{SG} q \left(\sum_{k=1}^{j_k} c_k \phi_k\right) v = \int_{SG} fv $

By the weak definition of the Laplacian (equation 5), we get

$\displaystyle \sum_{i=1}^{j_k} c_k \mathcal{E}(\phi_i,v) + \sum c_k \int_{SG} q \phi_i v = \int_{SG} f v$ (7)

Now, we can take $ v = \phi_\ell$ , a basis element for $ \mathcal{S}_k$ , and sum equation 7 $ j_k$ -times, once for each $ \ell = 1, \ldots, j_k$ , resulting in:

$\displaystyle \sum_{\ell=1}^{j_k} \sum_{i = 1}^{j_k} c_i \mathcal{E}(\phi_i, \p... ...i=1}^{j_k} \int_{SG} q c_i \phi_i \phi_\ell = \sum_{\ell = 1}^{j_k} f \phi_\ell$ (8)

The equations are clumsy in this notation. To fix this, define the constants

$\displaystyle e_{i,\ell} = \mathcal{E}(\phi_i, \phi_\ell) $

$\displaystyle q_{i,\ell} = \int_{SG} q \phi_i \phi_\ell $

$\displaystyle f_\ell = \int_{SG} f \phi_\ell $

and the matrices

$\displaystyle E = [ E_{i,\ell} ] $

$\displaystyle Q = [ Q_{i,\ell} ] $

$\displaystyle F = [ f_\ell ] $

$\displaystyle c = [ c_\ell ], $

where $ i$ is row index and $ \ell$ is the column index. Then, equation 8 becomes the following matrix equation:

$\displaystyle Ec + Qc = f. $

The solution is straightforward:

$\displaystyle c = (E + Q)^{-1} f. $

Of course, while this is the solution, in the actual numeric computation, there are more accurate and quicker methods to solve for $ c$ than finding the inverse of $ (E + Q)^{-1}$ .
Now that we know the value of $ c$ , we can rebuild our approximation to the solution, $ u$ :

$\displaystyle u \approx \tilde u = \sum_{i=1}^{j_k} c_i \phi_i $

This completes the solution procedure. Notice that $ E$ is easy to compute; exact values are known for the energy of the product of two harmonic functions. These need to be computed once, and reused afterward. The matrix $ Q$ is a bit more difficult, as it requires $ j_k^2$ integrations over all of $ SG$ . However, for many problems, $ q$ is identically zero, so we don't have to worry about it. In the case that $ Q = 1$ , we get that each entry of $ Q$ is

$\displaystyle \int_{SG} \phi_i \phi_\ell, $

or simply the inner product of the two spline basis functions. For this reason, $ Q$ is sometimes called the inner product matrix. Finally, the vector $ f$ requires $ j_k$ integrations. So, if a very efficient integration method is used, and $ q(x) \equiv 0$ (or identically a constant), the finite element method can be extremely quick.

Now that you've progressed through the basics of solving dynamic equations on the Sierpinski Gasket, there should be nothing stopping you from being able to implement the methods yourself. No need to reinvent the wheel though; in the next lesson, I walk you through working with a set of tools made specially for this tutorial that will allow you to use the finite element method (FEM) on your own.

Home - Lesson 1 - Lesson 2 - Lesson 3 - Lesson 4 - References.
The copyright on all content and design of this website belongs to Brian Bockelman. The content is licensed under the Creative Commons license, and code under the GPL.