Numerical Quadrature for Singular Integrals on Fractals

We present and analyse numerical quadrature rules for evaluating regular and singular integrals on self-similar fractal sets. The integration domain $\mathbb{R}^n$ is assumed to be the compact attractor of an iterated function system of contracting similarities satisfying the open set condition. Integration is with respect to any ``invariant'' (also known as ``balanced'' or ``self-similar'') measure supported on $\Gamma$, including in particular the Hausdorff measure $\mathcal{H}^d$ restricted to $\Gamma$, where $d$ is the Hausdorff dimension of $\Gamma$. Both single and double integrals are considered. Our focus is on composite quadrature rules in which integrals over $\Gamma$ are decomposed into sums of integrals over suitable partitions of $\Gamma$ into self-similar subsets. For certain singular integrands of logarithmic or algebraic type we show how in the context of such a partitioning the invariance property of the measure can be exploited to express the singular integral exactly in terms of regular integrals. For the evaluation of these regular integrals we adopt a composite barycentre rule, which for sufficiently regular integrands exhibits second-order convergence with respect to the maximum diameter of the subsets. As an application we show how this approach, combined with a singularity-subtraction technique, can be used to accurately evaluate the singular double integrals that arise in Hausdorff-measure Galerkin boundary element methods for acoustic wave scattering by fractal screens.


Introduction
In this paper we study numerical quadrature rules for the evaluation of integrals of the form where Γ and Γ are compact fractal subsets of R n -more precisely, the attractors of iterated function systems (IFSs) satisfying the open set condition (OSC) (see Subsection 2.2) -and µ and µ are "invariant" (also known as "balanced" or "self-similar") measures on Γ and Γ respectively (see Subsection 2.5) Our particular interest is in singular integrals, where, in the case of (1), f is singular at some point η ∈ Γ, and, in the case of (2), f is singular on x = y.
One context in which such integrals arise is in the discretization of certain boundary integral equation and volume integral equation formulations of boundary value problems for elliptic PDEs (such as the Laplace or Helmholtz equation) posed on domains with fractal boundary, for instance in the scattering of electromagnetic and acoustic waves by fractal obstacles [11,12], applications of which include antenna design in electrical engineering [35,38] and the quantification of the scattering effect of atmospheric ice crystals in climate modelling [37].Our main motivating example is the "Hausdorff boundary element method (BEM)" introduced in [9] for the solution of time-harmonic acoustic scattering in R n+1 (n = 1, 2) by a sound-soft fractal screen Γ screen ⊂ R n ×{0}, assumed to be the attractor of an IFS satisfying the OSC.The BEM proposed in [9] discretizes the associated single layer boundary integral equation on Γ screen using an approximation space consisting of products of the relevant Hausdorff measure with piecewise-constant functions on a "mesh" of Γ screen comprising self-similar fractal "elements", which are subsets of Γ screen obtained as scaled, rotated and translated copies of Γ screen via the IFS structure.The entries of the right-hand side vector in the Galerkin BEM system then involve integrals of the form (1) (with µ = H d | Γ ), where Γ is an element of the mesh and f (x) depends on the incident wave.The Galerkin BEM system matrix entries involve integrals of the form (2) (with µ = H d | Γ and µ = H d | Γ ), where Γ and Γ are elements of the mesh and f (x, y) = Φ(x, y), where Φ(x, y) is the fundamental solution of the Helmholtz equation in R n+1 , viz.
(1) 0 (k|x − y|), n = 1, e ik|x−y| 4π|x − y| , n = 2, where k > 0 is the wavenumber and H (1) ν denotes the Hankel function of the first kind of order ν.This choice of f (x, y) makes (2) singular when Γ = Γ .Although not studied in [9], one could also consider a collocation (as opposed to Galerkin) method for the same integral equation and approximation space, in which case the collocation matrix entries would involve integrals of the form (1) (with µ = H d | Γ ), where Γ is an element of the mesh and f (x) = Φ(x, η), with η denoting a collocation point, giving a singular integral when η ∈ Γ.A key goal of the current paper is to present a detailed derivation and rigorous error analysis of the quadrature rules used in the implementation of the Hausdorff BEM in [9].But we expect that the techniques we present will be of wider interest, since they apply to general invariant measures, to general regular integrands, and to a quite general class of singular integrands with logarithmic or algebraic singularities.
Our quadrature rules for (1) and ( 2) are based on decomposing integrals over Γ into sums of integrals over suitable partitions of Γ into self-similar subsets, generated using the IFS structure, in the same way that the Hausdorff BEM meshes are constructed in [9].After reviewing some preliminaries in Section 2, we start in Section 3 by considering regular integrands.Applying a one-point quadrature rule on each subset leads to a composite quadrature rule, which, when the quadrature nodes are chosen as the barycentres (with respect to µ or µ ) of the subsets, can achieve second-order convergence with respect to the maximum diameter of the subsets (Theorems 3.6 and 3.7).In Section 4 we then consider a special class of singular integrands, indexed by t ≥ 0, namely f (x) = Φ t (x, η) (in the case of (1)), and f (x, y) = Φ t (x, y) with Γ = Γ (in the case of (2)), where For these particular choices of f (and for a particular choice of η in the case of (1)), we use the fact (as noted in e.g.[7] for the case of Cantor sets) that the invariance property of µ and the homogeneity property of Φ t can be exploited to express the singular integral exactly in terms of regular integrals (Theorems 4.3 and 4.6), which can be evaluated using our composite barycentre rule, again with second-order accuracy (Corollaries 4.4 and 4.7).The results on Φ t are directly relevant to Hausdorff-BEM formulations of Laplace problems analogous to the Helmholtz ones described above.In Section 5 we combine this approach to integrating Φ t with a singularity-subtraction technique (cf.[21,3,36,2]) to propose and analyse a second-order accurate quadrature rule for the motivating example from [9] discussed above, exploiting the fact that the singular behaviour of Φ(x, y) matches that of Φ t (x, y) for t = n − 1, n = 1, 2.Here the main result is Theorem 5.7.In Section 6 we present some numerical results illustrating our theory.In the Appendix we collect some results concerning the integrability of singular functions with respect to invariant measures on IFS attractors.
For ease of reference, we list the quadrature rules that we propose: (25) for single regular integrals; • the barycentre rule Q Γ,Γ [f ] (32) for double regular integrals; • the rule Q h Γ,t,m (43) for single integrals of the singular integrand Φ t ; • the rule Q h Γ,Γ,t (48) for double integrals of the singular integral Φ t ; • the singularity-subtraction rule Q h Γ,Γ,Φ (61) for double integrals of the singular Helmholtz fundamental solution Φ; this rule reduces to (60) for small wavenumbers.
A link to our open-source implementation of these rules, and a pointer to an interactive notebook showing example usage, is provided in Section 6.
To put our results in the context of related work, we note that in the case n = 1, i.e. when Γ ⊂ R, Gauss quadrature rules can be derived for (1) (and applied to (2) iteratively), as discussed e.g. in [27,28].For sufficiently regular integrands these rules offer superior convergence rates when compared to our composite barycentre rules.However, one advantage of our low-order composite approach is that it is also generically applicable for Γ ⊂ R n , n > 1.By contrast, stable Gauss rules are not in general available for the case n > 1 (except when Γ is a subset of a line, in which case the n = 1 results apply).Moreover, the case n > 1 cannot in general be treated by taking Cartesian products of Gauss rules, since IFS attractors in R n , n > 1, are not in general the Cartesian product of IFS attractors in R. Furthermore, even when Γ ⊂ R n has such a Cartesian product structure, as is the case for the Cantor dust in Figure 1(I), the corresponding invariant measure is not the tensor-product measure of the respective lower-dimensional invariant measures (see e.g.[18,Proposition 7.1]).We stress, however, that if a stable Gauss rule (or any other quadrature rule for regular integrals) is available, it can be used in place of our composite barycentre rule within the context of our singularity-subtraction and invariance techniques for singular integrals, with corresponding analogues of the convergence results in Corollaries 4.4 and 4.7 and Theorem 5.7.
We also note that other quadrature rules for regular integrals on IFS attractors have been investigated in the abstract framework of uniform distributions and discrepancies in [23,13], and in the context of so-called "chaos games" (i.e.Monte-Carlo-type algorithms) in e.g.[19].In both settings, convergence is typically proved for general continuous integrands, but convergence rates for smoother integrands are not provided.We present a numerical comparison between our quadrature rules and a simple chaos game approach in Section 6.
We end this introduction with a comment relating to the practical evaluation of the quadrature rules presented in this paper.All our rules require knowledge of µ(Γ), since this quantity appears as a multiplicative factor in the formula (27) for the weights in the barycentre rule on which all our other quadrature rules are based.Somewhat surprisingly, even in the special case µ = H d | Γ the exact value of H d (Γ) is known only for a small number of IFS attractors, including the middle third Cantor set in R but not including the middle third Cantor dust in R 2 -we provide a more detailed commentary on the current state of knowledge regarding H d (Γ) in Remark 3.4.This means that, in practice, it is in general not possible to compute (1) and ( 2) even for f = 1!However, this does not compromise the utility of our quadrature rules for the main motivating application of this paper, namely the implementation of the Galerkin Hausdorff BEM for acoustic scattering by a fractal screen Γ screen in [9], and for similar possible applications to other differential and integral equations.This is because when one uses a BEM to solve a wave scattering problem, the physically relevant quantities such as the scattered wave field and its far-field pattern are unaffected by the choice of normalisation of the surface measure used in the BEM calculations.Working with the normalised measure H d (•) := H d (•)/H d (Γ screen ) in the BEM application leads to integrals of the form (1) and ( 2) with H d replaced by H d .Our quadrature rules apply mutatis mutandis to such integrals, requiring the value of H d (Γ), but this can be computed for any subcomponent Γ of Γ screen using the IFS structure, because H d has the same self-similarity scaling properties as H d on subcomponents of Γ screen (cf.(18) below), and H d (Γ screen ) is known (H d (Γ screen ) = 1 by the definition of H d ).For details see [9].

