\section{Introduction}

Over the last decades, fractal interpolation and approximation have been extensively research. This research originated with \cite{B2} where a special set-up was used to define so-called affine fractal interpolation functions. The graphs of these affine fractal interpolation functions are the attractors of a class of iterated function systems and thus geometrically motivated. An analytic construction of general fractal functions originated in \cite{bedford,dubuc1,dubuc2} where the concept of a Read-Bajractarevi\'c operator is first encountered. Numerous constructions of fractal functions based on Read-Bajractarevi\'c operators satisfying given interpolation and approximation conditions followed. Some of these constructions are introduced and summarized in \cite{massopust,massopust1}. The number of publications in fractal interpolation theory is enormous and the interested reader may want to search for fractal functions using terms such as ``hidden variable,'' ``$V$-variable,'' ``coalescent,'' ``super,'' and ``$\alpha$-fractal functions'' to name just a few. 

The intend of this chapter is to introduce the reader to the concept of fractal interpolation and its extensions from a global to a local setting, then to non-stationarity and finally to a quaternionic setting. In a certain sense, these are the main set-ups with the possible exclusion of \emph{unbounded} fractal interpolation \cite{m7}. It is understandable that such an endeavor must necessarily restrict itself to the main points of each construction and setting. However, the exposition will give the reader an overall perspective of the issues involved and the techniques used, and can be used as a starting point for a deeper investigation into each the topics.

The outline of this chapter is as follows. Section 2 introduces the global setting of fractal interpolation and exhibits a relationship between the Read-Bajractarevi\'c operator and the solution to a canonically associated system of functional equations. In the next section, we extend global interpolation to a local set-up giving more flexibility to the construction. This type of interpolation found deep applications to fractal imaging and fractal compression \cite{BH}. In Section 4, the recently introduced concept of non-stationary fractal interpolation is presented and it is shown that backward trajectories allow distinct features to be delineated at different interpolation scales. The final Section 5, describes the novel setting of fractal interpolation in the theory of quaternions and shows that the non-commutative character of quaternions introduces even more intricate fractal patterns.

\section{Global Fractal Interpolation}\label{GFI}

The purpose of the current section is to introduce global fractal interpolation and to relate the global fractal interpolant to the solution of a system of functional equations. We see that this system of functional equations defines in a canonical way a Read-Bajractarevi\'c (RB) operator and vice versa. This relationship will be encountered several times in the subsequent sections as well.

In the following, $(\mathsf{E},d_\mathsf{E})$ denotes a normed space and $(\mathsf{F},d_\mathsf{F})$ a Banach space. For $n\in \mathbb{N}$, we write $\mathbb{N}_n :=\{1,\ldots, n\}$ for the initial segment of the natural numbers $\mathbb{N}$ of length $n$. 

For a given normed space $(\mathsf{E}, \n{\cdot}_\mathsf{E})$ and a map $f: \mathsf{E} \to\mathsf{E}$, we define the Lipschitz constant associated with $f$ by
\[
\Lip (f) := \sup_{x,y \in \mathsf{E}, x \neq y} \frac{\n{f(x)-f(y)}_\mathsf{E}}{\n{x-y}_\mathsf{E}}.
\]
The map $f$ is called Lipschitz if $\Lip (f) < + \infty$ and a contraction (on $\mathsf{E}$) if $\Lip (f) < 1$.

\subsection{Bounded Solutions}
Let $\mathsf{X}$ be a nonempty bounded subset of $\mathsf{E}$. Suppose we are given a finite family $\{l_i\}_{i = 1}^{n}$ of injective contractions $\mathsf{X}\to \mathsf{X}$ generating a partition of $\mathsf{X}$ in the sense that
\begin{align}
&\forall\;i, j\in \mathbb{N}_n, i\neq j: l_i(\mathsf{X})\cap l_j(\mathsf{X}) = \emptyset;\label{c1}\\
&\mathsf{X} = \bigcup_{i=1}^n l_i(\mathsf{X}).\label{c2}
\end{align}
For simplicity, we write $\mathsf{X}_i := l_i(\mathsf{X})$.

Given the above set-up, we are looking for a global function $\psi:\mathsf{X} = \bigcup\limits_{i=1}^n \mathsf{X}_i\to\mathsf{F}$ satisfying $n$ functional equations of the form
\begin{equation}\label{psieq}
\psi (l_i (x)) = q_i (x) + s_i (x) \psi (x), \quad\text{on $\mathsf{X}$ and for $i\in \mathbb{N}_n$},
\end{equation}
where for each $i\in\mathbb{N}_n$, $s_i$ is a given bounded function $\mathsf{X}\to\mathbb{R}$ and $q_i$ a bounded function $\mathsf{X}\to\mathsf{F}$. Recall that a function $f:\mathsf{X}\to\mathsf{F}$ is called \emph{bounded} if there exists a finite $M>0$ such that $\n{f(x)}_\mathsf{F} \leq M$, for all $x\in \mathsf{X}$.

The idea is to consider \eqref{psieq} as the fixed point equation for an associated affine operator acting on an appropriately defined function space.

To this end, let $\mathcal{B} (\mathsf{X},\mathsf{F}) := \{f:\mathsf{X}\to\mathsf{F} : \text{$f$ is bounded}\}$ denote the the Banach space of bounded functions equipped with the supremums norm $\n{f} := \sup\limits_{x\in \mathsf{X}} \n{f(x)}_\mathsf{F}$. 

On the Banach space $\mathcal{B} (\mathsf{X},\mathsf{F})$, we define an affine operator $T: \mathcal{B} (\mathsf{X},\mathsf{F})\to \mathcal{B} (\mathsf{X},\mathsf{F})$, called a Read-Bajractarevi\'c (RB) operator, by 
\begin{equation}\label{eq3.17}
T f (x) = (q_i\circ l_i^{-1})(x) + (s_i\circ l_i^{-1})(x)\cdot (f\circ l_i^{-1})(x), 
\end{equation} 
for $x\in \mathsf{X}_i$ and $i\in \mathbb{N}_n$, or, equivalently, by
\begin{align*}
T f (x) &= \sum_{i=1}^n (q_i\circ l_i^{-1})(x)\, \chi_{\mathsf{X}_i}(x) + \sum_{i=1}^n (s_i\circ l_i^{-1})(x)\cdot (f\circ l_i^{-1})(x)\, \chi_{\mathsf{X}_i}(x)\\
&= T(0) + \sum_{i=1}^n (s_i\circ l_i^{-1})(x)\cdot (f\circ l_i^{-1})(x)\, \chi_{\mathsf{X}_i}(x),\quad x\in \mathsf{X},
\end{align*}
where $\chi_S$ denotes the characteristic function of a set $S$: $\chi_S(x) = 1$, if $x\in S$, and $\chi_S(x) = 0$, otherwise.

The following result is well-known (see, for instance, \cite{B2,massopust1}) but for the sake of completeness we reproduce the proof. We also refer the interested reader to \cite{SB} where a similar set-up is considered.

\begin{theorem}\label{sol}
The system of functional equations \eqref{psieq} has a unique bounded solution $\psi: \mathsf{X}\to\mathsf{F}$ provided that
\begin{enumerate}
\item conditions \eqref{c1} and \eqref{c2} are satisfied, and
\item $s:= \max\limits_{i\in \mathbb{N}_n} \sup\limits_{x\in \mathsf{X}} |s_i(x)| < 1$.
\end{enumerate} 
\end{theorem}

\begin{proof}
First note that, as the mappings $l_i$ are injective, the right-hand side of \eqref{psieq} can be written as the right-hand side of \eqref{eq3.17}.

As the functions $l_i$, $q_i$, and $s_i$ are all assumed to be bounded, $T$ maps $\mathcal{B} (\mathsf{X},\mathsf{F})$ into itself. For all $f,g\in \mathcal{B} (\mathsf{X},\mathsf{F})$, we have that
\begin{align*}
\sup_{x\in\mathsf{X}}\n{Tf(x) - Tg(x)}_\mathsf{F} &= \max_{i\in \mathbb{N}_n}\sup_{x\in\mathsf{X}_i} \n{(s_i\circ l_i^{-1})(x)\cdot (f-g)\circ l_i^{-1}(x)}_\mathsf{F}\\
&= \max_{i\in \mathbb{N}_n}\sup_{\xi\in\mathsf{X}} \n{s_i(\xi)\cdot (f-g)(\xi)}_\mathsf{F}\\
&\leq \max_{i\in \mathbb{N}_n}\sup\limits_{x\in \mathsf{X}} |s_i(x)| \sup_{x\in\mathsf{X}}\n{(f-g)(x)}_\mathsf{F},
\end{align*}
from which it follows that
\[
\n{Tf-Tg} \leq s \n{f-g}.
\]
Hence, $T$ is contractive on the Banach space $\mathcal{B} (\mathsf{X},\mathsf{F})$ and therefore, by the Banach Fixed Point Theorem, has a unique fixed point $\psi\in\mathcal{B} (\mathsf{X},\mathsf{F})$. This fixed point solves the functional equations \eqref{psieq}.
\end{proof}

\begin{remarks}\hfill
\begin{enumerate}
\item The fixed point $\psi\in\mathcal{B} (\mathsf{X},\mathsf{F})$ of the RB operator $T$ is also called a \emph{bounded fractal function}. In this context, Eqn. \eqref{psieq} is also referred to as a \emph{self-referential equation} for $\psi$. 
\item The self-referential equation $T\psi = \psi$ expresses the fractal nature of the $\gr\psi$: It is made up of a finite number of copies of itself with each copy being supported on the partitioning sets $\mathsf{X}_i$. Hence, the terminology \emph{fractal function} for $\psi$.
\item The proof of Banach's Fixed Point Theorem also provides an algorithm for the construction of $\psi$: Choose \emph{any} function $\psi_0\in \mathcal{B} (\mathsf{X},\mathsf{F})$ and iteratively define the following sequence of functions: 
$$
\psi_k := T \psi_{k-1}, \quad k\in \mathbb{N}. 
$$
Then, $\psi$ is given by $\psi = \lim\limits_{k\to\infty}\psi_k$ where the limit is taking with respect to the norm $\n{\cdot}$ on  $\mathcal{B} (\mathsf{X},\mathsf{F})$.
\item The afore-mentioned algorithm for the construction of $\psi$ together with the proof of the Banach Fixed Point Theorem gives an error estimate as well, namely,
\[
\n{\psi - \psi_{k}} \leq \frac{s^k}{1-s} \n{\psi_1 - \psi_0}, \quad k\in \mathbb{N}.
\]
\item The fixed point $\psi$ depends on $n$, the partition $(\mathsf{X}_i : i\in\mathbb{N}_n)$, and the functions $s_i$ and $q_i$ with different choices yielding different fractal functions.
\item Emphasizing the dependence of $\psi$ on the functions $s_i$, the expression $s$-fractal function can be found in the literature. (See, for instance, \cite{N}.) In this context, one considers a fractal function as the image under an operator $\mathcal{F}^s$ associating with a given (non-fractal) function its fractal analogue.
\item Functional equations such as \eqref{psieq} exhibit connections to so-called \emph{fractels} \cite{bhm2,massopust1} and also to the approximation of rough functions \cite{BHVV}.
\end{enumerate}
\end{remarks}

