1. Introduction

In recent years we have seen many applications of fractal sets in modeling sciences. Especially, to study several processes that can be modeled using fractals and specially self-similar structures. We mention processes related to diffusion on fractal sets witch have several possible applications, including diffusion of substances in biological structures and flow inside fractures in modeling fluid flow in fractured porous media, among other models. See [6, 9]. In this paper, we consider fractals defined as self-similar sets with the additional properties of being post-critically finite; see [6]. These self-similar sets can be approximated (in the Hausdorff metric) by a finite union of sets generated by removing a finite number of vertices from a graph approximation of the fractal set.

We recall the definition of the Laplace operator, a standard model for fractal diffusion. Then, introduce numerical approximation procedures of the presented model. It is im-portant to stress that our approximation procedure consists of renormalizing standard approximation methods on two- and three-dimensions such as the finite difference and the finite element method. See for example [3] where a finite element method is de-signed and analyzed. Note that in [3], the authors assume the renormalization constant is known. We do not need to compute the renormalization constant analytically and in-stead of of that, we approximate it numerically. We also mention that the finite element formulation implemented here is based on weak forms computed using standard length and area measures restricted to approximations of the fractal sets. In particular, in the case of the Sierpinski triangle, the approximation by the finite element method can be summarized as follows:

In the last step, the finite element procedure above is then renormalized with a pre-computation of the renormalization constant to obtain the correct approximation of the Laplace operator in the Sierpinski triangle. The computation of the renormalization constant involves the comparison of two (not-normalized) solutions at consecutive (levels of refinement) approximations of the fractal set. We call this procedure the renormalized FEM, or rFEM for short. We illustrate the performance of the rFEM numerically for the solution of a Dirichlet problem on the Sierpinski triangle and a realization of the Hata tree. We also present a renormalized Finite Difference (rFD) procedure using the same idea.

The rest of the paper is organized as follows. In Section 2 we recall some examples of self-similar sets. Section 3 is dedicated to reviewing the definition of the graph laplacian operator.

In Section 4 we present several formulations of the Dirichlet problem as well as the proposed procedure for the computation of the renormalization constant. We finish this section with numerical experiments and illustrations of the correctness of our method. In Section 5 we present some conclusions.

2. Examples of self-similar sets

This section reviews some facts related to self-similar sets, its constructions and also related to the approximation of fractal sets. We follow [⁵, ⁶, ⁹].

A self similar set is obtained by applying a fixed point functional iteration. Let (X, d) be a Hausdorff metric space and denote by C(X) the metric space of all compact subsets of a metric space equipped in the Hausdorff metric. Assume you have_Scontractions f _i : X → X, i = 1, 2, 3, ..., N and define F: C(X) → C(X) by F (A) = _1≤i≤N f _i (A) for all A ∈ C(X). Then F have a unique fixed point K. Also for any A ∈ C(X), F ⁿ (A) converge to K when n → ∞ with respect to Hausdorff metric. We then have that there exists a unique non-empty compact K ⊂ X such that

The set K is the self-similar set associated with {f ₁ , f ₂ , . . . , f _N } and this is a fractal set. The term self-similar is given to K because K is the union of images of itself by the contractions. See [⁶] and reference therein.

In order to obtain particular examples, we need only the initial metric and the finite set of contractions. In particular, we can consider the following examples of subsets of R ^d , d = 1, 2, 3 with the Euclidean distance together with a finite family of contractions. See [⁵, ⁶, ⁹].

Here r(θ) is the rotation matrix with angle θ. We define the set {W _n : n ∈ ℕ ₀ } ⊆ ℝ² inductively by,

We can define K = lim _n→∞ W _n with the limit in the Hausdorff metric. The set K is called the Kosh curve. In particular, W _n approximates the Kosh curve where we construct discrete differential operators that are approximate to differential operators defined on the Kosh curve K.

Figure 1 The first four iterations of the construction off the Kosh curve: W ₁ (Up right), W ₂ (Up left), W ₃ (down right), W ₄ (down left). 

