{'fraction_non_alphanumeric': 0.04576659038901602, 'fraction_numerical': 0.009010297482837529, 'mean_word_length': 3.802312271062271, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 32, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Numerics and Fractals

Local iterated function systems are an important generalisation of the standard (global) iterated function systems (IFSs). For a particular class of mappings, their fixed points are the graphs of local fractal functions and these functions themselves are known to be the fixed points of an associated Read-Bajactarevi\'c operator. This paper establishes existence and properties of local fractal functions and discusses how they are computed. In particular, it is shown that piecewise polynomials are a special case of local fractal functions. Finally, we develop a method to compute the components of a local IFS from data or (partial differential) equations.


Introduction.
Contractive operators on function spaces are important for the development of both the theory and algorithms for the solution of integral and differential equations. They are used in the theory of elliptic partial differential equations, Fredholm integral equations of the second kind, Volterra integral equations, and ordinary differential equations. This is just a small selection of instances were they appear in mathematics. Contractive operators are fundamental for the development of iterative solvers in general and wavelet-based solvers for elliptic problems [11] in particular.
One class of contractive operators is defined on the graphs of functions using a special kind of iterated function system (IFS). The fixed point of such an IFS is the graph of a fractal function. There is a vast literature on IFSs, see for example the recent review by the first author [1]. Computationally, IFSs are used in Computer Graphics in refinement methods which effectively compute points on curves and surfaces [8]. They are also used to compute function values of piecewise polynomial functions and wavelets. In fact, it can be shown that these applications use a variant of IFSs where the iterated functions are defined locally [6]. These local IFSs and, in particular, their computational application are the topic of the following discussion. In this first manuscript we will mostly consider functions of one real variable in the examples. Functions of multiple variables are planned to be covered in a future paper.
The remaining part of this introduction will provide some further background and motivation for our approach to utilise IFSs or local IFSs in computations. In the second section we introduce and review local IFSs. The third section applies local IFSs to graphs of functions to define local fractal functions. It will be seen that these functions are the fixed points of a Read-Bajactarević (RB) operator. (See also [20] for the use of such operators in the theory of (global) fractal functions.) Section 4 provides a reformulation of the RB operator in terms of matrices acting on vectors of function values over grids. Several examples of local fractal functions are then displayed. In Section 5 we discuss the important case of polynomials and their RB operators. In a penultimate section we discuss the determination of (approximate) iterated function systems both from data and from functional equations such as partial differential equations. We conclude this discussion with some general remarks and in particular with a connection between fractals and the active research area of tensor approximation.
1.1. Fractals and numerics. One can show that graphs of piecewise polynomial functions can be written as the fixed points of local IFSs. Thus the popular finite element method approximates solutions of PDEs with particular fractal functions. However, numerical methods do not usually use IFSs directly. Exceptions are the subdivision schemes used in computer graphics where (local) IFSs are employed mostly for the representation of smooth curves and surfaces.
We suggest the construction and use of IFSs for the solution of PDEs. This is done by choosing an initial IFS and then changing it iteratively until it approximates a desired function given by either data or functional (e.g. partial differential) equations. We use ideas based on the collage theorem to fit a given function class and refine the domains of the IFS if necessary.
In the following we will discuss the numerical application of local IFSs which is not based on a basis of a linear space but on the IFS itself. We anticipate that this approach has the following advantages over approaches that are based on a linear basis: • The same approach can be used to approximate and solve PDEs on very general grids defined by IFSs including fractal sets. • Visualisation and numerical solutions are computed simultaneously and can be done on the same or on neighboring processors of a multiprocessor system such that communication overhead may be reduced. • Dimensionality is handled much more flexibly in fractals -for example, one may use 1D solvers for higher-dimensional problems. • We can at the same time adapt the basis functions (or frames) as well as solving the problem. Searches over large collections of dictionaries of finite dimensional approximation spaces can be done locally during the solution. • The computational complexity is bounded by the resolution one requires.
• Adaptivity is naturally included as in wavelet-based methods and is a consequence of the iteration -one application of the IFS reduces the finest scale. • Convergence of the method can be controlled with few parameters and is driven by the convergence of the IFS. • The theory is based on the theory for fractals and IFSs which is well established. In addition, there has been a lot of work on wavelets and subdivision schemes which provides further firm foundations.

The Collage Theorem.
While it is usually assumed that the iterated function system (IFS) is given, a very important class of methods used in image processing determines the IFS from its fixed point. An important result used here is the Collage Theorem [2]. For the purposes of self-containment, we state this theorem below.
Theorem 1.1. Let (X, d X ) be a complete metric space. Denote by (H(X), d H ) the associated complete metric space based on the hyperspace of nonempty compact subsets of X endowed with the Hausdorff metric d H . Let M ∈ H(X) and ε > 0 be given. Suppose that F := {X; f 1 , . . . , f N } is a contractive IFS such that where A is the attractor of the IFS and and s : It has been demonstrated that approaches that are based on the Collage Theorem lead to very efficient image compression methods. The interested reader is referred to [6,13] for methodologies and to [18] for a summery of fractal-type approaches in an analytical setting. Note that the application of an IFS starts with points on a large scale and then moves to finer scales. This is very similar to some multigrid methods and wavelet methods.

Local Iterated Function
Systems. The concept of local iterated function system is a generalization of an IFS as defined in [2] and was first introduced in [6].
In the following, (X, d X ) denotes a complete metric space with metric d X and N := {1, 2, 3, . . .} the set of positive integers.
is a family of nonempty subsets of X. Further assume that for each X i there exists a continuous mapping f i : Note that if each X i = X, then Definition 2.1 coincides with the usual definition of a standard (global) IFS on a complete metric space. However, the possibility of choosing the domain for each continuous mapping f i different from the entire space X adds additional flexibility as will be recognized in the sequel.
A mapping f : U ⊂ X → X is called contractive on U or a contraction on U if there exists a constant s ∈ [0, 1) so that Definition 2.2. A local IFS F loc is called contractive if there exists a metric d equivalent to d X with respect to which all functions f ∈ F loc are contractive (on their respective domains).
Let 2 X := {S | S ⊆ X} be the power set of X. With a local IFS we associate a set-valued operator F loc : 2 X → 2 X by setting Here By a slight abuse of notation, we use the same symbol for a local IFS and its associated operator.
In (2.2) we allow for A ∩ X i to be the empty set. Thus, every local IFS has at least one local attractor, namely A = ∅. However, it may also have many distinct ones. In the latter case, if A 1 and A 2 are distinct local attractors, then A 1 ∪ A 2 is also a local attractor. Hence, there exists a largest local attractor for F loc , namely the union of all distinct local attractors. We refer to this largest local attractor as the local attractor of a local IFS F loc . We remark that there exists an alternative definition for (2.1). For given functions f i which are only defined on X i one could introduce set functions (which will also be called f i ) which are defined on 2 X by On the left-hand side f i (S ∩ X i ) is the set of values of the original f i as in the previous definition. This extension of a given function f i to sets S which include elements which are not in the domain of f i basically just ignores those elements. In the following we will assume this definition of the set function f i to be used.
In the case where X is compact and the X i , i ∈ N N closed, i.e., compact in X, and where the local IFS {X; (X i , f i ) | i ∈ N N } is contractive, the local attractor may be computed as follows. Let K 0 := X and set Then {K n | n ∈ N 0 } is a decreasing nested sequence of compact sets. If each K n is nonempty, then by the Cantor Intersection Theorem, Using [19,Proposition 3 (vii)], we see that where the limit is taken with respect to the Hausdorff metric on H(X). This implies that Thus, K = A loc . A (mild) condition guaranteeing that each K n is nonempty is that [6].) In the above setting where the f i have been extended to 2 X , one can derive a relation between the local attractor A loc of a contractive local IFS {X; (X i , f i ) | i ∈ N N } and the (global) attractor A of the associated (global) IFS {X; f i | i ∈ N N } where the extensions of f i to all sets are defined as above. To this end, let the sequence {K n | n ∈ N 0 } be defined as above. The unique attractor A of the IFS F := {X; f i | i ∈ N N } is obtained as the fixed point of the set-valued map F : H(X) → H(X), where B ∈ H(X). If the IFS F is contractive, then the set-valued mapping (2.3) is contractive on H(X) (with respect to the Hausdorff metric) and its fixed point can be obtained as the limit of the sequence of sets {A n | n ∈ N 0 } with A 0 := X and A n := F(A n−1 ), n ∈ N.
Note that K 0 = A 0 = X and, assuming that K n−1 ⊆ A n−1 , n ∈ N, it follows by induction that Hence, upon taking the limit with respect to the Hausdorff metric as n → ∞, we obtain A loc ⊆ A. This proves the next result. Proposition 2.4. Let X be a compact metric space and let X i , i ∈ N N , be closed, i.e., compact in X. Suppose that the local IFS F loc := {X; (X i , f i ) | i ∈ N N } and the IFS F := {X; f i | i ∈ N N } are both contractive. Then the local attractor A loc of F loc is a subset of the attractor A of F.
Contractive local IFSs are point-fibered if X is compact and the X i , i ∈ N N , are closed. To show this, define the code space of a local IFS by Ω := n∈N N N and endow it with the product topology T. It is known that Ω is metrizable and that T is induced by the Fréchet metric d F : Ω × Ω → R, where σ = (σ 1 . . . σ n . . .) and τ = (τ 1 . . . τ n . . .). (As a reference, see for instance [12], Theorem 4.2.2.) The elements of Ω are called codes.