Conditions \eqref{c1} and \eqref{c2} cannot be relaxed without adding some compatibility conditions to guarantee that the RB operator $T$ has the form given by Eqn. \eqref{eq3.17}. Should Eqn. \eqref{c1} not be satisfied, one would have to impose in our current setting the following compatibility conditions:
\begin{align}\label{wd}
\forall x_1,&x_2\in\mathsf{X}:\nonumber\\
& \; l_i (x_1) = l_j (x_2)\;\;\Longrightarrow\;\; q_i(x_1) + s_i (x_1) \psi (x_1) = q_j(x_2) + s_ j(x_2) \psi (x_2).
\end{align}
We refer to \cite{SB,SB1} for more details regarding this issue.\\

As an application of the above approach to obtain solutions to functional equations of the form \eqref{psieq} or equivalently finding the unique fixed points of the associated RB operator \eqref{eq3.17}, we provide the following example.
\begin{example}\label{ex1}
Let $\mathsf{E} := \mathbb{R} =: \mathsf{F}$ together with the Euclidean norm $\abs{\cdot}$. Further. let $\mathsf{X} := [0,1)\subset\mathsf{E}$. Assume that we are given two injective contractions $l_i:[0,1)\to [0,1)$, $i=1,2$, with $l_1(x) := \frac13 x$ and $l_2(x) := \frac23 x + \frac13$. Hence, $\mathsf{X}_1 = [0,\frac13)$ and $\mathsf{X}_2 = [\frac13,1)$. Clearly, $\mathsf{X} = \mathsf{X}_1\cup\mathsf{X}_2$ and $\mathsf{X}_1\cap\mathsf{X}_2 = \emptyset$. 

Now choose $q_1(x) = -1$, $q_2(x) = {x}$, $s_1(x) = \frac12 \sin(x)$, and $s_2(x):=-\frac23 \cos(x)$. The system of functional equations and the associated RB operator read then
\[
\psi(\tfrac13 x) = -1 + \tfrac12 \sin(x)\psi(x)\quad\text{and}\quad\psi(\tfrac23 x + \tfrac13) = {x} - \tfrac23\cos(x)\psi(x),
\]
and
\[
Tf(x) = \begin{cases} -1 + \frac12\sin(3x) f(3x), & 0\leq x < \tfrac13;\\
{3x-1} - \frac23\cos(\tfrac12(3x-1))f(\tfrac12(3x-1)), & \tfrac13 \leq x < 1,
\end{cases}
\]
respectively. 

As $s = \frac23 < 1$, $T$ is contractive. A plot of the solution, respectively, fixed point $\psi$, is shown in Fig. 1.
\begin{figure}[h!]
\begin{center}
\includegraphics[width=7cm, height= 4cm]{ex1.pdf}
\caption{The solution/fixed point $\psi$.}\label{fig1}
\end{center}
\end{figure}
\end{example}
\subsection{$L^p$ solutions}
In the following, we set $\mathsf{X}\subset\mathsf{E}:=\mathbb{R}^m$ and $\mathsf{Y}:=\mathbb{R}^k$ where the Euclidean spaces $\mathbb{R}^m$ and $\mathbb{R}^k$ are endowed with the corresponding canonical Euclidean norms.

Recall that the (real) Lebesgue spaces $L^p (\mathsf{X}, \mathbb{R}^k)$, where $\mathsf{X}\subset \mathbb{R}^m$ is nonempty, are defined as consisting of (equivalence classes of) functions $f:\mathsf{X}\to\mathbb{R}^k$ for which
\[
\n{f}_p :=\begin{cases}\left(\displaystyle{\int_\mathsf{X} \n{f(x)}^p dx} \right)^{1/p}, & 1\leq p < \infty;\\ \\
\esssup_{x\in\mathsf{X}} \n{f(x)}, & p = \infty.
\end{cases}
\]
is finite. Here, $\n{f(x)} := \sqrt{\abs{f_1(x)}^2 + \cdots + \abs{f_k(x)}^2}$ with $f := (f_1, \ldots, f_k)$.

We ask under what conditions on the functions $q_i$ and $s_i$ the solution $\psi$ is an element of $L^p(\mathsf{X}, \mathbb{R}^k)$, for $1\leq p < \infty$ and a bounded nonempty set $\mathsf{X}\subset\mathbb{R}^m$. 

To this end, note that in order for $\psi$ to be in $L^p(\mathsf{X}, \mathbb{R}^k)$, the RB operator $T$ must map $L^p(\mathsf{X}, \mathbb{R}^k)$ into itself. Therefore, the functions $q_i$ and $s_i$ must also be in $L^p(\mathsf{X}, \mathbb{R}^k)$. Moreover, $s_i$ needs to be in $L^\infty(\mathsf{X}, \mathbb{R}^k)$ for the product $s_i\cdot f$ to be in $L^p(\mathsf{X}, \mathbb{R}^k)$. Thus, as $\mathsf{X}$ is bounded it has finite measure and therefore $s_i\in L^\infty(\mathsf{X}, \mathbb{R}^k)$ implies that $s_i\in L^p(\mathsf{X}, \mathbb{R}^k)$ for all $1\leq p \leq \infty$.