● The Sierpinski triangle: The Sierpinski triangle is one of the most known examples of self-similar, see [⁴]. It can be consider a benchmark fractal where several ques-tions and problems can be test out. A construction of the triangle goes as follows. Let α ₀ = (0, 0), α ₁ = (1, 0), α ₂ = ( ) the vertices of the equilateral triangle ⊆ ℝ². We consider the set X with the Euclidean distance. For each i = 0, 1, 2 we define the affine mapping

Let W ₀ = [α ₀ , α ₁] ∪ [α ₁ , α ₂] ∪ [α ₂ , α ₀] and we define {W _n : n ∈ ℕ₀ } by

As before, we see that K = W _∞ = lim _n→∞ W _n (where the limit is taken in the Hausdorff metric), and W _n can be viewed as an approximation of K where differ-ential operators can be computed to approximation differential operators defined on K.

The Hata tree in the plane: W ₀ = [p ₁ , p ₂]. Define

Figure 2 The first four iterations of the Sierpinski triangle: W ₁ (up right), W ₂ (up left), W ₃ (up right), W ₄ (up left). 

Where r(θ) denotes the θ-rotation matrix. Define

The first four iterations are illustrated in Figure 3.

Figure 3 The first four iteration of the Hata tree W ₁ (up right), W ₂ (up left), W ₃ (down right) and W ₄ (down left). 

Where and is and orthonormal set.

Figure 4 The firstfour iteration of the Hata tree. See W ₁ (Up right), W ₂ (Up left), W ₃ (Down right) and W ₄ (Down Left). 

The previous fractals can be built from the set of initial vertices V ₀ and use the sequence V _n+1 = f ₁(V _n ) ∪ ... ∪ f _N (V _n ). We will use this notation later on. For example, for the Sierpinski triangle the initial vertices are V0 = {α1, α 2, α 3} and we will have the sequence V _n+1 = f ₁(V _n ) ∪ ... ∪ f ₃(V _n ), and Sierpinski triangle is V _∞ =∪^∞ _i=1 f _i (V _n ). From now on, the development of the theory will focus on the Sierpinski triangle.

3. Laplacian on a graph

In this section, we review the construction of the Laplace operator and the energy on a graph; in particular, we introduce the renormalization constant, which is important to get a finite limit of the energies associated with a family of graphs that approximates a fractal. See [⁵, ⁶, ⁹].

Let G(V, W ) a finite graph, where V = {v ₁ , v ₂ , ..., v _n } determines the set of vertices and W the set of edge (without orientation) of V . If v, w ∈ V and exist an edge between v and w we write v ∼ w ∈ W. define the adjacency matrix A _G associate to a graph G as the n × n matrix A _G = [α _ij ] _i,j ⁿ ₌₁, where

The weight matrix P _G of G is the diagonal matrix of dimension n × n defined by P _G = [p _ij ] ⁿ _i,j=1 with p _ij = 0 when i # j and p _ii is the number of adjacent vertices to v _i , i = 1, . . . , n. Therefore,

The Laplacian matrix associate to G(V, W ) is given by

If u: V → ℝ we define the energy of u by

The bilinear form associated to the energy is

We denote ε(u) = ε (u, u). If we introduce the vectors ⃗ given by ⃗ and ⃗ , then

We observe that matrix ∆ _G is the matrix representation of the energy ε _G . Now we consider the approximations to a fractal set K. We denote the energy associated to V _n by

The renormalized energy of level n = 2, 3, . . . can be computed for u _n : V _n → ℝ as

Here r ⁿ is a renormalization constant needed in order to obtain a non-increasing sequence of renormalized energies ε _n , n = 0, 1, . . . . This step is necessary to obtain well-defined energy defined on K that we introduce as a limit of the renormalized energies above. For more details see [⁶] and references therein.

3.1. The case of the Sierpinski triangle K

For the Sierpinski triangle, see [³], we obtain that the energy for f: defined for n ∈ ℕ ₀ by

where x ∼ _n y is the edge that connects to x with y and contained in W _n . That is, the renormalization constant is r = (5/3). Then, we can compute ε(u, u) := lim _n ε _n (u| _Vn , u| _Vn ). Introduce the renormalized Laplace operators by