Preliminaries
We begin by reviewing a number of basic results about IFS attractors and integration on them, and introduce the notation and terminology we will use throughout the paper.

Hausdorff measure and dimension
For E ⊂ R n and α ≥ 0 we recall (e.g.[18,Section 3]) the definition of the Hausdorff α-measure of E, where the infimum is over all countable covers of E by sets where the supremum of the empty set is taken to be zero.Given 0 < d ≤ n, we call a non-empty closed set Γ ⊂ R n a d-set if there exist c 2 > c 1 > 0 such that

Iterated function systems
Throughout the paper we assume that Γ is the attractor of an iterated function system (IFS) of contracting similarities (see e.g.[22,18]), by which we mean a collection {s 1 , s 2 , . . ., s M }, for some M ∈ N, M ≥ 2, where, for each m = 1, . . ., M , s m : R n → R n satisfies, 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 .We denote by η m := (I − ρ m A m ) −1 δ m (I being the n × n identity matrix) the fixed point of the contracting similarity s m , i.e. the unique point η m ∈ R n such that s m (η m ) = η m .Saying that Γ is the attractor of the IFS means that Γ is the unique non-empty compact set satisfying where We shall also assume throughout that the open set condition (OSC) [22, Section 5.2] holds, meaning that there exists a non-empty bounded open set O ⊂ R n such that Then Γ is a d-set (e.g.[40,Thm. 4.7]), where d ∈ (0, n] is the unique solution of , with ρ = 0.41.(III) 5 s 1 , . . ., s 4 as in (I) and s 5 (x) = 1 27 x + ( 4 27 , 4  27 ).

Further assumptions on the IFS
We say that the IFS is homogeneous (as in, e.g., [15]) if ρ m = ρ ∈ (0, 1) for m = 1, . . ., M .In this case (10) becomes We say that the IFS is disjoint (as in, e.g., [5,Defn 7.1]) if which holds if and only if the open set O in the OSC can be taken such that Γ ⊂ O (e.g.[9]).
We say that the IFS is hull-disjoint if where Hull(E) denotes the convex hull of a set E ⊂ R n .Clearly R Γ,Hull ≤ R Γ , so hull-disjointness implies disjointness.But the converse is not true -see (II) and (III) in Figure 1 for counterexamples. (I) Figure 1: Four IFS attractors in R 2 with different degrees of "disjointness".The formulas of all contractions are in Table 1.Homogeneity, disjointness and hull-disjointness were defined in Subsection 2.

Vector index notation
Our quadrature rules will be based on partitioning Γ into self-similar subsets, via the IFS structure.
For m ∈ I , Γ m is itself the attractor of an IFS, namely This implies that the "elements" of the "Hausdorff-BEM" proposed in [9] are themselves IFS attractors, so the quadrature rules developed here can be used in the implementation of that method.
When considering singular integrands it will be important to estimate the distance between subsets of Γ.To that end, given m = n ∈ ∪ ∈N I we define (cf. the partial ordering in [22, Section 2.1]) (m 1 , . . ., m ) ∈ I , and for any In particular, taking f ≡ 1 gives The measure H d | Γ is just one member of a general class of finite measures on Γ for which similar scaling results apply.Given a collection (p 1 , . . ., p M ) of positive weights (or "probabilities") satisfying 0 < p m < 1, m = 1, . . ., M, and there exists (see, e.g., [22,Sections 4 & 5]) a Borel regular finite measure µ supported on Γ, unique up to normalisation, called an "invariant" [22] (also known as "balanced" [6] or "self-similar" [34]) measure associated to Γ and (p 1 , . . ., p M ), such that µ(A) = M m=1 p m µ(s −1 m (A)) for every measurable set A ⊂ R n .By [34,Thm. 2.1] the OSC implies that µ(s m (Γ) ∩ s m (Γ)) = 0 for each m = m , and as a consequence we find that for ∈ N, m = (m 1 , . . ., m ) ∈ I , and any µ-measurable function f , and We shall assume henceforth that µ is an invariant measure on Γ in this sense, for some collection of associated weights (p 1 , . . ., p M ).The case µ = H d | Γ corresponds to choosing p m = ρ d m for m = 1, . . ., M , with (19) holding by (10).

Partitioning Γ
To define our composite quadrature rules we need to specify an index set One approach is to choose subsets with a fixed level of refinement, i.e. to take I = I for some fixed ∈ N.However, in the case of a non-homogeneous IFS, subsets chosen in this way may  differ significantly in size.An alternative approach is to choose subsets with approximately equal diameter, taking I = L h (Γ) for some fixed h > 0, where with Γ (m 1 ,....,m −1 ) replaced by Γ when = 1 and L h (Γ) := {0} when h ≥ diam(Γ) (recall our convention that Γ 0 = Γ).See Figures 3 and 4 for illustrations of the decomposition L h (Γ) for the Koch snowflake and a non-homogeneous Cantor set.If Γ is homogeneous then for 0 < h ≤ diam(Γ) we have L h (Γ) = I , where 3 Barycentre rule for regular integrals We now present our composite barycentre rules for the evaluation of regular integrals of the form (1) and (2).

