Numerical evaluation of singular integrals on non-disjoint self-similar fractal sets

We consider the numerical evaluation of a class of double integrals with respect to a pair of self-similar measures over a self-similar fractal set (the attractor of an iterated function system), with a weakly singular integrand of logarithmic or algebraic type. In a recent paper [Gibbs, Hewett and Moiola, Numer. Alg., 2023] it was shown that when the fractal set is"disjoint"in a certain sense (an example being the Cantor set), the self-similarity of the measures, combined with the homogeneity properties of the integrand, can be exploited to express the singular integral exactly in terms of regular integrals, which can be readily approximated numerically. In this paper we present a methodology for extending these results to cases where the fractal is non-disjoint but non-overlapping (in the sense that the open set condition holds). Our approach applies to many well-known examples including the Sierpinski triangle, the Vicsek fractal, the Sierpinski carpet, and the Koch snowflake.


Introduction
In this paper we consider the numerical evaluation of integrals of the form where (see §2 for details) Γ ⊂ R n is the attractor of an iterated function system (IFS) of contracting similarities satisfying the open set condition, µ and µ ′ are self-similar (also known as "invariant", or "balanced") measures on Γ, and Φ t (x, y) := |x − y| −t , t > 0, log |x − y|, t = 0, x, y ∈ R n . (2) In the case where t > 0 and µ ′ = µ the integral (1) is known in the fractal analysis literature as the "t-energy", or "generalised electrostatic energy", of the measure µ (see e.g.[9, §4.3], [20, §2.5] and [19,Defn 4]).Integrals of the form (1) also arise as the diagonal entries in the system matrix in Galerkin integral equation methods for the solution of PDE boundary value problems in domains with fractal boundaries, for instance in the scattering of acoustic waves by fractal screens [6].In such contexts the accurate numerical evaluation of these matrix entries is crucial for the practical implementation of the methods in question.
Numerical quadrature rules for (1) were presented recently in [13] for the case where Γ is disjoint (see §2 for our definition of disjointness), e.g. a Cantor set in R or a Cantor dust in R n , n ≥ 2. The approach of [13] is to decompose Γ into a finite union of self-similar subsets Γ 1 , . . ., Γ M (each similar to Γ) using the IFS structure, and to write I Γ,Γ as a sum of integrals over all possible pairs Γ m × Γ n .
In the non-disjoint case, however, distinct self-similar subsets of Γ may be non-disjoint, intersecting at discrete points (such as for the Sierpinski triangle, see §5.1) or at higher dimensional sets (such as for the Sierpinski carpet, see §5.3, or the Koch snowflake, see §5. 4).This means that some of the integrals over Γ m × Γ n , for m ̸ = n, are singular, reducing the accuracy of quadrature rules based on the representation formula of [13,Thm 4.6] (which assumes they are regular).In this paper we remedy this, showing that in many non-disjoint cases (including those mentioned above), by decomposing Γ further into smaller self-similar subsets it is possible to find a finite number of "fundamental" singular integrals (including I Γ,Γ itself) that satisfy a square system of linear equations that can be solved to express I Γ,Γ purely in terms of regular integrals that are amenable to accurate numerical evaluation.
We will describe our methodology in generality in §4, but for the benefit of the reader seeking intuition we illustrate here the basic idea for the simple case where Γ = [0, 1] 2 ⊂ R 2 is the unit square, µ = µ ′ is the Lebesgue measure on R 2 restricted to Γ, and t = 1.Despite not generally being regarded as a "fractal", the square Γ can be viewed as the self-similar attractor of an iterated function system comprising four contracting similarities, and can accordingly be split into a "level one" decomposition of 4 squares of side length 1/2, or a "level two" decomposition of 16 squares of side length 1/4, as illustrated in Figure 1.Let I m,n denote the integral Γm Γn |x − y| −1 dydx where Γ m and Γ n are any of the level one squares in Figure 1.One can then express (1) as the sum I Γ,Γ = 4 m=1 4 n=1 I m,n of 16 singular integrals over all the pairs of level one squares, which can be categorised as follows: 4 self interactions (I 1,1 , I 2,2 etc.), 8 edge interactions (I 1,2 , I 2,1 , I 2,3 etc.) and 4 vertex interactions (I 1,3 , I 3,1 , I 2,4 , I 4,2 ).By symmetry, each of the integrals in each category is equal to all the others, so that Figure 1: Level 1 (left) and level 2 (right) subsets of the unit square.Here the labels "1" and "13" indicate the subsets Γ 1 and Γ 13 etc.referred to in the text.
Furthermore, by (3) (with t = 1), combined with a change of variables, the self-interaction integral I 1,1 can be expressed in terms of the original integral, as Combining this with (4) we obtain the equation 1 2 I Γ,Γ = 8I 1,2 +4I 1,3 .To derive two further equations connecting I Γ,Γ , I 1,2 and I 1,3 we move to the level two decomposition, extending our notation in the obvious way, so that e.g.I 12,23 denotes the integral Γ 12 Γ 23 |x − y| −1 dydx, where Γ 12 and Γ 23 are the level two squares labelled "12" and "23" in Figure 1.Then the edge interaction integral I 1,2 can be written as a sum of 16 integrals over pairs of level two squares, which, after applying symmetry simplifications, gives where R 1,2 = 4I 11,21 + 4I 11,24 + 2I 11,22 + 2I 11,23 is a sum of regular integrals.Similarly, the vertex interaction integral I 1,3 can be written as where R 1,3 = 4I 13,32 + 4I and inserting these identities into ( 6) and (7) gives our two sought-after equations, namely To summarize, we have shown that the triple (I Γ,Γ , I 1,2 , I 1,3 ) T satisfies the linear system and solving the system gives which is an exact formula for I Γ,Γ in terms of the seven regular integrals I , which are all amenable to accurate numerical evaluation, for instance with a product Gauss rule.
Our goal in this paper is to derive formulas analogous to ( 9) and ( 10) for more general Γ, t, µ and µ ′ .The structure of the paper is as follows.In §2 we review some preliminaries concerning self-similar fractal sets and measures.In §3 we introduce the notion of "similarity" for integrals over pairs of subsets of Γ, and provide sufficient conditions under which it holds.In §4 we describe a general algorithm for generating linear systems of the form (9) using our notion of similarity.In §5 we apply our algorithm to a number of examples including the Sierpinski triangle, the Sierpinski carpet, and the Koch snowflake.Finally, in §7 we show how our results can be combined with numerical quadrature to compute accurate numerical approximations to the integral (1) in these and other cases.As an application, we show how our algorithm can be used in the context of the "Hausdorff boundary element method" of [6] to compute acoustic scattering by non-disjoint fractal screens.
Regarding related literature, we note that a three-dimensional version of the approach described above for integration over the unit square was used to compute the gravitational force between two cubes sharing a common face in [28].More generally, this sort of approach forms the basis of the "hierarchical quadrature" developed for singular integrals over cubical and simplicial domains by Börm and Hackbusch [5] and Meszmer [21,22].In the context of integration over fractals, we mention the work of Mantica [17] and Strichartz [25], where self-similarity techniques were used to derive exact formulas for integrals of polynomials over fractals.Our previous paper [13] and the current paper can be viewed as extensions of the results of [17] and [25] to singular integrands.