If we define the measure µ ⁿ on V _n by assigning full measure 1 to V _n and stating that each cell has the same measure (3 ⁻ⁿ ), we see that we have

Using this it is defined the Laplace operator in K by

Here µ is the standard self-similar measure associated to K that can computed as the limit of the measure ν ⁿ in the sense that

where we note that,

From here on, we will focus our study on the Sierpinski triangle.

4. Formulations for the Dirichlet problem

Given g: V _∞ → ℝ and u ₀: V ₀ → ℝ, we seek for u: V _∞ → ℝ such that

where ∆u(x) is defined in (9) and u ₀ is the values of u in x ∈ V ₀. We call this the strong formulation of the Dirichlet problem. Using the integration by parts formula we have

Therefore, we can write the problem as seeking for u with bounded energy, u ∈ H ¹, such that

Here, H ¹ = {u: V _∞ → ℝ, ε(v, v) < +∞} and H ₀ ¹ = {v ∈ H ¹: v(x) = 0, x ∈ V ₀ }. We refer to these formulation as the weak form of the Dirichlet problem.

4.2. Finite difference approximation

To approximate the solution of (11) in V _n we consider ∆ _n , the renormalized Laplace operator (defined for the Sierpinski triangle in (8)). Let the approximation be defined by u ^{F D} : V _n → ℝ that can be written as u ^FD = {u ^F _n ^D (x)} _x∈Vn and we partitioned it as follows

Note that u ₀ is know and corresponds to the boundary values. Analogously for b ^FD = {g(x)} _x∈Vn = [b ₀ , b _I ]. We obtain the block structure

We compute u _I as the solution of

The main issue with this discrete formulation is that we need to know the renormalization constant: the value 5 ⁿ in the case of the Sierpinski triangle, see (8). The renormalization constant can be viewed as a scaling of the forcing term in the linear system. This scaling can be approximated as explained next.

4.3. Computation of the renormalization constant

First we consider the finite difference method renormalization constant. The idea is to use the problem,

Denote by q ⁿ and approximation of the renormalization factor that we want to compute. The space of study is V _n /V ₀. We can approximate (13) by

We want to compute the value q ⁿ . Numerically we can compute the solution of the problem, see [7],

We must have that z _n (x) = q ⁿ u _n (x), x ∈ V _n since is nonsingular. For Vn+1 we will have z _n+1 = q ⁿ⁺¹ u _n+1 . Note that z _n and z _n+1 can be calculated without knowing q ⁿ . For n large enough we should have,

where u is the exact solution of (13). Therefore we should be able to use the approximation,

where x ∈ V _n \V ₀. See some numerical illustration in Table 1 for the case of the Sierpinski triangle. In this case we see hat lim _n→∞ q _n,n+1 = 5.

Table 1 Values of q _n,n+1 for the Sierpinski triangle, given for (16). 

A similar procedure can be implemented for the energy renormalization constants. The weak form of (13) is written as,

Compute the solution of the graph energy problem (without renormalization),

We can then define the approximation,

Here x ∈ V _n \ V ₀. See a numerical verification in Table 2. In this case we should have .

Table 2 Values of r _n,n+1 for the Sierpinski triangle, given by (17). 

4.4. Renormalized finite elements methods- rFEMs

Now we construct approximations for the weak form (12). That is, we propose approxi-mations for the computations of renormalized energy bilinear forms by rescaling standard approximations as introduced in Section 4.3.

We defined P ⁿ = {u: V _n → ℝ}. We then project the weak formulation into the space P ⁿ . In the Galerkin formulation we seek to find u ∈ P ⁿ such that

Both u and v are defined on V _n , then this formulation is equivalent to

We propose several classical finite element methods procedures in one and two dimensions to approximate the renormalized energy. Note that we can use the energy ε _n rather than the limit ε. The renormalization constant is then computed by a procedure similar to the one explained previously. We recall that in this section we consider the case of the Sierpinski triangle.

Integrals along the edges

Recall that . Introduce the bilinear form,