Single integrals
We first consider the single integral (1).Given a partitioning (22) of Γ into self-similar subsets, our quadrature nodes are the barycentres of the subsets (for an illustration in the case of the Koch snowflake see Figure 3), computed with respect to the measure µ, and the weights are the measures of the subsets.Definition 3.1 (Barycentre rule for single integrals).Let Γ and µ be as in Subsections 2.2 and 2.5, let I ⊂ ∪ ∈N I be an index set satisfying (22), and let f : Γ → C be continuous.Then for the approximation of the integral we define the barycentre rule where, for m= (m 1 , . . ., m l ) ∈ I, and The number of weights and nodes in this approximation is |I|.
Remark 3.2.While it always holds that x m ∈ Hull(Γ m ) (by the supporting hyperplane theorem), it does not in general hold that x m ∈ Γ m (the middle-third Cantor set provides a counterexample).
The weights w m in ( 27) can be computed using (21) as While the barycentres x m are defined a priori in terms of integrals with respect to µ, the following result shows how they can be computed using only information about the similarities (s 1 , . . ., s M ).This result coincides with the first step in the recursive procedure described in [26, 2.5.2] and [27,Section 2] for the calculation of moments of invariant measures.
Proposition 3.3.The barycentres x m defined by (26) can be evaluated as where x Γ := Γ x dµ(x) Γ dµ(x) = Γ x dµ(x) is the barycentre of Γ, which can be evaluated as where I is the n × n identity matrix and ρ m , A m and δ m , m = 1, . . ., M , are as in (7).
Proof.To prove (29) we first consider the case where s m = s m for some m ∈ {1, . . ., M }, for which by ( 7), ( 20) and ( 21) we have The general result follows by induction on the length of the vector index m.
To prove (30) we note that, by ( 26) and ( 21), Using (31), this can be written as from which (30) follows by matrix inversion.The invertibility of the matrix I − M m=1 p m ρ m A m follows from the fact that A m 2 = 1 and 0 < ρ m < 1 for each m = 1, . . ., M , so that, recalling (19), Remark 3.4.Evaluation of the quadrature weights w m defined in (27) requires knowledge of µ(Γ).If µ(Γ) is unknown then the quadrature rule (25) can only be evaluated up to the unknown factor µ(Γ).As mentioned in Section 1, even in the special case µ = H d | Γ the exact value of H d (Γ) is known only in certain special cases, and, to our knowledge, only for examples where d ≤ 1.In the Hausdorff BEM application that motivates this paper, this is unproblematic as one can simply work with an appropriately normalised measure (see the discussion in Section 1 and [9]).However, for completeness we comment briefly on the current state of knowledge regarding H d (Γ).The beststudied examples are Cantor-type sets in R (n = 1).Important early work in this area includes that of Marion [30,31] and Falconer [17], where it was proved that H d (Γ) = 1 for a large class of Cantor-type sets including the classical Cantor sets defined by (11) [17, Thm.1.14-1.15];for more recent related results see e.g.[4] and [43].For n > 1 it appears that the exact value of the Hausdorff measure of even the simplest IFS attractors is known only for d ≤ 1, see e.g.[41], where it is proved that for the Cantor dust defined in Example (I) of Table 1 [41,Cor. 1]; see also the earlier paper [42] where the case ρ = 1/4 was considered.For n > 1 and d > 1 it appears that only approximate results are available.For instance, when Γ ⊂ R 2 is the Sierpiński triangle, it is known that 0.77 ≤ H d (Γ) ≤ 0.81794 [33].One complication in the case n > 1 is that even if Γ is a Cartesian product of lower-dimensional IFS attractors, as is the case e.g. for the Cantor dust in Example (I) of Table 1, the measure of Γ cannot be computed as the product of lower-dimensional measures, since for sets Γ 1 and Γ 2 of dimension d 1 and d 2 respectively, in general we do not have

Double integrals
Double integrals of the form (2) can be treated by iterating the barycentre rule in the obvious way.Definition 3.5 (Barycentre rule for double integrals).Let Γ ⊂ R n and Γ ⊂ R n be as in Subsection 2.2 (possibly with different Hausdorff dimensions), and let µ and µ be invariant measures on Γ and Γ respectively, as in Subsection 2.5.Let I and I be index sets satisfying (22) for (Γ, µ) and (Γ , µ ) respectively, and let f : Γ × Γ → C be continuous.Then for the approximation of the iterated integral we define the iterated barycentre rule where, for m ∈ I, x m and w m are defined by (26) and (27), and, for m ∈ I , x m and w m are defined by the analogous formulas involving Γ and µ .