Preliminaries
Throughout the paper we assume that Γ is the attractor of an iterated function system (IFS) of contracting similarities satisfying the open set condition (OSC), meaning that (see e.g.[15]) (i) there exists M ∈ N, M ≥ 2, and a collection of maps {s 1 , s 2 , . . ., s M }, such that, for each for some ρ m ∈ (0, 1).Explicitly, for each m = 1, . . ., M we can write for some orthogonal matrix A m ∈ R n×n and some translation δ m ∈ R n ; (ii) Γ is the unique non-empty compact set such that where Then Γ has Hausdorff dimension dim H (Γ) = d, where d ∈ (0, n] is the solution of the equation We say that the IFS is disjoint which holds if and only if the open set O in ( 13) can be taken such that Γ ⊂ O [6].
For each m = (m 1 , . . ., m ℓ ) ∈ Σ ℓ , s m is a contracting similarity of the form where (with the convention that an empty product equals 1) The inverse of s m then has the representation Setting Σ ∅ := {∅}, Γ ∅ := Γ and s ∅ (x) := x, we define Σ : Given such a Γ, and a collection (p 1 , . . ., p M ) of positive weights (or "probabilities") satisfying 0 < p m < 1, m = 1, . . ., M, and there exists [15,Sections 4 & 5] a positive Borel-regular finite measure µ supported on Γ, unique up to normalisation, called a self-similar [23] (also known as invariant [15] or balanced [2]) measure, such that µ(E) = M m=1 p m µ(s −1 m (E)) for every measurable set E ⊂ R n .For such a measure, by [23,Thm. 2.1] the OSC (13) implies that µ(s m (Γ) ∩ s m ′ (Γ)) = 0 for each m ̸ = m ′ , and as a consequence we find that for m = (m 1 , . . ., m ℓ ) ∈ Σ, and any µ-measurable function f , where (again with the convention that an empty product equals 1) In particular, Given an IFS attractor Γ and two self-similar measures µ and µ ′ , with associated weights (p 1 , . . ., p M ) and (p ′ 1 , . . ., p ′ M ), we define t * > 0, if it exists, to be the largest positive real number such that the integral I Γ,Γ converges for 0 ≤ t < t * .In [13,Lem. A.4] we showed that if Γ is disjoint then t * exists and is the unique positive real solution of the equation Our conjecture is that the same holds for non-disjoint Γ, under the assumption of the OSC (13).
As yet, we have not been able to prove this conjecture in its full generality.However, we shall proceed under the assumption that it holds, noting that it is well known to hold, with t * = d, in the special case where µ = µ ′ = H d | Γ for d = dim H (Γ) (see e.g.[13, Corollary A.2]), which is the case of relevance for the integral equation application from [6] that we study in §7.We comment that when µ = µ ′ the quantity t * is sometimes referred to as the "electrostatic correlation dimension" of µ [19,Defn 6].This and related notions of the "dimension" of a measure µ give, amongst other things, lower bounds on the Hausdorff dimension of the support of µ (which may be strictly smaller than that of Γ), and important information about the asymptotic behaviour of the Fourier transform of µ -see e.g.[24].
We note that if the IFS is homogeneous then (21) can be solved analytically to give 1 A key step in the proof of this result is showing that µ(Γ m ′ ) = p m ′ µ(Γ) for m ′ = 1, . . ., M .To prove the latter, we first note that, by (18) and the positivity and self-similarity of µ, if E is measurable and µ(E) = 0 then µ(s −1 m (E)) = 0 for m = 1, . . ., M .In particular, since µ(sm(Γ) ∩ s m ′ (Γ)) = 0 for m ̸ = m ′ , we have that µ(s −1 m (Γm ∩ Γ m ′ )) = 0 for m ̸ = m ′ .Then, using the fact that µ is supported in Γ, we have µ( which reduces to t * = d in the case Self-similar measures sometimes possess useful symmetry properties.Let T : R n → R n be an isometry, with T (x) = A T x + δ T for some orthogonal matrix A T and some translation vector δ T .We say that µ is invariant under T if µ(T (E)) = µ(E) for all measurable E ⊂ R n .If µ is invariant under T then µ is also invariant under T −1 , the push-forward measure µ • T −1 : E → µ(T −1 (E)) coincides with µ, and T (Γ) = Γ, so that by [4,Thm. 3.6.1]we have that, for all measurable f ,
As mentioned in §1, our approach to deriving representation formulas for (1) will be based on decomposing the integral I Γ,Γ into sums of integrals over pairs of subsets of Γ.For any two vector indices m, m ′ ∈ Σ (possibly of different lengths) we define the sub-integral Note that the original integral I Γ,Γ = I ∅,∅ is included in this definition.Central to our approach will be identifying when two sub-integrals I m,m ′ and I n,n ′ are similar, in the sense that for some a > 0 and b ∈ R that we can determine explicitly in terms of the parameters of the IFS and the measures µ and µ ′ .
The simplest instance of similarity occurs when µ = µ ′ , in which case the symmetry of the integrand (i.e. the fact that Φ t (x, y) = Φ t (y, x)), combined with Fubini's theorem, provides the following elementary result (which was used in the example in §1, for which I 1,2 = I 2,1 etc.).
Other instances of similarity may be associated with the IFS structure (as for the derivation of ( 5) and (8), in the example in §1), and/or with other geometrical symmetries (as for the observation that I 1,2 = I 2,3 etc., in the example in §1).The following result provides sufficient conditions under which a given pair of sub-integrals I m,m ′ and I n,n ′ are similar in this manner.We recall that the notion of a self-similar measure being invariant under an isometry T was defined in §2, and that the question of determining for which isometries this holds was discussed in Remark 2.2.If no non-trivial isometries can be determined for µ or µ ′ one can always take T and T ′ to be the identity in the following.
Let T and T ′ be isometries of R n such that µ and µ ′ are invariant under T and T ′ respectively.Suppose there exists ϱ > 0 such that Then Proof.By (19), and the invariance of µ and µ ′ with respect to T and T ′ , we have that Condition (25) implies that which in turn, by (3), implies that and combining this with ( 27)-( 28) gives the result.
The following result provides an equivalent characterization of the sufficient condition (25) in terms of the scaling factors, orthogonal matrices and translations associated with the maps s m , s m ′ , s n , s n ′ , T and T ′ (see (16) for the definition of the notation ρ m , A m , δ m etc.).We remark that a necessary and sufficient condition for (30 Then condition (25) holds if and only if the following three conditions are satisfied: Proof.We first note that ( 25) is equivalent to It is easy to check that (30)-(32) are sufficient for (33) (and hence for ( 25)).To see that they are also necessary, suppose that (33) holds.Then taking x = y = 0 in (33) gives (32).Combining this with (33) gives and taking first x = 0 and y ̸ = 0, then x ̸ = 0 and y = 0 in this equation gives from which we deduce (30).Finally, combining (30) with (34) gives Now we note that if B is an n × n matrix and |x − By| = |x − y| for all x, y ∈ R n then B is the identity matrix.To prove this, suppose that B is not the identity matrix.Then there exists y ̸ = 0 such that By ̸ = y, and setting x = By gives |x − By| = 0 and is the identity matrix, which is equivalent to (31).
Remark 3.4.If A m , A m ′ , A n and A n ′ are all equal to the identity matrix, and A T = A T ′ , then condition (31) is automatically satisfied.If, further, δ T and δ T ′ are both zero, and condition (30) holds, then condition (32) reduces to

Algorithm for deriving representation formulas
We present our algorithm for deriving representation formulas for the integral I Γ,Γ in Algorithm 1 below.The output of the algorithm, when it terminates, is a linear system of the form where x = [I m s,1 ,m ′ s,1 , . . ., I m s,ns ,m ′ s,ns ] T ∈ R ns is a vector of "fundamental" singular sub-integrals (the subscript s standing for "singular"), with the original integral I m s,1 ,m ′ s,1 = I Γ,Γ as its first entry, r = [I m r,1 ,m ′ r,1 , . . ., I m r,nr ,m ′ r,nr ] T ∈ R nr is a vector of "fundamental" regular sub-integrals (the subscript r standing for "regular"), A ∈ R ns×ns and B ∈ R ns×nr are matrices, and b ∈ R ns is a vector of logarithmic terms present only in the case t = 0.The algorithm is based on repeated subdivision of the integration domain and the identification of similarities between the resulting sub-integrals (in the sense of ( 24)), and the word "fundamental" refers to a sub-integral which, when encountered in the subdivision algorithm, is not found to be similar to any other sub-integral previously encountered.Whether the algorithm terminates, and the resulting lengths n s and n r of the vectors x and r, depends on the measures µ and µ ′ , as we shall demonstrate in §5.
If the algorithm terminates, one can obtain numerical approximations to the values of the integrals in x, including the original integral I Γ,Γ , by applying a suitable quadrature rule to the regular integrals in r, solving the system (38), and extracting the relevant entry from the solution vector x.We discuss this in more detail in §6.Remark 4.2.Our algorithm (in line 9) requires us to specify a subdivision strategy.With reference to the notation in Algorithm 1, we have considered two such strategies: • Strategy 1: always subdivide both Γ ms,n and Γ m ′ s,n , i.e. take (with (∅, m) interpreted as (m)) • Strategy 2: subdivide only the larger of Γ ms,n and Γ m ′ s,n , i.e. take If the IFS is homogeneous then the two subdivision strategies coincide, since using Strategy 2 we never encounter pairs of subsets with different diameters.
Remark 4.3.Our algorithm (in line 11) requires the user to be able to determine whether an integral I n,n ′ is singular or regular, i.e. whether Γ n intersects Γ n ′ non-trivially or not.Deriving a criterion for this based solely on the indices n and n ′ and the IFS parameters appears to be an open problem.However, for the examples considered in §5 we were able to determine this by inspection on a case-by-case basis.We emphasize that one does not need to specify the type of singularity, i.e. the dimension of Γ n ∩ Γ n ′ , merely whether Γ n ∩ Γ n ′ is empty or not.
Remark 4.4.Our algorithm (in lines 12 and 20) requires a way of checking for "similarity" of pairs of subintegrals.For this we use Propositions 3.1-3.3,combined with a user-provided list of isometries T and T ′ under which µ and µ ′ are respectively invariant.Then, to verify (25) in Proposition 3.2, we use Proposition 3.3: we first check (30), then (31), then (32).As noted in Remark 3.4, in the special case where A m , A m ′ , A n and A n ′ are all equal to the identity matrix, A T = A T ′ , and δ T and δ T ′ are both zero, it is sufficient to check (30) and then (37).
The question of how to determine the permitted isometries T and T ′ was discussed in Remark 2.2.If the user is not able to provide the full list of isometries for µ and µ ′ , it may be that the algorithm still terminates, but does so with a larger number n s of fundamental singular integrals than would be obtained with the full list of isometries.However, in §5.5 we provide an example where failing to specify a non-trivial isometry would lead to non-termination of the algorithm.
Remark 4.5.The matrix A (when the algorithm terminates) depends not only on the IFS {s 1 , . . ., s M } and the weights p 1 , . . ., p M , but also on the subdivision strategy used in line 9 of the algorithm.For both subdivision strategies described in Remark 4.2, the first column of A is guaranteed to be of the form (α, 0 . . ., 0) T for some α > 0, because of the fact that µ(Γ m ∩ Γ m ′ ) = 0 for m ̸ = m ′ .For all the examples in §5.1- §5.4 the matrix A is upper triangular, with diagonal entries that are all non-zero when t < t * , so that A is invertible when t < t * .However, upper-triangularity of A is not guaranteed in general, as §5.5 illustrates (see in particular (47)), and proving that A is invertible whenever the algorithm terminates remains an open problem.
Augment A by an n × 1 column of zeros.
Set A n,ns = −1.end else if I n,n ′ is similar to an existing entry of r then Let i(n, n ′ ) and a(n, n ′ ) and b(n, n ′ ) be such that

Examples
We now report the results of applying Algorithm 1 to some standard examples.

Sierpinski triangle, Hausdorff measure
We first consider the case where Γ ⊂ R 2 is the Sierpinski triangle, the attractor of the homogeneous IFS with M = 3 and for which d = dim H (Γ) = log 3/ log 2 ≈ 1.59.The first two levels of subdivision of Γ are illustrated in Figure 2. We shall assume for simplicity that µ = µ ′ = H d | Γ , so that t * = d.Since we are working with Hausdorff measure, as discussed in Remark 2.2 the isometries T under which µ is invariant are precisely those for which T (Γ) = Γ, which in this case are the elements of the dihedral group D 3 corresponding to the symmetries of the equilateral triangle.The two subdivision strategies (39) and (40) coincide, since the IFS is homogeneous, and Algorithm 1 terminates after finding two fundamental singular sub-integrals, with x = (I Γ,Γ , I 1,2 ) T .The integral I 1,2 captures the interaction between neighbouring subsets of Γ of the same size, intersecting at a point.The linear system (38) satisfied by these unknowns is: where and, solving the system, we obtain the representation formulas

Vicsek fractal, Hausdorff measure
Next, we consider the case where Γ ⊂ R 2 is the Vicsek fractal (shown in Figure 3), the attractor of the homogeneous IFS with M = 5 and , In this case the isometries T under which µ is invariant are the elements of the dihedral group D 4 corresponding to the symmetries of the square.The first two levels of subdivision of Γ are illustrated in Figure 3, from which it is clear that the situation is similar to that for the Sierpinski triangle of §5.1, as the only new singularities at level one are point singularities, which are similar (in the sense of ( 25)) to those arising at level two.Again, our two subdivision strategies coincide because the IFS is homogeneous, and Algorithm 1 terminates after finding two fundamental singular subintegrals, with x = (I Γ,Γ , I 1,5 ) T .For brevity we do not present the full linear system satisfied by these unknowns, but rather just report the matrix A of (38), which is where , t ∈ [0, d).

Sierpinski carpet, Hausdorff measure
Next, we consider the case where Γ ⊂ R 2 is the Sierpinski carpet, the attractor of the homogeneous IFS with M = 8 and Figure 5: The level one, two and three subdivisions of the Koch snowflake.The integrals I 1,2 and I 2,4 capture the interaction between neighbouring subsets of the same size, intersecting along a line segment and at a point, respectively.In this case, the matrix of (38) is: , where, for t ∈ [0, d),

Koch snowflake, Lebesgue measure
Next, we consider the case where Γ ⊂ R 2 is the Koch snowflake, the attractor of the non-homogeneous IFS with M = 7 and for which d = dim H (Γ) = 2.The first three levels of subdivision of Γ are illustrated in Figure 5.We assume that µ and µ ′ are both equal to the Lebesgue measure on R where are linear combinations of regular integrals.Hence , and solving the system gives

Γ = [0, 1], including non-terminating examples
We now consider a class of simple one-dimensional examples that illustrates the dependence of the output of Algorithm 1 on the choice of measures µ and µ ′ , and the fact that it does not always terminate.
If ρ ∈ (0, 1/2) then the IFS is inhomogeneous, so Strategy 1 and Strategy 2 differ, and the outcome of Algorithm 1 depends on the choice of strategy, and on the measure µ.We consider four cases: • none of which is found to be similar to any other.To see this, note that for a sub-integral Γ m,m ′ in this sequence, with m = (1, 2, . . ., 2) (k 2's) and m ′ = (2, 1, . . ., 1) (k 1's), we have Then since 0 < ρ < 1 − ρ < 1, the sequence (R k ) ∞ k=0 is monotonically increasing, with R k ≥ 1 for k ≥ 1.This implies that (29) (and hence (30)) is not satisfied by any pair of elements of the sequence (48), except for I 1,2 and I 122,211 .However, since µ is not Lebesgue measure, I 1,2 and I 122,211 are not found to be similar, because µ is not invariant under T ref (so one cannot use it in Proposition 3.2), and (32) fails with T = T ′ the identity.
• Case 4: Non-Lebesgue measure, Strategy 2: 2 The matrix A arising in (38) in this case, with x = (IΓ,Γ, I1,2, I12,21) T , is given by where ω1 = ρ 2−t and ω2 In this case, Algorithm 1 can only terminate if ρ is a solution of a polynomial equation for some j, k ∈ N 0 with either j > 0 or k > 0. In particular, Algorithm 1 does not terminate if ρ is transcendental.To see that (50) is necessary for termination of the algorithm, we note that, in the subdivision of I 1,2 , Strategy 2 will produce pairs of subsets of Γ (intervals) that intersect at the point x = ρ, and the lengths of the intervals in each pair will be in the ratio ρ(1 − ρ) j : ρ k (1 − ρ) for some j, k ∈ N 0 .For the sub-integrals associated to any two distinct pairs of such intervals to be found to be similar, the ratio of their lengths must coincide (by ( 29)), implying that

Numerical quadrature and error estimates
Once Algorithm 1 has been applied, and the system (38) has been solved, producing a representation formula for the singular integral I Γ,Γ in terms of regular sub-integrals, a numerical approximation of I Γ,Γ can be obtained by applying a suitable quadrature rule to the regular sub-integrals.Let us call the resulting approximation Q Γ,Γ .We discuss some possible choices of quadrature rule below.But first we make a general comment on the error analysis of such approximations.Suppose that the quadrature rule chosen can compute each of the regular sub-integrals in the vector r with absolute error ≤ E for some E ≥ 0. Then the absolute quadrature error in computing I Γ,Γ using the representation formula (38) can be bounded by The constant ∥A −1 ∥ ∞ ∥B∥ ∞ depends on the problem at hand, and is expected to blow up as t → t * .Indeed, for the examples in §5.1- §5.4 (for which for some constant C > 0, independent of t.This follows from the fact that in all these examples the constant σ 1 = O(d − t) as t → d, while σ 2 , σ 3 etc.remain bounded away from zero in this limit.
We now return to the choice of quadrature rule for the approximation of the regular sub-integrals in r, which are double integrals of smooth functions over pairs of self-similar subsets of Γ with respect to a pair of invariant measures µ and µ ′ .We shall restrict our attention to tensor product quadrature rules, so that it suffices to consider methods for evaluating a single integral of a smooth function over a single self-similar subset Γ m of Γ with respect to a single invariant measure µ.In fact, it is enough to consider the case Γ m = Γ, since the more general case can then be treated using (19).Hence we consider quadrature rules for the evaluation of the integral for an integrand f that is smooth in a neighbourhood of Γ.We consider three types of quadrature: • Gauss rules: Highly accurate, but currently only practically applicable for the case n = 1, i.e.Γ ⊂ R. 3• Composite barycentre rules: Less accurate than Gauss rules, but can be applied to Γ ⊂ R n for n > 1.
• Chaos game rules: Monte-Carlo type rules which converge (in expectation) at a relatively slow but dimension-independent rate, which makes them well-suited to high-dimensional problems (large d = dim H (Γ)).
In the following three sections we provide further details of these methods, and any theory supporting them, before comparing their performance numerically in §7.

Gauss rules in the case Γ ⊂ R
In general, N -point Gauss rules require the existence of a set of polynomials {p j } N j=0 , orthogonal with respect to the measure µ.A sufficient condition for the existence of such polynomials is positivity of the Hankel determinant, which in the case of self-similar invariant measures, is implied by supp µ = Γ having infinitely many points [11, §1.1].We then define the N -point Gauss rule on Γ as where x j , j = 1, . . ., N , are the zeros of p N .Gauss rules are interpolatory, so the weights (also known as Christoffel numbers) may be defined by w j := Γ ℓ j (x)dµ(x), where ℓ j is the jth corresponding Lagrange polynomial (see e.g [29, (5.3)]).The weights are positive -see for example [12,Theorem 1.46], which generalises to any positive measure µ.
For classical µ, a range of algorithms (see e.g.[14,26]) exist for efficient O(N ) computation of the weights and nodes in (53).However, standard approaches involving polynomial sampling break down for singular measures [17,18].This presents an obstacle for the evaluation of (52) in our context of self-similar invariant measures, which are in general singular when dim H (Γ) ̸ = n.However, in the special case n = 1, where Γ ⊂ R, this issue can be overcome by applying the stable Stieltjes technique proposed in [17, §5] 4 .It seems that a stable and efficient algorithm for the evaluation of Gauss rules for the case where Γ ⊂ R n , n > 1, has not yet been developed.Hence in this paper we only consider Gauss rules for the case where Γ ⊂ R.
The error analysis for the Gauss rule follows the standard approach, giving the usual exponential convergence as N → ∞.In the following, Hull(Γ) denotes the convex hull of Γ.
Theorem 6.1.Let Γ ⊂ R be an IFS attractor and let µ be a self-similar measure supported on Γ.
If f is analytic in a neighbourhood of Hull(Γ) ⊂ R, then for some constants C > 0 and c > 0, independent of N .

The composite barycentre rule
The basic idea of the composite barycentre rule (for more detail see [13]) is to partition Γ into a union of self-similar subsets of approximately equal diameter, then to approximate f on each subset by its (constant) value at the barycentre of each subset.Given a maximum mesh width h > 0, we define a partition of Γ using the following set of indices: The composite barycentre rule is then defined as where the weights and nodes are defined by w m := µ(Γ m ) and x m := Γm x dµ(x)/µ(Γ m ) respectively.The weights and nodes can be computed using simple formulas involving the IFS parameters of (11), as (see [13, (28)-(30)], and recall (20)) 4 There is an error in [17,Equation (32)]: Γ n i,m Γ n i,m+1 δi(rm + rm+1) should be replaced by 2 with where I is the n × n identity matrix and ρ m , A m and δ m , m = 1, . . ., M , are as in (11).
The error analysis of the composite barycentre rule follows a standard Taylor series approximation argument.The following is a simplified version of results in [13].
(i) If f is Lipschitz continuous on Hull(Γ) then for some C > 0 independent of h.
(ii) If f is differentiable in a neighbourhood of Hull(Γ), and its gradient is Lipschitz continuous on Hull(Γ) then for some C > 0 independent of h.
If the IFS defining Γ is homogeneous then h and h 2 on the right-hand sides of the above estimates can be replaced by N −1/d and N −2/d respectively, where N := |L h (Γ)|.

Chaos game quadrature
Chaos game quadrature, described, e.g., in [10, (3.22)-(3.23)]and [16, § 6.3.1], is a Monte-Carlo type approach, defined by the following procedure: (i) Choose some x 0 ∈ R n , e.g.x 0 = x Γ , the barycentre of Γ; (ii) Select a realisation of the sequence {m j } j∈N of i.i.d.random variables taking values in {1, . . ., M } with probabilities {p 1 , . . ., p M }; (iii) Construct the stochastic sequence x j = s m j (x j−1 ) for j ∈ N; (iv) For a given N ∈ N, define the chaos game quadrature approximation by For continuous f , the chaos game rule (59) will converge to (52) with probability one (see the arguments in the appendix of [10]).While no error estimates were provided in [10] or [16], in the numerical experiments of [13, §6] and §7 below, convergence in expectation was observed at a rate consistent with an estimate of the form

Numerical results and applications
Algorithm 1, and the quadrature approximations described in §6, have been implemented in the open-source Julia code IFSIntegrals, available at www.github.com/AndrewGibbs/IFSintegrals.In this section we present numerical results illustrating the accuracy of our approximations, comparing different quadrature approaches, and applying our method in the context of a boundary element method for acoustic scattering by fractal screens.

Sierpinski triangle, Vicsek fractal, Sierpinski carpet, Koch snowflake
We first consider the application of our approach to the attractors considered in §5.1- §5.4,namely the Sierpinski triangle, Vicsek fractal, Sierpinski carpet and Koch snowflake.However, in contrast to §5.1- §5.4,where representation formulas were presented for the standard case where µ = µ ′ = H d | Γ (Lebesgue measure in the case of the Koch snowflake), to demonstrate the generality of our approach we present numerical results for completely generic self-similar measures µ ̸ = µ ′ , with randomly chosen probability weights p m and p ′ m , as detailed in Table 1.Table 1 also documents the resulting values of t * , as computed by solving (21), as well as the numbers n s , n r , of fundamental singular and regular sub-integrals discovered by our algorithm.In all cases our algorithm terminated, using the same subdivision strategies as in §5.1- §5.4,producing an invertible matrix A. However, for these non-standard examples there are no nontrivial isometries under which the measures are invariant, so our algorithm took T = T ′ to be the identity throughout.As a result (cf. the related discussion in Remark 4.4), the linear systems (38) are larger than those obtained in the standard case documented in §5.1- §5.4,where additional symmetries of the measures could be exploited.
In Figure 7 (solid curves) we plot the relative error in our quadrature approximation for I Γ,Γ , for three values of t = 0, 0.5, 1, obtained by solving the linear system (38) obtained by Algorithm 1 (using subdivision strategy 2 for the Koch snowflake), combined with composite barycentre rule quadrature for the evaluation of the regular sub-integrals, for different values of the maximum mesh width h.In more detail, we plot errors for h = diam(Γ)ρ ℓ , for ℓ = 0, . . ., ℓ ref − 1, where ℓ ref is the value of ℓ used for the reference solution (which is computed using the same method).For the three homogeneous attractors, we take ρ = ρ 1 = . . .= ρ M (the common contraction factor), while for the Koch snowflake, we take ρ = 1/ √ 3 (the largest contraction factor).For the Sierpinski triangle ℓ ref = 10, for the Vicsek fractal ℓ ref = 7, and for the Sierpinski Carpet and Koch snowflake ℓ ref = 6.According to our theory, we expect our method to give O(h 2 ) error, by Theorem 6.2(ii) and (51), and this is exactly the rate we observe in our numerical results in Figure 7.
In Figure 7 (dashed curves) we also show results obtained using the method of our previous paper [13, (48)], which we refer to as the "old method".This method is accurate for disjoint IFS attractors, but is expected to perform less well for non-disjoint attractors, because it only applies self-similarity to deal with the self-interaction integrals, and treats all other sub-integrals as being regular.Precisely, the old method corresponds to taking the equation corresponding to the first row in the linear system (38) obtained by Algorithm 1, solving this equation for I Γ,Γ , then applying the composite barycentre rule not just to the regular sub-integrals coming from the right-hand side of (38), but also to the fundamental singular sub-integrals I m s,i ,m ′ s,i , i = 2, . . ., n s .We expect that the resulting quadrature approximation should converge to I Γ,Γ as h → 0, but at a slower rate than our new method, because of the inaccurate treatment of the singular sub-integrals.This is borne out in our numerical results in Figure 7, with the errors for the old method being significantly larger than those for the new method.To quantify these observations, we present in Table 2 the empirical convergence rates (computed from the errors for the two smallest h values) observed for the old method for each of the three t values considered.The deviation from O(h 2 ) convergence is different for each example, but clearly increases as t, the strength of the singularity, increases, as one would expect.
For all the experiments in Figure 7, the total number of quadrature points N tot grows like N tot ≈ Ch −2d as h → 0, with the value of C depending on the number of fundamental regular sub-integrals that need to be evaluated.(Recall from §6.2 that for each regular sub-integral we use a tensor product rule with N 2 points, where N ≈ C ′ h −d for some C ′ .)For each choice of attractor, the value of C for the new method is slightly smaller than that for the old method, because the new method takes greater advantage of similarities between regular sub-integrals.The value of N tot used for the reference solutions is 1,291,401,630 for the Sierpinski triangle, 2,382,812,500 for the Vicsek fractal, 18,790,481,920 for the Sierpinski carpet, and 379,046,894,100 for the Koch snowflake.

Unit interval experiments
We now consider an attractor Γ ⊂ R, so that we can investigate the performance of the Gauss quadrature discussed in §6.1.The classic example of an IFS attractor Γ ⊂ R is the Cantor set, but since this is disjoint (in the sense of ( 15)), it can already be treated by our old method (of [13]).To demonstrate the efficacy of our new method for dealing with non-disjoint attractors we consider the case where Γ = [0, 1] ⊂ R, which, as discussed in §5.5 (taking ρ = 1/2), is the attractor of the homogeneous IFS with M = 2, s 1 (x) = x/2 and s 2 (x) = x/2 + 1/2.We consider the case where µ = µ ′ , with p 1 = p ′ 1 = 1/3 and p 2 = p ′ 2 = 2/3, so that, by (22), t * = log(9/5)/ log 2 ≈ 0.848.As discussed in §5.5, for this problem Algorithm 1 finds just two fundamental singular sub-integrals, I Γ,Γ and I 1,2 , and two fundamental regular sub-integrals, I 11,21 and I 11,22 .
In Figure 8 we report relative errors for the computation of I Γ,Γ with t = 0 and t = 1/2, using Algorithm 1 combined with Gauss, composite barycentre, and chaos game quadrature for the evaluation of the regular sub-integrals.For each method the total number of quadrature points satisfies N tot = 2N 2 , where N is the number of points used for each of the two iterated integrals in each of the two fundamental regular sub-integrals (recall that we are using tensor product rules).For the composite barycentre rule we have N ≈ 1 4h in this case.As the reference solution we use the result obtained using the Gauss rule with N = 100, which corresponds to N tot = 20, 000.For the Gauss rule we see the expected root-exponential O(e −c √ Ntot ) convergence predicted by Theorem 6.1 and (51), with c ≈ 1.77 and c ≈ 2.31 (for t = 0 and t = 1/2 respectively) in this case.As a result, the singular double integral I Γ,Γ can be computed to machine precision using N ≈ 10, which corresponds to N tot ≈ 200 quadrature points.The barycentre rule is significantly less accurate, converging like

Application to Hausdorff BEM for acoustic scattering by fractal screens
We conclude by demonstrating how our new representation formulas and resulting quadrature rules can be used to compute the scattering of acoustic waves by fractal screens using the "Hausdorff boundary element method" (BEM) of [6].For full details of the scattering problem and the Hausdorff BEM we refer the reader to [6] and the references therein; here we merely provide a brief overview.
The underlying scattering problem under consideration is the three-dimensional time-harmonic acoustic scattering of an incident plane wave e ikx•ϑ (for x = (x 1 , x 2 , x 3 ) T ∈ R 3 , wavenumber k > 0 and unit direction vector ϑ ∈ R 3 ) by a fractal planar screen Γ × {0} ⊂ R 3 , where Γ ⊂ R 2 is the attractor of an IFS satisfying the open set condition.Assuming that the total wave field u (which is a solution of the Helmholtz equation ∆u + k 2 u = 0 in R 3 \ (Γ × {0})) satisfies homogeneous Dirichlet ("sound soft") boundary conditions on the screen, it was shown in [6] that the scattering problem can be reduced to the solution of the integral equation Here S is the single layer boundary integral operator defined by Sϕ(x) = Γ Φ(x, y)ϕ(y) dy (with the integral interpreted in a suitable distributional sense), where Φ(x, y) := e ik|x−y| /(4π|x − y|) is the fundamental solution of the Helmholtz equation in three dimensions, ϕ is the unknown jump in the x 3 -derivative of u across the screen, and f is a known function depending on the incident wave.
The Hausdorff BEM in [6] discretises (60) using a Galerkin method with a numerical approximation space of piecewise constant functions multiplied by the Hausdorff measure H d | Γ .The mesh used for the piecewise constant functions is of the same form as that used in the composite barycentre rule in §6.2 -having chosen a maximum BEM mesh width h BEM we partition Γ using the index When Γ m and Γ n are disjoint, the integral (61) has a smooth integrand and can be evaluated using the composite barycentre rule with some maximum mesh width h ≤ h BEM , with error O(h 2 ) (by Theorem 6.2(ii)).When Γ m ∩ Γ n is non-empty the integral (61) is singular, and to evaluate it we adopt a singularity subtraction approach, writing (62) where Φ * (x, y) := Φ(x, y) − (4π|x − y|) −1 = (e ik|x−y| − 1)/(4π|x − y|).The first integral on the righthand side of (62) can be evaluated using the methods of this paper with t = 1.In more detail, if Γ m = Γ n this first integral will be similar to I Γ,Γ , and if Γ m ̸ = Γ n it will be similar to one of the other fundamental singular sub-integrals encountered in Algorithm 1.In both cases it can be evaluated by combining Algorithm 1 with the composite barycentre rule, again with mesh width h ≤ h BEM and error O(h 2 ).The second integral on the right-hand side of (62) has a Lipschitz continuous integrand, and hence can be evaluated using the composite barycentre rule directly.According to Theorem 6.2(i), the error in this approximation is guaranteed to be O(h).In fact, for disjoint homogeneous attractors the error in evaluating this second term was proved in [13, Proposition 5.5] to be O(h 2 ), and experiments in [13, Figure 8(a)] suggest that the same may be true for certain non-homogeneous disjoint attractors.In Figure 9 we present numerical results suggesting, furthermore, that the same may also be true for certain non-disjoint attractors.The plots in Figure 9 show the relative error (against a high order reference solution) in computing the second term in (62) using the composite barycentre rule, for the four attractors from §5.1- §5.4 and a range of wavenumbers, in the case where Γ m = Γ n = Γ.This case is chosen since it represents the most difficult case, in which Γ m and Γ n have full overlap.For all four examples we clearly observe O(h 2 ) error in the numerical results.However, we leave theoretical justification of this observation to future work.

Remark 4 . 1 .
If Γ is disjoint then the algorithm terminates with n s = 1 and recovers the result of[13, Thm 4.6].

Algorithm 1 :
n ′ as a new entry in r.Increment n r = n r + 1. Augment B by an n × 1 column of zeros.Set B n,nr = 1.end end end end Algorithm for deriving representation formulas for the integral (1).

Figure 2 :
Figure 2: Level 1 and 2 subsets of the Sierpinski triangle.

Figure 4 :
Figure 4: Level 1 and 2 subsets of the Sierpinski carpet.
,for which d = dim H (Γ) = log 8/ log 3 ≈ 1.89.The first two levels of subdivision of Γ are illustrated in Figure2.We again assume that µ = µ ′ = H d | Γ , so that t * = d.As for the previous example, the isometries T under which µ is invariant are the elements of the dihedral group D 4 .Again the two subdivision strategies (39) and (40) coincide, since the IFS is homogeneous, and now Algorithm 1 terminates after finding three fundamental singular sub-integrals, with x = (I Γ,Γ , I 1,2 , I 2,4 ) T .

Figure 7 :
Figure 7: Convergence of composite barycentre rule quadrature approximations for the four examples of §7.1, using the full linear system (38) (labelled "new method") and using only the first row of (38) (labelled "old method").

Figure 8 :
Figure 8: Comparison of the performance of the three quadrature rules presented in §6.1- §6.3, applied in the context of Algorithm 1, for the example considered in §7.2, which concerns the evaluation of singular double integrals on the unit interval Γ = [0, 1] with respect to a self-similar measure.

Figure 9 :
Figure 9: Convergence of composite barycentre rule approximations to the second integral on the right-hand side of (62) in the case Γ m = Γ n = Γ, as considered in §7.3.

)
Example 2.1.By choosing p m = ρ d m for m = 1, . . ., M and normalising appropriately, we can obtain µ = H d | Γ , where H d is the d-dimensional Hausdorff measure on R n (note that in this case (18) holds by (14)).We recall that H n is proportional to n-dimensional Lebesgue measure [9, §3.1].
Remark 2.2.Determining a complete list of isometries T under which a given self-similar measure µ is invariant, directly from the IFS {s 1 , ..., s M } and weights p 1 , ..., p M , appears to be an open problem.However, for specific examples it is often straightforward to determine the admissible T , as we demonstrate in §5 and §7.In the case µ = H d | Γ (see §5.1- §5.4), a necessary and sufficient condition for µ to be invariant under T is thatT (Γ) = Γ, because H d is invariant under isometries of R n .For µ ̸ = H d | Γ itis still necessary that T (Γ) = Γ, but no longer sufficient (see §5.5).Generically, a self-similar measure µ may not be invariant under any non-trivial isometries T (see §7.1).
Subdivide Γ ms,n = n∈In Γ n and Γ m ′ s,n = n ′ ∈I ′ n Γ n ′ according to some subdivision strategy producing index sets I n , I ′ n ⊂ Σ.
1, and n = 0. 2 while n < n s do 3Increment n = n + 1.4Append a 1 × n s row of zeros to A. 5 Set A n,n = 1.6 if B is non-empty then 7 Append a 1 × n r row of zeros to B. end 8 Append a single zero entry to b.9 10 for each pair (n, n ′ ) ∈ I n × I ′ n do 11 if I n,n ′ is singular then if I n,n ′ is similar to an existing entry of x then Let i(n, n ′ ), a(n, n ′ ) and b(n, n ′ ) be such that 2, restricted to Γ, so, again, t * = d.(Asmentioned in Example 2.1, µ is proportional to H 2 | Γ .)Theisometries T under which µ is invariant are the elements of the dihedral group D 6 corresponding to the symmetries of the hexagon.In this case both subdivision strategies produce a terminating algorithm, but with different results.Γ,Γ , I 1,2 , I 2,3 , I 11,25 ) T .The integral I 1,2 captures the interaction between neighbouring subsets of Γ, intersecting along a Koch curve.The integrals I 2,3 and I 11,25 both capture point interactions, but of different types: in I 2,3 the two interacting subsets are the same size, while in I 11,25 one is three times the diameter of the other.I 11,25 arises as a new fundamental sub-integral in the subdivision of I 1,2 , but in the subdivision of I 11,25 one obtains just one singular sub-integral, I 117,255 , which is similar to I 11,25 , so the algorithm terminates.For brevity we do not report the resulting linear system satisfied by (I Γ,Γ , I 1,2 , I 2,3 , I 11,25 ) T , but instead present the simpler result obtained with Strategy 2 (subdividing the subset with the largest diameter), for which Algorithm 1 terminates after finding only three fundamental singular subintegrals, with x = (I Γ,Γ , I 1,2 , I 2,3 ) T .With this strategy, in the subdivision of I 1,2 we subdivide only Γ 1 , leaving Γ 2 intact, obtaining the edge interaction sub-integrals I 12,2 and I 17,2 , both of which are similar to I 1,2 , and the point interaction sub-integral I 11,2 , which is similar to I 2,3 .The resulting linear system is giving (50).If ρ is the solution of a polynomial (50) then Algorithm 1 may terminate, but the number of fundamental singular sub-integrals encountered will depend on ρ.For instance, if ρ = ρ * the algorithm terminates with four fundamental singular sub-integrals I Γ,Γ , I 1,2 , I 1,21 and I 12,21 (since in this case I 122,211 is similar to I 1,21

Table 1 :
(21)om weights used for the self-similar measures considered in §7.1, along with the resulting value of t * solving(21), and the numbers n s and n r of fundamental singular and regular integrals encountered in Algorithm 1.

Table 2 :
Empirical convergence rates for the old method.