where u ^′ denotes the one dimensional derivative along [x, y] of the linear interpolation of the vertex values u(x) and u(y). The measure ν _n ¹ is defined as the length measure along the edges of W _n (rescaled to obtain total length 1). Note that for each n there are total 3 ⁿ⁺¹ edges (3 for each cell). Each edge of V _n has length 2 ⁻ⁿ . Given a total length of 3(3/2) ⁿ before length rescaling.

We consider the following discrete problem,

As before we can write u = u _I + u _G where u _G (x) = g(x), x ∈ V ₀; u _G (x) = 0, x ∈ V _n \ V ₀. Analogously u _I ∈ P ⁿ tal such that u _I (x) = 0 for x ∈ V ₀. We then have,

This is equivalent to the linear system,

Let V _n \ V ₀ = {x ₁ , x ₂ , ..., x _p } the set of interior vertices and we have,

where φ _xi is the linear interpolation of the characteristic functions of {x _i } in V _n . We also have

It is easy to see that

Recall that, on the left hand side above we use the piecewise liner interpolation of the nodal value of g. The previous formulation have to be renomarlized to

We note the following relation between the renormalized measure and the measure in-duced by the length measure, for any n, we have

that follows by computing the length integrals using the trapezoidal rule.

The renormalization constant r¹ can be approximated using the procedure described in Section 4.3 by solving consecutive refinement level approximations with the constant function 1 as the right hand side and Dirichlet boundary conditions. See (17). In Table 3 we show the results of computing the renormalization constant.

In this case the renormalization constant can also be computed analytically. Integration by parts and the fact that we use piecewise linear interpolation (u ^′′ = 0 inside edges) yields

Having into account that the length of the edges of the n approximation of K is 1/2 ⁿ , we get,

Therefor

Due to (23) and (7) we see that renormalization constant for the family . In Table 3 we show the results of computing the renormalization constant using the procedure explained in Section 4.3. This a numerical verification that for the case of the Sierpinski triangle the computation agrees with the exact value of the renormalization constant just derived.

Table 3 Values of r _n,n+1 for the Sierpiski triangle using integrals along Edges. 

Area integrals

Introduce the bilinear form,

where ∇u denotes the two-dimensional gradient of the two-dimensional linear interpo-lation of the nodal value of u in the triangle τ. The measure ν _n ² is the area measure restricted to Kn and normalized such that the total area of (all the triangles of) Vn is one. Note that, for each n, there is 3 ⁿ each of them of área for a total area of before rescaling.

We formulate the following discrete problem,

This time the previous formulation is equivalent to the linear system,

Where

and

As before, a renormalization is needed, that is,

The renormalization constant r ⁽² can be approximated using the procedure described in Section 4.3 by solving consecutive refinement level approximations with the constant function 1 as the right hand side and Dirichlet boundary conditions. See (17). In Table 4 we show the results of computing the renormalization constant.

In order to verify our computations we compute the renormalization constant analytically. This is possible in this case. Recall that,

We use standard finite element analysis. Introduce the reference basis functions

defined in the reference triangle τ^ with vertices (0, 0), (1, 0) and (0, 1). This reference triangle can be mapped into the triangles of K _n by an affine mapping in two dimensions. If we consider the triangle τ of K _n with vertices and ⃗ , this mapping is given by F _τ : τ^→ τ defined by

Define . Any linear function on τ is a linear combination of the basis functions P ₁ , P ₂ , P 3, in particular if u is a linear function on τ we have u(ψ) = u(x)P ₁(ψ)+u(y)P ₂(ψ)+u(z)P ₃(ψ). From the definition of Pi is easy to see that

Moreover, we also have,

We also recall that det where l = 2^-n then diameter of the triangle. We can compute

Table 4 Values of r _n,n+1 for the Sierpiski triangle using area measures. 

Analogously we have,

Having into account that each interior node belongs only to two-triangles and that two distinct nodes share an edge in at most one triangle, we conclude that the assembled global matrix is given by

where i and j corresponds to the index of interior nodes. We see that,

Recall that ∆~ _n was defined as the graph laplacian of the n approximation of the Sierpinski triangle.

We note the following relation between the renormalized measure and the measure in-duced by the area measure, for any n, we have