Now it remains to be shown that the RB operator $T$ is contractive on $L^p(\mathsf{X}, \mathbb{R}^k)$. For this purpose, let $f,g\in L^p(\mathsf{X}, \mathbb{R}^k)$. Then, with $\mathsf{X}_i:= l_i(\mathsf{X})$,
\begin{align*}
\n{Tf - Tg}_p^p &= \int_\mathsf{X} \n{Tf (x) - Tg(x)}^p dx\\
&= \int_\mathsf{X} \n{\sum_{i=1}^n s_i(l_i^{-1}(x))\cdot (f-g)(l_i^{-1}(x)) \chi_{\mathsf{X}_i}(x)}^p dx\\
&\leq \sum_{i=1}^n \int_{\mathsf{X}_i} \abs{s_i(l_i^{-1}(x))}^p \n{(f-g)(l_i^{-1}(x))}^p dx\\
&= \sum_{i=1}^n \int_{\mathsf{X}} \abs{(l_i^{-1})'(x)} \abs{s_i(x)}^p \n{(f-g)((x)}^p dx.
\end{align*}

If, for all $i\in \mathbb{N}_n$, $\n{l_i^{-1})'}_\infty =: \lambda_i<\infty$ and $\n{s_i}_\infty = :s_i<\infty$, then 
\[
\n{Tf - Tg}_p^p \leq \left(\sum_{i=1}^n \lambda_i s_i^p\right) \n{(f-g)((x)}_p^p,
\]
and $T$ is contractive provided that
\[
\sum_{i=1}^n \lambda_i s_i^p < 1.
\]
Hence, we arrived at the following result
\begin{theorem}\label{th3.2}
The system of functional equations
\[
\psi (l_i (x)) = q_i (x) + s_i (x) \psi (x), \quad\text{on $\mathsf{X}\subset\mathbb{R}^m$ and for $i\in \mathbb{N}_n$},
\]
has a unique solution $\psi\in L^p(\mathsf{X}, \mathbb{R}^k)$, $1\leq p \leq \infty$, respectively, the RB operator
\[
T f (x) = (q_i\circ l_i^{-1})(x) + (s_i\circ l_i^{-1})(x)\cdot (f\circ l_i^{-1})(x),\quad{x\in\mathsf{X}_i},\;i\in \mathbb{N}_n,
\]
a unique fixed point $\psi\in L^p(\mathsf{X}, \mathbb{R}^k)$ provided that
\begin{enumerate}
\item $q_i\in L^p(\mathsf{X}, \mathbb{R}^k)$, $s_i\in L^\infty (\mathsf{X}, \mathbb{R}^k)$ and 
\item $\sum\limits_{i=1}^n \lambda_i s_i^p < 1$, where $\lambda_i = \n{(l_i^{-1})'}_\infty $ and $s_i = \n{s_i}_\infty$.
\end{enumerate}
\end{theorem}

\begin{remark}
In a similar fashion, one can derive conditions such that the unique solutions/fixed points $\psi$ are elements of H\"older or Sobolev spaces. See, for instance, \cite{PRM,massopust,massopust1}.
\end{remark}
\subsection{Continuous Solutions}
So far, we only considered bounded solution/fixed points $\psi$ for \eqref{psieq} and \eqref{eq3.17}. However, in some instances, a continuous or even differentiable solution is required. We only present a result for continuous $\psi$ and make some remarks about how to obtain differential solutions. 

\begin{theorem}\label{th3.3}
The system of functional equations \eqref{psieq} has a unique continuous solution $\psi: \mathsf{X}\to\mathsf{F}$ provided that
\begin{enumerate}
\item $\mathsf{X} = \bigcup\limits_{i=1}^n l_i(\mathsf{X})$,
\item the functions $l_i$, $q_i$, and $s_i$ are continuous,
\item and for all $i,j\in \mathbb{N}_n$ and $x_1, x_2\in X$:
\begin{equation}\label{joinupc}
\lim_{x\to x_1} f_j(x) = f_i(x_2)\;\;\Longrightarrow\;\; \lim_{x\to x_1} q_j(x) + s_j(x) \psi(x) = q_i(x_2) + s_i(x_2) \psi(x_2).
\end{equation}
\end{enumerate} 
\end{theorem}
\begin{proof}
We refer the interested reader to \cite{SB} or \cite{massopust}. In the former reference, the proof follows the functional equation setting and in the latter the RB operator setting.
\end{proof}

\begin{example}

\end{example}
We connect up with the previous Example 1 but choose as $\mathsf{X} := [0,1]$. We modify the functions $q_i$ to be $q_1(x) := x$ and $q_2(x) := 1-x$. but keep $s_1$ and $s_2$ unchanged. Note that here we have $\mathsf{X}_1 \cap \mathsf{X}_2 = \{\frac13\}$ and we need to ensure that conditions \eqref{wd} are satisfied. In particular, as we have $l_1(1) = \frac13 = l_2(0)$, the following equality has to hold:
\begin{equation}\label{3.7}
q_1(1) + s_1(1)\psi(1) = q_2(0) + s_2(0)\psi(0).
\end{equation}
The functional equations \eqref{psieq} imply for $x\in \{0,1\}$
\[
\psi(0) = q_1 (0) + s_1(0) \psi(0)\quad\text{and}\quad \psi(1) = q_2 (1) + s_2(1)\psi(0),
\]
which gives the values of $\psi$ at the endpoints of $\mathsf{X}$:
\[
\psi(0) = \frac{q_1(0)}{1-s_1(0)}\quad\text{and}\quad \psi(1) = \frac{q_2(1)}{1-s_2(1)}.
\]
The validity of \eqref{3.7} guarantees the existence of a \emph{bounded} solution $\psi$ (since $s = \frac23 < 1$). As $q_1(0) = 0 = q_2(1)$, the solution $\psi$ also vanishes on the boundary of $\mathsf{X}$: $\psi(0) = 0 = \psi(1)$.

In order to obtain a \emph{continuous} solution, equations \eqref {joinup} must be satisfied. In our current setting, as all the functions involved are continuous on $\mathsf{X}$, we obtain
\[
\lim_{x\to 0-} q_2(x) + s_2(x)\psi(x) = q_1(1) + s_1(1)\psi(1),
\]
which is identical to \eqref{3.7}. 

The solution/fixed point $\psi$ is therefore continuous. The graph of $\psi$ is depicted in Figure \ref{fig2}.
\begin{figure}[h!]
\begin{center}
\includegraphics[width=7cm, height= 4cm]{ex2.pdf}
\caption{A continuous solution/fixed point $\psi$.}\label{fig2}
\end{center}
\end{figure}

In \cite{SB}, a similar setting is considered.
\begin{remark}
For many applications, in particular the setting where $\mathsf{E} := \mathbb{R}$, the partition $(\mathsf{X}_i : i\in \mathbb{N}_n)$ is induced by a given finite set of data points $\{(x_i, y_i)\in \mathbb{R}\times\mathsf{F} : x_0 < x_1 < \cdots < x_n\}$ with $\mathsf{X}_i := [x_{i-1},x_i]$. The injective contractions $l_i$ are then affine mappings $[x_0,x_n]\to [x_{i-1},x_i]$. For $\mathsf{F} := \mathbb{R}$, this type of continuous fractal interpolation was first introduced in \cite{B2} and the fixed point $\psi$ was termed a \emph{fractal interpolation function} as it was also required that $y_j = \psi(x_j)$, $j\in\{0,1,\ldots, n\}$.
\end{remark}
\section{Local Fractal Interpolation}
In this section, we introduce a generalization of global fractal interpolation. Previously, we considered a fixed nonempty bounded subset $\mathsf{X}$ of a normed space $\mathsf{E}$ generating a partition of $\mathsf{X}$ via a finite family of injective contractions. Now, we replace the single subset $\mathsf{X}$ by a finite family of subsets $\mathsf{X}_i$ of $\mathsf{X}$.

More precisely, let $\{\mathsf{X}_i : i \in \mathbb{N}_n\}$ be a family of nonempty subsets of a fixed subset $\mathsf{X}$ of a normed space $\mathsf{E}$. Suppose $\{l_i\}_{i = 1}^{n}$ is a collections of injective mappings from $\mathsf{X}_i\to \mathsf{X}$ generating a partition of $\mathsf{X}$ in the sense that
\begin{align}
&\forall\;i, j\in \mathbb{N}_n, i\neq j: l_i(\mathsf{X}_i)\cap l_j(\mathsf{X}_j) = \emptyset;\label{cl1}\\
&\mathsf{X} = \bigcup_{i=1}^n l_i(\mathsf{X}_i).\label{cl2}
\end{align}
Note that the $l_i$ need not be contractive mappings here. 
\begin{remark}
One can actually have $m<n$ subsets of $\mathsf{X}$ and still be able to define $n$ injective mappings satisfying \eqref{cl1} and \eqref{cl2}. In this case some of the injections $l_i$ share the same domain $X_j$, i.e., $\mathsf{X}_j$ is repeated a certain number of times $\mathsf{X}_j = \mathsf{X}_{j+1} = \cdots = \mathsf{X}_{j+r}$. This situation occurs in Section \ref{sect4.3}.
\end{remark}
Local fractal interpolation looks for local solutions $\psi: \mathsf{X} = \bigcup\limits_{i\in \mathbb{N}_n} l_i (\mathsf{X}_i)\to \mathsf{F}$ of functional equations or for fixed points of RB operators of the form
\begin{equation}\label{localpsi}
\psi (l_i (x)) = q_i (x) + s_i (x) \psi (x), \quad\text{$x\in\mathsf{X}_i$,\;\;$i\in \mathbb{N}_n$},
\end{equation}
respectively,
\begin{equation}\label{Tlocal}
T f (x) = (q_i\circ l_i^{-1})(x) + (s_i\circ l_i^{-1})(x)\cdot (f_i\circ l_i^{-1})(x), \quad\text{$x\in l_i(\mathsf{X}_i)$,\;\; $i\in \mathbb{N}_n$},
\end{equation} 
where $f_i := f\vert_{\mathsf{X}_i}$, on appropriate function spaces.
\subsection{Bounded Local Solutions.}
The extension of the results presented in Section \ref{GFI} carry over to the setting of local fractal interpolation. However, special care must be taken when considering the domains of the functions involved. 

Besides the function space $\mathcal{B} (\mathsf{X},\mathsf{F})$ already introduced in the previous section, we also need the local version of this space, namely, $\mathcal{B} (\mathsf{X}_i,\mathsf{F})$, $i\in \mathbb{N}_n$. To this end and in view of the result below, we now assume that
\begin{enumerate}
\item $s_i\in \mathcal{B}(\mathsf{X}_i,\mathbb{R})$ and
\item $q_i\in \mathcal{B}(\mathsf{X}_i,\mathsf{F})$.
\end{enumerate}
Then the RB operator $T$ defined in \eqref{Tlocal} maps $\mathcal{B} (\mathsf{X},\mathsf{F})$ into itself. Hence, we arrive at the local version of Theorem \ref{sol}.
\begin{theorem}
The system of functional equations \eqref{localpsi} has a unique bounded solution $\psi: \mathsf{X}\to\mathsf{F}$, respectively, the RB operator \eqref{Tlocal} has a unique bounded fixed point $\psi: \mathsf{X}\to\mathsf{F}$ provided that
\begin{enumerate}
\item conditions \eqref{cl1} and \eqref{cl2} are satisfied, and
\item $s:= \max\limits_{i\in \mathbb{N}_n} \sup\limits_{x\in \mathsf{X}_i} |s_i(x)| < 1$.
\end{enumerate} 
\end{theorem}

\begin{proof}
The proof follows the same arguments as those given in the proof of Theorem \ref{sol}; replace $\mathsf{X}_i$ and $\mathsf{X}$ there by $l_i(\mathsf{X}_i)$ and $\mathsf{X}_i$ here. We leave the details to the reader who may also consult \cite{SB}.
\end{proof}

As in the global setting of fractal interpolation, a few remarks are in order.
\begin{remarks}\hfill
\begin{enumerate}
\item The solution/fixed point $\psi$ is referred to as a \emph{bounded local fractal function}.
\item The function $\psi$ depends on $n$, the family of subsets $\{\mathsf{X}_i : i\in\mathbb{N}_n\}$, the partition $(l_i(\mathsf{X}) : i\in\mathbb{N}_n)$ induced by the injective mappings $l_i$, and the now locally defined functions $s_i$ and $q_i$.
\end{enumerate}
\end{remarks}

As observed above, conditions \eqref{cl1} and \eqref{cl2} cannot be relaxed without adding some compatibility conditions to guarantee that the RB operator $T$ has the form given by Eqn.~ \eqref{Tlocal}. Should Eqn. ~\eqref{cl1} not be satisfied, one would have to impose in our current setting the following, now local, compatibility conditions:
\begin{align}
\forall i,j\in &\mathbb{N}_n\,\forall x_1\in\mathsf{X}_i\,\forall x_2\in\mathsf{X}_j:\nonumber\\
& \; l_i (x_1) = l_j (x_2)\;\;\Longrightarrow\;\; q_i(x_1) + s_i (x_1) \psi (x_1) = q_j(x_2) + s_ j(x_2) \psi (x_2).
\end{align}
We again refer to \cite{SB,SB1} for more details regarding this issue.
\subsection{$L^p$ Local Solutions}
As an application of local fractal functions for a particular setting, we again consider local fractal functions in $L^p$ spaces for $p\in[0,\infty]$.

For this purpose, we choose the following set up. Let $\mathsf{E} := \mathbb{R} := \mathsf{Y}$ with the canonical Euclidean norm and let $\mathsf{X}:= [0,1]$. Suppose that we are given a partition of $\mathsf{X}$ of the form $\Delta := (0 =: x_0 < x_1 < \cdots < x_{n-1} < x_n := 1)$, for some integer $n>1$. Furthermore, suppose that $\{\mathsf{X}_i  :  i \in \mathbb{N}_n\}$ is a family of half-open intervals of $[0,1]$. 

We define affine mappings $l_i:\mathsf{X}_i$ onto $[x_{i-1}, x_i)$, $i = 1, \ldots, n-1$, and from $\mathsf{X}_n^+ := \mathsf{X}_n\cup l_n^{-1}(1-)$ onto $[x_{n-1},x_n]$, where $l_n$ maps $\mathsf{X}_n$ onto $[x_{n-1}, x_n)$. 


We have the following result for RB-operators defined on the Lebesgue spaces $L^p[0,1]$, $1\leq p \leq \infty$. (See, also, \cite{bhm}.)

\begin{theorem}\label{th4.2}
Assume that $q_i\in L^p (\mathsf{X}_i, [0,1])$ and $s_i\in L^\infty (\mathsf{X}_i, \mathbb{R})$, $i \in\mathbb{N}_n$. 
The system of functional equations
\[
\psi (l_i (x)) = q_i (x) + s_i (x) \psi (x), \quad x\in [0,1],\; i\in \mathbb{N}_n,
\]
has a unique solution $\psi\in L^p[0,1]$, $1\leq p \leq \infty$, respectively, the RB operator
\begin{equation}\label{Top}
T f (x) = (q_i\circ l_i^{-1})(x) + (s_i\circ l_i^{-1})(x)\cdot (f_i\circ l_i^{-1})(x),\quad{x\in\mathsf{X}_i},\;i\in \mathbb{N}_n,
\end{equation}
has a unique fixed point $\psi\in L^p(\Omega, \mathbb{R}^k)$, where $f_i = f\vert_{\mathbb{X}_i}$, provided that
\begin{equation}\label{condition}
\begin{cases}
\left(\displaystyle{\sum_{i=1}^n}\, a_i \,\|s_i\|_{\infty, \mathsf{X}_i}^p\right)^{1/p} < 1, & p\in[1,\infty);\\ 
\max\limits_{i\in\mathbb{N}_n}\|s_i\|_{\infty,\mathsf{X}_i} < 1, & p = \infty,
\end{cases}
\end{equation}
where $a_i$ denotes the Lipschitz constant of $(l_i^{-1})'$. Here, we wrote $\|s_i\|_{\infty,\mathsf{X}_i}$ for $\sup\limits_{x\in \mathsf{X}_i} |s_i(x)|$.
\end{theorem}

\begin{proof}
Note that under the hypotheses on the functions $q_i$ and $s_i$ as well as the mappings $l_i$, $T f$ is well-defined and an element of $L^p[0,1]$. It remains to be shown that under condition \eqref{condition}, $T$ is contractive on $L^p[0,1]
$. 

To this end, let $g,h \in L^p [0,1]$ and let $p\in [0,\infty)$. Then
\begin{align*}
\|T g - T h\|^p_{p} & = \int\limits_{[0,1]} |T g (x) - T h (x)|^p dx\\
& = \int\limits_{[0,1]} \left|\sum_{i=1}^{n} (s_i\circ l_i^{-1})(x) [(g_i\circ l_i^{-1})(x) - (h_i\circ l_i^{-1})(x)]\,\chi_{l_i(\mathbb{X}_i)}(x)\right|^p\, dx\\
& = \sum_{i=1}^{n}\,\int\limits_{[x_{i-1},x_i]}\left| (s_i\circ l_i^{-1})(x) [(g_i\circ l_i^{-1})(x) - (h_i\circ l_i^{-1})(x)]\right|^p\,dx\\
& = \sum_{i=1}^{n}\,a_i\,\int\limits_{\mathbb{X}_i} \left| s_i (x) [g_i(x)- h_i(x)]\right|^p\,dx\\
&  \leq \sum_{i=1}^{n}\,a_i\,\|s_i\|^p_{\infty, \mathbb{X}_i}\,\int\limits_{\mathbb{X}_i} \left| g_i(x) - h_i(x)\right|^p\,dx = \sum_{i=1}^{n}\,a_i\,\|s_i\|^p_{\infty, \mathbb{X}_i}\,\|f_i - g_i\|^p_{p,\mathbb{X}_i}\\
& = \sum_{i=1}^{n}\,a_i\,\|s_i\|^p_{\infty, \mathbb{X}_i}\,\|g_i - h_i\|^p_{p} \leq \left(\sum_{i=1}^{n}\,a_i\,\|s_i\|^p_{\infty, \mathbb{X}_i}\right) \|g - h\|^p_{p}.
\end{align*}
Now let $p= \infty$. Then,
\begin{align*}
\|T g - T h\|_{\infty} & = \left\|\sum_{i=1}^{n} (s_i\circ l_i^{-1})(x) [(g_i\circ l_i^{-1})(x) - (h_i\circ l_i^{-1})(x)]\,\chi_{l_i(\mathbb{X}_i)}(x)\right\|_\infty\\
& \leq \max_{i\in\mathbb{N}_n}\,\left\| (s_i\circ l_i^{-1})(x) [(g_i\circ l_i^{-1})(x) - (h_i\circ l_i^{-1})(x)]\right\|_{\infty,\mathbb{X}_i}\\
& \leq \max_{i\in\mathbb{N}_n}\|s_i\|_{\infty,\mathbb{X}_i} \left\|g_i - h_i]\right\|_{\infty,\mathbb{X}_i} = \max_{i\in\mathbb{N}_n}\|s_i\|_{\infty,\mathbb{X}_i} \left\|g_i - h_i]\right\|_{\infty}\\
&  \leq \left(\max_{{i\in\mathbb{N}_n}}\,\|s_i\|_{\infty,\mathbb{X}_i}\right) \left\|g - h]\right\|_{\infty}
\end{align*}
These calculations prove the claims.
\end{proof}
\begin{remark}
Although conditions \eqref{condition} resembles those presented in Theorem \ref{th3.2}, they are -- because of the local nature of the estimates -- more subtle than the former ones.
\end{remark}
\subsection{Continuous Local Solutions}\label{sect4.3}
As in the global case, we can consider continuous local solutions to the given set of functional equations, respectively, an RB operator. Here, we present the analog to Theorem \ref{th3.3} but refer to the literature for the proof. A good reference in the functional equation setting is again \cite{SB}. In the RB operator setting, we refer to \cite{m2}.