Error estimates
When the integrands are sufficiently smooth, error estimates for the quadrature rules in Definitions 3.1 and 3.5 can be derived by standard Taylor series arguments.The result for single integrals (Definition 3.1) is presented in Theorem 3.6 below.Before stating the theorem, we introduce some notation.Given a set E ⊂ R n and a function f : E → C we define If f is differentiable in an open set Ω ⊃ E, we denote its gradient by ∇f : Ω → C n , and define Note that we are allowing the possibility that L 0,E [f ] and L 1,E [f ] are infinite.If f is twice differentiable in Ω, we denote its Hessian by Hf : Ω → C n×n .For α ∈ N n 0 , D α denotes standard multi-index notation for partial derivatives.Finally, given h > 0 we define Γ Hull,h := Theorem 3.6.Let Γ and µ be as in Subsections 2.2 and 2.5.Let h > 0, and let Γ Hull,h be as in (33).Suppose that f : Γ Hull,h → C, and let Q h Γ denote the barycentre rule of Definition 3.1 with Proof.(i) Elementary estimation, combined with (22), gives (ii) The first bound follows from part (i) and the fact that, by the mean value theorem, L 0,Hull(Γm For the second bound we also apply the mean value theorem, noting for each m ∈ L h (Γ) and x ∈ Γ m there exists a point ξ m (x) on the segment between x and The first integral on the right-hand side vanishes by the definition of x m , with the result that, again using (22), (iii) By Taylor's theorem, for each m ∈ L h (Γ) and x ∈ Γ m there exists a point ξ m (x) on the segment between x and x m such that f and the result follows, by bounding diam(Γ m ) ≤ h, summing over m and using (22).
We now consider the double integral case (Definition 3.5).Higher order iterated integrals over the product of arbitrarily many IFS attractors can be analysed similarly, but are not considered here.
Theorem 3.7.Let Γ ⊂ R n and Γ ⊂ R n be as in Subsection 2.2, and let µ and µ be invariant measures on Γ and Γ respectively, as in Subsection 2.5.Let h > 0, and let Γ Hull,h be as in (33) and Γ Hull,h be as in (33) with Γ replaced by Γ .Suppose that f : Γ Hull,h × Γ Hull,h → C, and let denote the barycentre rule of Definition 3.5 with Proof.Follows similar arguments to those used to prove Theorem 3.6, in the setting of R n+n .The extra factors of √ 2 and 2 arise because diam Remark 3.9.The error bounds in Theorem 3.6 are written in terms of the discretization parameter h, which measures the diameter of the portion of Γ on which the integrand is approximated by a constant value.In practice it is useful to estimate the error in terms of the computational effort of the quadrature rule, writing the error bounds in terms of the number of quadrature points (and thus of the evaluations of the integrand) N := |L h (Γ)|.
We shall do this in the special case where Γ is homogeneous.Then L h (Γ) = I for ∈ N as in (24), so that N = M , and from (24) and (10), which gives ρ d = 1/M , we have We can substitute either of the last two expressions in place of h in the right-hand sides of the error bounds in Theorem 3.6, obtaining O(N −1/d ) and O(N −2/d ) convergence rates, with respect to increasing N .We observe that the lower the Hausdorff dimension d of Γ, the faster the convergence of the quadrature rule with respect to increasing N , all other factors being equal.
In the double-integral case, we can substitute either of the expressions , where N := |L h (Γ )|, in place of h in the error estimates of Theorem 3.7.In this case, the computational cost, measured as the number of integrand evaluations, is the product N N .
The barycentre rule (25) depends on the choice of the index set I. The choice I = L h (Γ) stipulated in Section 3.3, with L h (Γ) as in (23), allocates the quadrature nodes x m only according to the diameter of the associated subsets Γ m ⊂ Γ, and is motivated by the use of a Taylor-polynomial technique to bound the quadrature error.This choice is most effective for homogeneous IFSs and Hausdorff measures, where the fraction of the measure associated to each quadrature node is the same.For general IFSs and invariant measures, we expect that more sophisticated choices of the index set I, taking into account the measure of the subsets in the induced partition as well as their diameter, may be more efficient.However, we leave further discussion of this issue to future work.

Evaluation of singular integrals
We now turn to the case where the integrands in (1) and ( 2) are singular.In this section, we show how two classes of singular integrals involving the function Φ t (defined in (4)) can be expressed in terms of regular integrals, using the scaling properties from Subsection 2.5 and certain homogeneity properties of Φ t .This will allow us to derive and analyse quadrature rules for the singular integrals based on the barycentre rule described in the previous section.

Homogeneity and bounds on derivatives of Φ t
Key to our analysis will be the fact that, for any x = y ∈ R n and ρ > 0, We note also that Φ t (x, y) is smooth as a function of both x and y, away from the diagonal x = y, on which it is singular.More precisely, we shall need the following result concerning derivatives of Φ t .
Here the assumption that the multi-index α ∈ N 2n 0 means that α could correspond to differentiation with respect to the components of either x or y, and From |x i − y i | ≤ r(x, y) and the values of the following partial derivatives for i, j = 1, . . ., n: Inserting these in (36) gives, for all α, β ∈ N 2n 0 with |β| = 1 and |α| = 2, Now recall that Φ t (x, y) = Φt (r(x, y)), where Elementary computations show that Inserting these results into (37) (with F = Φt ) gives the claimed result.

Single integrals
We consider first the evaluation of the single integral for η ∈ R n .When η ∈ Γ the integral is regular, and combining Theorem 3.6(iii) and Lemma 4.1 gives the following error estimate for the barycentre rule.
When η ∈ Γ the integral is singular, but convergent for sufficiently small t.In the case µ = H d | Γ we have convergence for 0 ≤ t < d ≤ n, for any η ∈ Γ (see Corollary A.2).For a general invariant measure the situation is more complicated and the integrability threshold depends on η (see the discussion in Section A.2).In the special case where η is the fixed point of one of the contracting similarities s m defining Γ, we have convergence for 0 ≤ t < t m (see Lemma A.3), where Furthermore, in this case the singular integral can be written in terms of regular integrals, as the following result shows.We remind the reader that if (This holds, for instance if Γ is disjoint in the sense of (13).)Then the singular integral (39) is finite for 0 ≤ t < t m , where t m is defined in (40), and it can be represented in terms of regular integrals, as: Proof.The integrability result is proved in Lemma A.3.To prove the claimed decomposition, we first split the integral, writing Focusing on the singular term, by (20) and the fact that s m (η m ) = η m we can write using the fact that Φ t is translation and rotation invariant.Then, applying (35) gives Substituting ( 42) into (41) and solving for I Γ [Φ t (•, η m )], we obtain the result.
Theorem 4.3 can be combined with any quadrature rule capable of evaluating the regular integrals , m = m, to produce a quadrature rule for evaluating the singular integral (39) when η = η m .In particular, given h > 0, applying the barycentre rule of Definition 3.1 with Γ replaced by Γ m and I = L h (Γ m ), for each m = m, produces the following quadrature rule: This quadrature formula could be used for the implementation of a collocation-type discretisation of the integral equations in [9], with the collocation nodes chosen as fixed points of the self-similar subsets of the IFS attractor used as the BEM elements in [9].This discretisation is not investigated in [9].
Proof.For 0 ≤ t < t m we have from Theorem 4.3 and equation ( 43) that and the result follows by applying Proposition 4.2 to each term in the sum, and recalling (22) to see that The

Double integrals
We now consider the evaluation of the double integral When Γ and Γ are disjoint the integral is regular, and combining Theorem 3.7(iii) and Lemma 4.1 gives the following error estimate for the barycentre rule.
Proposition 4.5.Let Γ, Γ ⊂ R n and µ, µ be as in Subsections 2.2 and 2.5.Let h > 0 and let Γ Hull,h be as in (33) and Γ Hull,h be as in (33) with Γ replaced by Γ .Suppose that Γ Hull,h ∩Γ Hull,h = ∅.Then When Γ and Γ are not disjoint (45) is a singular integral, which converges only for sufficiently small t.Suppose for simplicity that Γ = Γ.Then in the case µ = µ = H d | Γ we have convergence for 0 ≤ t < d ≤ n (see Corollary A.2).For a more general pair of invariant measures µ and µ on Γ, with respective (possibly different) weights/probabilities (p 1 , . . ., p M ) and (p 1 , . . ., p M ), if Γ is disjoint then the integral converges for 0 ≤ t < t * (see Lemma A.4), where t * is the unique positive solution of In the disjoint case the singular integral (when it converges) can be written purely in terms of regular integrals, as the following result shows.This was noted previously for the case of Cantor sets in e.g.[7].We remind the reader that if Theorem 4.6.Let Γ ⊂ R n be as in Subsection 2.2, and let µ, µ be as in Subsection 2.5.Suppose that Γ is disjoint in the sense of (13).Then the singular double integral (45) with Γ = Γ converges for 0 ≤ t < t * , where t * is the unique positive solution of (46), and it can be represented in terms of regular integrals, as: Proof.The integrability result is proved in Lemma A.4.To prove the claimed decomposition, as in the single integral case we begin by splitting the integral, writing By applying (20) (for both µ and µ ), (35), and the fact that Φ t is translation and rotation invariant, the singular integrals in (47) can be written as Substituting this expression into (47) and solving for I Γ,Γ [Φ t ] gives the claimed result.
By combining Theorem 4.6 with a suitable quadrature rule for evaluating the regular integrals I Γm,Γ m [Φ t ], m = m, we can obtain a quadrature rule for evaluating the singular integral (45).In particular, given h > 0, applying the barycentre rule of Definition 3.5 with Γ, Γ replaced by Γ m , Γ m and with I = L h (Γ m ) and I = L h (Γ m ) for each pair (m, m ) such that m = m produces the following quadrature rule: Corollary 4.7.Let Γ, µ, µ and t be as in Theorem 4.6.Let 0 < h < diam(Γ), and assume that Hull(Γ m ) ∩ Hull(Γ m ) = ∅ for all m, m ∈ L h (Γ) such that m 1 = m 1 .Then the quadrature rule defined by (48) for the integral (45) satisfies the error estimate , where Proof.The proof is similar to that of Corollary 4.4.For 0 ≤ t < t * we have and the result follows by applying Proposition 4.5 to each term in the sum.
Remark 4.8.So far we have introduced four different parameters quantifying the distance between self-similar subsets of an IFS attractor Γ: • R Γ in (13), which measures the minimum distance between level-1 subsets, • R Γ,Hull in (14), which measures the minimum distance between the convex hulls of level-1 subsets, • R m,h in (44), which measures the minimum distance between a fixed point η m of s m and the convex hulls of subsets of Γ \ Γ m of diameter approximately h, • R Γ,Hull,h in (49), which measures the minimum distance between the convex hulls of pairs of subsets, taken from different level-1 subsets, of approximate diameter h.
Recall that R m,h and R Γ,Hull,h are defined only for 0 < h < diam(Γ).They satisfy the inequalities In particular, Corollaries 4.4 and 4.7 show how R m,h and R Γ,Hull,h quantify the expected deterioration of the quadrature accuracy due to the vicinity of the integrand singularity, for single and double integrals (I Γ [Φ t (•, η m )] and I Γ,Γ [Φ t ]), respectively.
Table 2 shows the values of these parameters for the four examples in Figure 1.The values of R m,h and R Γ,Hull,h are valid for all sufficiently small h (e.g.R 5,h = 0 for (III) and h ≥ √ 2/3).
We are not aware of any IFS attractor satisfying the open set condition with R m,h = 0 for some m, i.e. with η m ∈ Hull(Γ m ) for m ∈ L h (Γ) and The values of the separation parameters described in Remark 4.8 for the IFS attractors of Figure 1, for sufficiently small h.

Relative errors and dependence on N for homogeneous IFSs and Hausdorff measure
In this section we show how the error bounds we derived for the singular integrals in Subsection 4.2 and Subsection 4.3 can be written in terms of the quantity N := |L h (Γ)|, as was discussed for the regular integrals of Subsection 3.3 in Remark 3.9.As well as allowing us to determine the dependence of our error bounds on the computational cost of the quadrature rules, this also allows us to clarify the limiting behaviour of our bounds as (d − t) → 0. To this purpose, we now consider relative errors, and restrict our attention in this section to homogeneous IFSs and the case µ For a homogeneous IFS, we have N = |I | = M for as in (24).The number of evaluations of the integrand Φ t required for the computation of the quadrature formulas is (M − 1)M −1 = M −1 M N for the single-integral formula ( 43) and (M − 1)M 2 −1 = M −1 M N 2 for the double-integral formula (48).Let us consider first the case t > 0. From (71) we have the lower bounds: Recall that c1 > 0, defined in (6), is an intrinsic parameter of the d-set Γ, independent of its characterization as an IFS attractor.Then, using that for a homogeneous IFS and the case µ (12), Corollaries 4.4 and 4.7 imply the following relative error estimates: where To bound E t we first note that, for 0 < z, ρ < 1, we have ρ z ≤ 1 − z(1 − ρ) (by comparison of an affine and a convex function of z that coincide for z = 0 and z = 1), so that Moreover, from Then, using also d = log M | log ρ| , ρ = M −1/d and h ≤ diam(Γ) ρ N 1/d (see (34)), we can bound Combined with (50), this reveals the dependence of the relative error on the computational cost, through the parameter N ; specifically, the errors are O(N −2/d ) as N → ∞.The above bound can also be written in terms of the refinement level using N M = M −1 , giving exponential convergence at the rate M −2 /d as → ∞.
While the bounds on the absolute errors in Corollaries 4.4 and 4.7 blow up in the limit t d (with N fixed), the bounds (50) and (53) show that the corresponding relative errors are bounded in this limit, because the integrals being approximated also blow up at the same rate.Similarly, for a sequence of IFSs with d t > 0 (and with M constant, c1 and R Γ,Hull uniformly bounded away from zero, and H d (Γ) and diam(Γ) uniformly bounded above), the absolute errors blow up while the relative errors are uniformly bounded.Furthermore, the same is true (again for fixed N ) in the case where d 0 and t 0 with d > t, since the algebraic growth of the 1 d term in (53) is controlled by the exponential decay of the factor ( N M ) −2/d (provided ≥ 2).In the case t = 0 the logarithmic function Φ 0 changes sign, so it is not in general possible to bound its integrals from below.Thus we assume that diam(Γ) ≤ 1.Under this assumption, for all η ∈ Γ, (71) gives Proceeding as above and using a 0 = 2, log diam(Γ) ≤ 0, and 1 1−ρ d = M M −1 ≤ 2, the bound (50) on the relative error extends to the case t = 0 with , which tends to zero as d t = 0.

Application to Galerkin Hausdorff BEM for acoustic scattering
We now apply our previous results to derive and analyse quadrature rules for the evaluation of where Γ, Γ ⊂ R n are as in Subsection 2.2 and Φ(x, y) is the fundamental solution of the Helmholtz equation in R n+1 , defined in (3).As (54) suggests, our focus on this section is on the case As explained in Section 1, integrals of the form (54) arise as the elements of the Galerkin matrix in the "Hausdorff BEM" described in [9], for acoustic scattering by fractal screens.We first consider (54) in the non-singular case where Γ and Γ are disjoint, corresponding to the off-diagonal matrix entries in [9].Our quadrature rule in this case is the composite barycentre rule, and the main result is Proposition 5.2.We then consider (54) in the singular case where Γ = Γ , corresponding to the diagonal matrix entries in [9].Our quadrature rule for this case is defined in (60) and (61), and the main result is Theorem 5.7.
Before proceeding with the analysis we note the following regularity estimate on Φ.Here and henceforth a b means a ≤ Cb for some constant C > 0, independent of Γ, Γ , h and k, which may change from occurrence to occurrence.Proof.The proof is analogous to that of Lemma 4.1.We first note that Φ(x, y) = Φ(r(x, y)), where and by standard calculations (e.g.[1, 10.6.2-3])we find that 1 (kr), (ikr − 1)e ikr 4πr 2 , and Φ (r) = 1 kr H (1) Inserting these results into (37) (with F = Φ) gives the claimed result, after application of the following standard bounds (which follow from results in [1, Section 10]): We now consider (54) in the non-singular case where Γ and Γ are disjoint.The following result follows from Theorem 3.7(iii) and Lemma 5.1, and the fact that (1 + (kz) n/2+1 )/z n+1 is a positive decreasing function on (0, ∞) for n = 1, 2.
We now turn to the singular case where Γ = Γ and (54) becomes For fixed k > 0, we have (by [1, (10.8.2)] in the case n = 1) that where Furthermore, the function Φ is continuous across x = y (in fact, Lipschitz continuous), with (see [1, 10.8.2] for the case n = 1) This motivates a singularity-subtraction approach for evaluating (57), using the splitting integrating Φ n−1 using the quadrature rules from Section 4, and Φ * using the composite barycentre rule.However, the application of ( 59) is complicated by the fact that while (59) is designed to deal efficiently with the singular behaviour (58), it is not well adapted to the oscillatory behaviour of Φ as k|x − y| → ∞.To deal with this systematically, we introduce a parameter c osc > 0, proportional to the maximum number of wavelengths there can be across Γ for us to consider the integral non-oscillatory.We then define our quadrature rule for (57) differently depending on whether k diam(Γ) ≤ c osc or k diam(Γ) > c osc .Note that the non-oscillatory regime k diam(Γ) ≤ c osc is the one relevant for the BEM application in [9], since the diameter of the BEM elements in [9] needs to be small compared to the wavelength in order to achieve acceptable approximation error.
When k diam(Γ) > c osc we first partition Γ into subsets of diameter at most h * := c osc /k before applying (59) only to the singular terms in the resulting decomposition of (57), approximating where Q h Γm,Γm,n−1 is defined as in (48) with Γ replaced by Γ m .
For the error analysis of (60), we recall that an error estimate for Q h Γ,Γ,n−1 was presented in Corollary 4.7, so it remains to derive an error estimate for because while Φ * is Lipschitz continuous across x = y, its derivative is not Lipschitz.An O(h 2 ) estimate (matching that for Q h Γ,Γ,n−1 provided by Corollary 4.7) can be obtained via a first-principles analysis, which we present in Proposition 5.5.We then apply this result to give a full error analysis of both ( 60) and (61) in Theorem 5.7.Our analysis is restricted to the case of a homogeneous IFS, but we expect that with further non-trivial work a similar analysis could be carried out for the non-homogeneous case -see the discussion before Theorem 5.11, where a weaker O(h) estimate is proved for the non-homogeneous case.
In what follows we shall also make use of the fact that and that, in the case µ = H d | Γ , for m ∈ L h (Γ) it follows from ( 18) that We now consider the approximation of Suppose that Γ is homogeneous in the sense of Subsection 2.3, with contraction factor ρ ∈ (0, 1).Suppose also that Γ is hull-disjoint in the sense of (14).Let k > 0 and suppose that k diam(Γ) ≤ c osc .Let 0 < h ≤ diam(Γ).Then where the constant implied in depends only on c osc and Proof.We first note that For the analysis of S 1 we note that for m ∈ L h (Γ) By (62) we have that, when kr ≤ c osc , Since both z 2 (1 + | log z|) and z are increasing functions of z ∈ (0, ∞), applying a uniform upper bound on the integrand in (66) (with r ≤ h), then summing over m and using ( 22) and (64) gives For the analysis of S 2 we note that if m ∈ L h (Γ) then (Γ m ) Hull,h = Hull(Γ m ), so that for m = m we have, by (16) and the assumption that Γ is hull-disjoint, that where * (m, m ) was defined in (15).Now, since Γ is homogeneous we have that L h (Γ) = I , where satisfies (24).Also, In the case n = 2, given m ∈ L h (Γ) = I , by Lemma 5.4, Theorem 3.7(iii) and Lemma 2.1 we have where we used the fact that M ρ d = 1 (by (10)) and d > n − 1 = 1 (an assumption of the theorem) to deduce that ρM > 1, which means we can take the summation limit to infinity in the geometric series.Finally, summing over m and bounding In the case n = 1, given m ∈ L h (Γ) = I , by Lemma 5.4, Theorem 3.7(iii) and Lemma 2.1 we have and then summing over m, and noting that 1 ≤ 1/d and that 12) and the fact that M − 1 > log M , gives Combining the estimates for S 1 and S 2 then gives The quantity in braces here is h-dependent, but can be bounded uniformly in h to give (65).
Remark 5.6.The error estimate of Proposition 5.5 for n = 2 blows up as ρ 1 M , equivalently as d 1, assuming H d (Γ), diam(Γ), M and N are fixed and R Γ,Hull is bounded away from zero, because of the factor ρM ρM −1 = For example, for a family of Cantor dusts as in Figure 1(I), parametrised by ρ and with M = 4, this corresponds to the limit ρ 1 4 .Differently from the setting in Subsection 4.4, the relative error is also predicted to blow up in this limit, because the integral being approximated in this case is bounded, since Φ * ∈ L ∞ (Γ × Γ).By contrast, for n = 1 the estimate in Proposition 5.5 tends to zero as d 0, because for N > M the algebraic growth of the 1/d term in c is beaten by the exponential decay of the (N/M ) −2/d factor.
Finally, we can state and prove our main result for the approximation of I Γ,Γ [Φ].
Theorem 5.7.Let Γ ⊂ R n be as in Subsection 2.2, with d = dim H (Γ) > n − 1. Suppose that Γ is homogeneous in the sense of Subsection 2.3, with contraction factor ρ ∈ (0, 1).Suppose also that Γ is hull-disjoint in the sense of (14).Let k > 0. For the approximation of the integral (57) define the quadrature rule , where the constant implied in depends only on c osc , and Proof.By redefining h * := min{c osc /k, diam(Γ)}, (60) can be viewed as a special case of (61), since with h * = diam(Γ) we have L h * (Γ) = {0}, in which case the first sum in (61) reduces to (60) and the second sum is absent.We therefore present the proof of (61) and specialise to (60) at the end.
By the triangle inequality and the splitting Φ = C n Φ n−1 + Φ * , we have: To bound T 1 above, we first note that by Corollary 4.7 (with Γ replaced by Γ m and with t = n − 1 for n = 1, 2) for any m ∈ L h * (Γ) we have Let (h * ) be given by ( 24) with h replaced by h * .Noting that ρ (h * ) > ρh * / diam(Γ) we have that Hence, using (22) with I = L h * (Γ) and (64) with h replaced by h * , For T 2 , by Proposition 5.5 we have for any m ∈ L h * (Γ) that Using (67), and for n = 1 the facts that | log ρ| log M | log ρ| ≤ 1 ρ , and that kR Γm,Hull ≤ k diam Γ m ≤ kh * ≤ c osc , so that we can apply (63) with x 1 = ρkh * R Γ,Hull / diam(Γ m ), x 2 = kR Γm,Hull and x 3 = c osc , to find that For T 3 , by Proposition 5.2 we have for m, m ∈ L h * (Γ) with m = m that and by ( 16) and (67) it holds that R Γ,Hull , so that, since (1 + (kz) n/2+1 )/z n+1 is positive and decreasing on (0, ∞), and Recalling our redefinition of h * := min{c osc /k, diam(Γ)}, the result for the case k diam(Γ) > c osc then follows by combining the estimates for T 1 , T 2 and T 3 , and the result for the case k diam(Γ) ≤ c osc follows by noting that in that case T 3 is absent.We describe in more detail the four possible cases.
• For k diam(Γ) ≤ c osc (and thus h * = diam(Γ)) and n = 1 we have, since M/(M − 1) ≤ 2, , where the final bound follows from the fact that k ≤ cosc and if ρkR Γ,Hull > 1 then • For k diam(Γ) ≤ c osc and n = 2 we have, again using that k 1/R Γ,Hull , • For k diam(Γ) > c osc (and thus h * = c osc /k) and n = 1, • Finally, for k diam(Γ) > c osc and n = 2, Remark 5.8 (Number of function evaluations).For a homogeneous IFS, recall that N = |L h (Γ)| = |I | = M for as in (24).A priori, the quadrature rule (60) requires (M −1)M 2 −1 = M −1 M N 2 evaluations of Φ n−1 (see Subsection 4.4) and M 2 = N 2 evaluations of Φ * .If * is defined by (24) with h replaced by h * , i.e. * is the level of the partition whose elements have diameter approximately h * = c osc /k, then the quadrature rule (61) requires However, the number of function evaluations can be reduced (by a factor of a half in the limit h → 0) by exploiting the symmetry of Φ, Φ * , Φ n−1 , all of which satisfy f (x, y) = f (y, x).Remark 5.9 (Limit behaviour for d n − 1).We consider the behaviour of the estimates of Theorem 5.7 as d n − 1, assuming H d (Γ), diam(Γ), M and N are fixed, and R Γ,Hull is bounded away from zero.This limit corresponds to ρ 0 for n = 1 and ρ 1/M for n = 2.For n = 1 the absolute error tends to zero like O((N/M 2 ) −2/d ) as d 0, provided N > M 2 .Since in this limit the integral is dominated by the contribution from the singular function Φ 0 , the integral of which grows like 1/d as d 0 (as shown in Subsection 4.4), the relative error tends to zero like For n = 2 the estimate for the absolute error grows as d 1, being asymptotically proportional to 1/(1 − M 1/d−1 ) ∼ 1/((d − 1) log M ) as d 1.However, the relative error is bounded in this limit because, again, the contribution of the singular function (Φ 1 in this case) grows in proportion to 1/(d − 1) as d 1 (as shown in Subsection 4.4).We validate these statements numerically in Section 6 below.
Remark 5.10 (Behaviour for vanishing distance between subsets).The estimates of Theorem 5.7 also blow up in the limit R Γ,Hull 0, specifically like R −(n+1) Γ,Hull .However, our numerical investigations in Section 6 suggest that, at least in certain cases, this is overly pessimistic.Proposition 5.5 and Theorem 5.7 are stated only for the case of a homogeneous IFS.We expect that with non-trivial further work one should be able to extend the O(h 2 ) estimates in these results to the non-homogeneous case, but we defer this to future studies.The main difficulty is obtaining sharp estimates for the sum S 2 in the proof of Proposition 5.5.While we cannot currently prove O(h 2 ) estimates for the non-homogeneous case, we can at least prove weaker O(h) estimates.The following is the O(h) analogue of Theorem 5.7.The proof, which we do not provide here, essentially follows that of Theorem 5.7, but applies lower order estimates, and estimates the sum S 2 in the proof of Proposition 5.5 more simply (but less sharply) using a uniform bound over all summands.Theorem 5.11.Let Γ ⊂ R n satisfy the assumptions of Theorem 5.7, except that we no longer assume Γ is homogeneous in the sense of Subsection 2.3.Then where the constant implied by depends only on c osc and where ρ min = min m∈{1,...,M } ρ m .

