# High-order numerical integration on self-affine sets

Patrick Joly, Maryna Kachanovska, Zoïs Moitier

► To cite this version:

Patrick Joly, Maryna Kachanovska, Zoïs Moitier. High-order numerical integration on self-affine sets.
2024. hal-04711030

HAL Id: hal-04711030

https://inria.hal.science/hal-04711030v1

Preprint submitted on 26 Sep 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# High-order numerical integration on self-affine sets

Patrick JOLY¹

Maryna KACHANOVSKA¹

Zoïs MOITIER²

¹POEMS, CNRS, Inria, ENSTA Paris, Institut Polytechnique de Paris, 91120 Palaiseau, France.

²IDEFIX, Inria, ENSTA Paris, Institut Polytechnique de Paris, 91120 Palaiseau, France.

September 26, 2024

## Abstract

We construct an interpolatory high-order cubature rule to compute integrals of smooth functions over self-affine sets with respect to an invariant measure. The main difficulty is the computation of the cubature weights, which we characterize algebraically, by exploiting a self-similarity property of the integral. We propose an *h*-version and a *p*-version of the cubature, present an error analysis and conduct numerical experiments.

## 1 Introduction

A significant portion of research in modern numerical analysis is dedicated to numerical approximation of solutions to the problems posed on rough domains or of irregular functions, see [AST06; Ban07; Cap10; RV13; DR22; DRT22; Cef+21; BB21; RS24] and references therein. In particular, a series of recent articles [CHM17; CH18; Cha+21; BGH22; Cae+24] deals with approximation of wave scattering by fractal screens, a problem that has applications in antenna engineering. The underlying numerical method, suggested in [Cae+24], relies on approximating integrals on fractal screens. Very few works have so far addressed the question of constructing cubature on fractal sets, which is a goal of the present work.

It seems that one of the first works treating cubature on fractal sets is due to G. Montica [Man96], who constructs a Gaussian quadrature for 1D fractals. The chaos-game cubatures, based on Monte-Carlo method, were suggested in [FMV98]. In [GHM23], the authors design and analyze a barycentric rule for integration of regular functions on fractals, and propose a clever method in [GHM24] to compute singular integrals of BEM. Up to our knowledge, existing methods are either uni-dimensional, or of low order.

Our approach is inspired by the work of Strichartz [Str99], who suggests an elegant method to evaluate integrals of polynomials on fractals using self-similarity (see Appendix A). Adapting this idea to our setting yields a cubature method of arbitrary high order for integrating functions on self-similar sets in $\mathbb{R}^n$.

This article is organized as follows. In Section 2 we introduce the notions related to the geometry and measure of self-similar sets, as well as polynomial spaces. Section 3 and Section 4 describe the new method to construct cubature on self-affine sets. In Section 5 we propose two versions of the cubature: an *h*-version and a *p*-version, and discuss the related error estimates. Finally, Section 6 contains numerical experiments.

## 2 Problem setting

### 2.1 Geometry, measure and integration

#### 2.1.1 Self-affine sets

The exposition below follows [Fal14, Ch. 9]. We will denote by $|\cdot|_2$ the Euclidean norm on $\mathbb{R}^n$.
**Definition 2.1.** An Iterated Function System (IFS) on $\mathbb{R}^n$ is a finite family of contractive maps $\mathcal{S} = \{S_\ell: \mathbb{R}^n \to \mathbb{R}^n: \ell = 1, 2, \dots, L\}$, with $L \ge 2$. Namely, for $\ell \in \mathbb{L} := \{1, 2, \dots, L\}$, there exists $0 \le \rho_\ell < 1$ such that the maps $S_\ell$ satisfy the estimation $|S_\ell(x) - S_\ell(y)|_2 \le \rho_\ell |x - y|_2$, for all $x, y \in \mathbb{R}^n$. The fixed point of a contractive mapping $S_\ell$ is denoted by $c_\ell$.

**Definition 2.2.** Let $K$ denote the set of all nonempty compact sets of $\mathbb{R}^n$. For an IFS $\mathcal{S} = \{S_\ell: \ell \in \mathbb{L}\}$, we define the Hutchinson operator by

$$ \mathcal{H}: K \to K \quad \text{such that} \quad \mathcal{H}(E) := S_1(E) \cup \dots \cup S_L(E). $$

The above definition is justified by the following result from [Fal14, Thm. 9.1].

**Theorem 2.3.** For an IFS $\mathcal{S} = \{S_\ell: \ell \in \mathbb{L}\}$, there exists a unique nonempty compact set $\Gamma \in \mathcal{K}$, which is a fixed point of $\mathcal{H}$, meaning that $\Gamma = \mathcal{H}(\Gamma)$. This set is called an attractor of the IFS $\mathcal{S}$.

An important class of IFS are the *affine IFS*, for which all the maps $S_\ell$ are affine:

$$ S_{\ell}: x \mapsto A_{\ell}x + b_{\ell} \quad \text{and} \quad c_{\ell} = (I - A_{\ell})^{-1}b_{\ell}. \qquad (2.1) $$

Then the contraction $\rho_\ell$ is equal to the spectral norm of the matrix $A_\ell$. A classic subclass of affine IFS are the *similar IFS* where the maps $S_\ell$ are contractive similarities, meaning that the maps $S_\ell$ satisfy $|S_\ell(x) - S_\ell(y)|_2 = \rho_\ell |x - y|_2$, for all $x, y \in \mathbb{R}^n$ which is equivalent to saying that

$$ A_{\ell} = \rho_{\ell} T_{\ell}, \quad \text{where } T_{\ell} \text{ is an orthogonal matrix.} \qquad (2.2) $$

In what follows, we will refer to an attractor of an IFS in $\mathbb{R}^n$ as a *fractal set*. For the particular case when the IFS is affine (resp. similar), the corresponding attractor will be referred to as a *self-affine set* (resp. *self-similar set*). From the above definition, we have immediately the following result, which will be of importance later.

**Lemma 2.4.** The pre-image of $\Gamma$ by the maps $S_\ell$, $\ell \in \mathbb{L}$, satisfies $S_\ell^{-1}(\Gamma) \cap \Gamma = \Gamma$.

*Proof.* Since $\mathcal{H}(\Gamma) = \Gamma$, we have, with $S_\ell^{-1}(\Gamma)$ denoting the pre-image of $\Gamma$,

$$ S_{\ell}^{-1}(\Gamma) = S_{\ell}^{-1}(\mathcal{H}(\Gamma)) = S_{\ell}^{-1}(S_{1}(\Gamma)) \cup \dots \cup S_{\ell}^{-1}(S_{L-1}(\Gamma)) \supset S_{\ell}^{-1}(S_{\ell}(\Gamma)). $$

The desired result follows from $\Gamma \subset S_\ell^{-1}(S_\ell(\Gamma))$. $\square$

In general, fractal sets are non-Lipschitz, but can be approximated by sets of a simpler structure, called *pre-fractals* or *pre-attractors*. A sequence of pre-fractals can be constructed as $\mathcal{H}^p(F)$, with $F$ being an arbitrary compact set. In a well-chosen topology, see the proof of **Theorem 2.3** in [Fal14, Thm. 9.1], $\mathcal{H}^p(F) \to \Gamma$ as $p \to +\infty$.

*Example 2.5.* On Fig. 1, we have plotted classic examples of self-affine sets. The IFS corresponding to the Fat Sierpiński triangle, see [BMS04; Jor06], is composed of the three maps $x \to \rho x + (1 - \rho)c_\ell$ where $\rho = (\sqrt{5} - 1)/2$ and the $c_\ell$ are the vertices of an equilateral triangle. The one for the Koch snowflake is composed of seven contractive similarities, see [GHM23, Fig. 3], and the IFS for the Barnsley fern is composed of four maps where the coefficients are given by [Bar93, Tb. 3.8.3].

**Lemma 2.6.** For an IFS $\mathcal{S} = \{S_\ell: \mathbb{R}^n \to \mathbb{R}^n: \ell \in \mathbb{L}\}$, there exists a nonempty compact set $K \in \mathcal{K}$ such that $S_\ell(K) \subset K$, for all $\ell \in \mathbb{L}$.

*Proof.* We define $K = \bar{B}(z, r)$ as a closed ball of radius $r > 0$ centered at $z \in \mathbb{R}^n$. Let us show that $\mathcal{H}(K) \subset K$ for $r > 0$ sufficiently large. Since the maps $S_\ell$ are $\rho_\ell$-Lipschitz, we have $S_\ell(K) \subset \bar{B}(S_\ell(z), \rho_\ell r)$. In order to have $\mathcal{H}(K) \subset K$, it suffices that $r \ge |z - S_\ell(z)|_2 + \rho_\ell r$, for all $\ell \in \mathbb{L}$. Therefore, for any $r \ge \max_{\ell \in \mathbb{L}}|z - S_\ell(z)|_2/(1 - \rho_\ell)$, the nonempty compact set $K = \bar{B}(0, r)$ satisfies the lemma. $\square$
<figure><img src="image_2.png" /><figcaption>Figure 1: Examples of self-affine sets.</figcaption></figure>

**Remark 2.7.** *From the proof of Lemma 2.6, we see that we can take any ball with a large enough radius as pre-fractal of level 0 to get a decreasing sequence of pre-fractals. Naturally, to approximate the smallest ball, we can minimize the function*

$$z \in \mathbb{R}^n \mapsto \max_{\ell \in \mathbb{L}} \frac{|z - S_\ell(z)|_2}{1 - \rho_\ell}. \qquad (2.3)$$

### 2.1.2 Invariant measure