\begin{theorem}\label{th4.3}
The system of functional equations \eqref{localpsi} has a unique continuous solution $\psi: \mathsf{X}\to\mathsf{F}$ provided that
\begin{enumerate}
\item $\mathsf{X} = \bigcup\limits_{i=1}^n l_i(\mathsf{X}_i)$,
\item the functions $l_i$, $q_i$, and $s_i$ are continuous,
\item and $\forall i,j\in \mathbb{N}_n$, $i\neq j$, $\forall x_1\in\mathsf{X}$, $\forall x_2\in \mathsf{X}_i$:
\begin{equation}\label{joinup}
\lim_{\substack{x\to x_1\\ x\in\mathsf{X}_j}} f_j(x) = f_i(x_2)\;\;\Longrightarrow\;\; \lim_{\substack{x\to x_1\\ x\in\mathsf{X}_j}} q_j(x) + s_j(x) \psi(x) = q_i(x_2) + s_i(x_2) \psi(x_2).
\end{equation}
\end{enumerate} 
\end{theorem}
As an application and an example for a continuous local fractal interpolation, we present the following set up which plays an important role in fractal-based numerical analysis as discussed in \cite{bhm}. 

Suppose that $\mathsf{E} := \mathbb{R} =: \mathsf{F}$ and $\mathsf{X}:= [0,1]$. For an even integer $n\in\mathbb{N}$, define subsets
\begin{equation}\label{Xs}
\mathsf{X}_{2j-1} := \mathsf{X}_{2j} := \left[\tfrac{2j-2}{n},\tfrac{2j}{n}\right], \quad j \in \{1,\ldots, \tfrac{n}{2}\}
\end{equation}
and affine mappings $l_i: \mathsf{X}_i\to [0,1]$ by 
\begin{equation}\label{ls}
l_{2j-1} (x) := \tfrac{x}{2} + \tfrac{j-1}{n}\quad\text{and}\quad l_{2j} (x) := \tfrac{x}{2} + \tfrac{j}{n}, \quad x\in \mathsf{X}_{2j-1} = \mathsf{X}_{2j}.
\end{equation}
Note that
\[
l_i (\mathsf{X}_i) = \left[\tfrac{i-1}{n},\tfrac{i}{n}\right], \quad i\in\mathbb{N}_n.
\]
Further, let $x_i := l_i(\mathsf{X}_i)\cap l_{i+1}(\mathsf{X}_{i+1}) = \{\frac{i}{n}\}$, $i\in\mathbb{N}_{n-1}$, $x_0:= 0$, and $x_n := 1$. In the terminology of \cite{SB}, the elements of $\{x_i : i\in \mathbb{N}_{n-1}\}$ are called \emph{contact points}.

We denote the distinct endpoints of the partitioning intervals $\{l_i(\mathbb{X}_{i})\}$ by $\{x_0 < x_1 <\ldots < x_{N}\}$ where $x_0 = 0$ and $x_N = 1$ and refer to them as {knots}.

Furthermore, we assume that we are given interpolation values at the endpoints of the intervals $\mathbb{X}_{2j-1} = \mathbb{X}_{2j}$:
\begin{equation}\label{Delta}
\Delta := \left\{(x_{2j}, y_j) :  j = 0, 1, \ldots, n/2\right\}.
\end{equation}
Let 
$$
\mathcal{C}_{\Delta}(\mathsf{X}) := \{f\in \mathcal{C}(\mathsf{X}) : f(x_{2j}) = y_j, \,\forall\, j = 0, 1, \ldots, n/2\}.
$$
Here, $\mathcal{C}(\mathsf{X}) := \mathcal{C}(\mathsf{X},\mathbb{R})$ denotes the Banach space of all continuous functions $\mathsf{X}\to\mathbb{R}$ endowed with the supremum norm $\n{\cdot}$. Note that $\mathcal{C}_{\Delta}(\mathsf{X})$ is a closed, hence complete, \emph{metric} subspace of $\mathcal{C}(\mathsf{X})$ to which we can apply the Banach fixed point theorem. To this end, consider an RB operator $T$ of the form \eqref{Top} acting on $\mathcal{C}_{\Delta}(\mathsf{X})$. 

In order for $T$ to map $\mathcal{C}_{\Delta}(\mathsf{X})$ into itself we require that 
\begin{equation}\label{qsandss}
q_i, s_i \in \mathcal{C}(\mathsf{X}_i) := \mathcal{C}(\mathsf{X}_i,\mathbb{R}) := \{f: \mathsf{X}_i\to \mathbb{R} :  \text{$f$ continuous}\}
\end{equation}
and that
\begin{equation}\label{intcon}
y_{j-1} = (T f) (x_{2(j-1)}) \quad \wedge \quad y_{j} = (T f) (x_{2j}), \quad j = 1, \ldots, n/2,
\end{equation}
where $x_{2j} := \frac{2j}{n}$. 

We remark that the preimages of the knots $x_{2(j-1)}$ and $x_{2j}$ are the endpoints of $\mathsf{X}_{2j-1} = \mathsf{X}_{2j}$. Substitution of $T$, as given in \eqref{Top}, into \eqref{intcon} and simplification results in
\begin{equation}\label{intcons}
\begin{split}
q_{2j-1} (x_{2(j-1)}) + \left(s_{2j-1}(x_{2(j-1)}) - 1\right) y_{j-1} & = 0,\\
q_{2j} (x_{2j}) + \left(s_{2j}(x_{2j}) - 1\right) y_{j} & = 0.
\end{split}
\end{equation}

To ensure global continuity of $T f$ on $\mathsf{X}=[0,1]$, we also have to impose the following join-up conditions at the oddly indexed knots. (Note that these oddly indexed knots are the images of the midpoints of the intervals $\mathsf{X}_{2j-1} = \mathsf{X}_{2j}$.)
\begin{equation}\label{prejoin}
(T f) (x_{2j-1}-)  = (T f) (x_{2j-1}+) , \quad j = 1, \ldots, n/2.
\end{equation}
These join-up conditions imply that
\begin{equation}\label{4.16}
q_{2j} (x_{2(j-1)}) + s_{2j} (x_{2(j-1)}) y_{j-1} = q_{2j-1} (x_{2j}) + s_{2j-1} (x_{2j})  y_j.
\end{equation}
In the case that all functions $q_i$ and $s_i$ are constant, \eqref{4.16} reduces to the condition given in \cite[Example 2]{bhm}. 