Numerical experiments
In this section we present numerical results complementing our theoretical analysis.The code used for our numerical experiments is available at https://github.com/AndrewGibbs/IFSintegrals,where we provide a Julia-based [8] implementation of all the quadrature rules presented in this paper.Within this repository, the interactive notebook QuadratureExample.jpynb provides an overview of the main steps in our algorithm and examples of usage.The pseudocode for a simple recursive implementation of the quadrature rule Q Γ [f ] (25) for regular single integrals is shown in Algorithms 1-2 below.
Estimation of diam(Γ) is a key step in Algorithm 1.In the numerical experiments which follow, diam(Γ) can be derived analytically.But for more general cases where an analytic derivation is not possible, our implementation estimates diam(Γ) using an algorithm that follows from [14,Proposition 6].
We shall focus mainly on the validation of the quadrature rule Q h Γ,Γ,Φ defined by (60) and (61) for the calculation of the singular double integral I Γ,Γ [Φ] defined in (57), since this is the most challenging integral we consider in the paper, and since it is important for the Hausdorff BEM application of [9].However, our numerical results for Q h Γ,Γ,Φ also implicitly validate the quadrature rule Q h Γ,Γ,t defined in (48) for the integration of the singular function Φ t , and, more fundamentally, the barycentre rule of Definition 3.5 for regular integrands.
The definition of Q h Γ,Γ,Φ in Section 5 involves a parameter c osc > 0, which governs whether the integral I Γ,Γ [Φ] is treated as non-oscillatory (k diam(Γ) ≤ c osc ), in which case (60) is applied, or oscillatory (k diam(Γ) > c osc ), in which case (61) is used.If c osc is too small, accuracy will deteriorate because the singularity will not be properly captured, while if c osc is too large, accuracy will also deteriorate because the splitting (59) is being used outside of its range of applicability.Our experience, following a detailed numerical investigation, suggests that a value of c osc = 2π gives acceptable performance across all the examples we considered, and this is the value of c osc we use throughout this section.This means we classify the integral I Γ,Γ [Φ] to be oscillatory (and use (61) rather than (60)) whenever the diameter of Γ is larger than one wavelength.5.9.
Vanishing separation limit.Next we consider the behaviour of our quadrature rules as the parameter R Γ,Hull tends to zero.In Figure 7(a) we show absolute errors for Q h Γ,Γ,Φ with with k = 5 at three values of N , for Γ ⊂ R a Cantor set, defined by (11), with ρ = (1 − R Γ )/2, where R Γ = 1 − 2ρ ∈ {0.1, 0.01, 0.001, 0.0001, 0.00001}, and µ = H d | Γ .The reference solution is as for Figure 5.In Remark 5.10 we observed that as R Γ,Hull → 0 our theory predicts blow-up of the error like R n+1 Γ,Hull , i.e. like R 2 Γ in this case.However, the numerical results suggest that, at least in this case, the theoretical prediction is overly pessimistic, since the error appears to be bounded as R Γ,Hull → 0. In fact, the integral for the non-disjoint case ρ = 1/2 (so Γ = [0, 1], d = 1 and R Γ = 0) can be computed using our method, and the corresponding errors appear to follow the same N −2/d behaviour (in this case, N −2/d = N −2 since d = 1) with respect to increasing N as for the case 0 < ρ < 1/2 (see the dashed lines in the figure).To further investigate the non-disjoint case ρ = 1/2 (with Γ = [0, 1], d = 1 and R Γ = 0), in Figure 7(b) we plot the absolute error in the quadrature rule