Local Fractal Functions.
In this section, we exhibit a class of special attractors of local IFSs, namely local attractors that are the graphs of bounded functions. These functions will be called local fractal functions. We prove that the set of discontinuities of these bounded functions is countably infinite and we derive conditions under which local fractal functions are elements of the Lebesgue spaces L p .
To this end, we assume that 1 < N ∈ N and set N N := {1, . . . , N }. Let X be a nonempty connected set and {X i | i ∈ N N } a family of nonempty connected subsets of X. Suppose {u i : X i → X | i ∈ N N } is a family of bijective mappings with the property that Endowed with the metric Note that under the usual addition and scalar multiplication of functions, the spaces B(X i , Y) and B(X, Y) become metric linear spaces. A metric linear space is a vector space endowed with a metric under which the operations of vector addition and scalar multiplication are continuous.
For i ∈ N N , let v i : X i × Y → Y be a mapping that is uniformly contractive in the second variable, i.e., there exists an ∈ [0, 1) so that for all y 1 , y 2 ∈ Y Define a Read-Bajactarević (RB) operator Φ : where f i := f | Xi and Note that Φ is well-defined and since f is bounded and each v i contractive in the second variable, Φf ∈ B(X, Y). Moreover, by (3.1), we obtain for all f, g ∈ B(X, Y) the following inequality: To simplify notation, we set v(x, y) := N i=1 v i (x, y) χ Xi (x) in the above equation. In other words, Φ is a contraction on the complete metric space B(X, Y) and, by the Banach Fixed Point Theorem, has therefore a unique fixed point f * in B(X, Y). This unique fixed point will be called a local fractal function f * = f * Φ (generated by Φ). Next, we would like to consider a special choice for mappings v i . To this end, we require the concept of an F -space. We recall that a metric d : If in addition we require that the metric d Y is homogeneous, that is, then v i given by (3.4) satisfies condition (3.1) provided that the functions S i are bounded on X i with bounds in [0, 1) for then Here, • ∞,Xi denotes the supremum norm with respect to X i and s : Thus, for a fixed set of functions {λ 1 , . . . , λ N } and {S 1 , . . . , S N }, the associated RB operator (3.2) has now the form or, equivalently, is a family of bijective mappings satisfying property (P).
Define a mapping Φ : If max{ S i ∞,Xi | i ∈ N N } < 1 then the operator Φ(λ)(S) is contractive on the complete metric space B(X, Y) and its unique fixed point f * satisfies the self-referential equation or, equivalently Proof. The statements follow directly from the considerations preceding the theorem.
Remark 2. Note that the local fractal function f * generated by the operator defined by (3.5) does not only depend on the family of subsets The fixed point f * should therefore be written more precisely as f * (λ)(S). However, for the sake of notational simplicity, we usually suppress this dependence for both f * and Φ.
The following result found in [15] and in more general form in [21] is the extension to the setting of local fractal functions.
Proof. Let α, β ∈ R and let λ, µ ∈ . Injectivity follows immediately from the fixed point equation (3.6) and the uniqueness of the fixed point: Linearity follows from (3.6), the uniqueness of the fixed point and injectivity: and .
For surjectivity, we define We may construct local fractal functions on spaces other than B(X, Y). To this end, we assume again that the functions v i are given by (3.4) and that X := [0, 1] and Y := R. We consider the metric on R and [0, 1] as being induced by the L 1 -norm. Note that endowed with this norm B([0, 1], R) becomes a Banach space.
We have the following result for RB-operators defined on the Lebesgue spaces To this end, let g, h ∈ L p [0, 1] and let p ∈ [0, ∞). Then These calculations prove the claims. Remark 3. The proof of the theorem shows that the conclusions also hold under the assumption that the family of mappings {u i : . In this case, the a i 's are given by The a i 's are given as above in item (i). Next we investigate the set of discontinuities of the fixed point f * of the RBoperator (3.8).
Theorem 3.5. Let Φ be given as in (3.8). Assume that for all i ∈ N N the u i are contractive and the λ i are continuous on X i . Further assume that condition (3.9) is satisfied for p = ∞ and that the fixed point f * is bounded everywhere. Then the set of discontinuities of f * is at most countably infinite.
Proof. Let f be a real-valued function and U a nonempty open interval contained in its domain. The oscillation of f on U is defined as and the oscillation of a function f at a point x 0 inside an open interval contained in its domain is defined by The Banach Fixed Point Theorem implies that we may start with any bounded function, say f 0 = χ [0,1] , to construct a sequence of iterates f n := Φf n−1 , n ∈ N, which under the given hypotheses, converge in the L ∞ -norm to the fixed point f * .
Each iterate f n may have finite jump discontinuities at the interior knots {x j | j = 1, . . . , N − 1} of the partition P and also at the images u i1 • u i2 • · · · • u i−1 (x j ) of the interior knots. The number of possible discontinuities at level n is bounded above by N n−1 (N −1) since the sets X i may only contain a subset of the interior knots. Denote by E n the finite set of all finite jump discontinuities at level n and let E := n∈N E n . Note that E is at most countably infinite.
Let x ∈ [0, 1] \ E and let ε > 0. The fixed point equation for f * , Given any K ∈ N there exists a nonempty compact interval I K such that Since f * is bounded on [0, 1], |J| ≤ 1, and a K → 0 as K → ∞, we can choose a K large enough so that s K ω(f * ; J) < ε/2 and Λ |J| a K /|a − s| < ε/2. Thus, ω(f * ; I K ) < ε, which proves the continuity of f * at all points in [0, 1] \ E and completes the proof. Proof. This is a direct consequence of the above theorem and, for instance, Theorem 7.5 in [26].
Next, we exhibit the relation between the graph G of the fixed point f * of the operator Φ given by (3.2) and the local attractor of an associated contractive local IFS. To this end, we need to require that X is a closed subset of a complete metric space. Consider the complete metric space X × Y and define mappings w i : Assume that the mappings v i : X i × Y → Y in addition to being uniformly contractive in the second variable are also uniformly Lipschitz continuous in the first variable, i.e., that there exists a constant L > 0 so that for all y ∈ Y, Denote by a := max{a i | i ∈ N N } the largest of the Lipschitz constants of the mappings u i : X i → X and let θ : is then a metric for X × Y which is compatible with the product topology on X × Y. Theorem 3.7. The family W loc := {X × Y; (X i × Y, w i ) | i ∈ N N } is a contractive local IFS in the metric d θ and the graph G(f * ) of the local fractal function f * associated with the operator Φ given by (3.8) is an attractor of W loc . Moreover, where W loc denotes the set-valued operator (2.1) associated with the local IFS W loc .
The graph G(f * ) of f * is an attractor for the contractive local IFS W loc , for That (3.11) holds follows from the above computation and the fixed point equation for f * written in the form

Computation and Examples.
4.1. Computational remarks. The main step in the computation of a fractal function relates in one way or the other to the evaluation of the RB operator. We will discuss a discretisation of the RB operator here. Note that this discretisation does not involve any numerical approximations but is an exact restriction of the full RB operator and will thus (in exact arithmetic) deliver values of the full RB operator applied to a function.
For computational and visualation purposes we introduce a grid X g ⊂ X which is a finite subset. The numerical computations will then be done for functions f g : X g → Y. We introduce a restriction Φ g of the RB operator Φ by Due to the occurrence of f i (u −1 i (x)), this defines a mapping Φ g : Y X g → Y X g if the grid has the property that u −1 i (x) ∈ X g whenever x ∈ u i (X i ) ∩ X g for some X i . If a grid X g satisfies this property, we call it admissible. We then call Φ g the discrete RB operator corresponding to the RB operator Φ and the grid X g .
We will now rewrite the discrete RB operator slightly for the case where Y = R. Note that in this case f g is an element of the finite dimensional vector space R X g := R |X g | . First, we define the (potentially nonlinear) maps is then a sampling operator and we have in particular As the sets u i (X i )∩X g form a partition of X g one then has for the discrete RB operator For the special case where v i (x, y) = λ i (x)+S i (x)y, one introduces the restriction operator E i : R X g → R X g i defined by E i f (x) := f (x) for x ∈ X g i . The RB operator then is an affine mapping of the form where S i is the multiplication operator (diagonal matrix) with elements S i (x). Thus, one has Φ g f g = λ g + M f g where the matrix M is factorised in the following way: Both matrices U i and E i are sampling matrices, i.e., they contain at most one nonzero element (with value one) in each column. As the matrices S i are diagonal, one can further simplify the factorisation as Here the matrices U i S i U T i are square so that the first factor is a diagonal matrix and the factors U i E i are sampling matrices.
One sees that the discrete RB operator can be applied in parallel. However, a difficulty is still that in general the evaluation of the sampling operators U i E i may require substantial communication between the processors. This needs to be analysed for each particular case. In some (practically important) cases, however, one can reduce the amount of communication. This happens when the X i are uniquely partitioned by some u j (X g j ) in the sense that there exists a partition From the factorisation above one can derive that in this case the operator M has a block diagonal structure with M blocks. Furthermore, each block has a factorisation similar to the one above. This leads to highly efficient parallel algorithms which will be discussed elsewhere. We will refer to this case as having a local refinement. Typically, to each block belongs a standard (global) IFS so that the local IFS consists of M standard ones. The connection between the various IFSs is obtained through the choice of the λ i and S i .

Example 1:
The one-dimensional case with constant λ i and S i . For this example let X = [0, 1) and Y = R. Furthermore, let the number N of functions in the local IFS be even and let X 2j−1 = X 2j = [(j − 1)h, jh) for j = 1, . . . , N/2 where h = 2/N . Furthermore, let This choice for the mappings u i implies that The discrete grid is chosen to be uniform with h g = 1/N g and where N g is a multiple of N .
One sees that we have here a block structure as discussed at the end of the previous section with M = 2. Using vector notation, one gets with e = (1, . . . , 1) ∈ R Ng/N the vector λ = (λ 1 e, . . . , λ N e) T where F is the sampling matrix selecting every second element in the X i 's. The fractal function is defined on each domain and there it obeys the fixed point equation From these equations one can see that solving this iteratively using the fixed point iteration gives an error of the order of O((S 2 2j−1 + S 2 2j ) k/2 ) for k iterations. We selected the S i 's and the λ i 's randomly and iterated with the RB operator. The result is displayed in Figure 4   If in addition one would like to have continuity at the midpoint then one needs to require that The constants S i with odd index, S 2j−1 , were chosen randomly and those with even index as This particular choice implies that the convergence rate is at not any faster than 1/2. If one selects S i = 0.5 for all i, a piecewise linear interpolant is obtained. In Figure 4.2 we have displayed a couple of interpolants for (x(1 − x)) 0.2 . This shows that some of the interpolants have similar errors as the piecewise linear interpolant. However, it also shows that at the boundaries some of the interpolants perform substantially better than the piecewise linear interpolant.
The evaluation of the RB operator for the interpolation problem converges with the same rate as if one begins the iteration at a random point. If, however, one starts the iteration at zero one obtains finite termination for a finite grid. The number of iterations is of order O(log 2 (N g )) where N g is the number of numerical grid points.
This shows that for any function f and S i there exists a λ i . (See also Theorem 3.3.) But this λ i is as complex as the original function and thus there is no gain in representing f by λ i . In some cases, howevever, f cannot be simplified. In this case one might choose S i = 0 and thus A simple choice for the λ i is: λ i (x) = α i +β i x, for some constants α i and β i . As in the case of constants this simple model can also lead to rather complicated functions. (See Figure 4.3.) Again one observes a different behaviour on the four different ranges of the u i . Note, however, that the fractal function f is a linear function of the λ i so  [21] and [22] where this dimension is explicitly computed.) One can also determine λ i such that the resulting fractal function is interpolatory. If one chooses all S i = 0.25, one can select the λ i such that the resulting fractal function is differentiable at the boundary points between the domains of the u i and is therefore a Hermite interpolant at these points. (Cf. also [21,22].) In general, this function does have discontinuities, in particular, at the midpoints of the domains. Experiments suggest that the approximation order of this (discontinuous) interpolant is of third order in the size of the domains. This is the same order as one would expect from a piecewise quadratic function. An example of the error curve for the function exp(4x) can be seen in Figure 4.4. One can clearly observe that the error is differentiable at the grid points but has some large discontinuities within the domains X i .
While linear spaces of λ define linear function spaces of fractal functions, not every linear function space for f consists of fractal functions. For this to be the case the function space itself must be self-referential. We define a linear function space spanned by finitely many functions ψ 1 (x), . . . , ψ d (x), d ∈ N, to be self-referential if there exist matrices A i and vectors b i such that where ψ : X → R d is defined by In this case there exists λ i (x) = d j=1 c i,j ψ j (x) such that the fractal function defined by the λ i is an element of the function space. Prominent examples of such function spaces include polynomials and scaling functions. More generally, the condition of self-referantiabilty for bases is found in the subdivision schemes of computer graphics.