We summarize these results in the next theorem.
\begin{theorem}\label{thm8}
Let $\mathsf{X}:= [0,1]$ and let $n\in 2\mathbb{N}$. Suppose that subsets of $\mathsf{X}$ are given by \eqref{Xs} and the associated mappings $l_i$ by \eqref{ls}. Further suppose that the functions $q_i$ and $s_i$ satisfy \eqref{qsandss} and that the join-up conditions \eqref{intcon} and \eqref{prejoin} hold. Then, an RB operator $T$ of the form \eqref{Top} maps $\mathcal{C}_\Delta (\mathsf{X})$ into itself and is well-defined. 

If, in addition, 
\begin{equation}\label{condfors}
\max\left\{\|s_i\|_{\infty,\mathsf{X}_i} : i\in\mathbb{N}_n\right\} < 1,
\end{equation}
then $T$ is a contraction on $\mathcal{C}_\Delta (\mathsf{X})$ and thus possesses a unique continuous fixed point $\psi:[0,1]\to \mathbb{R}$ satisfying $\psi\in \mathcal{C}_\Delta (\mathsf{X})$.
\end{theorem}

This unique fixed point $\psi$ is called a \emph{continuous local fractal function}.

\begin{proof}
It remains to show that under the condition \eqref{condfors}, $T$ is contractive on $\mathcal{C}_\Delta (\mathsf{X})$. This, however, follows immediately from the case $p=\infty$ in the proof of Theorem \ref{th4.2}.
\end{proof}

For RB operators of the form \eqref{Top} mapping function spaces such as H\"older, Sobolev, Besov or Triebel Lizorkin into themselves, we refer the interested reader to \cite{m2,m3}.
\section{Non-stationary Fractal Interpolation}\label{sect5}
In this section, we extend the notion of global fractal interpolation to a non-stationary setting. In other words, we no longer assume that we keep the functions $q_i$ and $s_i$ the same at each level of iteration in the construction of the fixed point $\psi$.
\subsection{The non-stationary setting}
To this end, consider a doubly-indexed family of injective contractions $\{l_{i_k,k} : i_k\in \mathbb{N}_{n_k}, \, k\in \mathbb{N}\}$ from $\mathsf{X}\to \mathsf{X}$ where $\mathsf{X}$ is nonempty bounded subset of a normed space $\mathsf{E}$ generating a partition of $\mathsf{X}$ for each $k\in \mathbb{N}$ in the sense of \eqref{cl1} and \eqref{cl2}. 

Suppose that $\mathsf{F}$ is a Banach space, $\{q_{i_k,k}: i_k\in \mathbb{N}_{n_k}, \, k\in \mathbb{N}\} \subset \mathcal{B}(\mathsf{X},\mathsf{F})$, and $\{s_{i_k,k}: i_k\in \mathbb{N}_{n_k}, \, k\in \mathbb{N}\}\subset \mathcal{B}(\mathsf{X},\mathbb{R})$ are such that 
\[
s := \sup\limits_{k\in \mathbb{N}} \max\limits_{i_k\in \mathbb{N}_k} \|s_{i_k,k}\|_\infty < 1.
\] 
For each $k\in \mathbb{N}$, we define an RB operator $T_k: \mathcal{B}(\mathsf{X},\mathsf{F})\to \mathcal{B}(\mathsf{X},\mathsf{F})$ by
\begin{align}
(T_k f)(l_{i_k,k}(x)) &:= q_{i_k,k}(x) + s_{i_k,k}(x) \cdot f(x), \quad\forall x\in\mathsf{X}.\label{nonstat1}
\end{align}
It is not difficult to verify that each $T_k$ is a contraction on $\mathcal{B}(\mathsf{X},\mathsf{F})$ with Lipschitz constant 
\begin{equation}\label{Lip}
\Lip (T_k) = \max\limits_{i_k\in \mathbb{N}_k} \|s_{i_k,k}\|_\infty \leq s < 1.
\end{equation}
In order to continue, we require the following definition and result from \cite{LDV} adapted to our current setting.

\begin{definition}{\cite[Definition 3.6]{LDV}}
Let $\{T_k\}_{k\in \mathbb{N}}$ be a sequence of transformations $T_k:\mathcal{B}(\mathsf{X},\mathsf{F})\to \mathcal{B}(\mathsf{X},\mathsf{F})$. A subset $\mathscr{I}$ of $\mathcal{B}(\mathsf{X},\mathsf{F})$ is called an invariant set of the sequence $\{T_k\}_{k\in \mathbb{N}}$ if
\[
\forall\,k\in \mathbb{N}\;\forall\,x\in \mathscr{I}: T_k (x)\in \mathscr{I}.
\]
\end{definition}
A criterion for obtaining an invariant domain for a sequence $\{T_k\}_{k\in \mathbb{N}}$ of transformations on $\mathcal{B}(\mathsf{X},\mathsf{F})$ is also given in \cite{LDV}.

\begin{proposition}{ \cite[Lemma 3.7]{LDV}}\label{prop2.1}
Let $\{T_k\}_{k\in \mathbb{N}}$ be a sequence of transformations on $\mathcal{B}(\mathsf{X},\mathsf{F})$. Suppose there exists a $g\in \mathcal{B}(\mathsf{X},\mathsf{F})$ such that for all $f\in\mathcal{B}(\mathsf{X},\mathsf{F})$
\[
\n{T_k(x) - g} \leq \mu\,\n{f - g} + M,
\]
for some $\mu\in [0,1)$ and $M> 0$. Then the ball $B_r(g)$ of radius $r = M/(1-\mu)$ centered at $g$ is an invariant set for $\{T_k\}_{k\in \mathbb{N}}$.
\end{proposition}
\begin{proof}
The proof, although in a more general setting, is found in \cite{LDV}.
\end{proof}

\begin{proposition}\label{prop5.2}
Let $\{T_k\}_{k\in \mathbb{N}}$ be a sequence of RB operators of the form \eqref{nonstat1} on $(\mathcal{B}(\mathsf{X},\mathsf{F}),\n{\cdot})$. Suppose that the elements of $\{q_{i_k,k}: i_k\in \mathbb{N}_{n_k}, \, k\in \mathbb{N}\}$ satisfy
\begin{equation}\label{condq}
\sup\limits_{k\in \mathbb{N}} \max\limits_{i_k\in \mathbb{N}_k} \n{q_{i_k,k}} \leq M,
\end{equation}
for some $M > 0$. Then the ball $B_r(0)$ of radius $r=M/(1-s)$ centered at $0\in \mathcal{B}(\mathsf{X},\mathsf{F})$ is an invariant set for $\{T_k\}_{k\in \mathbb{N}}$.
\end{proposition}
\begin{proof}
Let $x\in \mathsf{X}$. Then there exists an $i_k\in \mathbb{N}_{n_k}$ with $x\in l_{i_k,k}(\mathsf{X})$. Thus, for any $f\in \mathcal{B}(\mathsf{X},\mathsf{F})$, 
\begin{align*}
\n{T_k f (x)}_\mathsf{F} &\leq \n{s_{i_k,k}\circ l_{i_k,k}^{-1}(x) \cdot f\circ l_{i_k,k}^{-1}(x)}_\mathsf{F} + \n{T_k(0)}_\mathsf{F} 
\end{align*}
By \eqref{condq},  $T_k(0)$ is uniformly bounded on $\mathcal{B}(\mathsf{X},\mathsf{F})$ by $M > 0$. Hence, 
\[
\n{s_{i_k,k}\circ l_{i_k,k}^{-1}(x) \cdot f\circ l_{i_k,k}^{-1}(x)}_\mathsf{F} \leq s \n{f\circ l_{i_k,k}^{-1}(x)}_\mathsf{F},
\]
which gives, after taking the supremum over $x\in\mathsf{X}$, $\n{T_k f} \leq s\,\n{f} + M$. Proposition \eqref{prop2.1} now yields the statement. 
\end{proof}

To arrive at the main result of this section, we require yet another definition. (See also \cite[Section 4]{LDV}.)
\begin{definition}
Suppose that $f_0\in\mathcal{B}(\mathsf{X},\mathsf{F})$ and that $\{T_k\}_{k\in \mathbb{N}}$ be a sequence of RB operators on $\mathcal{B}(\mathsf{X},\mathsf{F})$. The sequences
\begin{equation}
\Phi_k (f_0) := T_k\circ T_{k-1} \circ \cdots \circ T_1 (f_0)
\end{equation}
and
\begin{equation}
\Psi_k (f_0) := T_1\circ T_2 \circ \cdots \circ T_k (f_0)
\end{equation}
are called the forward, respectively, backward trajectory of $f_0$.
\end{definition}

In connection to the above definition, we need the following result from \cite[Corollary 11]{LDV} adapted to our setting.

\begin{theorem}\label{th5.1}
Suppose that $\{T_k\}_{k\in \mathbb{N}}$ is a sequence of RB operators of the form \eqref{nonstat1} on $\mathcal{B}(\mathsf{X},\mathsf{F})$. Further suppose that
\begin{enumerate}
\item there exists a nonempty closed invariant set $\mathscr{I}\subseteq\mathcal{B}(\mathsf{X},\mathsf{F})$ for $\{T_k\}_{k\in \mathbb{N}}$;
\item and
\begin{equation}\label{eq5.6}
\sum_{k=1}^\infty\prod_{j=1}^k \Lip (T_j) < \infty.
\end{equation}
\end{enumerate}
Then the backward trajectories $\Psi_k(f_0)$ converge for any initial $f_0\in\mathscr{I}$ to a unique function $\psi\in\mathscr{I}$.
\end{theorem}

With these preliminaries, we obtain the main result of this section.

\begin{theorem}\label{thm4.1}
The backwards trajectories $\{\Psi_k\}_{k\in \mathbb{N}}$ converge for any initial $f_0\in \mathscr{I}$ to a unique function $\psi \in \mathscr{I}$, where $\mathscr{I}$ is the closed ball in $\mathcal{B}(\mathsf{X},\mathsf{F})$ of radius $M/(1-s)$ centered at $0$.
\end{theorem}
\begin{proof}
By Proposition \ref{prop5.2} and Theorem \ref{th5.1} it remains to be shown that \eqref{eq5.6} is satisfied. This, however, follows directly from \eqref{Lip}:
\[
\prod_{j=1}^k \Lip (T_j) \leq s^k\quad\text{and}\quad\sum_{k=1}^\infty s^k = \frac{s}{1-s} < \infty.\qedhere
\]
\end{proof}

A fixed point $\psi$ generated by a sequence $\{T_k\}$ of non-stationary RB operators will be called a \emph{non-stationary fractal function (of class $\mathcal{B}(X,Y)$).}

Employing a non-stationary sequence of RB operators $\{T_k\}$ allows the construction of more general fractal functions exhibiting different local behavior at different scales. This is illustrated by the following example which is taken from \cite{m4}.