that follows by computing the length integrals using the trapezoidal rule.

From (28) and (29), we then have that the exact value of the renormalization constant is the same as the one for the construction based on length measures. This result verifies the numerical computations obtained in Table 4.

The Sierpinski triangle

For problems posed on the Sierpinski triangle recall that V ₀ = {α ₀ , α ₁ , α ₂ }. We want to approximate the solution of

In Figure 5 we illustrate some results.

Figure 5 Approximated solution with rFD (left), rFEM1D (center) and rFEM2D (right). Here g(x, y) = sin(x + y), u(α ₀) = 1, u(α ₁) = u(α ₂) = 0 and we consider V ₅ the fifth level approxima-tion of the Sierpinski triangle. 

The Kosh curve

For problems posed on the Kosh curve, recall that V ₀ = {(0, 0), (1, 0)}. We want to approximate the solution of

Table 5 Values of r _n,n+1 that the Kosh curve. 

We use the method introduced in section 4.4; if we associate the Laplacian with the energy defined by integrations over the edges, we have

where µ _n is the length measure restricted to edges and was computed in Table (6) using the procedure described in Section 4.3.

Figure 6 Approximated solution with rFD (left) and rFEM1D (right). Here g(x, y) = sin(x+y), u(α ₀) = 1, u(α ₁) = 0 and we consider V ₅ the fifth level approximation of the Kosh curve. 

Table 6 Values of r _n,n+1 for the Hata tree. 

Figure 7 Approximated solution with rFD (left) and rFEM1D (right). Here g(x, y) = 0, u(α ₀) = 1, u(α ₁) = 0 and we consider V ₃ the third level approximation of the Kosh curve. 

The Hata tree

In the find method, the solution of the following equation is the same as that studied before, where V ₀ = {(0, 0), (1, 0)}.

if we associate the Laplacian with the energy on edge, we have

Where µ _n is the self-similar measure that the Kosh curve and the renormalization con-stant is given for , see Table 6.

Figure 8 Approximated solution with rFD (left) and rFEM1D (right). Here g(x, y) = sin(x+y), u(α ₀) = 1, u(α ₁) = 0 and we consider V ₃ the third level approximation of the Hata tree. 

Figure 9 Approximated solution with rFD (left) and rFEM1D (right). Here g(x, y) = 0, u(α ₀) = 1, u(α ₁) = 0 and we consider V ₃ the third level approximation of the Hata tree. 

Hata tree in the space

For this fractal we use only rFD method (4.2), In Figure 10 we present a computed solution.

Figure 10 Approximated solution with rFD. Here g(x, y) = 0, u(α ₀) = 1, u(α ₁) = 0 and we consider V ₃ the third level approximation of the Hata tree in the space. 

5. Conclusions

This paper designed a numerical procedure to approximate solutions to diffusion problems on self-similar fractal sets. We start with a discrete approximation of the fractal and the derivatives in standard non-renormalized formulations. We can then precompute the renormalization constant needed to approximate the actual differential operators on the fractal set. In particular, we present examples with the Sierpinski triangle using standard graph weights and adjacency matrices (Finite Difference method) or using week forms with length or area measures (Finite Element method). In the Finite Element method with the length measure, the derivatives in the weak forms are classical derivatives along the edges with integration concerning the length measure. In the finite element method with the area measure, we use partial derivatives with integration in two dimensions on triangles of the approximation of the Sierpinski triangle. We also present additional illustrations with the Kosh curve and the Hata tree.

It is also important to mention that the implementation of proposed finite element methods are simple and do not have significant changes with respect to the finite element method for differential equations in open domains. Also, the renormalization constant does not need to be known a priori. We can use finite element codes that work on tri-angulations in general, and only the “ triangulation” or graph that approximates the fractal must be used as input for these codes. The renormalization constant can be pre-computed as proposed in this paper. We observe that diffusion processes on these fractal sets can be approximated by classical diffusion processes (involving classical derivatives) on fractal approximations. These must be rescaled by the scale parameter that can be precomputed. The authors will explore this idea and the related numerical analysis in future works.