Polynomial Fractals.
5.1. The Taylor series. In the following we will investigate the fractal nature of the graph of polynomials This research is done with a view to the development of efficient numerical algorithms. In the future we will consider complex-valued polynomials and also real-analytic analytic functions.
For our purposes, we denote by 0 the space of all real-valued sequences having only finitely many terms not equal to zero. As is common practice, we endow 0 with the "norm" a 0 := |a 0 | 0 + · · · |a n | 0 , a ∈ 0 . Here, we defined 0 0 := 0. Furthermore, we denote by p , p > 0, the space of all real-valued sequences Note that for p ≥ 1, |x n − y n | p defines a metric making p into a complete topological vector space which is not normable. In this setting, 0 may be thought of as p where p → 0+.
Let a ∈ 0 be a finite sequence and define a function v : Then the function p given by is a polynomial and a k is the value of its k-th derivative at zero. One can see that any derivative of p satisfies This motivates the introduction of a function f : [0, 1] → 0 with components As p is a polynomial only a finite number of components f k are not equal to zero.
One can now reformulate the Taylor series of p at any point x as Similar formulas for all the derivatives of p may be obtained in a similar fashion. This can all be stated using the matrices A(x) and V (t) defined by respectively. In the following choose indices i, j = 0, 1, . . . to always start at zero. Note that A(x) is a Hankel matrix and V (t) a Toeplitz matrix. From the above one can show that the Taylor series for all the derivatives takes the form The infinite matrix A(x) is a Hankel matrix and the anti-diagonals take the values f i (x). As f (x) only has a finite number of elements (corresponding to p being a polynomial) there exists some nonnegative integer M so that As P M A is regular lower triangular Toeplitz with diagonal elements a M the matrix L = I − 1 a M P M A is also lower triangular Toeplitz with zero diagonal. It follows that L M = 0 and thus one can obtain the inverse of [A] M using the geometric series for I + L. This leads to the stated formula.
One can get an explicit formula for the inverse. For simplicity we omit the x. As M is the degree of the polynomial and P M is the reversal permutation of the first M + 1 elements, the geometric series converges as any term with k ≥ M is zero.
The matrix V (t) also has a nice structure and one can see that where σ is the forward shift matrix given by [σ] i,j = δ i,j−1 . Consequently one obtains and thus Hence, the operators V (t) form a group with In addition, we also obtain that V (s) T v(t) = v(t + s).