\begin{example}
Let $\mathsf{X}:= [0,1]\subset\mathbb{R}$ and $\mathsf{F}:= \mathbb{R}$. Consider the two RB operators
\[
T_1 f (x) := \begin{cases} 2x + \frac12 f(2x), & x\in [0,\frac12),\\
2 - 2x + \frac12 f(2x-1), & x\in [\frac12,1],
\end{cases}
\]
and 
\[
T_2 f (x) := \begin{cases} 2x + \frac14 f(2x), & x\in [0,\frac12),\\
2 - 2x + \frac14 f(2x-1), & x\in [\frac12,1].
\end{cases}
\]
For both operators, $l_i (x) := \frac12 (x + i-1)$, $i=1,2$.

The sequence $T_1^k f \to \tau$, where $\tau$ denotes the Takagi function \cite{Tak} and $T_2^k\to q$, where $q(x)= 4x(1-x)$.

Consider the alternating sequence $\{T_k\}_{k\in \mathbb{N}}$ of RB operators given by
\[
T_k := \begin{cases} T_1, & 10(j-1) < k \leq 10j-5,\\
T_2, & 10j-5 < k \leq 10j,
\end{cases}\quad j\in \mathbb{N}.
\]
Two images of this hybrid attractor of the backward trajectory $\Psi_k$ starting with $f_0 \equiv 0$ are shown in Figure \ref{fig11}.
\begin{figure}[h!]
\begin{center}
\includegraphics[width=5cm, height=3cm]{frac1.pdf}\hspace*{1cm}\includegraphics[width=5cm, height=3cm]{frac2.pdf}
\end{center}
\caption{The hybrid $\tau-q$ attractor. It is smooth at one scale but fractal at another.}\label{fig11}
\end{figure}
\end{example}

\subsection{Non-stationary Fractal Interpolation}

In this subsection, we consider the case $\mathsf{X} := [0,1]$ and $\mathsf{F} := \mathbb{R}$. In the following, we use a sequence $\{T_k\}$ of RB operators of a particular form to obtain a continuous fixed point $\psi$. 

To this end, consider an RB operator defined by \eqref{eq3.17}, choose two arbitrary functions $f,b\in \mathcal{B}(\mathsf{X}) := \mathcal{B}(\mathsf{X}, \mathbb{R})$ and set
\begin{equation}\label{qs}
q_i := f\circ l_i - s_i\cdot b.
\end{equation}
The (stationary) RB operator associated with this particular choice is then given by
\begin{equation}
(Tg)(x) = f + s_i(l_i^{-1}(x))\cdot(g-b)(l_i (x)), \quad x\in l_i(\mathsf{X})
\end{equation}
\begin{remark}
The fixed point of the above defined RB operator may be thought of as the ``fractalization" of a given function $f$.
\end{remark}

To work in a non-stationary setting, let $k\in \mathbb{N}$ and let $\{l_{i_k,k} : i_k\in \mathbb{N}_{n_k}, \, k\in \mathbb{N}\}$ be family of injections from $[0,1]\to [0,1]$ generating a partition of $[0,1]$ in the sense of \eqref{cl1} and \eqref{cl2}. We assume w.l.o.g. that $l_{1,k}(0) = 0$ and $l_{n_k,k}(1) = 1$ and define
\begin{align*}
x_{i_k-1,k} := l_{i_k,k}(0), \quad x_{i_k,k} := l_{i_k,k}(1), \quad i_k\in \mathbb{N}_{n_k}
\end{align*}
where $x_{0,k} :=0$ and $x_{n_k,k} := 1$. By relabelling -- if necessary -- we may further assume that $0 = x_{0,k}  < \cdots < x_{i_k-1,k} < x_{i_k,k} < \cdots x_{n_k,k} = 1$.

Let $f\in \mathcal{C}(\mathsf{X})$ be arbitrary. Define a metric subspace of $\mathcal{C}(\mathsf{X})$ by 
\[
\mathcal{C}_*(\mathsf{X}) := \{g\in \mathcal{C}(\mathsf{X}) : g(0) = f(0) \wedge g(1) = f(1)\}
\] 
and note that $\mathcal{C}_*(\mathsf{X})$ becomes a complete linear metric space when endowed with the metric $d$ that is induced by the sup-norm from $\mathcal{C}(\mathsf{X})$. 