The exposition on invariant measure of an IFS is taken from [Hut81]. Given a compact subset $K \subset \mathbb{R}^n$, equipped with the usual Euclidean topology, one can define an outer measure on $K$ as a map $\nu: 2^K \to [0, +\infty]$ such that $\nu(\emptyset) = 0$ (the null empty set condition) and, for any $A, B_0, B_1, \dots \subset K$, $A \subset \bigcup_{i \in \mathbb{N}} B_i$ implies $\nu(A) \le \sum_{i \in \mathbb{N}} \nu(B_i)$ (the countably subadditive condition). We say that a set $A \subset K$ is $\nu$-measurable, if $\nu(E) = \nu(E \cap A) + \nu(E \setminus A)$ for all $E \subset K$ (the Carathéodory's criterion). The family of measurable sets forms a $\sigma$-algebra. When restricted to $\nu$-measurable sets, the outer measure $\nu$ produces a measure. In the sequel, for a reason that will appear in Theorem 2.8, we shall restrict ourselves to Borel regular measures: an outer measure $\nu$ is Borel regular if, and only if, all Borel sets are measurable, and, for all $A \subset K$ there exists a Borel set $B$, s.t. $A \subset B$ and $\nu(A) = \nu(B)$. The support of the measure $\nu$ is defined by $\text{supp} \nu := K \setminus \bigcup \{V \text{ open } | \nu(V) = 0\}$. Let us further introduce

$$M_K^1 := \{\nu \mid \nu \text{ is Borel regular and } \nu(K) = 1\}.$$

Given $g: K \to K$, measurable, the push-forward measure $\nu \circ g^{-1} \in M_K^1$ is defined by

$$\forall A \subset K, \quad \nu \circ g^{-1}(A) := \nu(g^{-1}(A)), \text{ with } g^{-1}(A) := \{x \in K \mid g(x) \in A\}. \qquad (2.4)$$

Let us now construct an invariant measure associated to an affine IFS $\mathcal{S}$ given by the maps $\{S_\ell : \ell \in \mathbb{L}\}$. First, let $K$ satisfy the assumption of Lemma 2.6 (i.e. $\mathcal{H}(K) \subset K$).

Given $\boldsymbol{\mu} = (\mu_0, \dots, \mu_{L-1}) \in (0, 1)^L$ such that $|\boldsymbol{\mu}|_1 := \sum_{\ell \in \mathbb{L}} \mu_\ell = 1$, let us define the map $\mathcal{M}: M_K^1 \to M_K^1$ via the following identity<fnref n="1" />:

$$\mathcal{M}: \nu \mapsto \sum_{\ell \in \mathbb{L}} \mu_{\ell} \nu \circ S_{\ell}^{-1}. \qquad (2.5)$$

A measure satisfying $\nu = \mathcal{M}\nu$ is called invariant with respect to $(\mathcal{S}, \boldsymbol{\mu})$.

**Theorem 2.8.** Let $\mathcal{S} = \{S_\ell : \ell \in \mathbb{L}\}$ and $\boldsymbol{\mu}$, $\mathcal{M}$ be as in Eq. (2.5). Then there exists a unique invariant measure $\mu \in M_\Gamma^1$ on $\Gamma$, the attractor of $\mathcal{S}$, such that $\mu = \mathcal{M}\mu$. Furthermore, the support of the measure $\mu$ is $\Gamma$ (which implies that for any measurable $A$, $\mu(A) = \mu(A \cap \Gamma)$).

<footnote n="1">¹The interested reader can verify that $\mathcal{M}$ maps $M_K^1$ into itself because $|\boldsymbol{\mu}|_1 = 1$.</footnote>
Proof. This proof is due to [Hut81]; here we fill in some arguments omitted in the above reference. We equip $M_K^1$ with the Monge-Kantorovich distance $d(\mu, \nu)$:

$$ \forall (\mu, \nu) \in M_K^1, \quad d(\mu, \nu) = \sup \left\{ \int_K f(x) \, d\mu - \int_K f(x) \, d\nu \mid f: K \to \mathbb{R}, \text{ Lip}(f) \le 1 \right\}. $$

Then $(M_K^1, d)$ is a complete metric space, see the discussion after [Bog07, Thm. 8.10.43, Vol. II]; in the latter reference, the author restricts his attention to the Radon measures, however, since $K$ is complete and separable, measures in $M_K^1$ are Radon measures, see [Bog07, Thm. 7.1.7, Vol. II]. The operator $\mathcal{M}$ is a contractive map on this metric space, see the proof of [Hut81, Thm. 4.4(1)]. Existence and uniqueness of $\mu$ follows by the fixed point theorem. The property $\text{supp} \mu = \Gamma$ is proven in [Hut81, Thm. 4.4(4)]. □

**Remark 2.9** (Hausdorff measure and application to the wave scattering by fractal screens). Given an IFS, one can construct infinitely many invariant measures, parametrized by $\mu$. Nonetheless, in some cases, this measure coincides with the (normalized) $d$-dimensional Hausdorff measure $\mathcal{H}^d(\Gamma)$ (see [Bog07, 3.10(iii)]). Assume that the IFS $\mathcal{S} = \{S_\ell : \ell \in \mathbb{L}\}$ satisfy the Open Set Condition (OSC). More precisely, there exists a nonempty bounded open set $U$ such that $\mathcal{H}(U) \subset U$ and $S_i(U) \cap S_j(U) = \emptyset$ for $i \neq j \in \mathbb{L}$. Moreover, let the IFS be a similar IFS, cf. Eq. (2.2). In this case one can show that the Hausdorff dimension $d$ of the $\mathcal{S}$ is a unique solution $d$ to $\sum_{\ell \in \mathbb{L}} \rho_\ell^d = 1$, see [Fal14, Thm. 9.3]. Choosing $\mu_\ell = \rho_\ell^d$ yields the invariant measure $\mu(E) = \mathcal{H}^d(E \cap \Gamma)/\mathcal{H}^d(\Gamma)$, see [Hut81, Thm. 5.3(1)(iii)].

This example is of importance for applications in the computational scattering theory. Recently, Hausdorff measure based boundary element methods (BEM) were developed for scattering by fractal screens, see [Cae+24]. A non-exhaustive list of self-similar sets satisfying the OSC includes Sierpinski gasket, Koch snowflake and Cantor sets. Let us remark that one of favorable properties of the Hausdorff measure compared is its invariance with respect to orthogonal transformations. This can be used for efficient computation of singular integrals arising in the Hausdorff BEM, cf. [GHM23].

**Remark 2.10.** The above construction of the invariant measure remains valid and yields the same measure if the compact $K$ is replaced by $\Gamma$.

### 2.1.3 Integrals with respect to invariant measures and a Ruelle operator

Let us fix the sets $K$ as in Lemma 2.6, $\mu \in (0, 1)^L$, with $|\mu|_1 = 1$ and $\mu$ the associated invariant measure as defined in Theorem 2.8.

As discussed in the introduction, the goal of this work is to approximate numerically $\int_{\Gamma} f \, d\mu$, with $f$ being sufficiently regular and defined in the vicinity of $\Gamma$ (e.g. in $\mathcal{C}^0(K, \mathbb{C})$). This will be done by using a certain self-similarity property of the integral, which is introduced in the present section. The proof of this property is in turn based on the following change-of-variables formula, see [Bog07, Thm. 3.6.1.]:

$$ \int_{\Gamma} f(x) \, d(\mu \circ g^{-1}) = \int_{g^{-1}(\Gamma)} f \circ g(x) \, d\mu, \qquad (2.6) $$

In the above, $g: \Gamma \to \Gamma$ is measurable, and the measurable function $f$ is such that $f \circ g \in \mathbb{L}^1(\Gamma)$; the push-forward measure $\mu \circ g^{-1}$ is defined in Eq. (2.4).

**Motivation** By analogy with the Hutchinson operator $\mathcal{H}$ acting on sets, recall Definition 2.2, and the Hutchinson-like operator $\mathcal{M}$ acting on measures, see Eq. (2.5), one can introduce a so-called Ruelle operator [FL99] acting on complex-valued functions defined in $\mathbb{R}^n$, for instance in $\mathcal{C}(\mathbb{R}^n)$, by:

$$ \mathcal{F}: f \in \mathcal{C}(\mathbb{R}^n) \mapsto \sum_{\ell \in \mathbb{L}} \mu_\ell f \circ S_\ell \in \mathcal{C}(\mathbb{R}^n). \qquad (2.7) $$
The reader will note that, with respect to the operator $\mathcal{H}$, see **Definition 2.2**, compact sets are replaced by functions and the union is replaced by the weighted addition.

**Proposition 2.11.** Given an IFS $\mathcal{S} = \{S_\ell : \ell \in \mathbb{L}\}$, its attractor $\Gamma$, and an invariant measure $\mu$ associated to $\mathcal{S}$, for $f \in \mathcal{C}(\Gamma)$, we have

$$ \int_{\Gamma} f \, d\mu = \sum_{\ell \in \mathbb{L}} \mu_{\ell} \int_{\Gamma} f \circ S_{\ell} \, d\mu = \int_{\Gamma} \mathcal{F}[f] \, d\mu, $$

with $\mathcal{F}$ defined in Eq. (2.7).

*Proof.* Given $f \in \mathcal{C}(\Gamma)$, we compute, by definition of $\mathcal{M}$ in (2.5) and $\mathcal{M} \mu = \mu$,

$$ \int_{\Gamma} f \, d\mu = \int_{\Gamma} f \, d(\mathcal{M} \mu) = \sum_{\ell \in \mathbb{L}} \mu_{\ell} \int_{\Gamma} f \, d(\mu \circ S_{\ell}^{-1}). \quad (2.8) $$

From the change of variable formula Eq. (2.6) and $\text{supp}(\mu) = \Gamma$, we get

$$ \int_{\Gamma} f \, d(\mu \circ S_{\ell}^{-1}) = \int_{S_{\ell}^{-1}(\Gamma)} f \circ S_{\ell} \, d\mu = \int_{S_{\ell}^{-1}(\Gamma) \cap \Gamma} f \circ S_{\ell} \, d\mu = \int_{\Gamma} f \circ S_{\ell} \, d\mu, $$

where the latter identity follows from **Lemma 2.4**. We conclude with Eq. (2.8). $\square$

**Ruelle operator on spaces of polynomials for affine IFS** We shall be particularly interested in the action of $\mathcal{F}$ on polynomials. From now on, we will limit the discussion to *affine* IFS. We denote by $\mathbb{P}$ the space of polynomials of $n$ variables and, for any $k \in \mathbb{N}$, by $\mathbb{P}_k$ the subspace of polynomials of total degree less or equal to $k$. If case of ambiguity, we add the number of variables $n$ as a superscript to $\mathbb{P}^n$ and $\mathbb{P}_k^n$. The space of polynomials can be decomposed into a direct sum of homogeneous polynomials of degree $k$: for any multi-index $\alpha \in \mathbb{N}^n$, we define the monomial $x^\alpha := x_1^{\alpha_1} \cdots x_n^{\alpha_n} \in \mathbb{P}$ of total degree $|\alpha|_1 = \alpha_1 + \cdots + \alpha_n$ and

$$ \mathbb{H}_k := \operatorname{span}\{x^\alpha \mid |\alpha|_1 = k\}, \quad \text{so that} \quad \mathbb{P}_k = \bigoplus_{k' \le k} \mathbb{H}_{k'} \text{ and } \mathbb{P} = \bigoplus_{k \in \mathbb{N}} \mathbb{H}_k \quad (2.9) $$

Seen as a subspace of $\mathcal{C}(\mathbb{R}^n)$, $\mathbb{P}$ is kept invariant by $\mathcal{F}$ and the action of $\mathcal{F}$ does not increase the total degree, in other words, for all $p \in \mathbb{P}_k$, we have $\mathcal{F} p \in \mathbb{P}_k$. This means that $\mathcal{F}$ has an upper triangular structure with respect to the decomposition (2.9) of $\mathbb{P}$: there exists a family of linear operators $\mathcal{F}_{k,k'} \in \mathcal{L}(\mathbb{H}_k, \mathbb{H}_{k'})$, with $k' \le k$ such that

$$ \mathcal{F} p = \sum_{k' \le k} \mathcal{F}_{k,k'} p, \quad \text{for all } p \in \mathbb{H}_k. \quad (2.10) $$

In **Section 3**, we shall need some information about the eigenvalues of $\mathcal{F}$ restricted to the space of polynomials $\mathbb{P}$, which we denote by $\mathcal{F}_{\mathbb{P}}$. Let us denote by $\sigma(\mathcal{F}_{\mathbb{P}})$ the set of eigenvalues of $\mathcal{F}_{\mathbb{P}}$ (namely, values $\lambda \in \mathbb{C}$, s.t. there exists $p \in \mathbb{P}$, $p \ne 0$, for which $\mathcal{F}_{\mathbb{P}} p = \lambda p$). Due to the upper-triangular like structure of $\mathcal{F}_{\mathbb{P}}$, it holds that

$$ \sigma(\mathcal{F}_{\mathbb{P}}) = \bigcup_{k \in \mathbb{N}} \sigma(\mathcal{F}_{k,k}). \quad (2.11) $$

For $k = 0$, one easily sees that, as $\mu_1 + \cdots + \mu_L = 1$,

$$ \mathcal{F}_{0,0} = \mathrm{Id}, \quad \text{(the identity operator in } \mathbb{P}_0), \quad (2.12) $$

which means that 1 is an eigenvalue of $\mathcal{F}_{\mathbb{P}}$. Furthermore, one can localize the eigenvalues of the operators $\mathcal{F}_{k,k}$, for $k \ne 0$, as made precise in the following lemma.
**Lemma 2.12.** For any $k \ge 1$, $\sigma(\mathcal{F}_{k,k}) \subset \bar{D}(0, r_k)$ (the closed disk of center 0 and radius $r_k$), where

$$r_k := \sum_{\ell \in \mathbb{L}} \mu_\ell \rho_\ell^k \le \rho_{\max}^k < 1, \quad \rho_{\max} = \max_{\ell \in \mathbb{L}} \rho_\ell.$$

As a corollary, 1 is a simple eigenvalue of the operator $\mathcal{F}_{\mathbb{P}}$ (i.e. its algebraic, and thus, geometric, multiplicity equals to 1). Moreover, it is also the largest eigenvalue of $\mathcal{F}_{\mathbb{P}}$.

*Proof.* With $S_\ell(x) = A_\ell x + b_\ell$, one sees immediately that for all $p \in \mathbb{H}_k$, we have $p \circ S_\ell(x) = p(A_\ell x) + R_\ell(x)$ where $R_\ell \in \mathbb{P}_{k-1}$. From this observation, we infer the following identity, valid for all $p \in \mathbb{H}_k$:

$$\mathcal{F}_{k,k} p(x) = \sum_{\ell \in \mathbb{L}} \mu_\ell p(A_\ell x).$$

Let $B = \bar{B}(0, 1)$ be the unit closed ball of $\mathbb{R}^n$ and $p \in \mathbb{H}_k$ be an eigenvector of $\mathcal{F}_{k,k}$ associated to $\lambda \in \mathbb{C}$. Of course $\|p\|_{\mathbb{L}^\infty(B)} \ne 0$ and there exists $x_* \in B$ such that $|p(x_*)| = \|p\|_{\mathbb{L}^\infty(B)}$. However, $\mathcal{F}_{k,k} p = \lambda p$ implies $\mathcal{F}_{k,k} p(x_*) = \lambda p(x_*)$, thus

$$|\lambda p(x_*)| \le \sum_{\ell \in \mathbb{L}} \mu_\ell |p(A_\ell x_*)|$$

and since each matrix $A_\ell$ is $\rho_\ell$-contractive, $A_\ell x_* \in \rho_\ell B$, thus, by homogeneity,

$$|p(A_\ell x_*)| \le \max_{x \in \rho_\ell B} |p(x)| = \max_{y \in B} |p(\rho_\ell y)| = \rho_\ell^k \max_{y \in B} |p(y)| = \rho_\ell^k |p(x_*)|.$$

We then conclude the proof since, as $|\boldsymbol{\mu}|_1 = 1$,

$$|\lambda| \le \sum_{\ell \in \mathbb{L}} \mu_\ell \rho_\ell^k \le \rho_{\max}^k \left( \sum_{\ell \in \mathbb{L}} \mu_\ell \right) = \rho_{\max}^k.$$

□

**Remark 2.13.** In dimension $n = 1$, since all the spaces $\mathbb{H}_k$ are one-dimensional, the operators $\mathcal{F}_{k,k}$ are reduced to the multiplication operators. The eigenvalues $\sigma(\mathcal{F}_{\mathbb{P}})$ of the operator $\mathcal{F}_{\mathbb{P}}$ are known explicitly; they are real and positive. Indeed, from the proof of Lemma 2.12, one easily sees that

$$\sigma(\mathcal{F}_{\mathbb{P}}) = \left\{ \sum_{\ell \in \mathbb{L}} \mu_{\ell} \rho_{\ell}^{k}, k \in \mathbb{N} \right\}.$$

The above reasoning can be extended to a higher dimensional case, when all matrices $A_\ell$ are diagonal, $A_\ell = \operatorname{diag}(a_{\ell,1}, \dots, a_{\ell,n})$ with $|a_{\ell,i}| \le \rho_\ell < 1$. In this case, it is readily seen that $(S_\ell(x))^\alpha = a_\ell^\alpha x^\alpha + R_\alpha$, where $a_\ell^\alpha = a_{\ell,1}^{\alpha_1} \cdots a_{\ell,n}^{\alpha_n}$ and $R_\alpha \in \mathbb{P}_{k-1}$. Thus the operator $\mathcal{F}_{k,k}$ is diagonal in the basis $\{x^\alpha, \text{ s.t. } |\alpha|_1 = k\}$ of $\mathbb{H}_k$ and we have

$$\sigma(\mathcal{F}_{\mathbb{P}}) = \bigcup_{k \in \mathbb{N}} \left\{ \sum_{\ell \in \mathbb{L}} \mu_\ell a_\ell^\alpha, |\alpha|_1 = k \right\}.$$

## 2.2 Interpolatory cubature on self-affine sets

### 2.2.1 Definition and associated operators

Let us assume that $\Gamma \subset K$, where $K$ is defined as in Lemma 2.6. As discussed in the introduction, the goal of this article is to construct a cubature for an approximation of the integral of the function $f \in \mathcal{C}(K)$ over $\Gamma$ with the invariant measure $\mu$. One of the commonly used types of cubature formulae provides such an approximation of the form

$$\int_{\Gamma} f(x) \, d\mu \approx \mathcal{Q}[f] := \sum_{1 \le i \le M} w_i f(x_i), \qquad (2.13)$$
where $\mathcal{X} := \{x_i: 1 \le i \le M\} \subset K$ is the set of cubature points and the vector $\mathbf{w} = (w_1, \dots, w_M) \in \mathbb{R}^J$ is the vector of cubature weights.

**Remark 2.14.** As the cubature formula depends of on both $\mathcal{X}$ and $\mathbf{w}$, a better notation would be $Q_{\mathcal{X},\mathbf{w}}$. However, we have abandoned this option in order to avoid heavy notation, the role of $\mathcal{X}$ and $\mathbf{w}$ being implicit.