Self
Then the mapping w defined by and consequently w leaves the graph of v invariant. A similar observation has also been reported in a forthcoming publication Barnsley et al. [7].
Next, we like to find functions w under which the graph of f is invariant. To this end, consider Recall from above that A(x)v(t) = f (x + t). Therefore, one concludes that An argument similar to the one given in the previous example yields implying that w leaves the graph of f invariant. According to [7] we call the mappings w which leave a polynomial invariant fractels. For more details and fundamental properties of fractels, we refer the reader to the upcoming publication [7].

Affine IFSs for given polynomials.
Here we combine two fractels w from the previous section to form an IFS. The infinite matrix W s (x) = A(x)D s A(x) + has the following properties: • W s (x) is of rank M + 1; • Most eigenvalues of W s (x) are thus equal to zero. The nonzero eigenvalues are 1, s, s 2 , . . . , s M ; • W s (x) is lower triangular (and the eigenvalues are on the diagonal). It is possible to use the fractels introduced in the last section but due to the occurrence of the eigenvalue 1, the fixed point of the resulting IFS is not unique and typically depends on the starting point. Note that if linear mappings are used with all eigenvalues less than zero the only fixed point is the zero function. Thus in this case one needs eigenvalues of value 1. Such an approach may be suited for the case of projective spaces, here however we consider affine spaces. Therefore, we replace the linear w from last section by for some θ ∈ [0, 1]. In practice the choice θ = 0.5 was very stable but sometimes lead to slow convergence. Choosing θ = 0 was faster but less stable. This behaviour will be investigated further and the results reported elsewhere.
The particular choice of affine function also leaves the graph of f invariant as where we have a 0 = f 0 (0) and b 0 = f 0 (1).
6. Algorithms. In this section we present some algorithmic aspects which are mostly motivated by the Collage Theorem. We first consider convex optimisation, then grids and finally subdivision. Here we only provide a rough outline. A more detailed treatment is under development.
6.1. Collage fitting. In this section a new kind of approximant for the solution of elliptic problems is introduced. We call this approximant collage fit. Like the common Ritz method this approximation is shown to be quasi-optimal. Let in the following H be a Hilbert space and a(·, ·) be a symmetric H-elliptic form. We consider here the problem of determiningû It can be shown that the Ritz method minimises the energy norm of the errorû N −û, i.e., where v E = a(v, v). A consequence of the H-ellipticity is that the energy norm is equivalent to the H-norm, i.e., there exist c 1 , c 2 > 0 such that It follows directly that the Ritz approximation is quasi-optimal, and in particular We define V N as a set of fractal functions as follows. Let F (·; α) : H → H denote a family of RB operators (as defined in a previous section) parameterised by a parameter vector α ∈ R M . We will assume that the RB operators are contractive, i.e., that for some constant c ∈ (0, 1). We will also assume a stronger condition, namely that γ := c c 2 c 1 < 1.
Finally, we will assume that F (u; α) is a linear function of (u, α) ∈ V × R M . These assumptions hold for commonly used RB operators. The fixpoint u α of an RB operator F (·; α) is a fractal function. As approximation set for our elliptic problem we consider As F is linear in (u, α) the set V N is a finite-dimensional linear space, and, in addition that F can be decomposed as It follows that F (u; α) − F (v; α) = F (u; 0) − F (v; 0) and thus all the F (·; α) are contractive with a constant c independent of α.
For the following let W := {F (0; α) | α ∈ R M }. Note that W is a linear space and define the affine space We now introduce the operator G : H → H by where Ψ(v) is the quadratic form defined previously. Proposition 6.1.
• Let Ψ be an H-elliptic quadratic form which defines an energy norm · E for which there exist c 1 , c 2 > 0 such that c 1 v ≤ v E ≤ c 2 v for all v ∈ H. • Let F (u; α) = F (u; 0) + F (0; α) define an RB operator which is contractive with constant c such that c < c 1 /c 2 . • Let G(u) = argmin w∈W (u) Ψ(w). Then the so defined operator G is contractive and Proof. As G(u) is the best approximation in W (u) toû one can show thatû−G(u) is orthogonal to the space W with respect to the scalar product a(·, ·) and the same holds forû − G(v). Thus G(u) − G(v) is orthogonal to W in the same scalar product. It follows that G(u) − G(v) E is the distance between W (u) and W (v) in the energy norm. As this distance is the minimum distance between any point of W (u) and any point of W (v) one has in particular One then has: Corollary 6.2 (Existence of collage fitũ N ). Let G be as in Proposition 6.1. Then there exists a uniqueũ N ∈ V N such thatũ N = G(ũ N ).
Proof. As G is contractive there exists a uniqueũ N ∈ H such thatũ N = G(ũ N ).
Thus the collage fitũ N ∈ V N is defined to be the fixpoint of G. Note that this is an approximation ofû which is in V N , it is,however, in general different from the Ritz approximationû N . Nonetheless it is also a quasi-optimal approximation: Proposition 6.3 (quasi-optimality of collage fit). Letũ N be the collage fit for the quadratic form Ψ as defined in Corollary 6.2. If all the assumptions of this corollary hold and if γ and c are as defined in this corollary, then one has Proof. Let u N ∈ V N and α ∈ R M such that u N = F (u N ; α). Asũ N minimises the energy norm in W (ũ N ) one has By the triangle inequality and contractivity of F (·; α) one has The claimed inequality follows directly. We then compute the collage fitũ N iteratively using the fixpoint algorithm for G:
The algorithm converges because of the contractivity of the operator G. Applications of this algorithm include quasi-optimal approximations of L 2 functions by classes of fractal functions. In practice we found that these approximations are very close to the best L 2 approximations. Other applications are the computation of fractal approximations to the solution of Fredholm integral equations of the first kind using Tikhonov regularisation. Finally, this approach can also be used to solve elliptic PDEs numerically with fractal functions. More details on these applications will be provided in a forthcoming paper.