Q h
Γ,Γ,0 ≈ I Γ,Γ [Φ 0 ] for this case, for which we have the exact result where we used the fact that H 1 coincides with the Lebesgue measure on R.
Vanishing separation limit.The results in Figure 7 suggest that our assumption that Γ should be hull-disjoint, or even disjoint at all, may not be necessary in Theorem 5.7.To investigate this further we compute Q h Γ,Γ,Φ for two non-hull-disjoint examples and µ = H d | Γ : example (II) in Table 1, which is disjoint but not hull-disjoint, and example (IV) in Table 1, which is not disjoint.Both attractors are shown in Figure 1.Absolute errors (scaled by H d (Γ) 2 ) for these cases for k = 2 and a range of h values are presented in Figure 8(a), and for both examples it seems we obtain O(h 2 ) convergence, even though our theoretical error analysis does not cover these cases.Results for the middle-third Cantor dust (example (I) in Table 1) are included in the same figure for reference.
In Figure 8(a) we also include results for a hull-disjoint but non-homogeneous IFS with M = 4 and non-homogeneous hull-disjoint cases our current analysis only provides an O(h) convergence result (see Theorem 5.11).But the results in Figure 8 for the IFS (70) suggest that, at least in this case, our analysis may not be sharp in this respect, since we seem to obtain O(h 2 ) convergence in practice.We leave further investigation of this to future work.
For this non-Hausdorff invariant measure it is instructive to compare how the two quadrature rules deal with the non-uniform way in which the mass of the measure is distributed across Γ.In Figure 9 we plot the nodes and weights for the barycentre rule for the case h = 0.01 (which corresponds to N = 35839) alongside those for one realisation of the chaos game rule with the same N = 35839.Each node is represented by a small dot, coloured according to the corresponding quadrature weight.
For the chaos game, the weights are uniform, all being equal to 1/N , but the nodes are distributed non-uniformly, being concentrated in the regions where the measure has greatest mass.By contrast, the nodes for the barycentre rule are distributed approximately uniformly, but the weights vary according to the measure.
In Figure 10  Here the reference value I ref Γ [f ] was computed using the barycentre rule with N = 2876335 (which corresponds to h = 10 −3 ).For the chaos game rule, the plots show both the individual errors for each of 1000 random realisations (thin blue lines) and the average of these individual errors (thick red line), which represents an approximation to the statistical expectation of the error for the chaos game rule.The error for the barycentre rule clearly decays like h2 ∼ N −2/d = N −1 , consistently with Theorem 3.6(iii) and Remark 3.9, even if the latter does not directly apply to non-homogeneous IFSs, 1 while the average error for the chaos game rule decays like N −1/2 , as one expects from a Monte-Carlo-type stochastic method.So for this problem the barycentre rule clearly outperforms the chaos game rule.Comparing the convergence rates h 2 ∼ N −2/d (for the homogeneous case and the present one) and N −1/2 , we expect that the advantage provided by the barycentre rule over the chaos game rule is even stronger for the lower-dimensional attractors considered in the previous numerical experiments.
Figure 10: Convergence of the barycentre rule and the chaos game rule for integration of a smooth integrand with respect to a non-Hausdorff invariant measure on the Koch snowflake.For the chaos game rule we show both the individual errors for each of 1000 random realisations (thin blue lines) and the average of these individual errors (thick red line).
For higher-dimensional problems we might expect the stochastic approach to become more competitive.To investigate this we consider the case where Γ is a high-dimensional Cantor dust.Specifically, we take Γ = (Γ ρ ) 6 ⊂ [0, 1] 6 ⊂ R 6 , a 6-fold Cartesian product of a homogeneous dyadic Cantor set Γ ρ ⊂ R (with contractions s 1 (x) = ρx and s 2 (x) = 1 − ρ + ρx for some 0 < ρ < 1/2) with itself, so that Γ is the attractor of a homogeneous, disjoint IFS with M = 2 6 = 64 and d = dim H (Γ) = log M/ log(1/ρ) = 6 log 2/ log(1/ρ).Figure 11  For all three values of the dimension d, the barycentre rule converges like h 2 ∼ N −2/d , as predicted by Theorem 3.6(iii), while the average error for the chaos game rule converges consistently like N −1/2 .Hence for d = 3 the barycentre rule converges faster; for d = 4 the two methods converge at the same rate, and for d = 5 the chaos game rule converges faster (although we note that in this particular experiment the errors for the barycentre rule were smaller than the expected errors for the chaos game even for d = 5).