Let us remark that common cubature formulae (e.g. designed for $\Gamma$ being an interval and $\mu$ being the Lebesgue's measure) rely on a set of points $\mathcal{X} \subset \Gamma$. For self-affine sets, the condition $\mathcal{X} \subset \Gamma$ can be difficult to ensure, and therefore, we do not provide any constraints on the location of the points $\mathcal{X}$, but the ones described below. For constructing $\mathcal{Q}$, we shall follow the classic paradigm for constructing a cubature formula (interpolatory cubature), where one requires that "$\simeq$" becomes "=" in Eq. (2.13), when $f$ belongs to a finite dimensional space $\mathcal{P} \subset \mathbb{P}$ of polynomials. The idea behind being that any continuous function can be approximated uniformly by a sequence of polynomials. In the rest of the article, we assume that $\mathbb{P}_0 \subset \mathcal{P}$.

**Remark 2.15.** (Choice of the space $\mathcal{P}$). The usual practical choice for the space is $\mathcal{P} = \mathbb{P}_k$ or $\mathcal{P} = \mathbb{Q}_k$, both choices being particularly well-suited for the error analysis. Recall that the space $\mathbb{Q}_k$ is defined as $\operatorname{span}\{x^{\alpha} \mid |\alpha|_{\infty} \le k\}$ with $|\alpha|_{\infty} = \max_{1 \le i \le n} \alpha_i$. However, all theoretical arguments extend to a general space $\mathcal{P}$.

Given $\mathcal{P}$, the usual procedure is first to consider a set $\mathcal{X} := \{x_i: 1 \le i \le M\}$, where $M = \dim \mathcal{P}$, of cubature points which is $\mathcal{P}$-unisolvent.

**Definition 2.16.** The evaluation operator associates to any continuous functions the set of its values on $\mathcal{X}$:

$$ \mathcal{E}: \mathcal{C}(\mathbb{R}^n) \to \mathbb{R}^M, \quad f \mapsto \mathbf{f} := (f(x_1), \dots, f(x_M)). $$

Then, by definition the set of points $\mathcal{X}$ is $\mathcal{P}$-unisolvent if, and only if, the map $\mathcal{E}$ is injective from $\mathcal{P}$ into $\mathbb{R}^M$.

If $\mathcal{X}$ is $\mathcal{P}$-unisolvent, $\mathcal{E}_{\mathcal{P}} := \mathcal{E}|_{\mathcal{P}}$ is thus bijective from $\mathcal{P}$ to $\mathbb{R}^M$. The $\mathcal{P}$-unisolvency in $\mathcal{X}$ is a necessary and sufficient condition for constructing the set of Lagrange polynomials $\mathcal{L}_i$ in $\mathcal{P}$, associated to the set $\mathcal{X}$, defined by $\mathcal{L}_i(x_j) = \delta_{i,j}$, which is equivalent to $\mathcal{L}_i = \mathcal{E}_{\mathcal{P}}^{-1}(e_i)$, where $\{e_i\}$ is the canonical basis of $\mathbb{R}^M$. The set $(\mathcal{L}_i)_{i=1}^M$ is a basis of $\mathcal{P}$ and

$$ \forall p \in \mathcal{P}, \quad p(x) = \sum_{1 \le i \le M} p(x_i) \mathcal{L}_i(x). \qquad (2.14) $$

**Remark 2.17.** Even though this is not explicit in the notation, the Lagrange polynomials $\mathcal{L}_j$ do depend on $\mathcal{P}$ and $\mathcal{X}$.

**Remark 2.18.** The inverse of the operator $\mathcal{E}_{\mathcal{P}}$ is nothing but the interpolation operator that associates to a vector $\mathbf{a} = (a_1, \dots, a_M)$ the unique polynomial of $\mathcal{P}$ that takes the value $a_i$ at point $x_i$, namely

$$ \mathcal{I}_{\mathcal{P}}: \mathbb{R}^M \to \mathcal{P} \quad \text{such that} \quad \mathcal{I}_{\mathcal{P}}(\mathbf{a}) = \sum_{1 \le i \le M} a_i \mathcal{L}_i. \qquad (2.15) $$

We shall use later the adjoint operator $\mathcal{I}_{\mathcal{P}}^*$: $(\mathcal{P}; \|\cdot\|_{\mathbb{L}^2(\Gamma;\mu)}) \to \mathbb{R}^M$ (where $\mathbb{R}^M$ is equipped with the euclidean inner product). The reader will easily verify that

$$ \mathcal{I}_{\mathcal{P}}^* f = \mathbf{f}, \quad f_i = \int_{\Gamma} f \mathcal{L}_i \, d\mu, \quad i = 1, \dots, M. \qquad (2.16) $$

Note that the cubature formula Eq. (2.13) can be rewritten in algebraic form as

$$ \forall f \in \mathcal{C}(K), \quad \mathcal{Q}[f] = \mathbf{w} \cdot \mathcal{E}(f) \qquad (2.17) $$

where $x \cdot y$ the inner product in $\mathbb{R}^M$. Given $\mathcal{P}$ and $\mathcal{X}$ a $\mathcal{P}$-unisolvent set, we would like to find $\mathbf{w}$ in such a way that we integrate exactly all the functions in $\mathcal{P}$. The (trivial) answer is provided by the following lemma.
**Lemma 2.19.** Given $\mathcal{P} \subset \mathbb{P}$ and a $\mathcal{P}$-unisolvent set $\mathcal{X}$, there exists a unique choice of cubature weights $\{w_i : 1 \le i \le M\}$ such that $\int_{\Gamma} p(x) \mathrm{d}\mu = \mathcal{Q}[p]$, for all $p \in \mathcal{P}$. These cubature weights are given by

$$w_i = \int_{\Gamma} \mathcal{L}_i(x) \mathrm{d}\mu, \quad \forall 1 \le i \le M. \qquad (2.18)$$

Lemma 2.19 is essentially of theoretical interest since exploiting it would require to compute the integrals in Eq. (2.18). This justifies an alternative approach to define cubature weights that we develop in the next sections.

Before proceeding into the definition of the cubature weights, we will need to introduce the notion of $\mathcal{S}$-invariance of a polynomial space. We shall see, in Section 3, that cubature weights based on such polynomial spaces can be computed in a purely algebraic manner. On the other hand, as we shall see in Section 4, for polynomial spaces not satisfying this condition, the natural extension of the algebraic method will provide a set of practically useful cubature weights.

### 2.2.2 Definition of $\mathcal{S}$-invariant spaces

**Definition 2.20.** A finite-dimensional subspace $\mathcal{P}$ of $\mathbb{P}$ is said to be $\mathcal{S}$-invariant if and only if $\mathbb{P}_0 \subset \mathcal{P}$ and $\mathcal{P}$ is stable by the operator $\mathcal{F}$.

**Remark 2.21.** The reader will easily verify that for any $\mathcal{S}$ and any $k$, the space $\mathbb{P}_k$ is $\mathcal{S}$-invariant. If, for all $\ell$, the matrix $A_\ell$ has in each row and column only one non-zero element, the space $\mathbb{Q}_k$ is $\mathcal{S}$-invariant. However, in general, the space $\mathbb{Q}_k$ is not $\mathcal{S}$-invariant, as shown in the next example.

**Example 2.22.** Let $n = 2$, $L = 2$, and the matrices $A_k$ in Eq. (2.1) are defined as follows: $A_0 = \rho \mathrm{Id}$ and $A_1 = \rho T$, where

$$T = \frac{\sqrt{2}}{2} \begin{pmatrix} 1 & 1 \\ 1 & -1 \end{pmatrix} \iff Tx = \frac{\sqrt{2}}{2} (x_1 + x_2, x_1 - x_2)^\mathrm{T}.$$

We observe that $p(x) = x_1^k x_2^k \in \mathbb{Q}_k$ but $p(Tx) = \frac{1}{2}(x_1 + x_2)^k (x_1 - x_2)^k \notin \mathbb{Q}_k$. The space $\mathcal{F}(\mathbb{Q}_k)$ is included into $\mathbb{P}_{nk}$. This inclusion can be strict. Indeed, in the above case, for $k = 1$, we have that, the polynomials $p_{i,j}(x) = x_1^i x_2^j$, $i, j \in \{0, 1\}$, satisfy $p_{i,j} \circ T \in \mathbb{Q}_1$ if $ij = 0$, and $p_{1,1} \circ T = \frac{1}{2}(x_1^2 - x_2^2)$. The above computation implies that $\mathcal{F}(\mathbb{Q}_k) = \mathrm{span}\{1, x_1, x_2, x_1x_2, x_1^2 - x_2^2\}$, which is a strict subspace of $\mathbb{P}_2$.

# 3 Cubature based on $\mathcal{S}$-invariant polynomial spaces

The key assumption of this section, unless stated otherwise, is that the space $\mathcal{P}$ is $\mathcal{S}$-invariant.

## 3.1 Algebraic definition of the cubature weights

Recall that our goal is to construct a cubature rule Eq. (2.13), which would be exact for $p \in \mathcal{P}$. As shown in Lemma 2.19, this condition is ensured by the unique choice $w_i = \int_{\Gamma} \mathcal{L}_i \mathrm{d}\mu$, for $1 \le i \le M$. On the other hand, by Proposition 2.11, we can relate the integrals of the Lagrange polynomials to the integrals of their images by the Ruelle operator:

$$w_i = \int_{\Gamma} \mathcal{L}_i \mathrm{d}\mu = \int_{\Gamma} \mathcal{F}[\mathcal{L}_i] \mathrm{d}\mu. \qquad (3.1)$$

As $\mathcal{P}$ is $\mathcal{S}$-invariant, $\mathcal{F}[\mathcal{L}_i] \in \mathcal{P}$, and is thus integrated exactly by the cubature rule $\mathcal{Q}$. In other words, the cubature weights satisfy the following identity:

$$w_i = \sum_{1 \le j \le M} w_j \mathcal{F}[\mathcal{L}_i](x_j) = \sum_{1 \le j \le M} \sum_{\ell \in \mathbb{L}} \mu_\ell \mathcal{L}_i \circ S_\ell(x_j) w_j.$$
In order to write the above in the algebraic form, let us define the $M \times M$ matrix

$$ S_{i,j} := \sum_{\ell \in \mathbb{L}} \mu_{\ell} \mathcal{L}_j \circ S_{\ell}(x_i), \quad \text{for } 1 \le i, j \le M. \qquad (3.2) $$

With this definition, the above said is summarized in the following lemma, which provides an algebraic property of the cubature weights.

**Lemma 3.1.** Assume that $\mathcal{P}$ is $\mathcal{S}$-invariant and that $\mathcal{X}$ is $\mathcal{P}$-unisolvent, then the vector $\mathbf{w}$ given by Eq. (2.18) is an eigenvector of $\mathbf{S}^\mathrm{T}$, associated to the eigenvalue 1:

$$ \mathbf{S}^\mathrm{T}\mathbf{w} = \mathbf{w}. \qquad (3.3) $$

Obviously, Eq. (3.3) not sufficient to characterize $\mathbf{w} \in \mathrm{Ker}(I - \mathbf{S}^\mathrm{T})$. We must complete it by requiring that the cubature formula is exact for the constant function equal to 1, namely, using the algebraic version Eq. (2.17) of $\mathcal{Q}$,

$$ \mathbf{w} \cdot \mathbf{1} = \sum_{1 \le i \le M} w_i = 1, \quad \text{where } \mathbf{1} := \mathcal{E}(x \mapsto 1). \qquad (3.4) $$

Thus, the problem to be solved, to compute the vector $\mathbf{w}$ of cubature weights writes

$$ \text{Find } \mathbf{w} \in \mathrm{Ker}(I - \mathbf{S}^\mathrm{T}) \text{ satisfying Eq. (3.4).} \qquad (3.5) $$

By Lemma 3.1, we know that this problem admits at least one solution $\mathbf{w}$ given by Eq. (2.18). However, a priori it is unclear that this solution would be unique. The answer to this question is directly linked to the spectral structure of $\mathbf{S}^\mathrm{T}$ and more precisely to the multiplicity of the eigenvalue 1. Because the spectra of the matrices $\mathbf{S}^\mathrm{T}$ and $\mathbf{S}$ (together with multiplicities) coincide, we will work with the matrix $\mathbf{S}$, which appears to have a convenient structure for the analysis.

## 3.2 Properties of the matrix $\mathbf{S}^\mathrm{T}$ and well-posedness of Eq. (3.5)

The matrix Eq. (3.2) is strongly linked to the Ruelle operator $\mathcal{F}$, which we will use to characterize its spectrum. The first part of the following lemma will be used immediately, while the second part will be of use later.

**Lemma 3.2.** For all $p \in \mathcal{P}$, where the space $\mathcal{P}$ is not necessarily $\mathcal{S}$-invariant, it holds that $\mathcal{E}(\mathcal{F} p) = \mathbf{S}(\mathcal{E} p)$. We also have $\mathcal{E}_{\mathcal{P}} \mathcal{I}_{\mathcal{P}} \mathcal{E}(\mathcal{F} p) = \mathbf{S}(\mathcal{E} p)$.

*Proof.* From the definition of $\mathcal{F}$, see Eq. (2.7) and the definition of $\mathbf{S}$, see Eq. (3.2), the columns of $\mathbf{S}$ are given by $\mathbf{S}_{\cdot,j} = \mathcal{E}(\mathcal{F} \mathcal{L}_j)$, for $1 \le j \le M$. Next, with the Lagrange decomposition Eq. (2.14) of $p$ and by linearity of $\mathcal{E}$, it holds that

$$ \mathcal{E}(\mathcal{F} p) = \sum_{1 \le j \le M} p(x_j) \mathcal{E}(\mathcal{F} \mathcal{L}_j) = \sum_{1 \le j \le M} \mathbf{S}_{\cdot,j} p(x_j). $$

The first statement is obtained using the definition of $\mathcal{E}$. The second statement follows by recalling that $\mathcal{E}_{\mathcal{P}} \mathcal{I}_{\mathcal{P}} = \mathrm{Id}_{\mathbb{R}^M}$ and applying $\mathcal{E}_{\mathcal{P}} \mathcal{I}_{\mathcal{P}}$ to the first identity. $\square$

The $\mathcal{S}$-invariance of $\mathcal{P}$ implies in particular that $\mathcal{F}_{\mathcal{P}} = \mathcal{F}|_{\mathcal{P}} \in \mathcal{L}(\mathcal{P})$. Using the fact that $\mathcal{E}_{\mathcal{P}}$ is invertible (recall that its inverse is the interpolation operator $\mathcal{I}_{\mathcal{P}}$, see Remark 2.18), Lemma 3.2 implies that

$$ \mathcal{F}_{\mathcal{P}} = \mathcal{E}_{\mathcal{P}}^{-1} \mathbf{S} \mathcal{E}_{\mathcal{P}}, \qquad (3.6) $$

in other words that the spectra of $\mathcal{F}_{\mathcal{P}}$ and $\mathbf{S}$ coincide.
**Lemma 3.3.** The eigenvalue 1 is a simple eigenvalue of $\mathbf{S}$ (or, equivalently, $\mathcal{F}_P$ according to Eq. (3.6)) and all other eigenvalues are strictly less than 1 in modulus. Moreover, the corresponding eigenspace is given by $\text{span}\{\mathbf{1}\}$.

*Proof.* Given the $\mathcal{S}$-invariant space $\mathcal{P}$, let the smallest $k$ be such that $\mathcal{P} \subset \mathbb{P}_k$. Let $W$ be a complement of $\mathcal{P}$ in $\mathbb{P}_k$. Let $\mathcal{F}_k \in \mathcal{L}(\mathbb{P}_k)$ be the restriction of $\mathcal{F}$ to $\mathbb{P}_k$, the $\mathcal{S}$-invariance property means that, with respect to the decomposition $\mathbb{P}_k = \mathcal{P} \oplus W$, the operator $\mathcal{F}_k$ has a triangular block structure

$$ \mathcal{F}_k = \begin{pmatrix} \mathcal{F}_{\mathcal{P} \to \mathcal{P}} & \mathcal{F}_{W \to \mathcal{P}} \\ 0 & \mathcal{F}_{W \to W} \end{pmatrix}. $$

As a consequence, we have the inclusion $\sigma(\mathcal{F}_{\mathcal{P} \to \mathcal{P}}) \subset \sigma(\mathcal{F}_k)$. On the other hand, remark that $\sigma(\mathcal{F}_k) = \bigcup_{0 \le m \le k} \sigma(\mathcal{F}_{m,m})$, cf. (2.10). The announced result about the eigenvalues is easily deduced from Lemma 2.12 and (2.12).

It remains to show that the corresponding eigenspace is given by $\text{span}\{\mathbf{1}\}$. This is easily seen from (3.6): for $\mathbf{v} \in \mathbb{R}^M$, it holds that $\mathbf{S}\mathbf{v} = \mathbf{v}$ if and only if $\mathcal{F}_P \mathcal{E}_P^{-1}\mathbf{v} = \mathcal{E}_P^{-1}\mathbf{v}$. On the other hand, $\mathcal{F}_P \mathbf{1} = \mathbf{1}$, and $\mathcal{E}_P \mathbf{1} = \mathbf{1}$, and thus $\mathbf{v}$ satisfying $\mathbf{S}\mathbf{v} = \mathbf{v}$ necessarily belongs to $\text{span}\{\mathbf{1}\}$. $\square$

**Remark 3.4.** The proof of Lemma 3.3 provides a more precise characterization of the spectrum $\mathcal{F}_P$, namely, with $\mathcal{F}_{k,k}$ defined in the proof of Lemma 2.12, see Eq. (2.10):

$$ k := \min\{k \in \mathbb{N} | \mathcal{P} \subset \mathbb{P}_k\} \Rightarrow \sigma(\mathcal{F}_P) \setminus \{\mathbf{1}\} \subset \bigcup_{1 \le m \le k} \sigma(\mathcal{F}_{m,m}) $$

where each $\sigma(\mathcal{F}_{k,k})$ is localized in a “small” ball centered at the origin (see Lemma 2.12).

As a consequence of Lemma 3.3, we have the following result.

**Theorem 3.5.** The problem Eq. (3.5) is well-posed and characterizes the vector $\mathbf{w}$ of cubature weights.

The proof of this result relies partially on a standard linear algebra argument, which we repeat for the convenience of the reader and which will be used later in the paper.

**Lemma 3.6.** Assume that $\mathbf{A} \in \mathbb{R}^{M \times M}$, and let $\lambda \in \mathbb{R}$ be a simple eigenvalue of $\mathbf{A}$ (i.e. its algebraic, and thus, geometric multiplicity equals to 1). Then the corresponding left and right eigenvectors $\mathbf{v}_l$ and $\mathbf{v}_r$ are not orthogonal. I.e. if $\mathbf{v}_l, \mathbf{v}_r \in \mathbb{R}^M$, $\mathbf{v}_l \neq 0, \mathbf{v}_r \neq 0$, are s.t. $\mathbf{v}_l^\top \mathbf{A} = \lambda \mathbf{v}_l^\top$ and $\mathbf{A}\mathbf{v}_r = \lambda \mathbf{v}_r$, then, necessarily, $\mathbf{v}_l \cdot \mathbf{v}_r \neq 0$.

*Proof.* We reason by contradiction. Assume that $\mathbf{v}_l \cdot \mathbf{v}_r = 0$. Next, recall the decomposition of the space $\mathbb{R}^M = \text{Ker}(\lambda - \mathbf{A}) \oplus \text{Im}(\lambda - \mathbf{A}^\top)$. Since $\mathbf{v}_r \in \text{Ker}(\lambda - \mathbf{A})$, and $\mathbf{v}_l \cdot \mathbf{v}_r = 0$, then, necessarily, $\mathbf{v}_l \in \text{Im}(\lambda - \mathbf{A}^\top)$. Therefore, there exists $\mathbf{y} \in \mathbb{R}^M$, s.t. $(\lambda - \mathbf{A}^\top)\mathbf{y} = \mathbf{v}_l$. Remark that $\mathbf{y} \notin \text{span}\{\mathbf{v}_l\}$, and is thus a generalized eigenvector. This contradicts the fact that $\lambda$ is a simple eigenvalue of $\mathbf{A}^\top$ (equivalently, of $\mathbf{A}$). $\square$

*Proof of Theorem 3.5.* Because the eigenvalues of $\mathbf{S}$ and $\mathbf{S}^\top$ coincide (together with multiplicities), by Lemma 3.3 we conclude that there exists $\mathbf{v} \neq 0$, s.t. $\mathbf{S}^\top \mathbf{v} = \mathbf{v}$, and the corresponding eigenvalue is simple. Using Lemma 3.6, we also have that $\mathbf{S}\mathbf{1} = \mathbf{1}$. Applying Lemma 3.6 with $\mathbf{A} = \mathbf{S}$, $\lambda = 1$ shows that, necessarily, $\mathbf{v} \cdot \mathbf{1} \neq 0$.

It is easy to verify that $\tilde{\mathbf{w}} = (\mathbf{v} \cdot \mathbf{1})^{-1}\mathbf{v}$ is a unique solution to the problem (3.5), which is thus well-posed. Since the exact cubature weights satisfy (3.5), necessarily, $\tilde{\mathbf{w}}$ coincides with the vector $\mathbf{w}$. $\square$

In what follows, we will need a rewriting of Theorem 3.5, which will serve us in a sequel (see Theorem 4.3).
**Corollary 3.7.** Consider the following problem:

$$ \text{Find } \omega \in \mathcal{P}, \text{ s.t. } \int_{\Gamma} \omega p \, d\mu = \int_{\Gamma} \omega \mathcal{F} p \, d\mu, \quad \text{and} \quad \int_{\Gamma} \omega \, d\mu = 1. \qquad (3.7) $$

This problem is well-posed, and the unique solution is given by $\omega = 1$.

*Proof.* The result follows from the equivalence of the problem (3.7) and the problem (3.5), in the following sense:

(A) if $w \in \mathbb{R}^M$ solves (3.5), then $\omega = (\mathcal{I}_P^*)^{-1} w$ satisfies (3.7);

(B) if $\omega \in \mathcal{P}$ satisfies (3.7), then $w = \mathcal{I}_P^* \omega$ solves (3.7),

where we recall that $\mathcal{I}_P^*$ is given by (2.16). To prove (A) and (B), we remark that

$$ \text{for all } f, p \in \mathcal{P}, \quad \int_{\Gamma} f p \, d\mu = \sum_{1 \le i \le M} p(x_i) \int_{\Gamma} f \mathcal{L}_i \, d\mu = \mathcal{I}_P^* f \cdot \mathcal{E} p. $$

As a consequence, using the $\mathcal{S}$-invariance of the space $\mathcal{P}$, it holds that

$$ \int_{\Gamma} \omega p \, d\mu = \int_{\Gamma} \omega \mathcal{F} p \, d\mu \iff \mathcal{I}_P^* \omega \cdot \mathcal{E} p = \mathcal{I}_P^* \omega \cdot \mathcal{E} \mathcal{F} p. $$

By Lemma 3.2, we rewrite the above as

$$ \int_{\Gamma} \omega p \, d\mu = \int_{\Gamma} \omega \mathcal{F} p \, d\mu \iff \mathcal{I}_P^* \omega \cdot \mathcal{E} p = \mathcal{I}_P^* \omega \cdot \mathcal{S}(\mathcal{E} p). \qquad (3.8) $$

Since $\mathcal{E}$ is surjective from $\mathcal{P}$ to $\mathbb{R}^M$, (3.8) leads to the following equivalence:

$$ \int_{\Gamma} \omega p \, d\mu = \int_{\Gamma} \omega \mathcal{F} p \, d\mu, \quad \forall p \in \mathcal{P} \iff \mathcal{I}_P^* \omega \cdot \mathbf{p} = \mathcal{I}_P^* \omega \cdot \mathbf{S} \mathbf{p}, \quad \forall \mathbf{p} \in \mathbb{R}^M. \qquad (3.9) $$

In a similar manner,

$$ \int_{\Gamma} \omega \, d\mu = 1 \iff \mathcal{I}_P^* \omega \cdot \mathcal{E} 1 = \mathcal{I}_P^* \omega \cdot \mathbf{1} = 1. \qquad (3.10) $$

This proves (A) and (B). Finally, we recognize in the right-hand side of equivalences (3.9), (3.10) the problem (3.5); the desired conclusion follows from Theorem 3.5. $\square$

**Remark 3.8.** Remarkably, given two different sets of cubature points $\mathcal{X}$ and $\mathcal{Y}$ each of them being $\mathcal{P}$-unisolvent, the spectrum of the corresponding matrices $\mathcal{S}_\mathcal{X}$ and $\mathcal{S}_\mathcal{Y}$ (which differ a priori) is the same. Indeed, these matrices are representations of the same operator $\mathcal{F}_\mathcal{P}$, cf. Eq. (3.6).

# 4 Cubature based on non-$\mathcal{S}$-invariant polynomial spaces

One motivation for looking at the case of non-$\mathcal{S}$-invariant polynomial spaces stems from a practical difficulty of choosing cubature points that would be unisolvent in $\mathbb{P}_k$, and possess Lebesgue's constants of mild growth, with respect to the number of points (this property is favorable for convergence of quadratures formulae, as we will see in Section 5.2). Therefore, instead, we will look at the tensor product spaces $\mathbb{Q}_k$, which has an additional favorable property: for well-chosen unisolvent sets of points, the associated Lagrange polynomials are easy to compute using barycentric formulae, cf. [Tre13]. It appears that the spaces $\mathbb{Q}_k$ in general are not $\mathcal{S}$-invariant, as discussed in Section 2.2.2. Therefore, a priori, the cubature weights on the space $\mathbb{Q}_k$ no longer satisfy the algebraic property Eq. (3.5). This is illustrated numerically in Section 6.
## 4.1 An alternative definition of the cubature weights

Let $\mathcal{P} \subset \mathbb{P}$ with $\dim \mathcal{P} < +\infty$ and $\mathbb{P}_0 \subset \mathcal{P}$. According to Lemma 2.19, the cubature weights $\mathbf{w}$ given by Eq. (2.18) provide a cubature formula $\mathcal{Q}[f]$ which is exact in $\mathcal{P}$. In the previous section, under the assumption that $\mathcal{P}$ is $\mathcal{S}$-invariant, we have shown that the corresponding weights can be characterized as a unique solution to the purely algebraic problem Eq. (3.5). However, this is no longer true if the spaces $\mathcal{P}$ are not $\mathcal{S}$-invariant.

Our idea is to abandon the constraint imposed by the cubature formula of being exact in $\mathcal{P}$, but instead use Eq. (3.5) as a characterization of the cubature weights, in the same way as the property (3.5) characterized exact cubature weights for $\mathcal{S}$-invariant spaces, see Theorem 3.5. More precisely, instead of the exact weights $\mathbf{w}$, defined in Lemma 2.19, we will look for weights $\tilde{\mathbf{w}}$ as a solution to the following problem:

$$ \text{Find } \tilde{\mathbf{w}} \in \text{Ker}(I - \mathbf{S}^{\top}) \text{ satisfying } \tilde{\mathbf{w}} \cdot \mathbf{1} = 1. \quad (4.1) $$

Recall that the constraint $\tilde{\mathbf{w}} \cdot \mathbf{1} = 1$ in Eq. (4.1) ensures that $\mathbb{P}_0$ is integrated exactly. This new cubature rule, defined by

$$ \int_{\Gamma} f \, d\mu \approx \tilde{\mathcal{Q}}[f] := \sum_{1 \le i \le M} \tilde{w}_i f(x_i), $$

is consistent with the invariance property Proposition 2.11 of the integral, in sense of

**Lemma 4.1.** If $\tilde{\mathcal{Q}}[f] = \sum_{1 \le i \le M} \tilde{w}_i f(x_i)$, the following two properties are equivalent:

$$ (i) \text{ for all } p \in \mathcal{P}, \tilde{\mathcal{Q}}[p] = \tilde{\mathcal{Q}}[\mathcal{F} p]; \quad (ii) \tilde{\mathbf{w}} \in \text{Ker}(I - \mathbf{S}^{\top}). $$

We leave the proof of this result to the reader. Following the above discussion, we single out two important questions about the “new” cubature weights Eq. (4.1):

* the well-posedness of the problem Eq. (4.1);

* the error committed by replacing the cubature rule $\mathcal{Q}$ with its perturbed version $\tilde{\mathcal{Q}}$ when integrating polynomials $p \in \mathcal{P}$.

## 4.2 On the well-posedness of the problem Eq. (4.1)

Unfortunately, we were able to show only a partial well-posedness result.

**Theorem 4.2.** The value $\lambda = 1$ is an eigenvalue of the matrix $\mathbf{S}^{\top}$. If this eigenvalue is simple, then $\tilde{\mathbf{w}} \cdot \mathbf{1} \neq 0$.

*Proof.* The first statement follows by remarking that $\mathbf{S}1 = \mathbf{1}$ (the result follows by the same reasoning as in the proof of the second statement of Lemma 3.3). To prove the second statement, let us assume that $\tilde{\mathbf{w}} \neq 0$ is s.t. $\mathbf{S}^{\top}\tilde{\mathbf{w}} = \tilde{\mathbf{w}}$. Applying Lemma 3.6 with $\mathbf{A} = \mathbf{S}$ and $\lambda = 1$ implies that $\tilde{\mathbf{w}} \cdot \mathbf{1} \neq 0$. $\square$

According to the proof of the above result, the condition that $\lambda = 1$ is a simple eigenvalue ensures that (4.1) is well-posed, in particular, the constraint $\tilde{\mathbf{w}} \cdot \mathbf{1} = 1$ is satisfied. One can verify that since $\mathbf{S}$ is real, $\tilde{\mathbf{w}}$ satisfying (4.1) is necessarily real, too.
## 4.3 Subspaces of $\mathcal{P}$ integrated exactly by the inexact cubature

Since the error of the cubature is closely related to the largest polynomial space integrated exactly by the cubature, it is natural to ask a question of characterization of spaces of polynomials in $\mathcal{P}$ that are integrated exactly, i.e.

$$ \mathcal{P}_Q := \{ p \in \mathcal{P} \mid \int_{\Gamma} p \, d\mu = \tilde{Q}[p] \}. \qquad (4.2) $$

Since the integral and $\tilde{Q}$ are linear forms, this space has co-dimension 1 (thus dimension $\dim \mathcal{P} - 1$) in $\mathcal{P}$ independently of the vector $\tilde{w}$ (!). A more precise characterization of the space $\mathcal{P}_Q$ is given in the theorem below.

**Theorem 4.3.** Assume that $\mathcal{S}^\top \tilde{w} = \tilde{w}$ and $\tilde{w} \cdot \mathbf{1} = 1$. Let $\mathcal{P}_\mathcal{S}$ be the largest $\mathcal{S}$-invariant subspace of $\mathcal{P}$. Then the inexact cubature formula $\tilde{Q}$ integrates polynomials in the space $\mathcal{P}_\mathcal{S}$ exactly, in other words, $\mathcal{P}_\mathcal{S} \subset \mathcal{P}_Q$.

*Proof.* First, let us argue that we can rewrite the algebraic problem Eq. (4.1) in the variational form resembling the statement of Corollary 3.7. Recall Eq. (4.1) satisfied by $\tilde{w} \in \mathbb{R}^M$, rewritten in an equivalent form, see Lemma 3.2:

$$ \tilde{w}^{\top} = \tilde{w}^{\top} \mathcal{S} \iff \tilde{w} \cdot \mathcal{E} p = \tilde{w} \cdot \mathcal{E}_P \mathcal{I}_P \mathcal{E} \mathcal{F} p, \quad \forall p \in \mathcal{P}. \qquad (4.3) $$

Remark that in the above identity, it is the operator $\mathcal{I}_P \mathcal{E} \mathcal{F}$ that appears, rather than $\mathcal{F}$. Moreover, because of non-$\mathcal{S}$-invariance of $\mathcal{P}$, $\mathcal{I}_P \mathcal{E} \mathcal{F} \neq \mathcal{F}$ on $\mathcal{P}$. Denoting $\tilde{\omega} := (\mathcal{I}_P^*)^{-1} \tilde{w}$ (cf. the proof of the Corollary 3.7), we obtain the following counterpart of the identity in the right of Eq. (4.3):

$$ \int_{\Gamma} \tilde{\omega} p \, d\mu = \int_{\Gamma} \tilde{\omega} \mathcal{I}_P \mathcal{E} \mathcal{F} p \, d\mu, \quad \forall p \in \mathcal{P}. \qquad (4.4) $$

Next, let us decompose the space $\mathcal{P}$ into two subspaces orthogonal with respect to $L^2(\Gamma)$-inner-product: $\mathcal{P} = \mathcal{P}_\mathcal{S} \oplus \mathcal{P}_\mathcal{S}^\perp$. In particular, any $v \in \mathcal{P}$ can be written as orthogonal sum $v = v_\mathcal{S} + v_\perp$, where $v_\mathcal{S} \in \mathcal{P}_\mathcal{S}$ and $\int_\Gamma v_\perp v_\mathcal{S} \, d\mu = 0$. In a similar manner we rewrite $\tilde{\omega} = \tilde{\omega}_\mathcal{S} + \tilde{\omega}_\perp$. Applying this decomposition to Eq. (4.4), and restricting the space of test functions to $\mathcal{P}_\mathcal{S}$ yields the new identity for $\tilde{\omega}_\mathcal{S}$:

$$ \int_{\Gamma} \tilde{\omega}_{\mathcal{S}} p \, d\mu = \int_{\Gamma} \tilde{\omega}_{\mathcal{S}} \mathcal{F} p \, d\mu, \quad \forall p \in \mathcal{P}_{\mathcal{S}}, \qquad (4.5) $$

where we used the fact that $\mathcal{I}_P \mathcal{E} \mathcal{F} p = \mathcal{F} p$ for all $p \in \mathcal{P}_\mathcal{S}$: since $\mathcal{P}_\mathcal{S}$ is $\mathcal{S}$-invariant, $\mathcal{E} \mathcal{F} p \in \mathcal{P}$ when $p \in \mathcal{P}_\mathcal{S}$. Moreover, from $\tilde{w} \cdot \mathbf{1} = 1$, and $\mathbb{P}_0 \subset \mathcal{P}_\mathcal{S}$, we have

$$ \int_{\Gamma} \tilde{\omega}_{\mathcal{S}} \, d\mu = 1. \qquad (4.6) $$

In (4.5), (4.6), we recognize the problem of Corollary 3.7, where $\mathcal{P}$ is replaced by $\mathcal{P}_\mathcal{S}$. Therefore, $\tilde{\omega}_\mathcal{S} = 1$. since for all $p \in \mathcal{P}_\mathcal{S}$, it holds that

$$ \tilde{Q}[p] = \tilde{w} \cdot \mathcal{E} p = \int_{\Gamma} \tilde{\omega} p \, d\mu = \int_{\Gamma} \tilde{\omega}_{\mathcal{S}} p \, d\mu = \int_{\Gamma} p \, d\mu, $$

where the second equality above follows from the definition $\tilde{\omega} = (\mathcal{I}_P^*)^{-1} \tilde{w}$, and the third is a consequence of the orthogonality, we conclude that $\mathcal{P}_\mathcal{S} \subset \mathcal{P}_Q$. $\square$

The above result may seem somewhat surprising: indeed, even if the matrix $\mathcal{S}^\top$ admits several eigenfunctions satisfying the assumptions of Theorem 4.3, the cubature rule based on any such eigenfunction is exact in $\mathcal{P}_\mathcal{S}$.
**Remark 4.4.** The above result enables us to use the inexact cubature weights to compute the weights exact in the space $\mathcal{P}_\mathcal{S}$. For example, to integrate polynomials in the space $\mathbb{Q}_k$ exactly, one could have computed the inexact weights in $\mathbb{Q}_{nk}$. Such weights would have integrated the polynomials in the space $\mathbb{P}_{nk}$ (which is $\mathcal{S}$-invariant) exactly, in particular the Lagrange polynomials for $\mathbb{Q}_k$, which define the exact cubature weights. However, this can be potentially computationally expensive.

# 5 Error estimates

In this section we will quantify the approximation error of $\int_\Gamma f \, d\mu$ using the cubature rule introduced in the previous sections. Such estimates will be obtained from quite standard arguments for cubature's error estimates and cubature on regular sets. We suggest two versions of approximation of $\int_\Gamma f \, d\mu$:

* The h-version, where, very roughly speaking, we decompose $\Gamma$ into smaller subsets $\mu(\Gamma_j) \le h$, and approximate $\int_{\Gamma_j} f \, d\mu$ by a cubature of a fixed degree $M$. Convergence will be assured by taking $h \to 0$.

* The p-version, where we approximate directly $\int_\Gamma f \, d\mu$ by the cubature of the degree $M$. Convergence is then assured by taking $M \to +\infty$.

In both cases, the error analysis relies on the following (classic) estimate, which we recall here. Let $\mathcal{P}_Q$ be defined as in (4.2). Then

$$|e[f]| \le \inf_{p \in \mathcal{P}_Q} \|f - p\|_{L^\infty(K)} [1 + |w|_1], \quad \text{where } e[f] := \int_\Gamma f \, d\mu - Q[f] \quad (5.1)$$

which follows from the fact that, for any $p \in \mathcal{P}_Q$, $e[f] = \int_\Gamma (f - p) \, d\mu - Q[f - p]$.

## 5.1 Convergence by refining the mesh

### 5.1.1 Nested meshes on the IFS

To define a mesh on $\Gamma$, it is natural to decompose it into sub-fractals using the IFS contractive maps. Let us fix the mesh size $h > 0$. The first, simplest decomposition, can be obtained from the following observation. Since $\Gamma = \mathcal{H}(\Gamma)$ is a fixed point of the Hutchinson operator $\mathcal{H}$, see Theorem 2.3, it is also a fixed point of the p-times iterated operator $\mathcal{H}^p$. This gives the following decomposition of $\Gamma$, using $m = (m_1, \dots, m_p) \in \mathbb{L}^p$:

$$\Gamma = \bigcup_{m \in \mathbb{L}^p} \Gamma_m \quad \text{where } \Gamma_m = S_{m_1} \circ \dots \circ S_{m_p}(\Gamma). \quad (5.2)$$

The maximal refinement level $p$ is chosen so that $\text{diam} \Gamma_m < h$, for all $m \in \mathbb{L}^p$. Remark that in the above, a priori $\mu(\Gamma_m \cap \Gamma_{m'}) \ne 0$ for $m \ne m' \in \mathbb{L}^p$. The use of the mesh (5.2) for the integral evaluation can be computationally inefficient, since the diameters of the patches $\Gamma_m$ in the above decomposition can vary significantly (in other words, some regions of $\Gamma$ are refined too finely). Therefore, instead we will use an approach suggested in [GHM23]. Let us introduce some useful notation (with $m \in \mathbb{L}^p$)

$$\left| \begin{array}{l} S_m := S_{m_1} \circ \dots \circ S_{m_p}, \quad \rho_m := \rho_{m_1} \dots \rho_{m_p}, \quad \mu_m := \mu_{m_1} \dots \mu_{m_p}. \\ \mathbb{L}^p = \mathbb{L}^1 \cup \dots \cup \mathbb{L}^p, \quad \text{for any } p \in \mathbb{N}^* \cup \{\infty\}. \end{array} \right. \quad (5.3)$$

For $p = 0$, by convention, we define $\mathbb{L}^0 = \{\emptyset\}$, $\Gamma_\emptyset = \Gamma$, $S_\emptyset = \text{Id}$, $\rho_\emptyset = 1$, and $\mu_\emptyset = 1$. Given $h > 0$, let us define the set of multi-indices $\mathcal{L}_h \subset \mathbb{L}^\infty$ by the following:

$$m \in \mathcal{L}_h \iff \rho_m \text{ diam } \Gamma \le h, \quad \text{and } \rho_{\ell}^{-1} \rho_m \text{ diam } \Gamma > h, \text{ for some } \ell \in \mathbb{L}, \quad (5.4)$$
and the corresponding decomposition of $\Gamma$ via

$$ \Gamma = \bigcup_{m \in \mathcal{L}_h} \Gamma_m. \qquad (5.5) $$

We will see below that these definitions and decompositions make sense. Let us remark that in (5.4), we used $\rho_m$ diam $\Gamma$ instead of diam $\Gamma_m$, because for self-affine sets diam $\Gamma_m$ is not easily computable. On the other hand, diam $\Gamma_m \le \rho_m$ diam $\Gamma$. The first condition in (5.4) thus ensures that diam $\Gamma_m < h$ for all patches $\Gamma_m$, while the second condition is a requirement of minimality of the set $\mathcal{L}_h$, in the sense that the patches $\Gamma_m$ cannot become too small.

The set $\mathcal{L}_h$ can be obtained by the following algorithm that constructs a sequence $\mathcal{L}_h^k$, with $\mathcal{L}_h^k \subset \mathbb{L}^k$, which is stationary after a finite number of steps. Then, $\mathcal{L}_h$ is defined as the stationary point of this algorithm. The algorithm proceeds as follows. Starting from $\mathcal{L}_h^0 = \{\emptyset\}$, we deduce $\mathcal{L}_h^{k+1}$ from $\mathcal{L}_h^k$ by looking at the (possibly empty) set $\mathcal{N}_h^k = \{m \in \mathcal{L}_h^k | \rho_m \text{ diam } \Gamma > h\}$ and setting

$$ \mathcal{L}_h^{k+1} = (\mathcal{L}_h^k \setminus \mathcal{N}_h^k) \cup \{(m_1, \dots, m_k, \ell) : (m_1, \dots, m_k) \in \mathcal{N}_h^k \text{ and } \ell \in \mathbb{L}\}. \qquad (5.6) $$

The reader will easily verify that 1) the stationary point $\mathcal{L}_h$ is characterized by (5.4); 2) for all $k \ge 0$, $\Gamma = \bigcup_{m \in \mathcal{L}_h^k} \Gamma_m$ (by induction on $k$), which implies (5.5).

**Remark 5.1.** Setting $\rho_{\max} = \max_{\ell \in \mathbb{L}} \rho_{\ell}$, the algorithm stops after $k_{\max}$ iteration with

$$ k_{\max} \le k_* := \inf \{k \in \mathbb{N} \mid \rho_{\max}^k \text{ diam } \Gamma \le h\} \quad \text{i.e.} \quad k_* = \left\lceil \frac{\log(h/\text{diam } \Gamma)}{\log \rho_{\max}} \right\rceil. $$

**Remark 5.2.** For similar IFS, it holds that $\rho_{\min}h < \text{diam } \Gamma_m \le h$, for all $m \in \mathcal{L}_h$.

Just like in (5.2), $\Gamma_m$ and $\Gamma_{m'}$ in (5.5) may have non-trivial intersections. However, this decomposition defines a 'partition' of $\Gamma$ in the sense of the following result.

**Lemma 5.3.** For $h > 0$ and $f \in \mathbb{L}^1(\Gamma; \mu)$, we have, for $k \ge 0$,

$$ (i) \quad \int_{\Gamma} f(x) \, d\mu = \sum_{m \in \mathcal{L}_h^k} \mu_m \int_{\Gamma} f \circ S_m \, d\mu, \quad (ii) \quad \sum_{m \in \mathcal{L}_h^k} \mu_m = 1 \qquad (5.7) $$

*Proof.* We proceed by induction. The result is true for $k = 0$ (trivial); for $k = 1$, we remark that $\mathcal{L}_h^1 = \mathbb{L}$. Then (5.7) is direct from Proposition 2.11. Next, assuming Eq. (5.7), by definition (5.6) of $\mathcal{L}_h^{k+1}$, it holds that

$$ \sum_{m \in \mathcal{L}_h^{k+1}} \mu_m \int_{\Gamma} f \circ S_m \, d\mu = \sum_{m \in \mathcal{L}_h^k \setminus \mathcal{N}_h^k} \mu_m \int_{\Gamma} f \circ S_m \, d\mu + \sum_{m \in \mathcal{N}_h^k} \mu_m \sum_{\ell \in \mathbb{L}} \mu_{\ell} \int_{\Gamma} f \circ S_m \circ S_{\ell} \, d\mu. $$

By Proposition 2.11, $\sum_{\ell \in \mathbb{L}} \mu_{\ell} \int_{\Gamma} f \circ S_m \circ S_{\ell} \, d\mu = \int_{\Gamma} f \circ S_m \, d\mu$, thus

$$ \sum_{m \in \mathcal{L}_h^{k+1}} \mu_m \int_{\Gamma} f \circ S_m \, d\mu = \sum_{m \in \mathcal{L}_h^k \setminus \mathcal{N}_h^k} \mu_m \int_{\Gamma} f \circ S_m \, d\mu + \sum_{m \in \mathcal{N}_h^k} \mu_m \int_{\Gamma} f \circ S_m \, d\mu, \\ = \sum_{m \in \mathcal{L}_h^k} \mu_m \int_{\Gamma} f \circ S_m \, d\mu, $$

which proves (5.7)(i). The equality (5.7)(ii) follows from taking $f = 1$. $\square$

<figure><img src="image_3.png" /></figure>
### 5.1.2 Mesh refinement cubature

Let $K \supset \Gamma$ be a compact of $\mathbb{R}^n$ as in Lemma 2.6, and $\mathcal{X} := \{x_i\}_{i=1}^M \subset K$ be a set of points that are $\mathcal{P}$-unisolvent on $K$. Let us assume that the weights $\{w_i\}_{i=1}^M$ are such that the cubature rule $\mathcal{Q}$ is exact on $\mathbb{P}_k \subset \mathcal{P}$ for some $k \in \mathbb{N}$. Following Lemma 5.3, we define a cubature rule $\mathcal{Q}_h$:

$$ \mathcal{Q}_h[f] := \sum_{m \in \mathcal{L}_h} \mu_m \mathcal{Q}[f \circ S_m]. \qquad (5.8) $$

We define the set

$$ K_h = \bigcup_{m \in \mathcal{L}_h} S_m(K) \text{ which satisfies } \Gamma \subset K_h \subset K \quad \text{(by Lemma 2.6).} \qquad (5.9) $$

It can be shown that $K_h$ converges to $\Gamma$ as $h \to 0$, for the Hausdorff distance. By construction, the $h$-cubature formula (5.8) uses the values of the function $f$ in the set $K_h$, rather than the whole compact $K$ (the interest of this is discussed in Remark 5.5). Accordingly, we have the following result.

**Theorem 5.4.** Let $K$ be a convex set and satisfy Lemma 2.6. Let $H > 0$, $k \in \mathbb{N}$, and let a cubature $\mathcal{Q}$ be exact on $\mathbb{P}_k$. For all $f \in \mathcal{C}^{k+1}(K_H)$ and $h \le H$, we have

$$ |\mathcal{Q}_h[f] - \int_{\Gamma} f \, d\mu| \le C h^{k+1} \|\mathcal{D}^{k+1} f\|_{L^{\infty}(K_H)}, $$

where $\|\mathcal{D}^{k+1} f\|_{L^{\infty}(K_H)} = \max_{|\beta|=k+1} \|\partial^{\beta} f\|_{L^{\infty}(K_H)}$ and the constant $C$ is defined by

$$ C = \left( \frac{\mathrm{diam} K}{\mathrm{diam} \Gamma} \right)^{k+1} \left( |\mathbf{w}|_1 + 1 \right) \sum_{|\beta|=k+1} \frac{1}{\beta!}. $$

**Remark 5.5.** As seen from the statement of the above result, the convergence of the cubature depends on the regularity of the integrand in the set $K_H$. This allows to use efficiently such an $h$-quadrature to integrate functions that have singularities inside a convex set $K$, but regular in a small vicinity of $\Gamma$, cf. Fig. 1.

*Proof of Theorem 5.4.* Let $f \in \mathcal{C}^{k+1}(K_H)$. Using Lemma 5.3 and Eq. (5.8),

$$ \left| \mathcal{Q}_h[f] - \int_{\Gamma} f \, d\mu \right| = \left| \sum_{m \in \mathcal{L}_h} \mu_m \left[ \mathcal{Q}[f \circ S_m] - \int_{\Gamma} f \circ S_m \, d\mu \right] \right| \le \max_{m \in \mathcal{L}_h} e_m \qquad (5.10) $$

where

$$ e_m = \left| \mathcal{Q}[f \circ S_m] - \int_{\Gamma} f \circ S_m \, d\mu \right|, $$

and we used (5.7)(ii). To estimate each of the terms $e_m$, we use (5.1), where we recall that the polynomials $\mathbb{P}_k$ are integrated exactly:

$$ |e_m| \le \inf_{p \in \mathbb{P}_k} \|f \circ S_m - p\|_{L^{\infty}(K)} [1 + |\mathbf{w}|_1]. \qquad (5.11) $$

To estimate the right-hand side, we use the multivariate Taylor series expansion of $f$ in the point $z_m = S_m z$, with arbitrary $z \in K$. It holds that $f = t_m + r_m$, where

$$ t_m(y) = \sum_{|\alpha|=k} \frac{\partial^\alpha f(z_m)}{\alpha!} (y - z_m)^\alpha, \quad r_m(y) = \sum_{|\beta|=k+1} h_\beta(y) (y - z_m)^\beta, \qquad (5.12) $$

$$ h_\beta(y) = \frac{k+1}{\beta!} \int_0^1 (1-t)^k D^\beta f(z_m + t(y-z_m)) \, dt. \qquad (5.13) $$
Remark that $t_m \in \mathbb{P}_k$, and thus $t_m \circ S_m \in \mathbb{P}_k$; then (5.11) implies that

$$|e_m| \le \|r_m \circ S_m\|_{L^\infty(K)}[1 + |\mathbf{w}|_1]. \qquad (5.14)$$

From (5.12) it follows that

$$\|r_m \circ S_m\|_{L^\infty(K)} \le \sum_{|\beta|=k+1} \|h_\beta\|_{L^\infty(S_m(K))} |S_m y - S_m z|_2^{k+1}, \qquad (5.15)$$

where we used the bound $|(y - z_m)^\beta| \le |y - z_m|_2^{|\beta|_1}$. Using (5.13), convexity of K and the fact that $z_m = S_m z \in K_h$, we conclude that

$$\|h_\beta\|_{L^\infty(S_m(K))} \le \frac{1}{\beta!} \|\partial^\beta f\|_{L^\infty(K_h)}. \qquad (5.16)$$

On the other hand,

$$|S_m y - S_m z|_2^{k+1} \le \rho_m^{k+1} |y - z|^{k+1} \le \left( \frac{h \operatorname{diam} K}{\operatorname{diam} \Gamma} \right)^{k+1}, \qquad (5.17)$$

where we use the fact that $m \in \mathcal{L}_h$, and thus $\rho^m \operatorname{diam} \Gamma \le h$. Combining Eqs. (5.16) and (5.17) into Eq. (5.15), inserting the resulting bound into Eq. (5.14) and using (5.10) yields the desired result. $\square$

## 5.2 Convergence by increasing the quadrature order

Alternatively, we can evaluate the cubature by computing $Q[f]$ directly, without splitting $\Gamma$ into sub-sets. In this case the convergence is achieved by increasing the number of cubature points $M$ (we will denote the corresponding cubature by $Q^M[f]$). The goal of this section is to quantify the speed of convergence with respect to $M$. Recall the estimate Eq. (5.1), where we emphasize the dependence of all quantities on $M$, by using appropriate indices:

$$|e^M[f]| \le \inf_{p \in \mathcal{P}_Q^M} \|f - p\|_{\infty} [1 + |\mathbf{w}^M|_1]. \qquad (5.18)$$

Unlike in the previous section, both factors in the above estimate will play an important role, since, in general, $|\mathbf{w}^M|_1$ can be unbounded as $M \to +\infty$. We will start our discussion with an estimate on $|\mathbf{w}^M|_1$. Next, we present convergence estimates for $\mathcal{C}^k$-functions in Theorem 5.10 and finally estimates for analytic functions in Theorem 5.11.

Many of the estimates will be made more explicit for the tensor-product cubature rule satisfying the following assumption.

**Assumption 5.6.** Assume that $K = \Pi_{a,b}$, where $\Pi_{a,b} = [a_1, b_1] \times \cdots \times [a_\ell, b_\ell]$. Let additionally $\mathcal{X} = \mathcal{X}_1^N \times \cdots \times \mathcal{X}_\ell^N$, where $\mathcal{X}_i^N$ are $(N+1)$ Chebyshev (first or second kind) quadrature nodes on $[a_i, b_i]$ [Tre13, Chapter 2 and Exercise 2.4]. The corresponding unisolvent polynomial space is the tensor product space $\mathcal{P} = \mathbb{Q}_N$ with $M = (N+1)^n$.

We denote by $\mathbf{w}_C^{N,n}$ and $\mathcal{Q}_C^{N,n}$ (the index C stands for “Chebyshev”) the respective cubature weights and cubature rule, computed with the method of this paper.

**Remark 5.7.** Remark that the cubature points depend on the choice of the compact $K$, i.e. $K = \Pi_{a,b}$.

We choose the Chebyshev points in particular because, for this case, we are able to obtain the most optimal error estimates (this is related to the behavior of the bounds on $|\mathbf{w}^M|_1$, which we connect in Section 5.2.1 to the Lebesgue constants of the points $\mathcal{X}_i^N$, and that behave asymptotically optimally for Chebyshev points).
### 5.2.1 An estimate on $|w^M|_1$ for $\mathcal{S}$-invariant spaces

In this case, the weights provided with our algebraic approach satisfy $w_i^M = \int_{\Gamma} \mathcal{L}_i d\mu$ (cf. Section 3) and Theorem 3.5. These weights may not be positive (see the experiments in Section 6.2), and therefore, even though they sum to 1, $|w^M|_1$ can still be unbounded as $M \to +\infty$. Nonetheless, this quantity grows very mildly, provided a suitable choice of the cubature points. This is summarized in the following result, whose proof is left to the reader.

**Lemma 5.8.** The exact weights $w_j^M = \int_{\Gamma} \mathcal{L}_j d\mu$ satisfy the following bound:

$$|w^M|_1 = \sum_{1 \le j \le M} |w_j^M| \le \Lambda_X^1(\Gamma) := \sup_{x \in \Gamma} \sum_{1 \le i \le M} |\mathcal{L}_i(x)|,$$

where $\Lambda_X^1(\Gamma)$ is the Lebesgue's constant for $\mathcal{X}$.

In particular, for cubature based on Assumption 5.6, we have the following bound.

**Proposition 5.9.** For the cubature as in Assumption 5.6, there exists $C(n) > 0$, s.t.

$$|w_{C}^{N,n}|_{1} \leq C(n) \log^{n}(N+1), \quad \text{for all } N \geq 1.$$

*Proof.* We use the bound of Lemma 5.8. In our setting Lagrange polynomials are tensor products of univariate Lagrange polynomials, thus

$$\Lambda_{\mathcal{X}}^{\Gamma} \leq \Lambda_{\mathcal{X}}^{K} \leq \prod_{1 \leq i \leq n} \Lambda_{\mathcal{X}_i}^{[a_i, b_i]},$$

where $\Lambda_{\mathcal{X}_i}^{[a_i, b_i]}$ is a Lebesgue constant for interpolation in $\mathcal{X}_i$ on $[a_i, b_i]$, cf. [Tre13, Ch. 15]. It remains to apply to the above the exact asymptotic of Lebesgue's constant of the Chebyshev points, cf. [Bru97] and references therein. $\square$

### 5.2.2 Error analysis for $\mathcal{C}^k$-functions

The goal of this section is to provide error estimates for the quadrature rule $Q^M[f]$, provided that $f \in \mathcal{C}^k(K)$. Recall the error bound (5.18). To obtain convergence estimates, it remains to quantify

$$|e_{\mathcal{P}}^{M}[f]| := \inf_{p \in \mathcal{P}_{\mathcal{Q}}^{M}} \|f - p\|_{L^{\infty}(K)}, \qquad (5.19)$$

where $\mathcal{P}_{\mathcal{Q}}^M = \mathcal{P}_{\mathcal{Q}}$ is defined in (4.2). The following estimate is fairly well-known.

**Theorem 5.10.** Let $K$ be a convex compact set in $\mathbb{R}^n$. If $f \in \mathcal{C}^k(K)$, and if $\mathbb{P}_N \subset \mathcal{P}_{\mathcal{Q}}^M$, for some $N \ge 1$, then there exists a constant $C(K, n, k)$ independent of $f$ and $N$, such that

$$|e_{\mathcal{P}}^{M}[f]| \le \frac{C(K, n, k)}{N^k} \sum_{|\alpha|=k} \|\partial^{\alpha} f\|_{L^{\infty}(K)}. \qquad (5.20)$$

*Proof.* Please see the multivariate version of Jackson's theorem as stated in [Sch69, Thm. 4.10] (see p. 168 of [Sch69] for the justification that int $K$ is a regular set, as required by [Sch69, Thm. 4.10]), cf. as well [BBL02, Thm. 2] for a similar result. In the first reference, the result is formulated for $\mathbb{Q}_N$ instead of $\mathbb{P}_N$; nonetheless, it is still valid for $\mathbb{P}_N$ (with a different constant compared to [Sch69, Thm. 4.10]), since $\mathbb{Q}_{[N/n]} \subset \mathbb{P}_N$. $\square$

Combining the estimates of Theorem 5.10 and Lemma 5.8 yields the following result.
**Theorem 5.11.** Let the cubature rule be like in Assumption 5.6. Assume that $f \in \mathcal{C}^k(K)$. Then, there exists $C(f, K, n, k) > 0$, s.t. for all $N \ge 1$,

$$ \left| Q_{C}^{N,n}[f] - \int_{\Gamma} f \, d\mu \right| \le C(f, K, n, k) \left( 1 + |\mathbf{w}_{C}^{N,n}|_{1} \right) N^{-k}. \quad (5.21) $$

In the case when $\mathcal{P}_C^{N,n}$ are $\mathcal{S}$-invariant for all $N \ge 0$, the above estimate yields

$$ \left| Q_{C}^{N,n}[f] - \int_{\Gamma} f \, d\mu \right| \le c \cdot C(f, K, n, k) N^{-k} \log^n(N+1), \quad \text{with some } c > 0. \quad (5.22) $$

*Proof.* Both results follow from the expression Eq. (5.18) combined with Theorem 5.10. The estimate in the $\mathcal{S}$-invariant case is a corollary of Proposition 5.9. $\square$

Compared to the error estimates of the tensor-product cubatures on the domains $\Pi_{i=1}^n[a_i, b_i]$, we see that in the $\mathcal{S}$-invariant case our estimates are worse by a factor of $\log^n(N+1)$, due to a potential growth of cubature weights.

### 5.2.3 Error analysis for analytic functions

In this section, we fully restrict our attention to cubatures satisfying Assumption 5.6. We will also work with analytic functions. This is of interest when evaluating regular integrals in boundary element methods, e.g. $\int_{\Gamma_m} \int_{\Gamma_n} G(x, y) d\mu_x d\mu_y$, where G is the fundamental solution of the Helmholtz equation. The ideas and definitions that we present here are now standard, and we follow the exposition in [SS11, Sec. 5.3.2.2], which we simplify and adapt to our setting. We will work with componentwise analytic functions, as defined below.

**Definition 5.12.** The function $f: \Pi_{a,b} \to \mathbb{C}$ is componentwise analytic, if there exists $\zeta > 1$, such that, for each $i$, and for all $\mathbf{y}_i := (y_1, \dots, y_{i-1}, y_{i+1}, \dots, y_n)$ with $y_\ell \in [a_\ell, b_\ell]$, the function

$$ f_{\mathbf{y}_i}(t) := f(y_1, \dots, y_{i-1}, t, y_{i+1}, \dots, y_n): [a_i, b_i] \to \mathbb{C}, $$

admits an analytic extension $F_{\mathbf{y}_i}(z)$ in the Bernstein ellipse

$$ \mathcal{E}_{a_i, b_i}^{\zeta} = \left\{ \frac{a_i + b_i}{2} + \frac{a_i - b_i}{2} \frac{z + z^{-1}}{2} : z \in B_{\zeta}(0) \right\}, $$

and is additionally continuous in the closure of this ellipse.

Next, to estimate Eq. (5.1), in particular, the expression (5.19), we will choose a particular polynomial from $\mathcal{P}_{\mathbb{Q}}^M$. To construct it, let us define a Chebyshev weighted L$^2$-projection operator of degree $m$ in the direction $i$, which is helpful because we work with tensorized cubature rules. In particular, let $\{T_k(t)\}_{k=0}^\infty$, $T_k(t) = \cos(k \arccos t)$, $t \in [-1, 1]$, be the set of Chebyshev polynomials of the first kind, see [Tre13, Ch. 3 and Thm. 3.1] for the relevant discussion. For each given $\mathbf{y}_i$ as in Definition 5.12, we set

$$ \left| \begin{array}{ll} \Pi_i^{(m)} : (\mathcal{C}(K), \|\cdot\|_{\infty}) \to (\mathcal{C}(K), \|\cdot\|_{\infty}), & \Pi_i^{(m)} f(\mathbf{y}) := \sum_{k=0}^{m} a_{\mathbf{y}_i, k} T_k(\mathbf{y}_i), \\ a_{\mathbf{y}_i, k} := c_k \int_{-1}^{1} T_k(t) f_{\mathbf{y}_i} \left( \frac{b_i - a_i}{2} t + \frac{a_i + b_i}{2} \right) (1 - t^2)^{-1/2} \, dt, & \end{array} \right. $$

where $c_0 = 2/\pi$ and $c_k = 1/\pi$ for $k > 0$. Let us introduce the product operator

$$ \Pi^{(m)} f := \Pi_1^{(m)} \Pi_2^{(m)} \cdots \Pi_n^{(m)} f \in \mathbb{Q}_m. \quad (5.23) $$

We shall use this operator to prove the following approximation result.
**Proposition 5.13.** Assume that $f$ is as in Definition 5.12. Assume that $\mathbb{Q}_N \subset \mathcal{P}_Q^M$, $N \ge 1$. Then there exists $C_f > 0$, such that

$$e_{\mathcal{P}}^M = \inf_{p \in \mathcal{P}_Q^M} \|f - p\|_{L^\infty(K)} \le C_f \log^{n-1}(N+1) \zeta^{-N}/(\zeta-1).$$

*Proof.* We start by bounding $e_{\mathcal{P}}^M \le e_N := \|f - \Pi^{(N)}f\|_{L^\infty(K)}$. Since $f - \Pi^{(N)}f = f - \Pi_1^{(N)}f + \Pi_1^{(N)}f - \Pi_1^{(N)}\Pi_2^{(N)}\cdots\Pi_n^{(N)}f$, by the triangle inequality,

$$e_N \le \|f - \Pi_1^{(N)}f\|_{L^\infty(K)} + \|\Pi_1^{(N)}\| \|f - \Pi_2^{(N)}\cdots\Pi_n^{(N)}f\|_{L^\infty(K)}.$$

Repeating this procedure $n-2$ times, with different operators $\Pi_k^{(N)}$, yields the bound:

$$e_N \le \|f - \Pi_1^{(N)}f\|_{L^\infty(K)} + \sum_{k=1}^{n-1} \left( \prod_{\ell=1}^{k} \|\Pi_\ell^{(N)}\| \right) \|f - \Pi_{k+1}^{(N)}f\|_{L^\infty(K)}.$$

One concludes with the help of the following bounds

$$\|\Pi_i^{(N)}\| \le C \log(N+1) \quad \text{and} \quad \|f - \Pi_{k+1}^{(N)}f\|_{L^\infty(K)} \le C_f \zeta^{-N}/(\zeta-1),$$

easily deduced from existing 1D bounds in [Pow67; Mas80] and [Tre13, Thm. 8.2] respectively. $\square$

Combining Proposition 5.13 and Lemma 5.8 allows to quantify the cubature error.

**Theorem 5.14.** Assume that $f$ is like in Definition 5.12 with $\zeta > 1$. In the situation of Assumption 5.6, there exists a constant $C_f > 0$, such that the following holds true. If the spaces $\mathbb{Q}_N$, $N \ge 1$, are $\mathcal{S}$-invariant, then,

$$\left| \mathcal{Q}_C^{N,n}[f] - \int_{\Gamma} f \, d\mu \right| \le \frac{C_f \zeta^{-N}}{\zeta-1} \log^{2n-1}(N+1). \quad (5.24)$$

Otherwise,

$$\left| \mathcal{Q}_C^{N,n}[f] - \int_{\Gamma} f \, d\mu \right| \le \frac{C_f \zeta^{-\lfloor N/n \rfloor}}{\zeta-1} \left( 1 + \left| \boldsymbol{w}_C^{N,n} \right|_1 \right) \quad (5.25)$$

*Proof.* Assume that $\mathbb{Q}_N$, $N \ge 1$, are $\mathcal{S}$-invariant. Then (5.24) ifollows from Eq. (5.18), Proposition 5.9 and Proposition 5.13. If $\mathbb{Q}_N$ is not $\mathcal{S}$-invariant, by Theorem 4.3, the respective cubature rule integrates the polynomials in $\mathbb{P}_N$ exactly. Because $\mathbb{Q}_{\lfloor N/n \rfloor} \subset \mathbb{P}_N$, we can again combine Eq. (5.18) and Proposition 5.13 to get (5.25). $\square$

**Remark 5.15.** In the $\mathcal{S}$-invariant case, compared to the multivariate integration on product domains based on the tensor-product Gauss cubature (cf. [SS11, Thm. 5.3.13]), our cubature error is slightly worse, by a factor $\log^{2n-1}(N+1)$.

Moreover, according to the estimate (5.25), the non-$\mathcal{S}$-invariant case might suffer from a deteriorated convergence rate. We think that this estimate is non-optimal, due to the technique of proof based on $\mathbb{Q}_N$ spaces. We believe that one could recover, under stronger assumptions on $f$ than the ones in Definition 5.12, the same convergence rate as for the $\mathcal{S}$-invariant case by adapting the approach of [Tre17] to $\mathbb{P}_N$ spaces. This approach uses a more elaborate theory of analytic functions of several variables, which is out of scope of this paper. See also the numerical results of Section 6.3.

# 6 Numerical experiments

In this section we present several numerical experiments, illustrating statements of different results in our paper, as well as performance of the new cubature rule.
## 6.1 Algorithmic realization and implementation

We implemented the cubature rule following Assumption 5.6, as a Julia code, see [Moi24]. In all experiments we used Chebyshev points of the first kind.

**Finding K** Finding a bounding box for an IFS is a subject of research [CC03; Mar09]. The bounding rectangle $K = \prod_{i=1}^{n}[a_i, b_i]$ is either explicitly known (for “classic” fractals like Cantor dust), or can be computed as the smallest hypercube containing the fractal. To approximate the latter, we used a minimization procedure applied to Eq. (2.3), with $|\cdot|_2$ replaced by $|\cdot|_{\infty}$-norm.<fnref n="2" />

**Computation of cubature weights** To compute the cubature weights based on the algebraic characterization Eq. (3.5), we need to evaluate the entries of the matrix $S$ (Lagrange polynomials), and solve the corresponding eigenvalue problem. The Lagrange polynomials were evaluated using the barycentric interpolation formula from [Tre13, Ch. 5]. Recall solving the eigenvalue problem in the $\mathcal{S}$-invariant case amounts to finding the eigenvector of $S^T$ corresponding to the largest (in modulus) eigenvalue $\lambda = 1$, cf. Lemma 3.3. Thus we used the power iteration method, which converges geometrically fast. We used the same method in the non-$\mathcal{S}$-invariant case. Despite the absence of theory, in all our numerical experiments, $\lambda = 1$ appeared to be the largest in modulus and a simple eigenvalue.

**Reference solution and errors** Reference values for the integrals were always computed using a highly-refined $h$-version of the method of order $O(h^{15})$, as described in Section 5.1.1. The depicted errors are relative errors.

**Self-similar sets and measures** Many of our calculations were performed on different Vicsek-type fractals, defined on $\mathbb{R}^2$ using 5 contracting maps:

$$S_0x = R_\theta x, \quad S_\ell x = \rho x + (1 - \rho)c_\ell, \quad \ell = 1, \dots, 4, \quad b_\ell = (\pm 1, \pm 1),$$

where $\rho = 1/3$, $R_\theta$ is a rotation by angle $\theta$, see also Fig. 2. We use the self-similar measure with $\mu_\ell = \frac{1}{5}$, $\ell \in \mathbb{L}$. Remark that for $\theta \in \{0, \frac{\pi}{2}\}$, the polynomial space $\mathbb{Q}_N$ is $\mathcal{S}$-invariant, while this is

<figure><img src="image_4.png" /><figcaption>Figure 2: Three Vicsek prefractals.</figcaption></figure>

not the case for other values of $\theta \in (0, \frac{\pi}{2})$.

In some experiments, we worked with the 1D Cantor set, where we used the Hausdorff measure, see Remark 2.9.

<footnote n="2">²We used the Julia package Optim.jl [MR18] for this. This is justified by Lemma 2.6; remark that a ball in $|\cdot|_{\infty}$-norm is exactly a hypercube.</footnote>
## 6.2 Behavior of cubature weights

In our first experiment, we compute cubature weights for two examples of fractals: Cantor set ($S_1(x) = \rho x$ and $S_2(x) = \rho x + (1 - \rho)$) and Viscek fractal with $\theta = 0$, see Fig. 3. First, the

<figure><img src="image_5.png" /><figcaption>Figure 3: Left figure: the values of weights (y-axis) vs the quadrature points (x-axis). The intersection of vertical lines with $y = 0$ shows an approximate location of the Cantor set. Right figure: location of cubature points in the 2D-set $K$ used to construct a quadrature on the Viscek fractal. Values of corresponding quadrature weights are indicated in color.</figcaption></figure>

weights can be negative. Second, they seem to be larger in modulus for cubature points close to $\Gamma$ as one expects intuitively.

To illustrate that negative weights do not seem to affect the quadrature accuracy for large $M$, see (5.18), we plot $|w^M|_1$ in Fig. 4 as a function of $M$.

<figure><img src="image_6.png" /><figcaption>Figure 4: The dependence of the norm $|w^M|_1$ on $M$.</figcaption></figure>

## 6.3 Convergence of the h-version and p-version

We illustrate in Fig. 5 the statements of Theorem 5.4 and 5.14 by evaluating $Q[f]$ for

$$f(x) = e^{i\kappa|x-x_0|_2}/|x-x_0|_2, \text{ for } \kappa = 5 \text{ and } x_0 = (0.1, -2), \qquad (6.1)$$

which is used in wave scattering applications. Remark that the location of the singularity $x_0$ is chosen outside $K$.

As expected, cf. Theorem 5.4, since the cubature is exact in both $\mathcal{S}$-invariant and $\mathcal{S}$-non-invariant cases, the results given by h-version do not differ significantly in these two cases. For the p-version we also observe that the convergence rates for non-$\mathcal{S}$-invariant case seem to be the
same as in the $\mathcal{S}$-invariant case. This seems to confirm that our estimate (5.25) in Theorem 5.14 is not optimal (see Remark 5.15).

<figure><img src="image_7.png" /><figcaption>Figure 5: Convergence for Viscek without and with rotation and $\rho = 1/3$.</figcaption></figure>

## 6.4 Integrating almost singular functions

In this section, we perform the experiments for a 2D Cantor dust (a tensor product of the Cantor sets with $\rho = 1/3$), located inside $K = [-1, 1]^2$. We work with the Hausdorff measure, cf. Remark 2.9.

We compute integrals of the function $f$ for $\kappa = 5$, see (6.1), with different choices of $y$ in $x_0 = (0.1, y)$. For $y \in (-2, -1)$, $x_0$ is located outside $K$ and when $y$ approaches $-1$ from below, $x_0$ approaches $K$, meaning that the integrand $f$ becomes almost singular. For $y = -1$, $x_0$ lies on the boundary of $K$, and for $y = 0$, $x_0$ is inside $K$.

When analyzing the error of the $h$-version of the cubature, we remark that for the two latter cases, $f$ is not continuous in $K$, but for $H$ small enough, $f$ is smooth in the domain $K_H \subsetneq K$ defined in (5.9). We expect the cubature (5.8) to converge at the maximal rate, cf. Remark 5.5. This is confirmed by the numerical results in Fig. 6, right.

On the other hand, the convergence analysis of the $p$-version relies on the smoothness of $f$ inside $K$ (and its analyticity properties), cf. Theorem 5.11 and 5.14. Thus, we expect the convergence to deteriorate as $x_0$ approaches $K$. This is confirmed by numerical results in Fig. 6, left, where the worst results are observed for $y \ge -1.25$.

## 6.5 Integrating over more exotic attractors of IFS

This section is dedicated to computation of integrals on less classic fractals. We approximate $\int_{\Gamma} f d\mu$ for $f$ defined in (6.1), over four different self-affine sets: the fat Sierpiński triangle, the Barnsley fern, the Koch snowflake (Example 2.5), as well as a non-symmetric Cantor dust<fnref n="3" />. The fat Sierpiński triangle does not satisfy an open set condition. The IFS of the Barnsley fern is not similar, and, moreover, for one of the matrices $A_\ell$, Ker $A_\ell \ne \{0\}$. The IFS for the Koch snowflake and the non-symmetric Cantor dust contain transformations with non-trivial rotations; moreover, the Hausdorff dimension of the Koch snowflake is $d = 2$. The invariant measures for these examples were chosen as follows: for the fat Sierpiński triangle, $\mu_\ell = 1/3$, $\ell \in \mathbb{L}$; for the Koch snowflake we choose $\mu_\ell = \rho_\ell^2$, $\ell \in \mathbb{L}$; for the Barnsley fern $\mu_\ell$'s have been chosen according

<footnote n="3">³$A_\ell = \rho_\ell R \theta_\ell$ and $b_\ell = (1 - \rho_\ell)c_\ell$ with $(\rho_\ell, \theta_\ell, c_\ell) = (0.25, 0.4, (-1.4, -1.1))$, $(0.35, 0.2, (0.8, -0.7))$, $(0.3, 0.3, (1.2, 1.3))$, $(0.4, 0.1, (-1.3, 0.9))$.</footnote>
<figure><img src="image_8.png" /><figcaption>Figure 6: Comparison of the *p*- and *h*-versions. In the right plot, we mark with the same colors the curves corresponding to different $x_0$ but the same value of *p* in the space $Q_p$.</figcaption></figure>

to the Table [Bar93, Tb. 3.8.3]; finally, for the non-symmetric Cantor dust, we again choose $\mu_\ell = \rho_\ell^d$, $\ell \in \mathbb{L}$, with $d$ solving $\sum_{\ell \in \mathbb{L}} \rho_\ell^d = 1$, cf. Remark 2.9.

The results for the *h*- and *p*- versions of the cubature for the Barnsley fern are shown in Fig. 7, and for the rest of the fractals in Fig. 8. We see that the *h*-version converges quite neatly for all

<figure><img src="image_9.png" /><figcaption>Figure 7: Convergence for the Barnsley fern. Note the difference of the scales in both figures.</figcaption></figure>

the fractals in question. At the same time, the *p*-version performs slightly worse for the Cantor dust and the Barnsley fern. We do not know a precise reason for this (perhaps this is related to the size of the computed bounding boxes which are quite large in this case).

# A An extension of the Strichartz algorithm for polynomial integration

In this section we generalize the method of [Str99] for computing integrals of polynomials on self-affine sets, defined by the IFS $\mathcal{S} = \{S_\ell | \ell \in \mathbb{L}\}$ via Theorem 2.3, with respect to the invariant
<figure><img src="image_10.png" /><figcaption>Figure 8: Convergence for the fat Sierpinski gasket, Koch snowflake and Cantor dust.</figcaption></figure>

measure $\mu$, defined via Theorem 2.8. Recall the invariance property Proposition 2.11:

$$ \int_{\Gamma} f(x) d\mu = \sum_{\ell \in \mathbb{L}} \mu_{\ell} \int_{\Gamma} f \circ S_{\ell}(x) d\mu = \int_{\Gamma} \mathcal{F}[f](x) d\mu, \quad \forall f \in \mathcal{C}^{0}(\Gamma). $$

We want to compute the moments $m_{\alpha}$ defined as the integrals of the monomials $x^{\alpha}$. This is done by induction on the total degree. Let $k \ge 1$, assume that $m_{\beta}$ is known for $|\beta|_1 < k$. For $|\alpha|_1 = k$, to the triangular structure of $\mathcal{F}$, see (2.10), we have

$$ m_{\alpha} = \int_{\Gamma} x^{\alpha} d\mu = \int_{\Gamma} \mathcal{F}[x^{\alpha}] d\mu = \int_{\Gamma} \mathcal{F}_{k,k}[x^{\alpha}] d\mu + \sum_{k'<k} \int_{\Gamma} \mathcal{F}_{k,k'}[x^{\alpha}] d\mu. $$

This can be rewritten as a linear system with unknown vector $M_k = (m_{\alpha})_{|\alpha|=k}$ of the form $(I - F_k)M_k = R_k$ where $F_k$ is the matrix of the operator $\mathcal{F}_{k,k}$ on the monomials basis $(x^{\alpha})_{|\alpha|=k}$ and the right hand side $R_k$ is known since $\mathcal{F}_{k,k'}[x^{\alpha}]$ has total degree $k' < k$. Lemma 2.12 implies that 1 is not an eigenvalue of $F_k$ and therefore the matrix $I - F_k$ is invertible, and $U_k = (I - F_k)^{-1}R_k$.<fnref n="4" />

## Acknowledgements

We are grateful to David Hewett (UCL) and Andrew Gibbs (UCL) for fruitful discussions. Additionally, we thank Andrew Gibbs and Simon Chandler-Wilde (U. Reading) for suggesting the numerical experiment in Section 6.4.

We acknowledge the use of ChatGPT to correct occasional grammar and style errors (in particular, in the introduction).

## References

[AST06] Y. Achdou, C. Sabot, and N. Tchou. “A multiscale numerical method for Poisson problems in some ramified domains with a fractal boundary”. In: *Multiscale Model. Simul.* 5.3 (2006). ISSN: 1540-3459,1540-3467.

[BBL02] T. Bagby, L. Bos, and N. Levenberg. “Multivariate simultaneous approximation”. In: *Constr. Approx.* 18.4 (2002). ISSN: 0176-4276,1432-0940.

<footnote n="4">⁴This method has been implemented, see folder symbolic-polynomial-integration in [Moi24], using Python and SymPy [Meu+17].</footnote>
[Ban07] L. Banjai. "Eigenfrequencies of fractal drums". In: *J. Comput. Appl. Math.* 198.1 (2007). ISSN: 0377-0427,1879-1778.

[BB21] L. Banjai and L. Boulton. "Computation of sharp estimates of the Poincaré constant on planar domains with piecewise self-similar boundary". In: *J. Fractal Geom.* 8.2 (2021). ISSN: 2308-1309,2308-1317.

[BGH22] J. Bannister, A. Gibbs, and D. P. Hewett. "Acoustic scattering by impedance screens/cracks with fractal boundary: well-posedness analysis and boundary element approximation". In: *M³AS*. 32.2 (2022). ISSN: 0218-2025,1793-6314.

[Bar93] M. F. Barnsley. *Fractals everywhere*. 2nd ed. Revised with the assistance of and with a foreword by Hawley Rising, III. Academic Press Professional, Boston, MA, 1993. ISBN: 0-12-079061-0.

[Bog07] V. I. Bogachev. *Measure theory. Vol. I and II*. Berlin: Springer, 2007. ISBN: 9783540345145.

[BMS04] D. Broomhead, J. Montaldi, and N. Sidorov. "Golden gaskets: variations on the Sierpiński sieve". In: *Nonlinearity* 17.4 (2004). ISSN: 0951-7715.

[Bru97] L. Brutman. "Lebesgue functions for polynomial interpolation – a survey". In: *Ann. Numer. Math.* 4.1-4 (1997). ISSN: 1021-2655.

[Cae+24] A. M. Caetano et al. "A Hausdorff-measure boundary element method for acoustic scattering by fractal screens". In: *Numer. Math.* 156.2 (2024). ISSN: 0029-599X,0945-3245.

[Cap10] R. Capitanelli. "Asymptotics for mixed Dirichlet-Robin problems in irregular domains". In: *J. Math. Anal. Appl.* 362.2 (2010). ISSN: 0022-247X,1096-0813.

[Cef+21] M. Cefalo et al. "Approximation of 3D Stokes flows in fractal domains". In: *Fractals in Engineering: Theoretical Aspects and Numerical Approximations*. SEMA SIMAI Springer Ser. Springer International Publishing, 2021. ISBN: 9783030618032.

[CH18] S. N. Chandler-Wilde and D. P. Hewett. "Well-posed PDE and integral equation formulations for scattering by fractal screens". In: *SIAM J. Math. Anal.* 50.1 (2018). ISSN: 0036-1410,1095-7154.

[CHM17] S. N. Chandler-Wilde, D. P. Hewett, and A. Moiola. "Sobolev spaces on non-Lipschitz subsets of $\mathbb{R}^n$ with application to boundary integral equations on fractal screens". In: *Integral Equations and Operator Theory* 87.2 (2017). ISSN: 0378-620X,1420-8989.

[Cha+21] S. N. Chandler-Wilde et al. "Boundary element methods for acoustic scattering by fractal screens". In: *Numer. Math.* 147.4 (2021). ISSN: 0029-599X,0945-3245.

[CC03] H.-T. Chu and C.-C. Chen. "On bounding boxes of iterated function system attractors". In: *Comput. Graph.* 27.3 (2003). ISSN: 0097-8493.

[DR22] A. Dekkers and A. Rozanova-Pierrat. "Dirichlet boundary valued problems for linear and nonlinear wave equations on arbitrary and fractal domains". In: *J. Math. Anal. Appl.* 512.1 (2022). ISSN: 0022-247X,1096-0813.

[DRT22] A. Dekkers, A. Rozanova-Pierrat, and A. Teplyaev. "Mixed boundary valued problems for linear and nonlinear wave equations in domains with fractal boundaries". In: *Calc. Var. Partial Differential Equations* 61.2 (2022). ISSN: 0944-2669,1432-0835.

[Fal14] K. Falconer. *Fractal geometry. Mathematical foundations and applications*. 3rd ed. Hoboken, NJ: John Wiley & Sons, 2014. ISBN: 978-1-119-94239-9.

[FL99] A. H. Fan and K.-S. Lau. "Iterated function system and Ruelle operator". In: *J. Math. Anal. Appl.* 231.2 (1999). ISSN: 0022-247X,1096-0813.

[FMV98] B. Forte, F. Mendivil, and E. R. Vrscay. ""Chaos games" for iterated function systems with grey level maps". In: *SIMA* 29.4 (1998). ISSN: 0036-1410.
[GHM23] A. Gibbs, D. Hewett, and A. Moiola. "Numerical quadrature for singular integrals on fractals". In: *Numer. Algorithms* 92.4 (2023). ISSN: 1017-1398,1572-9265.

[GHM24] A. Gibbs, D. P. Hewett, and B. Major. "Numerical evaluation of singular integrals on non-disjoint self-similar fractal sets". In: *Numer. Algorithms* 97.1 (2024). ISSN: 1017-1398,1572-9265.

[Hut81] J. E. Hutchinson. "Fractals and self similarity". In: *Indiana Univ. Math. J.* 30.5 (1981). ISSN: 0022-2518.

[Jor06] T. Jordan. "Dimension of fat Sierpiński gaskets". In: *Real Anal. Exch.* 31.1 (2006). ISSN: 0147-1937.

[Man96] G. Mantica. "A stable Stieltjes technique for computing orthogonal polynomials and Jacobi matrices associated with a class of singular measures". In: *Constr. Approx.* 12.4 (1996). ISSN: 0176-4276,1432-0940.

[Mar09] T. Martyn. "The smallest enclosing disc of an affine IFS fractal". In: *Fractals* 17.3 (2009). ISSN: 0218-348X.

[Mas80] J. C. Mason. "Near-best multivariate approximation by Fourier series, Chebyshev series and Chebyshev interpolation". In: *J. Approx. Theory* 28.4 (1980). ISSN: 0021-9045,1096-0430.

[Meu+17] A. Meurer et al. "SymPy: symbolic computing in Python". In: *PeerJ Comput. Sci.* 3 (2017). ISSN: 2376-5992.

[MR18] P. K. Mogensen and A. N. Riseth. "Optim: A mathematical optimization package for Julia". In: *J. Open Source Softw.* 3.24 (2018).

[Moi24] Z Moitier. *IFSCubature*. Version v1.0.0. 2024. URL: https://github.com/zmoitier/IFSCubature.

[Pow67] M. J. D. Powell. "On the maximum errors of polynomial approximations defined by interpolation and by least squares criteria". In: *Comput. J.* 9 (1967). ISSN: 0010-4620.

[RV13] M. Rosaria Lancia and P. Vernole. "Semilinear fractal problems: approximation and regularity results". In: *Nonlinear Anal.* 80 (2013). ISSN: 0362-546X,1873-5215.

[RS24] F. Rösler and A. Stepanenko. "Computing eigenvalues of the Laplacian on rough domains". In: *Math. Comp.* 93.345 (2024). ISSN: 0025-5718,1088-6842.

[SS11] S. A. Sauter and C. Schwab. *Boundary element methods*. Vol. 39. Springer Ser. Comput. Math. Berlin: Springer, 2011. ISBN: 978-3-540-68092-5; 978-3-540-68093-2.

[Sch69] M. H. Schultz. "*L*∞-multivariate approximation theory". In: *SINUM* 6 (1969). ISSN: 0036-1429.

[Str99] R. S. Strichartz. "Analysis on fractals". In: *Not. Am. Math. Soc.* 46.10 (1999). ISSN: 0002-9920.

[Tre13] L. N. Trefethen. *Approximation theory and approximation practice*. Vol. 128. Society for Industrial and Applied Mathematics, 2013. ISBN: 978-1-611972-39-9.

[Tre17] L. N. Trefethen. "Multivariate polynomial approximation in the hypercube". In: *Proc. Amer. Math. Soc.* 145.11 (2017). ISSN: 0002-9939,1088-6826.