Evaluation of functions on grids.
Grids are very important objects for numerical computations. They are a collection of points where during the computations one needs to compute the unknown function in order to get the value at the points one is interested in.
In the simple case of local interpolation one requires just neighboring points. However, if one would like to solve a PDE one needs a whole field.
In the end, however, the values of interest are a function of a certain collection of values at other points. This is a type of self-referentiality and we now proceed to define self-referential grids. This approach is based on a upcoming paper by Barnsley et al. [7] on the computation of function values. Here we only consider the simple case discussed above where we have an IFS with two functions and our functions are defined over [0,1].
To this end, consider x ∈ [0, 1]. Then x has a numerical representation of the form Then x is defined by the recursion If f (x) is the vector of derivatives of some polynomial evaluated at x then one may use the recursion to obtain y = f (x). This is essentially the method of function evaluation discussed in [7].
Here we note that in order to obtain the value of f at the point x one requires the values of f on all x (k) . This is the "grid" required to determine f (x). This "grid" is nothing else but the path of the shift function σ starting at point x where More generally, we define a self-referential grid γ as a finite set of points We now have Proposition 6.4. A self-referential grid γ is invariant under σ, i.e., σ(γ) ⊂ γ.

Proof.
As γ is self-referential there exists for every x ∈ γ a z ∈ γ such that x = l i (z), for some i ∈ {0, 1}.
Hence, we now can define for any finite set M ⊂ [0, 1] a self-referential grid γ M = k σ k (γ). If we know the IFS we can then determine the values on all the points of γ M recursively (as outlined above). In particular, one then also obtains the values on the set M . One could use this for multiscale modelling where one models fine scale behaviour on just a small subset of a very fine grid and uses self-referentiality to get the overall solution.
6.3. Subdivision schemes. Subdivision schemes are widely used in computer graphics for modelling curves and surfaces. An introduction and survey of the mathematics of subdivision schemes can be found in [8,9,23,27,28]. A subdivision scheme is a collection of mappings (called refinement rules) R k : V k → V k+1 between linear spaces V k of real functions defined on nested meshes (at most countable sets of isolated points) N 0 ⊂ N 1 ⊂ · · · ⊂ R s . Iterated function systems (and LIFSs) provide a rich source of subdivision schemes. For example, consider an IFS with X = [0, 1), N = 2 and u 1 (x) = 1 2 x and u 2 (x) = 1 2 + x 2 for x ∈ [0, 1). Furthermore, let v i (x, y) be continuous for i = 1, 2 with v 1 (1, y) = v 2 (0, y) for y ∈ R. One then obtains a subdivision scheme with meshes N k = 2 −k Z 2 k and Z 2 k = {0, . . . , 2 k − 1} by choosing the refinement rules R k : R N k → R N k+1 to be Note that these rules are well defined as 2ξ ∈ N k in the first case and 2ξ − 1 ∈ N k in the second case. An important question regards the convergence of subdivision schemes to continuous functions. In the examples generated by LIFSs one obtains this convergence directly from the convergence of the LIFS itself. Subdivision algorithms -like LIFSs -are used to generate the values of graphs of functions. Generalising the concept of polyomial fractals discussed in the previous section, one now may obtain LIFSs from the common subdivision schemes, see also the book by Prautzsch et al [28] for a different angle of this discussion based on Bézier splines. More specifically, Micchelli and Prautzsch [27,23] discuss refinement algorithms which use a refined basis based on uniform subdivision. They present the unified structure of a large class of smoothing methods. In particular, they show that the obtained curves are uniformly refinable or self-referential in the sense that the curve may be patched together from scaled subcomponents of itself. This fundamentally defines a local IFS and, in particular, generalizes methods used for Bezier curves which are based on polynomials.