Furthermore, let $b\in \mathcal{C}_*[0,1]$ be the unique affine function whose graph connects the points $(0, f(0))$ and $(1, f(1))$:
\begin{equation}\label{b}
b(x) = (f(1) - f(0)) x + f(0).
\end{equation}
Let $\{\mathcal{P}_k\}_{k\in \mathbb{N}}$ be a family of sets of points in $\mathsf{X}\times \mathsf{F}$ where 
\[
\mathcal{P}_k :=\{(x_{j_k}, f(x_{j,k})\in \mathsf{X}\times \mathsf{F} : j = 0,1,\ldots, n\}.
\]
For $k\in \mathbb{N}$, define an RB operator $T_k: \mathcal{C}_*[0,1]\to \mathcal{C}_*[0,1]$ by
\begin{equation}\label{eq5.1}
T_k g = f + \sum_{i_k=1}^{n_k} s_{i_k,k}\circ l_{i_k,k}^{-1} \cdot (g-b)\circ l_{i_k,k}^{-1}\, \chi_{l_{i_k,k} [0,1]},
\end{equation}
where $\{s_{i_k,k}\}_{i_k=1}^{n_k}\subset \mathcal{C}[0,1]$ such that 
\[
\sup_{k\in \mathbb{N}} \max_{i_k\in \mathbb{N}_{i_k}} \|s_{i_k,k}\|_\infty < 1.
\]
Note that 
\begin{align*}
T_k g (x_{i_k,k}-) &= f(x_{i_k,k}-) + s_{i_k,k}\circ l_{i_k,k}^{-1}(x_{i_k,k}-) \cdot (g-b)\circ l_{i_k,k}^{-1}(x_{i_k,k}-) \\
&= f(x_{i_k,k}) + s_{i_k,k} (1)\cdot (f - b)(1) = f(x_{i_k,k})
\end{align*}
and
\begin{align*}
T_k g (x_{i_k,k}+) &= f(x_{i_k,k}+) + s_{i_k+1,k}\circ l_{i_k+1,k}^{-1}(x_{i_k,k}+) \cdot (g-b)\circ l_{i_k+1,k}^{-1}(x_{i_k,k}+) \\
&= f(x_{i_k,k}) + s_{i_k+1,k} (0)\cdot (f - b)(0) = f(x_{i_k,k}).
\end{align*}
implying that $T_k g$ is continuous at the points $x_{i_k,k}\in [0,1]$:
\[
T_k g (x_{i_k,k}-) = T_k g (x_{i_k,k}+), \quad\forall\,i_k\in\{1, \ldots, n-1\}.
\]
Hence,  $T_k g\in \mathcal{C}_*[0,1]$ and $T_k g$ interpolates $\mathcal{P}_k$ in the sense that
\[
T_k g (x_{i_k,k}) = f(x_{i_k,k}), \quad \forall\,i_k\in \mathbb{N}_{n_k}.
\]

\begin{proposition}
A nonempty closed invariant set for $\{T_k\}_{k\in \mathbb{N}}$ is given by the closed ball in $\mathcal{C}_*(\mathsf{X})$,
\begin{equation}\label{I}
\mathscr{I} = \left\{g\in \mathcal{C}_*[0,1] : \|g\| \leq \frac{\|f\| + s \|b\|}{1-s}\right\},
\end{equation}
where $s$ is defined by \eqref{Lip}.
\end{proposition}

\begin{proof}
From the form \eqref{qs} of the functions $q_{i_k,k}$, we obtain from \eqref{condq} the estimate $\|q_{i_k,k}\| \leq \|f\| + s \|b\|$, which by Proposition \ref{prop2.1} yields the result.
\end{proof}

Together with Theorem \ref{thm4.1}, the above arguments prove the next theorem.

\begin{theorem}
Let $\{T_k\}_{k\in \mathbb{N}}$ be a sequence of RB operators of the form \eqref{eq5.1} each of whose elements acts on the complete metric space $(\mathcal{C}_*(\mathsf{X}), d)$ where $f\in \mathcal{C}_*(\mathsf{X})$ is arbitrary and $b$ is given by \eqref{b}. Furthermore, let the family of functions $\{s_{i_k,k}\} \subset \mathcal{C}(\mathsf{X})$ satisfy \eqref{Lip}. Then, for any $f_0\in \mathscr{I}$, the backward trajectories $\Psi_k (f_0)$ converge to a function $\psi \in \mathscr{I}$ which interpolates $\mathcal{P}_k$.
\end{theorem}

We refer to the fixed point $\psi\in \mathcal{C}_*(\mathsf{X})$ as a \emph{continuous non-stationary fractal interpolation function}.

\begin{remark}
As $f_0$ one may choose $f$ or $b$.
\end{remark}

To illustrate the above results, we present the following example from \cite{m4}.

\begin{example}
Consider the two RB operators $T_i: C[0,1]\to C[0,1]$, $i=1,2$, defined by
\[
(T_1 f)(x) = \begin{cases} -\frac{1}{2}\,f(4x), & x\in[0,\frac{1}{4}),\\
 -\frac{1}{2} + \frac{1}{2}\,f(4x-1), & x\in[\frac{1}{4},\frac{1}{2}),\\
 \frac{1}{2}\,f(4x-2), & x\in[\frac{1}{2},\frac{3}{4}),\\
 \frac{1}{2} + \frac{1}{2}\,f(4x-3), & x\in[\frac{3}{4},1],\end{cases}
\]
and 
\[
(T_2 f)(x) := \begin{cases} \frac34 f(2x), & x \in [0,\frac12),\\
\frac34 + \frac14 f(2x-1), &  x \in [\frac12,1].
\end{cases}
\]
The RB operators $T_1$ and $T_2$ generate \emph{Kiesswetter's fractal function} \cite{Kiess}, respectively, a \emph{Casino function} \cite{DS}. 

Consider again the alternating sequence $\{T_i\}_{i\in \mathbb{N}}$ of RB operators given by
\[
T_k := \begin{cases} T_1, & 10(j-1) < k \leq 10j-5,\\
T_2, & 10j-5 < k \leq 10j,
\end{cases}\quad j\in \mathbb{N}.
\]
Two images of the hybrid attractor of the backward trajectory $\Psi_k$ starting with the function $f_0 (x) = x$, $x\in\mathsf{X}$, are shown below in Figure \ref{fig22}.
\begin{figure}[h!]
\begin{center}
\includegraphics[width=5cm, height=4cm]{frac4.pdf}\hspace*{1cm}\includegraphics[width=5cm, height=4cm]{frac3.pdf}
\end{center}
\caption{The hybrid Kiesswetter-Casino attractor.}\label{fig22}
\end{figure}

\end{example}

\section{Quaternionic Fractal Interpolation}
In this section, we extend fractal interpolation to a quaternionic setting. As quaternions from a non-commutative division algebra, the non-commutativity generates more intricate fractal patterns. 

First, we give a short introduction to quaternions and present some of those properties that are relevant for the remainder of this section. The interested reader is referred to the literature on the subject, a short and subjective list of which is \cite{Bour,GHS,GS,Krav,MGS}.


\subsection{A Brief Introduction to Quaternions}
Let $\{e_1, e_2, e_3\}$ be the canonical basis of the Euclidean vector space $\mathbb{R}^3$. We call $\{e_1, e_2, e_3\}$ imaginary units and require that the following multiplication rules hold:
\begin{gather}
e_1^2=e_2^2=e_3^2=-1,\label{6.1}\\ 
e_1e_2=e_3=-e_2e_1,\quad e_2e_3=e_1=-e_3e_2,\quad e_3e_1=e_2=-e_1e_3.\label{6.2}
\end{gather}
\begin{remark}
Note that \eqref{6.2} is equivalent to $e_1e_2e_3 = -1$.
\end{remark}
A \emph{real quaternion $q$} is then an expression of the form
\[
q = a+\sum_{i=1}^3v_ie_i, \quad a, v_1, v_2, v_3\in\mathbb{R}.
\]
The addition and multiplication of two quaternions $q_1 = a+\sum_{i=1}^3v_ie_i$ and $q_2 = b+\sum_{i=1}^3w_ie_i$ is defined by
\begin{align*}
q_1 + q_2 &:= (a+b) + \sum_{i=1}^3 (v_i+w_i) e_i\\
q_1 q_2&:= (ab - v_1w_1 - v_2w_2 - v_3w_3)\\
&\;\;+(a {w_1}+b {v_1}+{v_2} {w_3}-{v_3} {w_2})e_1\\
&\;\;+(a {w_2}+b {v_2}-{v_1} {w_3}+{v_3} {w_1})e_2\\
&\;\;+(a {w_3}+b {v_3}+{v-1}{w_2}-{v_2} {w_1})e_3.
\end{align*}
Each quaternion $q=a+\sum\limits_{i=1}^3v_ie_i$ may be decomposed as $q=\Sc (q)+\Ve (q)$ where $\Sc (q)=a$ is the {\it scalar part} of $q$ and $\Ve (q)=v=\sum\limits_{i=1}^3v_ie_i$ is the {\it vector part} of $q$. $v=\Ve (q)$ is also called a \emph{quaternionic vector}.

The {\it conjugate} $\overline{q}$ of the real quaternion $q=a+v$ is the quaternion $\overline{q}=a-v$. Note that $q\overline{q}=\overline{q}q=a^2+|v|^2=a^2+\sum\limits_{i=1}^3v_i^2$. Therefore, we can define a norm on $\mathbb{H}$ by setting
\[
\abs{q} := \sqrt{q\overline{q}}.
\]
The inverse of a quaternion $q$ is given by
\[
q^{-1} = \frac{\overline{q}}{\abs{q}^2}.
\]

It is straight-forward to establish that the collection of all real quaternions
\[
\mathbb{H} := \mathbb{H}_\mathbb{R} := \left\{a+\sum_{i=1}^3v_ie_i : a, v_1, v_2, v_3\in{\mathbb R}\right\},
\]
is a four-dimensional associative normed division algebra over $\mathbb{R}$. Due to the multiplication rules \eqref{6.1} and \eqref{6.2}, $\mathbb{H}$ is not commutative. We also note that $\mathbb{H}$, as already indicated above, is a four-dimensional vector space over $\mathbb{R}$ with basis $\{e_0,e_1,e_2,e_3\}$, where $e_0 := 1$.

\begin{remark}
We note that $\mathbb{H}$ and $\mathbb{R}^4$ are identical as \emph{point sets} but differ in their algebraic structure. The vector space $\mathbb{R}^4$ is not an algebra whereas $\mathbb{H}$ is one, albeit non-commutative.
\end{remark}

\begin{remark}
As we can define a norm and therefore a metric on $\mathbb{H}$, $\mathbb{H}$ becomes a topological space where open sets are defined via the metric. All other topological concepts such as limits, convergence, compactness etc. then follow. $\mathbb{H}$ thus becomes a topological space and, moreover, also a complete metric space.
\end{remark}

We also remark that if $v=\sum\limits_{j=1}^3v_je_j$ and $w=\sum\limits_{j=1}^3w_je_j$ are quaternionic vectors, then 
\begin{equation}
vw=-\langle v,w\rangle +v\wedge w\label{vec mult},
\end{equation}
where $\langle v,w\rangle :=\sum\limits_{j=1}^3v_jw_j$ is the scalar product of $v$ and $w$ and 
$$v\wedge w :=(v_2w_3-v_3w_2)e_1+(v_3w_1-v_1w_3)e_2+(v_1w_2-v_2w_1)e_3$$ 
is the vector (cross) product of $v$ and $w$.

For our purposes, we need to introduce the analog of a Banach space in the quaternionic setting. To this end, we begin with the following definitions. (See, also \cite{Bour}.)
\begin{definition}
A real vector space $\mathsf{V}$ is called a left quaternionic vector space if it is a left $\mathbb{H}$-module, i.e., if there exists a mapping $\mathbb{H}\times\mathsf{V}\to \mathsf{V}$, $(q,v)\mapsto q v$ which satisfies
\begin{enumerate}
\item $\forall v\in \mathsf{V}\,\forall q_1,q_2\in \mathbb{H}:$ $(q_1+q_2)v = q_1 v + q_2 v$.
\item $\forall v_1,v_2\in \mathsf{V}\,\forall q\in \mathbb{H}:$ $q (v_1+v_2) = q v_1 + q v_2$.
\item $\forall v\in \mathsf{V}\,\forall q_1,q_2\in \mathbb{H}:$ $q_1(q_2 v) = (q_1 q_2) v$.
\end{enumerate}
\end{definition}

\begin{remark}
In analogous fashion, one defines a \emph{right} quaternionic vector space as a right $\mathbb{H}$-module where the mapping is now $\mathsf{V}\times\mathbb{H}\to V$, $(v,q)\mapsto v q$.
\end{remark}

A \emph{two-sided quaternionic vector space} $\mathsf{V}$ is a left and right quaternionic vector space such that $\lambda v = v \lambda$, for all $\lambda\in\mathbb{R}$ and for all $v\in\mathsf{V}$. An example of a two-sided quaternionic vector space is given by $\mathbb{H}$ itself.

\begin{remark}
A quaternionic vector space becomes a real vector space when its scalar multiplication is restricted to $\mathbb{R}$.
\end{remark}

One can start with any real vector space $\mathsf{V}_\mathbb{R}$ and construct a two-sided quaternionic vector space by setting
\[
\mathsf{V}_\mathbb{H} := \left\{\sum_{i=0}^3 v_i\otimes e_i : v_i\in \mathsf{V}_\mathbb{R}\right\},
\]
where $\otimes$ denotes the algebraic tensor product. (See, for instance, \cite{Ng}.) 

On the other hand, given any two-sided quaternionic vector space and defining
\[
\mathsf{V}_\mathbb{R} := \{v\in \mathsf{V} : \lambda v = v \lambda,\;\forall\lambda\in\mathbb{H}\},
\]
then $\mathsf{V}_\mathbb{R}$ is a real vector space called the \emph{real part of $\mathsf{V}$}.

The proof of the next result can be found in, i.e., \cite{Ng}.

\begin{proposition}
Let $\mathsf{V}$ be a two-sided quaternionic vector space and let $\mathsf{V}_\mathbb{R}$ denote its real part. Then $V\cong \mathsf{V}_\mathbb{R}\otimes\mathbb{H}$.
\end{proposition}

Next, we introduce quaternionic normed spaces.

\begin{definition}
Let $\mathsf{V}$ be a left quaternionic vector space. A function $\n{\cdot}: \mathsf{V}\to \mathbb{R}_0^+$ is called a norm on $V$ if
\begin{enumerate}
\item $\n{v} = 0$ iff $v=0$.
\item $\n{q v} = \abs{q} \n{v}$, for all $v\in \mathsf{V}$ and $q\in \mathbb{H}$.
\item $\n{v+w} \leq \n{v} + \n{w}$, for all $v,w\in \mathsf{V}$.
\end{enumerate}
A left quaternionic vector space endowed with a norm will be called a left quaternionic normed space. 
\end{definition}

A left quaternionic normed space $\mathsf{E}$ is called \emph{complete} if it is a complete metric space with respect to the metric $d(x,y) = \n{x-y}$ induced by the norm $\n{\cdot}_\mathsf{E}$. In this case, we refer to $\mathsf{E}$ as \emph{left quaternionic Banach space}.

\begin{remark}
A left (or right) quaternionic Banach space becomes a real Banach space if the left (right) scalar multiplication is restricted to $\mathbb{R}$. (Cf., \cite[Section 2.3]{CGK}.)
\end{remark}

\begin{example}
The space $\mathbb{H}^k$ consisting of $k$-tuples of quaternions is both a left and a right quaternionic vector space. We represent elements $\xi\in\mathbb{H}^k$ as column vectors and define the quaternionic conjugate ${}^*$ of $\xi$ by
\[
\begin{pmatrix} \xi_1\\ \vdots\\ \xi_k
\end{pmatrix}^* := \begin{pmatrix} \overline{\xi_1} & \cdots & \overline{\xi_k}
\end{pmatrix},
\]
where each $\xi_j\in\mathbb{H}$. When endowed with the norm
\begin{equation}\label{dist}
\n{\xi}_k := \sqrt{\xi^*\xi} = \sqrt{\sum\limits_{l=1}^k \abs{\xi_j}^2},
\end{equation}
$\mathbb{H}^k$ becomes a two-sided quaternionic Banach space as $\lambda v = v \lambda$, for all $\lambda\in\mathbb{R}$ and for all $v\in\mathbb{H}^k$.

Note that under the norm $\n{\cdot}_k$, $\mathbb{H}^k$ becomes a topological space and also a complete metric space.
\end{example}

The final concept we need to introduce is that of left linear mapping between left quaternionic vector spaces.

\begin{definition}
Let $\mathsf{V}_1$ and $\mathsf{V}_2$ be left quaternionic vector spaces. A mapping $f:\mathsf{V}_1\to\mathsf{V}_2$ is called left linear if
\[
f(q v + w) = q f(v) + f(w), \quad\forall v,w\in \mathsf{V}, \forall q\in \mathbb{H}.
\]
A left linear mapping is called bounded if
\[
\n{f} := \sup_{x,y \in \mathsf{V}_1, x \neq y} \frac{\n{f(x)-f(y)}_{\mathsf{V}_2}}{\n{x-y}_{\mathsf{V}_1}} < \infty.
\]
\end{definition}
\subsection{Quaternionic Fractal Interpolation}
In this section, we introduce the novel concept of quaternionic fractal interpolation. For illustrative purposes, we do not choose the most general set up but restrict ourselves to the case where the left quaternionic Banach spaces are $\mathsf{E} := \mathbb{H}^k =:\mathsf{F}$, $k\in \mathbb{N}$.

For our purposes, we need the following function space. Let $\mathsf{X}\subset\mathbb{H}^k$ be compact (as defined via the norm $\n{\cdot}_k$) and let
\[
\mathcal{B}(\mathsf{X},\mathbb{H}^k) := \left\{f:\mathsf{X}\to\mathbb{H}^k : \text{$f$ is bounded}\right\}.
\]
A function $f:\mathsf{X}\to\mathbb{H}^k$ is called \emph{bounded} if there exists a real number $M>0$ such that $\n{f}_k \leq M$. If we define for $x\in\mathsf{X}$ and $\lambda\in\mathbb{H}$
\[
(f+g) (x) := f(x) + g(x)\quad\text{and}\quad (\lambda\cdot f)(x) := \lambda\cdot f(x)
\]
then $\mathcal{B}(\mathsf{X},\mathbb{H}^k)$ becomes a left quaternionic vector space. Setting for each $f\in \mathcal{B}(\mathsf{X},\mathbb{H}^k)$
\[
\n{f} := \sup_{x\in\mathsf{X}} \n{f(x)}_k,
\]
then $\mathcal{B}(\mathsf{X},\mathbb{H}^k)$ becomes a left Banach space. (See, for instance, \cite[Chapter IV.E.1.]{J} for the case $k=1$. The extension to $k > 1$ is straight-forward.)

For the sake of simplicity and the purpose of understanding the underlying issues, we concentrate on the special case that 
\[
\mathsf{X} := \left\{q\in \mathbb{H} : \max\limits_{i=0,1,2,3} \abs{q_i} = 1\right\} \cong [-1,1]^4
\]
and that $k:=1$. The interested reader is encouraged to consider the general set-up.

We take the four-dimensional cube $\mathsf{X}$ and divide it into $n:=2^4$ congruent four-dimensional subcubes $\mathsf{X}_i$ each similar to $\mathsf{X}$ and such that $\{\mathsf{X}_i\}_{i=1}^{n}$ forms a partition of $\mathsf{X}$ in the sense of \eqref{c1} and \eqref{c2}. We leave it to the diligent reader to derive closed expressions for the $n$ injections $l_i$.

On the left Banach space $\mathcal{B}(\mathsf{X},\mathbb{H})$, we consider the RB operator $T: \mathcal{B}(\mathsf{X},\mathbb{H})\to\mathcal{B}(\mathsf{X},\mathbb{H})$ given by
\begin{equation}\label{6.5}
Tf(l_i(x)) := q_i(x) + s_i(x) f(x), \quad x\in \mathsf{X}, \;\;i \in \mathbb{N}_n,
\end{equation}
where $q_i, s_i:\mathsf{X}\to\mathbb{H}$ are bounded functions. Clearly, $T f\in \mathcal{B}(\mathsf{X},\mathbb{H})$ under these assumptions. Note that we left-multiply $f$ by $s_i$. 

\begin{remark}
In the case of real Banach spaces, the RB operator is affine, i.e., $T - T(0)$ is a linear operator. In the quaternionic setting this is no longer true: $T - T(0)$ is not a left linear operator. It is only if $s_i (x) \in \mathbb{R}$.
\end{remark}

To simplify notation, we introduce the following abbreviation for the $m$-fold composition of functions from an IFS. Let $\mathcal{F}:=\{f_1, \ldots, f_n\}$. We write
\[
f_{i_m i_{m-1} \cdots i_1} := f_{i_m}\circ f_{i_{m-1}}\circ f_{i_1},
\]
where each $i_j\in \mathbb{N}_n$.

For each $m\in \mathbb{N}$, the $m$-fold application of $T$, $T^m f := T(T^{m-1} f)$, can be written as
\[
T^m f (l_{i_m i_{m-1} \cdots i_1}(x)) = \sum_{k=1}^{m} \prod_{j=1}^{k-1} s_{i_j} (x) q_{k} (x)+ \prod_{k=1}^m s_{i_k}(x) f(x),
\]
where the factors in the products $\prod$ are left-multiplied and where we set the empty product equal to 1. Notice the reverse order of the indices $i_m i_{m-1} \cdots i_1$ on the left- and right-hand side.

Continuing as in the previous section to ensure the existence of a fixed point $\psi$ for the RB operator \eqref{6.5} or a bounded solution of the associated system of functional equations, we obtain in summary the next result.
\begin{theorem}
For the above setting, the RB operator $T$ defined in \eqref{6.5} has a unique fixed point $\psi\in\mathcal{B}(\mathsf{X},\mathbb{H})$, i.e.,
\begin{equation}
T\psi = \psi\quad\Longleftrightarrow\quad\psi (l_i(x)) = q_i(x) + s_i(x) \psi(x), \quad x\in \mathsf{X}, \;\; i \in \mathbb{N}_n,
\end{equation}
provided that 
\[
\max\limits_{i\in \mathbb{N}_n} \sup\limits_{x\in \mathsf{X}_i} |s_i(x)| < 1.
\]
\end{theorem}
The fixed point $\psi$ is then called a \emph{bounded quaternionic fractal function}. 

\begin{remark}
A similar result exists of course for right Banach spaces and RB operators where the functions $s_i$ are right-multiplied onto $f$. Note that the fixed point $\psi$ depends on this right or left multiplication.
\end{remark}

\begin{example}
This example shows the versatility of quaternionic fractal interpolation. We choose as $\mathsf{X} := \{q\in\mathbb{H} : \Sc q \in [0,1) \wedge \Ve q = 0\}$. Note that $\mathsf{X} \cong [0,1] \subset\mathbb{R}$. Define injections $l_i:\mathsf{X}\to\mathsf{X}$ as follows:
\[
l_1 (x) : = \tfrac12 x\quad \text{and} \quad l_2 (x):= \tfrac12(x+1).
\]
Clearly, $\{l_1(\mathsf{X}), l_2(\mathsf{X})\}$ is a partition of $\mathsf{X}$ satisfying \eqref{c1} and \eqref{c2}.

Moreover, let $q_1 := e_0 + 2e_1-e_3+3e_4$ and $q_2 := -e_0 - 2e_1 + 2e_3+e_4$ be two quaternions. We set
\[
q_1 (x) := (1-q_1) x \quad \text{and} \quad q_2(x) := q_2 x^2. 
\]

Then, $q_1,q_2\in \mathcal{B}(\mathsf{X},\mathbb{H})$. We define an RB operator $T$ by
\[
Tf (x) := \begin{cases} 2(1-q_1) x + s_1 f(2x), & x\in [0,\frac12),\\
q_2(2x-1)^2 + s_2 f(2x-1), & x\in [\frac12,1),
\end{cases}
\]
where $s_1 := \frac{1}{10}e_0 + \frac12 e_1 - \frac15 e_2 - \frac{1}{10} e_3$ and $s_1 := -\frac{1}{5}e_0 + \frac15 e_1 - \frac35 e_2 + \frac{1}{10} e_3$. Note that $\abs{s_1} = \frac{1}{10}\sqrt{31}$ and $\abs{s_2} = \frac{1}{10}\sqrt{45}$. Thus, $\max\{s_1,s_2\} < 1$.

Therefore, the RB operator $T$ is a contraction on $\mathcal{B}(\mathsf{X},\mathbb{H})$ and has a unique fixed point $\psi\in\mathcal{B}(\mathsf{X},\mathbb{H})$. The following Figure \ref{fig5} shows the projection of the graph of $\psi$ onto the $(e_0,e_i)$-plane for $i=0,1,2,3$.

\begin{figure}[h!]
\begin{center}
\includegraphics[width=5cm, height=4cm]{qf1.pdf}\hspace*{1cm}\includegraphics[width=5cm, height=4cm]{qf2.pdf}\\
\includegraphics[width=5cm, height=4cm]{qf3.pdf}\hspace*{1cm}\includegraphics[width=5cm, height=4cm]{qf4.pdf}
\end{center}
\caption{The projections of $\psi$ onto the $(e_0,e_i)$-planes.}\label{fig5}
\end{figure}
As $\psi$ can be written as $\psi = \sum\limits_{i=0}^3 \psi_i e_i$, we display in Figure \ref{fig6} the parametric plots $(\psi_0, \psi_1, \psi_2)$ and $(\psi_0, \psi_2, \psi_4)$.
\begin{figure}[h!]
\begin{center}
\includegraphics[width=5cm, height=4cm]{qf123.pdf}\hspace*{1cm}\includegraphics[width=5cm, height=4cm]{qf124.pdf}\end{center}
\caption{Some parametric plots of the components of $\psi$.}\label{fig6}
\end{figure}
\end{example}
\section{Summary}
In this chapter, we considered several aspects of fractal interpolation and identified their commonalities. Fractal interpolation is based either on solving a system of functional equations \cite{SB,SB1} or on obtaining the fixed point of an operator \cite{massopust,massopust1}. Both approaches yield the same results and are interchangeable. 

In particular, we presented
\begin{itemize}
\item global fractal interpolation. This was the original approach undertaken in \cite{B2} which is based on a geometric setting. The analytic setting based on an RB operator commenced in \cite{bedford,dubuc1,dubuc2} and produced numerous generalizations and extensions. The main idea is to use a fixed bounded subset $\mathsf{X}$ of a normed space $\mathsf{E}$ to produce a partition of $\mathsf{E}$ over which copies of functions defined on the global set $\mathsf{X}$ reside. Each of the partitioning sets is then iteratively partitioned into subsets thus producing the fractal nature of the limiting object.
\item local fractal interpolation. Here, we allow subsets of a bounded subset $\mathsf{X}$ of a normed space $\mathsf{E}$ to produce a partition of $\mathsf{X}$. These subsets can be repeated in the construction of the partition and need no longer be defined by contractive injections. Here, functions defined on the local subsets reside over partitions of $\mathsf{X}$. This type of local interpolation was successfully employed in fractal compression and image analysis. (Cf., for instance, \cite{BH}.)
\item non-stationary fractal interpolation. In this approach, we do no longer require that we keep the same quantities at each level of iteration but allow them to vary. This set-up is very similar to non-stationary subdivision and was investigated in this respect in \cite{DLM,LDV,m4}. This new methodology shows great potential for future investigations into the subject. 
\item quaternionic fractal interpolation. In this novel approach, we kept the set-theoretic structure but allowed for a change in the algebraic structure of the underlying sets. By requiring that we are using a division algebra such as the quaternions or more generally hypercomplex numbers, we obtain even more flexibility in the construction of fractal interpolants. The investigation into this type of setting commenced in \cite{m5} and promises to be a fruitful field of investigation in years to come.
\end{itemize}