A Integrability of singular functions with respect to invariant measures
In this appendix we collect some results concerning the integrability of singular functions with respect to invariant measures of the type defined in Section 2.5.In particular, we study for which t ≥ 0 the single integral I Γ [Φ t (•, η)] (defined in (39)) and the double integral I Γ,Γ [Φ t ] (defined in (45)) are finite. A
From the above observations it follows that if dim loc µ(η) exists then t µ (η) = dim loc µ(η).By [20, Thm 2] we have that if Γ is disjoint (see also [39,Thm 7.4] for the general case) then dim loc µ(η) = t a.As a consequence, t µ (η) = t a.e. , for µ-a.e.η ∈ Γ.But t µ (η) is not in general equal to t a.e. on the whole of Γ.By [18,Thm 17.4], if Γ is disjoint then for all points η ∈ Γ where the local dimension exists (which we know from the above is µ-a.e.) we have: Suppose that η m ∈ Γ m for any m ∈ {1, . . ., M }, m = m.(This holds, for instance if Γ is disjoint in the sense of (13).)Then there exist C 2 > C 1 > 0 such that, for all sufficiently small r > 0, where t m := log p m log ρ m .
Furthermore, H d (s m (Γ)) ∩ s m (Γ)) = 0 for m = m , a property known as self-similarity [22, 5.1(4)(ii)].The best-known example of an IFS attractor is the Cantor set Γ ⊂ R, defined by Figure 1: Four IFS attractors in R 2 with different degrees of "disjointness".The formulas of all contractions are in Table 1.Homogeneity, disjointness and hull-disjointness were defined in Subsection 2.3.(I) Top left: a homogeneous hull-disjoint IFS (a Cantor dust with ρ = 1/3).(II) Top right: a homogeneous disjoint IFS that is not hull-disjoint and such that η m / ∈ Hull(Γ m ) for all m = m (recall that η m is the fixed point of s m ).The fixed points are the three vertices and the centre of Hull(Γ).This condition on the fixed points is relevant e.g. in Corollary 4.7.(III) Bottom left: a non-homogeneous disjoint IFS that is not hull-disjoint and such that η m ∈ Hull(Γ m ) for some m = m.(In particular, s 5 (Γ) ⊂ Hull(s 3 (Γ)).)(IV) Bottom right: a homogeneous IFS attractor that is not disjoint (the Vicsek fractal).All these examples satisfy the OSC (9): (I) and (IV) for the open square O = (0, 1) 2 , (II) and (III) (and (I) again) for a neighbourhood O = Γ + B (0) with sufficiently small .