Conclusions and Final Remarks.
We have demonstrated that fractal functions defined by local iterated function systems can be used for computations. In fact, many known methods including piecewise polynomial approximation and wavelets and more generally subdivision schemes can be described within the fractal framework because the underlying components (the polynomials and wavelets) have a fractal nature.
While this fractal nature has been observed in particular in the subdivision and wavelet literature, one observes that even some of the newest numerical approximation schemes do have a fractal nature. As an illustration thereof, we consider here the QTT (quantized tensor train) method. It considers functions which can be represented by matrix products of the form f (x) = α1,...,α d where x has the binary representation The approximation of functions using their binary digits in this way was motivated by the work on high-dimensional approximation and quantum mechanics. QTT was introduced by Oseledets in [24]. The summation ranges of the indices α k = 1, . . . , r k are defined by the tensor train ranks r k . For computational efficiency it is important that these ranks are small. Except for special cases (the exponential function, trigonometric functions and piecewise polynomials) little is known [10,14,25] about which functions can be approximated by QTT functions with low ranks. We briefly remark that fractals admit such a representation. This demonstrates that the fractal framework considered here is also useful for the analysis of the QTT method.
Consider in particular a fractal function defined by f (x/2 + 1/2) = λ 2 + S 2 f (x). Then the recursion for the fractal function can be rewritten as or in matrix form as If one now iterates this for f (y) one gets the factorisation This provides an explicit QTT representation for the fractal function f (x) and shows that these fractal functions have QTT rank 2. Note, however, that we have only considered a function class with 4 parameters λ i and S i . Rank 2 QTT functions allow the parameters to depend on the levels or position of digits. This can also be discussed in the fractal framework and will be considered in future works as will local IFSs.