Figure 3 :
Figure 3: Examples of the partitioning (22) corresponding to the index set I = L h (Γ) defined by (23) in the case where Γ ⊂ R 2 is the Koch snowflake.The barycentres (defined by (26)) of the self-similar subsets in the partitioning are indicated with red crosses.The Koch snowflake can be written as a non-homogeneous IFS attractor with M = 7, s m (x, y) = 1 3 (x, y) + 2 3 (cos α m , sin α m ) with α m = (2m−1)π 6

Remark 3 . 8 .
The fact that x m is the barycentre of Γ m only enters in the proof of the O(h 2 ) estimates in Theorems 3.6 and 3.7.The O(h) estimates remain true if the barycentre x m is replaced by any other point of Hull(Γ m ).

Corollary 4 . 4 .
Let Γ, m, η m and t m be as in Theorem 4.3.Let 0 < h < diam(Γ), and suppose that η m ∈ Hull(Γ m ) for each m ∈ L h (Γ) such that m 1 = m.Then the quadrature rule defined by(43) for the integral(39) with η = η m satisfies the error estimate

Figure 5 :
Figure 5: Convergence of Q h Γ,Γ,Φ for a collection of Cantor sets with parameters ρ approaching 0.

Figure 9 :
Figure 9: Visual representation of barycentre rule (left) and chaos game quadrature (right) on the Koch snowflake, for a randomly chosen invariant measure.Each quadrature node is represented by a small dot, coloured according to the corresponding quadrature weight.
we plot the relative quadrature errors|Q Γ [f ]−I ref Γ [f ]|/|I ref Γ [f ]| and |Q CG Γ [f ]−I ref Γ [f ]|/|I ref Γ [f ]| against the number N of point evaluations of f , for a range of values of N between 463 and 320503.

Figure 11 :
Figure 11: Convergence of the barycentre rule and the chaos game rule for the approximation of a single smooth integral on three Cantor dusts Γ ⊂ R 6 with d = dim H (Γ) = 3, 4, 5.
min m=1,...,M log p m log ρ m ≤ t µ (η) = dim loc µ(η) ≤ max m=1,...,M log p m log ρ m .The upper and the lower bounds coincide, i.e. log p m / log ρ m is the same for all m = 1, . . ., M , if and only if p m = ρ d m , i.e. µ = H d | Γ .Moreover, the extremal values are attained, as the following lemma shows.Lemma A.3.Let Γ and µ be as in Subsections 2.2 and 2.5.Fix m ∈ {1, . . ., M } and let η m denote the fixed point of the contracting similarity s m , i.e. the unique point η m ∈ Γ such that s m (η m ) = η m .
. A special case is where µ = H d | Γ and µ = H d | Γ , where H d and H d are Hausdorff measures, with d, d denoting the Hausdorff dimensions of Γ and Γ (see Subsection 2.1), and | Γ denotes the restriction to Γ in the sense that

Table 1 :
The contractions corresponding to the IFS attractors in Figure1.
separation parameter R m,h introduced in (44) is used only in the statement of Corollary 4.4, and is compared to other related parameters in Remark 4.8.The relative error and the behaviour of the bound in the limit (d − t) → 0 are analysed for the case of homogeneous IFSs with µ = H d | Γ in Subsection 4.4.