1 Introduction

The theory of fractal interpolation induces new insights into the fields of integration and approximation [16]. It is considered an analogue to the conventional interpolation schemes. The functions generated through fractal interpolation are typically used to analyse signals with complex characteristics [4]. The self similarity property of fractals enables them to analyse the naturally occurring signals, extract information from them, and process the obtained information more properly and accurately [15, 22].

Unlike traditional interpolation methods, the concept of fractal interpolation proceeds by creating the copies of images at different scales. Fractal interpolation functions are accompanied by a complete metric space and a finite collection of contractions defined on the complete metric space [6]. The set of finite contractions is known as the iterated function system (IFS) for the fractal interpolation functions. The coefficients of the IFS are crucial in determining the shape of the signal and also to evaluate the integral results. In order to exemplify the one-dimensional chaotic signals, Barnsley [5] introduced the idea of fractal interpolation functions. The concept is then extended to two dimensions by Massopust [19] for the first time, following which there have been several attempts by authors to construct fractal interpolation functions over two-dimensional regions. The major difficulty in defining a fractal interpolation function over a higher dimensional region is the lack of continuity of the related fractal interpolation operator. In order to tackle this problem, Dalla [10] used collinear interpolation points for the construction. Without restricting the interpolation points, continuous fractal surfaces are constructed in [13] over polygonal regions. Redefining the IFS, Drakopoulos et al. [12] solved the problem of continuity. Another attempt by Ruan et al. [25] includes modifying the endpoint conditions of the contraction mappings and adding suitable criteria to the IFS. In order to make the fractal interpolation operator continuous, [30] uses a continuous contraction function, instead of the freely chosen contractivity factors. The construction proposed in [20] uses the tensor product of single variable fractal interpolation functions. The transformations like rotations and reflections have been used in [23] to create continuous fractal surfaces and volumes over rectangular lattices. The stability of the fractal interpolation functions have been studied in [14]. Instead of using Banach contraction mapping theorem, [24] employs Rakotch fixed point theorem and Matkowski fixed point theorem for the construction. Malysz [17] rectified the problem of continuity by formulating a new IFS, along with changing the endpoint conditions of the contractions. The construction we present here is based on the approach made by Malysz [17]. Instead of considering a constant vertical scaling factor as in [17], a formula has been proposed for the same in our work.

The theory of fractal interpolation has been used in a variety of applications such as curve fitting [18], image compression [7, 11, 26], tumor perfusion [9] and Covid-19 data analysis [2, 28]. Since most naturally occurring signals have sudden peaks and valleys in their structure, the Euclidean measurements of specific parameters such as length, area, and perimeter become irrelevant. However, it is appropriate to use the fractal numerical integration method to evaluate such parameters. In [27], authors propose a new method of numerical integration using the trapezoidal rule wherein the signal has been split into fractal fragments. The observed integral value corresponds to the energy needed on the electrical load curve.

The classical integral of a linear fractal interpolation function is investigated in [1]. It is trivial that the integration of a function is associated with the area under the graph of the function. When the function is not known explicitly, it is difficult to plot the functions directly. But, several numerical methods are available to interpolate the data and find the approximation function that passes through the given points. Techniques are also available to numerically evaluate the integral results [3] from the given data values without explicitly drawing the graph and calculating the area under the graph. The most popular techniques include the trapezoidal rule, midpoint rule, and Simpson’s rule. Another method of definite integral calculation using triangles is proposed in [8]. However, it is true that the traditional numerical schemes may fail to portray and numerically evaluate the data precisely, especially when signals of chaotic nature are involved. For example, in the trapezoidal rule, the data set has been approximated with straight lines, and the integral is calculated as the area of the trapezoid under the straight lines. Approximating the non linear data sets with straight lines may fail to capture the sudden spikes in the signals. As a result, the area calculated under the straight lines may be far from the desired value. Hence, the motivation behind choosing fractal interpolation techniques in evaluating numerical integral results can be summarized as follows:

  • Since the fractal interpolation functions capture the irregularities in the signals, the area as well as volume calculated using the fractal numerical integration method will be more precise.

  • Fractal numerical integration provides the desired results in fewer iterations.

The method of fractal numerical integration for the data sets corresponding to single variable functions was introduced in [21]. The integration formula is based on the coefficients of the iterated function systems. These coefficients are derived using the endpoint conditions of the iterated function systems. The integral results are compared with the conventional numerical integration schemes in [21]. This paper aims to extend the fractal numerical integration method proposed in [21] to two-dimensional signals defined over rectangular interpolating regions.

The integration formula has been defined based on the construction proposed in [17] and considering the coefficients of the IFS therein. Instead of using a constant vertical scaling factor, this work suggests a formula for determining the same. The formula has been derived from the given data values. The computations are done using the newly obtained formula of the vertical scaling factor. This paper then illustrates the derived numerical integration formula, considering the data set of four test functions. The obtained results are tabulated and compared with their actual integral values.

The paper has been organized into eight sections, where the first section briefly introduces the theory of fractal interpolation functions and numerical integration. The second section deals with constructing bivariate fractal interpolation functions over rectangular interpolating domains. In the third section, we provide the proposed method of numerical integration. A recursive relation connecting the bivariate fractal interpolation function with the bilinear interpolation function is derived in the fourth section. The fifth section deals with the method of selection of the vertical scaling factor that minimizes the approximation errors. The convergence of the proposed method is verified through certain theorems in section six. The illustration of the technique through sufficient examples is given in the seventh section. The paper concludes in section eight by summarizing the results and observations obtained.

2 Construction of bivariate fractal interpolation functions over rectangular regions

Let \(\{x_{n}:n=0,1,...,N\}\) and \(\{y_{m}:m=0,1,...,M\}\) be two sets of real numbers such that \(0=x_{0}<x_{1}<...<x_{N}=1,\) and \(0=y_{0}<y_{1}<...<y_{M}=1.\) Let I be the closed interval containing the former set and J be the closed interval containing the latter set of real numbers. The data set considered be denoted by \(\{(x_{n}, y_{m}, z_{n,m}):n=0,1,...,N,\,\, m=0,1,...,M\}\) where \(z_{n,m}\) is the value of the function at the point \((x_{n},y_{m}).\) Set \(I_{n}=[x_{n-1}, x_{n}], \,\, J_{m}=[y_{m-1},y_{m}]\) for \(n=1,2,...,N\) and \(m=1,2,...,M.\) Define N number of contractive homeomorphisms \(L_{n}:I \rightarrow I_{n}\) such that

$$\begin{aligned}{} & {} L_{n}(x_{0}) = \left\{ \begin{array}{lr} x_{n-1}, &{} \text {if } n \,\,\text {is odd }\\ x_{n}, &{} \quad \text {if } n \,\,\text {is even } \end{array} \right. \end{aligned}$$
(1)
$$\begin{aligned}{} & {} L_{n}(x_{N}) = \left\{ \begin{array}{lr} x_{n}, &{} \text {if } n \,\,\text {is odd }\\ x_{n-1}, &{} \quad \text {if } n \,\,\text {is even } \end{array} \right. \end{aligned}$$
(2)

and

$$\begin{aligned}{} & {} \mid L_{n}(x)-L_{n}(x^{'}) \mid \le c_{1} \mid x-x^{'} \mid ,\,\,\\{} & {} \quad \text {for all} \,\, x,\,\, x^{'} \in [x_{0},x_{N}], \text {for some}\,\, c_{1} \in [0,1)\\ & \quad \text {and} \,\, n=1,2,...,N. \end{aligned}$$

Similarly, consider M number of contractive homeomorphisms \(K_{m}:J \rightarrow J_{m}\) with the conditions

$$\begin{aligned}{} & {} K_{m}(y_{0}) = \left\{ \begin{array}{lr} y_{m-1}, &{} \text {if } m \,\,\text {is odd }\\ y_{m}, &{} \quad \text {if } m \,\,\text {is even } \end{array} \right. \end{aligned}$$
(3)
$$\begin{aligned}{} & {} K_{m}(y_{M}) = \left\{ \begin{array}{lr} y_{m}, &{} \text {if } m \,\,\text {is odd }\\ y_{m-1}, &{} \quad \text {if } m \,\,\text {is even } \end{array} \right. \end{aligned}$$
(4)

and

$$\begin{aligned}{} & {} \mid K_{m}(y)-K_{m}(y^{'}) \mid \le c_{2} \mid y-y^{'} \mid , \,\, \\{} & {} \quad \text {for all} \,\, y, \,\,y^{'}\,\, \in [y_{0},y_{M}], \,\,\\{} & {} \quad \text {for some} \,\, c_{2} \in [0,1) \,\, \text {and} \,\, m=1,2,...,M. \end{aligned}$$

The functions \(L_{n}\) and \(K_{m}\) used here, satisfying the above stated conditions are defined as follows:

$$\begin{aligned} L_{n}(x)=r_{n}x+s_{n}, \,\, K_{m}(y)=q_{m}y+w_{m} \end{aligned}$$
(5)

where the coefficients \(r_{n}, s_{n}, q_{m}\) and \(w_{m},\) obtained by solving the endpoint conditions (1), (2), (3) and (4) are as follows:

$$\begin{aligned}&r_{n}= \frac{(-1)^{\sigma (n+1)}}{N}, \,\, s_{n}=\frac{n-\sigma (n)}{N} \nonumber \\&\quad q_{m}= \frac{(-1)^{\sigma (m+1)}}{M}, \,\, w_{m}=\frac{m-\sigma (m)}{M} \end{aligned}$$
(6)

where \(\sigma (k)=k \,\,mod \,\,2,\) for \(k\in {\mathbb {N}}.\) \(\,\,i.e,\) \(\sigma (k)\) is the remainder when k is divided by 2.

Choose a number \(d_{n,m}\) between -1 and 1. This number is known as the vertical scaling factor and the corresponding scale vector is denoted as \(\overline{d_{n,m}}.\) Define \(F=[0,1]\times [0,1]\times R\) and \(F_{n,m}:F \rightarrow R\) of the form

$$\begin{aligned} F_{n,m}(x,y,z)&= d_{n,m}z+q_{n,m}(x,y) \end{aligned}$$
(7)

where

$$\begin{aligned} q_{n,m}(x,y)&= a_{n,m}x+b_{n,m}y+c_{n,m}xy+f_{n,m} \end{aligned}$$
(8)

satisfying the conditions

$$\begin{aligned} F_{n,m}(\sigma (k), \sigma (l), z_{\Delta (k,l)}) = z_{k,l} \end{aligned}$$
(9)

for \((k,l) \in \{n-1,n\} \times \{m-1,m\}, n=1,2,...,N\) and \(m=1,2,...,M\) where \(\Delta (k,l)=(N\sigma (k), M\sigma (l)).\)

The coefficients in \(F_{n,m},\) except the vertical scaling factor \(d_{n,m},\) are obtained by solving the endpoint conditions (9). Since there are four different sets of endpoint conditions, based on the indices n, m are odd or even, there are four different expressions for each of the coefficients \(a_{n,m}, b_{n,m}, c_{n,m}\) and \(f_{n,m}.\) Therefore, the function \(F_{n,m}\) can be written as:

$$\begin{aligned} F_{n,m}(x,y,z) = \left\{ \begin{array}{lr} a_{n,m}^{(1)}x+b_{n,m}^{(1)}y+c_{n,m}^{(1)}xy+f_{n,m}^{(1)}+d_{n,m}z,&{} \text {if } n, m \,\,\text {are odd } \\ a_{n,m}^{(2)}x+b_{n,m}^{(2)}y+c_{n,m}^{(2)}xy+f_{n,m}^{(2)}+d_{n,m}z,&{} \text {if } n, m \,\,\text {are even }\\ a_{n,m}^{(3)}x+b_{n,m}^{(3)}y+c_{n,m}^{(3)}xy+f_{n,m}^{(3)}+d_{n,m}z,&{} \quad \text {if } n \,\,\text {is odd and} \,\,m \,\,\text {is even}\\ a_{n,m}^{(4)}x+b_{n,m}^{(4)}y+c_{n,m}^{(4)}xy+f_{n,m}^{(4)}+d_{n,m}z,&{} \quad \text {if } n \,\,\text {is even and} \,\,m \,\,\text {is odd}\\ \end{array} \right. \end{aligned}$$
(10)

with the following expressions for the coefficients:

$$\begin{aligned}& a_{n,m}^{(1)}=z_{n,m-1}-z_{n-1,m-1}-d_{n,m}(z_{N,0}-z_{0,0}) \nonumber \\&\quad b_{n,m}^{(1)}=z_{n-1,m}-z_{n-1,m-1}-d_{n,m}(z_{0,M}-z_{0,0}) \nonumber \\&\quad c_{n,m}^{(1)}=z_{n,m}-z_{n-1,m}-z_{n,m-1}+z_{n-1,m-1}\\ & \quad -d_{n,m}(z_{N,M}-z_{0,M}-z_{N,0}+z_{0,0}) \nonumber \\&\quad f_{n,m}^{(1)}=z_{n-1,m-1}-d_{n,m}z_{0,0} \end{aligned}$$
(11)
$$\begin{aligned}&a_{n,m}^{(2)}=z_{n-1,m}-z_{n,m}-d_{n,m}(z_{N,0}-z_{0,0}) \nonumber \\&\quad b_{n,m}^{(2)}=z_{n,m-1}-z_{n,m}-d_{n,m}(z_{0,M}-z_{0,0}) \nonumber \\&\quad c_{n,m}^{(2)}=z_{n,m}-z_{n-1,m}-z_{n,m-1}+z_{n-1,m-1}\\ & \quad -d_{n,m}(z_{N,M}-z_{0,M}-z_{N,0}+z_{0,0}) \nonumber \\&\quad f_{n,m}^{(2)}=z_{n,m}-d_{n,m}z_{0,0} \end{aligned}$$
(12)
$$\begin{aligned}&a_{n,m}^{(3)}=z_{n,m}-z_{n-1,m}-d_{n,m}(z_{N,0}-z_{0,0}) \nonumber \\&\quad b_{n,m}^{(3)}=z_{n-1,m-1}-z_{n-1,m}-d_{n,m}(z_{0,M}-z_{0,0}) \nonumber \\&\quad c_{n,m}^{(3)}=z_{n-1,m}-z_{n,m}-z_{n-1,m-1}+z_{n,m-1}\\ & \quad -d_{n,m}(z_{0,0}-z_{N,0}-z_{0,M}+z_{N,M}) \nonumber \\&\quad f_{n,m}^{(3)}=z_{n-1,m}-d_{n,m}z_{0,0} \end{aligned}$$
(13)
$$\begin{aligned}&a_{n,m}^{(4)}=z_{n-1,m-1}-z_{n,m-1}-d_{n,m}(z_{N,0}-z_{0,0}) \nonumber \\&\quad b_{n,m}^{(4)}=z_{n,m}-z_{n,m-1}-d_{n,m}(z_{0,M}-z_{0,0}) \nonumber \\&\quad c_{n,m}^{(4)}=z_{n-1,m}-z_{n,m}-z_{n-1,m-1}+z_{n,m-1}\\ & \quad -d_{n,m}(z_{0,0}-z_{N,0}-z_{0,M}+z_{N,M}) \nonumber \\&\quad f_{n,m}^{(4)}=z_{n,m-1}-d_{n,m}z_{0,0} \end{aligned}$$
(14)

Using the functions \(\phi _{\Delta (k,l)}:[0,1]^{2} \rightarrow [0,1]\) defined by

$$\begin{aligned} \phi _{N,M}(x,y)&=xy \nonumber \\ \phi _{0,0}(x,y)&=(1-x)(1-y) \nonumber \\ \phi _{0,M}(x,y)&=(1-x)y \nonumber \\ \phi _{N,0}(x,y)&=x(1-y), \end{aligned}$$
(15)

for \((k,l) \in \{n-1,n\} \times \{m-1,m\},\) the function \(F_{n,m}\) can be rewritten as follows:

$$\begin{aligned}&F_{n,m}(x,y,z)= dz\nonumber \\&\quad + \sum _{(k,l)\in \{n-1,n\} \times \{m-1,m\}} (z_{k,l}-dz_{\Delta (k,l)})\phi _{\Delta (k,l)}(x,y) \end{aligned}$$
(16)

Now define the IFS as:

$$\begin{aligned} W_{n,m}(x,y,z)=(L_{n}(x), K_{m}(y), F_{n,m}(x,y,z)) \end{aligned}$$
(17)

for \(n=1,2,...,N\) and \(m=1,2,...,M.\)

Note that when \(\sigma (n)=1,\)

$$\begin{aligned} F_{n,m}(1,y,z)&= dz+(z_{n,m}-dz_{N, M\sigma (m)})\nonumber \\&\quad \phi _{N,M\sigma (m)}(1,y)\nonumber \\&\quad +(z_{n,m-1}-dz_{N, M\sigma (m-1)})\nonumber \\&\quad \phi _{N,M\sigma (m-1)}(1,y) \,\,\,\text {and} \end{aligned}$$
(18)
$$\begin{aligned} F_{n+1,m}(1,y,z)&= dz+(z_{n,m}-dz_{N, M\sigma (m)})\nonumber \\&\quad \phi _{N,M\sigma (m)}(1,y) \nonumber \\&\quad +(z_{n,m-1}-dz_{N, M\sigma (m-1)})\nonumber \\&\quad \phi _{N,M\sigma (m-1)}(1,y) \end{aligned}$$
(19)

for \(n=1,2,...,N-1.\) Similar conditions are obtained for \(\sigma (n)=0, \,\, \sigma (m)=1\) and \(\sigma (m)=0.\)

The hyperbolicity of the IFS can be established as explained in proposition 2.1 in [10], for any values of n and m.

Theorem 1

Consider two positive integers N, M greater than 1. Let (17) be an IFS for the data set \(\{(x_{n}, y_{m}, z_{n,m}):n=0,1,...,N,\,\, m=0,1,...,M\}\) with the vertical scaling factor \(-1<d_{n,m}<1\) for \(n=1,2,...,N\) and \(m=1,2,...,M.\) Then, there exists a unique, continuous function \(f:[0,1]\times [0,1] \rightarrow R,\) interpolating the data set such that the graph of f coincides with the attractor of the IFS.

Proof

Consider the set \({\mathcal {F}}\) of all continuous functions \(f: [0,1]^{2} \rightarrow R\) such that f satisfies the following endpoint conditions,

$$\begin{aligned}&f(x_{N\sigma (k)},y_{M\sigma (l)}) \\&\quad =z_{N\sigma (k),M\sigma (l)}, \,\, (k,l)\in \{N-1,N\} \times \{M-1,M\}. \end{aligned}$$

Then, \({\mathcal {F}},\) with the supremum metric is a complete metric space.

Now, consider the operator \(T:{\mathcal {F}} \rightarrow {\mathcal {F}}\) such that

$$\begin{aligned}&(Tf)(x,y)= F_{n,m}(L_{n}^{-1}(x), K_{m}^{-1}(y), f(L_{n}^{-1}(x),\nonumber \\&\quad K_{m}^{-1}(y))), \,\, n=1,2,...,N, m=1,2,...,M. \end{aligned}$$
(20)

Clearly, T satisfies endpoint conditions on the complete metric space \({\mathcal {F}}.\) This can be established using the endpoint conditions on the functions \(L_{n},\) \(K_{m}\) and \(F_{n,m}.\) The operator T is continuous on the interior of each subrectangle \(I_{n}\times J_{m}, n=1,2,...,N, m=1,2,...,M.\) It remains to prove the continuity of T at the common boundaries of the two subrectangles, say, \(I_{n}\times J_{m}\) and \(I_{n+1}\times J_{m}.\)

Consider the common side \(x=x_{n}, \,\, y_{m-1}\le y \le y_{m}\) of the two subrectangles \(I_{n}\times J_{m}\) and \(I_{n+1}\times J_{m}.\) Along this common side,

$$\begin{aligned} (Tf)(x_{n},y)&= F_{n,m}(L_{n}^{-1}(x_{n}), K_{m}^{-1}(y), f(L_{n}^{-1}(x_{n}), K_{m}^{-1}(y))) \\&= F_{n,m}(x_{N}, K_{m}^{-1}(y), f(x_{N}, K_{m}^{-1}(y))) \end{aligned}$$

By the condition (18), this can be written as:

$$\begin{aligned}&F_{n,m}(x_{N}, K_{m}^{-1}(y), f(x_{N}, K_{m}^{-1}(y))) = F_{n+1,m}(x_{N}, K_{m}^{-1}(y),\\&\quad f(x_{N}, K_{m}^{-1}(y))), \,\, \text {if} \,\, n \,\, \text {is odd} \end{aligned}$$

When n is even,

$$\begin{aligned} (Tf)(x_{n},y)&= F_{n,m}(L_{n}^{-1}(x_{n}), K_{m}^{-1}(y), f(L_{n}^{-1}(x_{n}), K_{m}^{-1}(y))) \\&= F_{n,m}(x_{0}, K_{m}^{-1}(y), f(x_{0}, K_{m}^{-1}(y))) \end{aligned}$$

By the condition (18), this can be written as:

$$\begin{aligned}&F_{n,m}(x_{0}, K_{m}^{-1}(y), f(x_{0}, K_{m}^{-1}(y))) \\&\quad = F_{n+1,m}(x_{0}, K_{m}^{-1}(y), f(x_{0}, K_{m}^{-1}(y))), \\ \end{aligned}$$

Similarly, we can prove the continuity of T at all the common sides of the subrectangles. Therefore, T is continuous on each \(I_{n}\times J_{m}.\)

The contractivity of T is proved in Lemma 1 [17] as follows:

$$\begin{aligned} \mid \mid T(f)-T(g) \mid \mid _{\infty }&= sup_{n,m} \mid \mid F_{n,m}\Big (L_{n}^{-1}(x), \\&\quad K_{m}^{-1}(y), f(L_{n}^{-1}(x), K_{m}^{-1}(y))\Big ) \\&\quad -F_{n,m}\Big (L_{n}^{-1}(x), K_{m}^{-1}(y), g(L_{n}^{-1}(x), K_{m}^{-1}(y))\Big )\mid \mid _{\infty } \\&= sup_{n,m} \mid \mid d_{n,m} \Big (f(L_{n}^{-1}(x), K_{m}^{-1}(y))\\&\quad -g(L_{n}^{-1}(x), K_{m}^{-1}(y))\Big )\mid \mid _{\infty } \\&\le S \mid \mid f-g\mid \mid _{\infty } \end{aligned}$$

where \(S=sup_{n,m} \mid d_{n,m}\mid .\) Since T is a contraction, there exists a unique function f in \({\mathcal {F}}\) such that f is the fixed point of T and graph of f is the required fractal interpolation surface. Since f is the unique, fixed point of T,  it satisfies the recursive relation

$$\begin{aligned}&f(x,y) = d_{n,m}f(L_{n}^{-1}(x), K_{m}^{-1}(y)) \nonumber \\&\qquad \qquad \qquad + q_{n,m}(L_{n}^{-1}(x), K_{m}^{-1}(y)), \nonumber \\&\qquad \qquad \qquad \text {where} \,\, (x,y) \in I_{n}\times J_{m}, n=1,2,...,N,\nonumber \\&\qquad \qquad \qquad m=1,2,...,M. \end{aligned}$$
(21)

\(\hfill\square\)

The formula used for the vertical scaling factor is as provided below:

$$\begin{aligned} d_{n,m} = \frac{z_{n_{p},{m_{q}}}-0.25(z_{n,m}+z_{n-1,m}+z_{n,m-1}+z_{n-1,m-1})}{z_{\frac{N}{2}, \frac{M}{2}}-0.25(z_{N,M}+z_{N,0}+z_{0,M}+z_{0,0})} \end{aligned}$$
(22)

where \(z_{n_{p},m_{q}} =f(\frac{x_{n}+x_{n-1}}{2},\frac{y_{m}+y_{m-1}}{2}), \,\, z_{\frac{N}{2},\frac{M}{2}} = f(\frac{x_{N}+x_{0}}{2},\frac{y_{M}+y_{0}}{2})\)

This formula has been derived in Sect. 5 as a special case of the general expression for the vertical scaling factor.

3 Numerical integration over rectangular region

Let \(D_{0}\) denote the double integral of a bivariate fractal interpolation function f over the rectangular region \([0,1]\times [0,1].\) Then,

$$\begin{aligned}&D_{0}= \int \limits _{I}\int \limits _{J} f(x,y) dxdy \\&\quad = \sum _{n=1}^{N} \sum _{m=1}^{M} \int \limits _{I_{n}} \int \limits _{J_{m}} f(x,y) dxdy \\&\quad = \sum _{n=1}^{N} \sum _{m=1}^{M} \int \limits _{I_{n}} \int \limits _{J_{m}} (Tf)(x,y) dxdy \\&\quad = \sum _{n=1}^{N} \sum _{m=1}^{M} \int \limits _{I_{n}} \int \limits _{J_{m}} F_{n,m}\Big (L_{n}^{-1}(x),K_{m}^{-1}(y),\\&\qquad \quad f\big (L_{n}^{-1}(x),K_{m}^{-1}(y)\big )\Big ) dxdy \end{aligned}$$

Since \(F_{n,m}\) is defined differently based on the indices n, m are even or odd,

$$\begin{aligned} D_{0}&= \sum _{n=1}^{N/2} \sum _{m=1}^{M/2} \int \limits _{I_{2n-1}} \int \limits _{J_{2m-1}} F_{2n-1,2m-1}\Big (L_{2n-1}^{-1}(x),K_{2m-1}^{-1}(y),\nonumber \\&\quad f\big (L_{2n-1}^{-1}(x),K_{2m-1}^{-1}(y)\big )\Big ) dxdy \nonumber \\&\quad + \sum _{n=1}^{N/2} \sum _{m=1}^{M/2} \int \limits _{I_{2n}} \int \limits _{J_{2m}} F_{2n,2m}\Big (L_{2n}^{-1}(x),K_{2m}^{-1}(y),\nonumber \\&\quad f\big (L_{2n}^{-1}(x),K_{2m}^{-1}(y)\big )\Big ) dxdy \nonumber \\&\quad + \sum _{n=1}^{N/2} \sum _{m=1}^{M/2} \int \limits _{I_{2n-1}} \int \limits _{J_{2m}} F_{2n-1,2m}\Big (L_{2n-1}^{-1}(x),K_{2m}^{-1}(y),\nonumber \\&\quad f\big (L_{2n-1}^{-1}(x),K_{2m}^{-1}(y)\big )\Big ) dxdy \nonumber \\&\quad + \sum _{n=1}^{N/2} \sum _{m=1}^{M/2} \int \limits _{I_{2n}} \int \limits _{J_{2m-1}} F_{2n,2m-1}\Big (L_{2n}^{-1}(x),K_{2m-1}^{-1}(y),\nonumber \\&\quad f\big (L_{2n}^{-1}(x),K_{2m-1}^{-1}(y)\big )\Big ) dxdy \end{aligned}$$
(23)

Using \(u_{1}=L_{2n-1}^{-1}(x), \,\,\, v_{1}=L_{2\,m-1}^{-1}(y), \,\,\, u_{2}=L_{2n}^{-1}(x), \,\,\, v_{2}=L_{2\,m}^{-1}(y),\)

\(D_{0}\) can be written as

$$\begin{aligned} D_{0}&= \sum _{n=1}^{N/2} \sum _{m=1}^{M/2} \int \limits _{I} \int \limits _{J} F_{2n-1,2m-1}\big (u_{1},v_{1},\nonumber \\&\quad f(u_{1},v_{1})\big )r_{2n-1}q_{2m-1} du_{1}dv_{1} \nonumber \\&\quad + \sum _{n=1}^{N/2} \sum _{m=1}^{M/2} \int \limits _{I} \int \limits _{J} F_{2n,2m}\big (u_{2},v_{2},\nonumber \\&\quad f(u_{2},v_{2})\big )r_{2n}q_{2m} du_{2}dv_{2} \nonumber \\&\quad + \sum _{n=1}^{N/2} \sum _{m=1}^{M/2} \int \limits _{I} \int \limits _{J} F_{2n-1,2m}\big (u_{1},v_{2},\nonumber \\&\quad f(u_{1},v_{2})\big )r_{2n-1}q_{2m} du_{1}dv_{2} \nonumber \\&\quad + \sum _{n=1}^{N/2} \sum _{m=1}^{M/2} \int \limits _{I} \int \limits _{J} F_{2n,2m-1}\big (u_{2},v_{1},\nonumber \\&\quad f(u_{2},v_{1})\big )r_{2n}q_{2m-1} du_{2}dv_{1} \nonumber \\ \end{aligned}$$
(24)
$$\begin{aligned}&= \sum _{n=1}^{N/2} \sum _{m=1}^{M/2} \int \limits _{I} \int \limits _{J} \Big [a_{2n-1,2m-1}^{(1)}u_{1} +b_{2n-1,2m-1}^{(1)}v_{1}\\&\quad +c_{2n-1,2m-1}^{(1)}u_{1}v_{1}+f_{2n-1,2m-1}^{(1)} \\&\quad +d_{2n-1,2m-1}f(u_{1},v_{1})\Big ]r_{2n-1}q_{2m-1}du_{1}dv_{1}\\&\quad + \sum _{n=1}^{N/2} \sum _{m=1}^{M/2} \int \limits _{I} \int \limits _{J} \Big [a_{2n,2m}^{(2)}u_{2} +b_{2n,2m}^{(2)}v_{2}\\&\quad +c_{2n,2m}^{(2)}u_{2}v_{2}+f_{2n,2m}^{(2)}\\&\quad +d_{2n,2m}f(u_{2},v_{2})\Big ]r_{2n}q_{2m}du_{2}dv_{2}\\&\quad + \sum _{n=1}^{N/2} \sum _{m=1}^{M/2} \int \limits _{I} \int \limits _{J}\Big [a_{2n-1,2m}^{(3)}u_{1} +b_{2n-1,2m}^{(3)}v_{2}\\&\quad +c_{2n-1,2m}^{(3)}u_{1}v_{2}+f_{2n-1,2m}^{(3)}\\&\quad +d_{2n-1,2m}f(u_{1},v_{2})\Big ]r_{2n-1}q_{2m}du_{1}dv_{2}\\&\quad + \sum _{n=1}^{N/2} \sum _{m=1}^{M/2} \int \limits _{I} \int \limits _{J} \Big [a_{2n,2m-1}^{(4)}u_{2} +b_{2n,2m-1}^{(4)}v_{1}\\&\quad +c_{2n,2m-1}^{(4)}u_{2}v_{1}+f_{2n,2m-1}^{(4)}\\&\quad +d_{2n,2m-1}f(u_{2},v_{1})\Big ]r_{2n}q_{2m-1}du_{2}dv_{1}\\ \end{aligned}$$

which implies

$$\begin{aligned} D_{0}&= P_{01}+E_{01} \int \limits _{I}\int \limits _{J} f(u_{1},v_{1}) du_{1}dv_{1} \nonumber \\&+ P_{02}+E_{02} \int \limits _{I}\int \limits _{J} f(u_{2},v_{2}) du_{2}dv_{2} \nonumber \\&+ P_{03}+E_{03} \int \limits _{I}\int \limits _{J} f(u_{1},v_{2}) du_{1}dv_{2} \nonumber \\&+ P_{04}+E_{04} \int \limits _{I}\int \limits _{J} f(u_{2},v_{1}) du_{2}dv_{1} \end{aligned}$$
(25)

where

$$\begin{aligned} P_{01}&= \sum _{n=1}^{N/2} \sum _{m=1}^{M/2} \int \limits _{I} \int \limits _{J} \Big [a_{2n-1,2m-1} ^{(1)}u_{1}\\&\quad +b_{2n-1,2m-1}^{(1)}v_{1}+c_{2n-1,2m-1}^{(1)}u_{1}v_{1}\\&\quad +f_{2n-1,2m-1}^{(1)}\Big ]r_{2n-1}q_{2m-1}du_{1}dv_{1}\\ P_{02}&= \sum _{n=1}^{N/2} \sum _{m=1}^{M/2} \int \limits _{I} \int \limits _{J} \Big [a_{2n,2m}^{(2)} u_{2}+b_{2n,2m}^{(2)}v_{2}\\&\quad +c_{2n,2m}^{(2)}u_{2}v_{2}+f_{2n,2m}^{(2)}\Big ]r_{2n}q_{2m}du_{2}dv_{2}\\ P_{03}&= -\sum _{n=1}^{N/2} \sum _{m=1}^{M/2} \int \limits _{I} \int \limits _{J} \Big [a_{2n-1,2m} ^{(3)}u_{1}\\&\quad +b_{2n-1,2m}^{(3)}v_{2}+c_{2n-1,2m}^{(3)}u_{1}v_{2}\\&\quad +f_{2n-1,2m}^{(3)}\Big ]r_{2n-1}q_{2m}du_{1}dv_{2}\\ P_{04}&=- \sum _{n=1}^{N/2} \sum _{m=1}^{M/2} \int \limits _{I} \int \limits _{J} \Big [a_{2n,2m-1} ^{(4)}u_{2}+b_{2n,2m-1}^{(4)}v_{1}\\&\quad +c_{2n,2m-1}^{(4)}u_{2}v_{1}\\&\quad +f_{2n,2m-1}^{(4)}\Big ]r_{2n}q_{2m-1}du_{2}dv_{1}\\ \end{aligned}$$

and

$$\begin{aligned} E_{01}&= \sum _{n=1}^{N/2} \sum _{m=1}^{M/2} d_{2n-1,2m-1} r_{2n-1}q_{2m-1} \nonumber \\ E_{02}&= \sum _{n=1}^{N/2} \sum _{m=1}^{M/2} d_{2n,2m} r_{2n}q_{2m} \nonumber \\ E_{03}&= - \sum _{n=1}^{N/2} \sum _{m=1}^{M/2} d_{2n-1,2m} r_{2n-1}q_{2m} \nonumber \\ E_{04}&= -\sum _{n=1}^{N/2} \sum _{m=1}^{M/2} d_{2n,2m-1} r_{2n}q_{2m-1} \end{aligned}$$
(26)

Therefore, (25) can be written as

$$\begin{aligned}&D_{0}= (P_{01}+E_{01}D_{0}) + (P_{02}\nonumber \\&\quad +E_{02}D_{0}) + (P_{03}+E_{03}D_{0}) +(P_{04}+E_{04}D_{0}) \end{aligned}$$
(27)

which implies,

$$\begin{aligned} D_{0}&= \frac{P_{01}+P_{02}+P_{03}+P_{04}}{1-(E_{01}+E_{02}+E_{03}+E_{04})} \end{aligned}$$
(28)

4 Bivariate fractal interpolation functions and bilinear interpolation functions on the rectangular domain

Theorem 2

Consider the bivariate fractal interpolation function f constructed for the data set \(\{(x_{n},y_{m},z_{n,m}):n=0,1,...,N, m=0,1,...,M\}\) on the rectangular domain \(I\times J.\) Let h be the corresponding bilinear interpolation function defined on \(I_{n}\times J_{m}, n=1,2,...,N, \,\, m=1,2,...,M.\) Define r as the bilinear interpolation function defined on the rectangle \(I\times J.\) Then, \(f,\,\,h,\,\,r\) are related by the recursive relation

$$\begin{aligned} f(x,y)&=h(x,y)+d_{n,m}(f-r)o(L_{n}^{-1}(x),K_{m}^{-1}(y)) \end{aligned}$$
(29)

Proof

Let h be a bilinear interpolation function for the data set

$$\begin{aligned}{} & {} \big \{(x_{n-1},y_{m-1},z_{n-1,m-1}), (x_{n},y_{m-1},z_{n,m-1}), \\{} & {} \quad (x_{n-1},y_{m},z_{n-1,m}),(x_{n},y_{m},z_{n,m})\big \} \end{aligned}$$

of the form

$$\begin{aligned} h(x,y)&= p_{n,m}x+u_{n,m}y+v_{n,m}xy+t_{n,m}. \end{aligned}$$
(30)

Denote r as another bilinear interpolation function for the data set

$$\begin{aligned}{} & {} \big \{(x_{0},y_{0},z_{0,0}),(x_{N}, y_{0}, z_{N,0}), \\{} & {} \quad (x_{0}, y_{M}, z_{0,M}),(x_{N}, y_{M}, z_{N,M})\big \} \end{aligned}$$

of the form

$$\begin{aligned} r(x,y)&= P_{n,m}x+U_{n,m}y+V_{n,m}xy+T_{n,m}. \end{aligned}$$
(31)

Since the IFS (17) satisfies the condition (9), there are four different expressions for the coefficients involved in \(h, \,\, r\) in the following four cases:

  • case 1: n, m are odd

  • case 2: n, m are even

  • case 3: n is odd, m is even

  • case 4: n is even, m is odd

The proof of the theorem is similar in all the four cases. Hence the proof for the first case is only provided here. When n, m are odd integers the coefficients in \(h, \,\,r\) are obtained as follows:

$$\begin{aligned} p_{n,m}&= Nm(z_{n,m-1}-z_{n-1,m-1}) \\&\quad + N(m-1)(z_{n-1,m},z_{n,m}) \\ u_{n,m}&=Mn(z_{n-1,m}\\&\quad -z_{n-1,m-1})+M(n-1)(z_{n,m-1}-z_{n,m}) \\ v_{n,m}&= NM(z_{n,m}-z_{n-1,m}\\&\quad -z_{n,m-1}+z_{n-1,m-1}) \\ t_{n,m}&= mn(z_{n,m}-z_{n-1,m}\\&\quad -z_{n,m-1}+z_{n-1,m-1}) + n(z_{n-1,m}-z_{n,m}) \\&\quad + m(z_{n,m-1}-z_{n,m}) +z_{n,m} \\ P_{n,m}&= z_{N,0}-z_{0,0}\\ U_{n,m}&= z_{0,M}-z_{0,0} \\ V_{n,m}&= z_{N,M}-z_{0,M}-z_{N,0}+z_{0,0}\\ T_{n,m}&=z_{0,0}. \end{aligned}$$

In order to prove (29), one should establish the following relation:

$$\begin{aligned} h(L_{n}(x),K_{m}(y)) - d_{n,m} r(x,y)&= q_{n,m} (x,y) \end{aligned}$$
(32)

(32) is to be proved by equating the coefficients of \(x, \,\, y, \,\, xy\) and constant term in the LHS and RHS.

Note that when \(n, \,\, m\) are odd integers, the coefficients in \(q_{n,m}\) are as given in Eq. (11).

The coefficient of x in the LHS of (32) is:

$$\begin{aligned}&(p_{n,m}+v_{n,m}w_{m})r_{n} -d_{n,m}P_{n,m} \\&\quad = z_{n,m-1}-z_{n-1,m-1}-d_{n,m}(z_{N,0}-z_{0,0}) \end{aligned}$$

which is nothing but the coefficient of x in the RHS.

The coefficient of y in the LHS of (32):

$$\begin{aligned}&(u_{n,m}+v_{n,m}s_{n})q_{m} -d_{n,m}U_{n,m} \\&\quad = z_{n-1,m}-z_{n-1,m-1}-d_{n,m}(z_{0,M}-z_{0,0}) \end{aligned}$$

which is nothing but the coefficient of y in the RHS.

The coefficient of xy in the LHS of (32):

$$\begin{aligned}&v_{n,m}r_{n}q_{m}-d_{n,m}V_{n,m} = z_{n,m}\\&\quad -z_{n,m-1}-z_{n-1,m}+z_{n-1,m-1} \\&\quad -d_{n,m}(z_{N,M}-z_{N,0}-z_{0,M}+z_{0,0}) \end{aligned}$$

which is nothing but the coefficient of xy in the RHS.

The constant term in the LHS of (32):

$$\begin{aligned}&p_{n,m}s_{n}+u_{n,m}w_{m}+v_{n,m}s_{n}w_{m}\\&\quad +t_{n,m}-d_{n,m}T_{n,m} = z_{n-1,m-1}-d_{n,m}z_{0,0} \end{aligned}$$

which is nothing but the constant term in the RHS.

Therefore,

$$\begin{aligned} h(L_{n}(x),K_{m}(y)) - d_{n,m} r(x,y)&= q_{n,m} (x,y) \end{aligned}$$

Now, substituting (32) in (21), it is obtained that

$$\begin{aligned} f(x,y)&=h(x,y)+d_{n,m}(f-r)o(L_{n}^{-1}(x),K_{m}^{-1}(y)) \end{aligned}$$

Hence the proof. \(\hfill\square\)

5 Selection of vertical scaling factor

Consider the given data set \(\{(x_{n},y_{m},z_{n,m}):n=0,1,...,N, m=0,1,...,M\}\) where \(x_{n}=x_{0}+nh^{'}, \,\, y_{m}=y_{0}+mk^{'}.\) Divide each of the subintervals \(I_{n}, \,\, J_{m}\) into \(p, \,\, q\) number of equal parts respectively and consider the intermediate points \(x_{i}^{'} \in I_{n}, \,\, y_{j}^{'} \in j_{m}, \,\, i=1,2,...,p-1, \,\, j=1,2,...,q-1\) where \(p \ge 2, \,\, q\ge 2.\) The expressions for \(x_{i}^{'}, \,\, y_{j}^{'}\) can be written as

$$\begin{aligned} x_{i}^{'}&= \frac{(p-i)x_{n-1}+ix_{n}}{p}\\ y_{j}^{'}&= \frac{(q-j)y_{m-1}+jy_{m}}{q} \end{aligned}$$

If f is the fractal interpolation function to the data set, then

$$\begin{aligned}&z_{i,j}^{'} = f(x_{i}^{'}, y_{j}^{'}) \\&\quad = d_{n,m} f( L_{n}^{-1}(x_{i}^{'}), K_{m}^{-1}(y_{j}^{'})) \\&\qquad + q_{n,m} (L_{n}^{-1}(x_{i}^{'}), K_{m}^{-1}(y_{j}^{'}))\\&\quad = d_{n,m} f( L_{n}^{-1}(x_{i}^{'}), K_{m}^{-1}(y_{j}^{'})) + h(x_{i}^{'}, y_{j}^{'}) \\&\quad - d_{n,m} r(L_{n}^{-1}(x_{i}^{'}), K_{m}^{-1}(y_{j}^{'})) \end{aligned}$$

Approximating f by h, 

$$\begin{aligned} z_{i,j}^{'}&= h(x_{i}^{'}, y_{j}^{'}) - d_{n,m} (h-r)o(L_{n}^{-1}(x_{i}^{'}), K_{m}^{-1}(y_{j}^{'})) \end{aligned}$$

Now, considering the optimization problem,

$$\begin{aligned}&min (E(d_{n,m})) = \sum _{i=1}^{p-1} \sum _{j=1}^{q-1} \big [z_{i,j}^{'}-h(x_{i}^{'}, y_{j}^{'})\\&\quad -d_{n,m}(h-r)o(L_{n}^{-1}(x_{i}^{'}), K_{m}^{-1}(y_{j}^{'}))\big ]^{2} \end{aligned}$$

with \(-1< d_{n,m} <1.\)

For solving the problem, let \(\frac{d(E(d_{n,m}))}{d(d_{n,m})}=0,\) which implies

$$\begin{aligned}&-2 \sum _{i=1}^{p-1} \sum _{j=1}^{q-1} z_{i,j}^{'} G_{i,j}+ 2 \sum _{i=1}^{p-1} \sum _{j=1}^{q-1} h(x_{i}^{'}, y_{j}^{'}) G_{i,j} \\&\quad + 2 d_{n,m} \sum _{i=1}^{p-1} \sum _{j=1}^{q-1} G_{i,j}^{2} =0 \end{aligned}$$

where \(G_{i,j} = (h-r)o(L_{n}^{-1}(x_{i}^{'}), K_{m}^{-1}(y_{j}^{'})).\)

This implies

$$\begin{aligned} d_{n,m}&= \frac{\sum _{i=1}^{p-1} \sum _{j=1}^{q-1} \big [z_{i,j}^{'}-h(x_{i}^{'}, y_{j}^{'})\big ] G_{i,j}}{\sum _{i=1}^{p-1} \sum _{j=1}^{q-1} G_{i,j}^{2}} \end{aligned}$$
(33)

Hence, the error gets minimized for this particular value of \(d_{n,m}.\)

When \(p=2,\,\,q=2,\) the formula for the vertical scaling factor reduces to (22).

In order to find the upper bound of the vertical scaling factor \(d_{n,m}\) apply Cauchy-Schwartz inequality to (33). Then,

$$\begin{aligned} \mid d_{n,m} \mid&\le \frac{[\sum _{i=1}^{p-1} \sum _{j=1}^{q-1} \mid z_{i,j}^{'}-h(x_{i}^{'}, y_{j}^{'}) \mid ^{2}]^{1/2} \big [\sum _{i=1}^{p-1} \sum _{j=1}^{q-1} \mid G_{i,j}\mid ^{2}\big ]^{1/2}}{\sum _{i=1}^{p-1} \sum _{j=1}^{q-1} \mid G_{i,j}\mid ^{2}} \end{aligned}$$

Put \(T_{h} = [\sum _{i=1}^{p-1} \sum _{j=1}^{q-1} \mid G_{i,j}\mid ^{2}]^{1/2},\) which implies

$$\begin{aligned} \mid d_{n,m} \mid&\le \frac{\big [\sum _{i=1}^{p-1} \sum _{j=1}^{q-1} \mid z_{i,j}^{'}-h(x_{i}^{'}, y_{j}^{'}) \mid ^{2}\big ]^{1/2}}{T_{h}} \end{aligned}$$
(34)

It is to be noted that the vertical scaling factor \(d_{n,m}\) varies in the interval \((-1,1).\) However, to find the upper bound of \(d_{n,m}\) particularly for a given data set, (34) is used.

6 Error analysis

Consider a continuous function g defined on the rectangular region \(I\times J.\) Then, define \(\Vert g\Vert _{\infty }\) as

$$\begin{aligned} \Vert g\Vert _{\infty } = max\{\mid g(x,y) \mid : x \in I, \,\, y \in J\}. \end{aligned}$$

Let \(w_{g}(\delta )\) represents the modulus of continuity of g,  then

$$\begin{aligned}w_{g}(\delta )&= sup \{\mid g(x,y)-g(x^{'},y^{'})\mid : d((x,y),(x^{'},y^{'}))\le \delta , \,\, \\ &\quad (x,y), \,\, (x^{'},y^{'}) \in I\times J\} \end{aligned}$$

where the metric considered here is given by

$$\begin{aligned} d((x,y),(x^{'},y^{'}))=\mid x-x^{'}\mid +\mid y-y^{'} \mid . \end{aligned}$$

Lemma 3

If g is a continuous function providing the data \(\{(x_{n},y_{m},z_{n,m}):n=0,1,...,N, m=0,1,...,M\}\) with \(h^{'}=x_{n}-x_{n-1}, \,\, k^{'}=y_{m}-y_{m-1}.\) If f is the corresponding bivariate fractal interpolation function with scale vector \(\overline{d_{n,m}},\) then

$$\begin{aligned} \Vert g-f\Vert _{\infty } \le 4 w_{g}(h^{'}+k^{'})+ \frac{\mid d_{n,m} \mid _{\infty } \Vert h-r\Vert _{\infty }}{1-\mid d_{n,m} \mid _{\infty }}. \end{aligned}$$

Proof

Since the IFS satisfies the condition (9), the lemma has to be proved by considering the following cases:

  • case 1: n, m are odd

  • case 2: n, m are even

  • case 3: n is odd, m is even

  • case 4: n is even, m is odd

Since the proof is similar for all the four cases, only the first case is proved here.

Let n, m be odd integers.

Consider the bilinear interpolation function h defined as in (30) through the points

$$\begin{aligned}{} & {} \big \{(x_{n-1},y_{m-1},z_{n-1,m-1}), (x_{n},y_{m-1},z_{n,m-1}), \\{} & {} \quad (x_{n-1},y_{m},z_{n-1,m}),(x_{n},y_{m},z_{n,m})\big \}. \end{aligned}$$

By rearranging, the function h can also be written as

$$\begin{aligned} h(x,y)&= \frac{(x-x_{n-1})(y-y_{m-1})z_{n,m}}{(x_{n}-x_{n-1})(y_{m}-y_{m-1})}\\ & \quad + \frac{(x-x_{n})(y-y_{m-1})z_{n-1,m}}{(x_{n-1}-x_{n})(y_{m}-y_{m-1})} \\&\quad + \frac{(x-x_{n-1})(y-y_{m})z_{n,m-1}}{(x_{n}-x_{n-1})(y_{m-1}-y_{m})}\\ & \quad + \frac{(x-x_{n})(y-y_{m})z_{n-1,m-1}}{(x_{n-1}-x_{n})(y_{m-1}-y_{m})} \end{aligned}$$

Now, consider \(\Vert g-h\Vert\)

Let \(w_{g}(h^{'}+k^{'}) = \{\mid g(x,y)-g(x^{'},y^{'})\mid : d((x,y),(x^{'},y^{'}))\le h^{'}+k^{'}\}\)

be the modulus of continuity of g.

Then, \(\mid g(x,y)-h(x,y) \mid\)

$$\begin{aligned} &= \Bigg \vert \frac{(g(x,y)-z_{n,m}) (x-x_{n-1})(y-y_{m-1})}{(x_{n}-x_{n-1})(y_{m}-y_{m-1})} \\{} & {} \quad - \frac{(g(x,y)-z_{n-1,m})(x-x_{n})(y-y_{m-1})}{(x_{n}-x_{n-1})(y_{m}-y_{m-1})} \\{} & {} \quad -\frac{(g(x,y)-z_{n,m-1}) (x-x_{n-1})(y-y_{m})}{(x_{n}-x_{n-1})(y_{m}-y_{m-1})} \\{} & {} \quad + \frac{(g(x,y)-z_{n-1,m-1})(x-x_{n})(y-y_{m})}{(x_{n}-x_{n-1})(y_{m}-y_{m-1})} \Bigg \vert \\{} & \le \Bigg \vert \frac{(g(x,y)-z_{n,m}) (x-x_{n-1})(y-y_{m-1})}{(x_{n}-x_{n-1})(y_{m}-y_{m-1})} \Bigg \vert \\{} & {} \quad + \Bigg \vert \frac{(g(x,y)-z_{n-1,m}) (x-x_{n})(y-y_{m-1})}{(x_{n}-x_{n-1})(y_{m}-y_{m-1})} \Bigg \vert \\{} & {} \quad + \Bigg \vert \frac{(g(x,y)-z_{n,m-1}) (x-x_{n-1})(y-y_{m})}{(x_{n}-x_{n-1})(y_{m}-y_{m-1})} \Bigg \vert \\{} & {} \quad + \Bigg \vert \frac{(g(x,y)-z_{n-1,m-1}) (x-x_{n})(y-y_{m})}{(x_{n}-x_{n-1})(y_{m}-y_{m-1})} \Bigg \vert \\{} & \le \frac{\Bigg \vert (g(x,y)-z_{n,m}) \Bigg \vert \Bigg \vert (x-x_{n-1})\Bigg \vert \Bigg \vert (y-y_{m-1})\Bigg \vert }{\Bigg \vert (x_{n}-x_{n-1}) \Bigg \vert \Bigg \vert (y_{m}-y_{m-1})\Bigg \vert } \\{} & {} \quad + \frac{\Bigg \vert (g(x,y)-z_{n-1,m}) \Bigg \vert \Bigg \vert (x-x_{n})\Bigg \vert \Bigg \vert (y-y_{m-1})\Bigg \vert }{\Bigg \vert (x_{n}-x_{n-1})\Bigg \vert \Bigg \vert (y_{m}-y_{m-1})\Bigg \vert } \\{} & {} \quad + \frac{\Bigg \vert (g(x,y)-z_{n,m-1})\Bigg \vert \Bigg \vert (x-x_{n-1})\Bigg \vert \Bigg \vert (y-y_{m})\Bigg \vert }{\Bigg \vert (x_{n}-x_{n-1})\Bigg \vert \Bigg \vert (y_{m}-y_{m-1})\Bigg \vert } \\{} & {} \quad + \frac{\Bigg \vert (g(x,y)-z_{n-1,m-1}) \Bigg \vert \Bigg \vert (x-x_{n})\Bigg \vert \Bigg \vert (y-y_{m})\Bigg \vert }{\Bigg \vert (x_{n}-x_{n-1})\Bigg \vert \Bigg \vert (y_{m}-y_{m-1})\Bigg \vert } \\{} & \le sup_{\{x \in I, y \in J\}}\left\{\frac{\bigg \vert (g(x,y)-z_{n,m}) \bigg \vert \bigg \vert (x-x_{n-1})\Bigg \vert \Bigg \vert (y-y_{m-1})\Bigg \vert }{\Bigg \vert (x_{n}-x_{n-1}) \Bigg \vert \Bigg \vert (y_{m}-y_{m-1})\Bigg \vert }\right\}\\{} & {} \quad + sup_{\{x \in I, y \in J\}} \left\{\frac{\Bigg \vert (g(x,y)-z_{n-1,m}) \Bigg \vert \Bigg \vert (x-x_{n})\Bigg \vert \Bigg \vert (y-y_{m-1})\Bigg \vert }{\Bigg \vert (x_{n} -x_{n-1})\Bigg \vert \Bigg \vert (y_{m}-y_{m-1})\Bigg \vert }\right\} \\{} & {} \quad + sup_{\{x \in I, y \in J\}} \left\{\frac{\Bigg \vert (g(x,y)-z_{n,m-1})\Bigg \vert \Bigg \vert (x-x_{n-1}) \Bigg \vert \Bigg \vert (y-y_{m})\Bigg \vert }{\Bigg \vert (x_{n}-x_{n-1})\Bigg \vert \Bigg \vert (y_{m}-y_{m-1}) \Bigg \vert }\right\} \\{} & {} \quad + sup_{\{x \in I, y \in J\}} \left\{\frac{\Bigg \vert (g(x,y)-z_{n-1,m-1}) \Bigg \vert \Bigg \vert (x-x_{n})\Bigg \vert \Bigg \vert (y-y_{m})\Bigg \vert }{\Bigg \vert (x_{n}-x_{n-1})\Bigg \vert \Bigg \vert (y_{m}-y_{m-1})\Bigg \vert }\right\} \end{aligned}$$

Since \(x, x_{n}, x_{n-1},\) are in \(I_{n}\) and \(y, y_{m}, y_{m-1}\) are in \(J_{m}\) and \(\mid I_{n} \mid = h^{'},\)

\(\mid J_{m} \mid = k^{'},\) implies \(\mid x-x_{n-1} \mid \le h^{'}, \,\, \mid x-x_{n} \mid \le h^{'}, \,\, \mid y-y_{m-1} \mid \le k^{'}, \,\,\)

and \(\mid y-y_{m} \mid \le k^{'}.\)

Therefore, the above inequality becomes

$$\begin{aligned}&= sup_{\{x \in I, y \in J\}} \{\mid g(x,y)-g(x_{n},y_{m}) \mid \}\\&+ sup_{\{x \in I, y \in J\}} \{\mid g(x,y)-g(x_{n-1},y_{m}) \mid \} \\&+ sup_{\{x \in I, y \in J\}} \{\mid g(x,y)-g(x_{n},y_{m-1}) \mid \} \\&+ sup_{\{x \in I, y \in J\}} \{\mid g(x,y)-g(x_{n-1},y_{m-1}) \mid \}\\&\le 4w_{g}(h^{'}+k^{'}) \end{aligned}$$

Thus, \(\Vert g-h\Vert \le 4 w_{g}(h^{'}+k^{'}).\)

Now, in order to find \(\Vert h-f\Vert ,\) consider the relation

$$\begin{aligned} f(x,y)&= h(x,y)+d_{n,m}(f-r)o(L_{n}^{-1}(x), K_{m}^{-1}(y)) \end{aligned}$$

Thus,

$$\begin{aligned} \Vert f-h\Vert _{\infty }&\le \mid d_{n,m} \mid _{\infty } \Vert f-r\Vert _{\infty } \\&\le \mid d_{n,m} \mid _{\infty } \big [\Vert f-h\Vert _{\infty } + \Vert h-r\Vert _{\infty }\big ]\\ \end{aligned}$$

i.e,

$$\begin{aligned} \Vert f-h \Vert _{\infty } \le \frac{\mid d_{n,m} \mid _{\infty } \Vert h-r\Vert _{\infty }}{1-\mid d_{n,m} \mid _{\infty }} \end{aligned}$$

Therefore, \(\Vert g-f\Vert \le 4 w_{g}(h^{'}+k^{'})+ \frac{\mid d_{n,m} \mid _{\infty } \Vert h-r\Vert _{\infty }}{1-\mid d_{n,m} \mid _{\infty }}.\) Hence the proof. \(\hfill\square\)

Theorem 4

If g is a continuous interpolation function to the data set \(\{(x_{n},y_{m},z_{n,m}):n=0,1,...,N, m=0,1,...,M\}\) and f is the corresponding bivariate fractal interpolation function, then the double integral value calculated by the fractal numerical integral method converges to the exact integral value as \(h^{'},k^{'} \rightarrow 0,\) where \(h^{'}=x_{n}-x_{n-1}, \,\, k^{'}=y_{m}-y_{m-1}\) and

$$\begin{aligned}&\Bigg \vert \int _{I}\int _{J} g - \int _{I}\int _{J} f \Bigg \vert\le \big \vert I \big \vert \big \vert J \big \vert w_{g}(h^{'}+k^{'})\\&\quad \Big [4+\frac{T_{H}\Vert h-r\Vert _{\infty }}{1-w_{g}(h^{'}+k^{'})T_{h}}\Big ] \end{aligned}$$

Proof

Let \(\Big \vert E_{h^{'}+k^{'}} \Big \vert = \Big \vert \int _{I}\int _{J} g - \int _{I}\int _{J} f\Big \vert \le \big \vert I \big \vert \big \vert J \big \vert \Vert g-f\Vert _{\infty }\)

By lemma 3

$$\begin{aligned} \Vert g-f\Vert _{\infty }&\le 4 w_{g}(h^{'}+k^{'})+ \frac{\mid d_{n,m} \mid _{\infty } \Vert h-r\Vert _{\infty }}{1-\mid d_{n,m} \mid _{\infty }} \end{aligned}$$

If g is a function interpolating the data, then,

$$\begin{aligned} \mid z_{i,j}^{'} - h(x_{i}^{'}, y_{j}^{'}) \mid&= \mid g(x_{i}^{'}, y_{j}^{'}) - h(x_{i}^{'}, y_{j}^{'}) \mid \\&\le sup_{i,j} \{\mid g(x_{i}^{'}, y_{j}^{'}) - h(x_{i}^{'}, y_{j}^{'}) \mid \} \\&= \Vert g-h\Vert _{\infty } \\&\le 4w_{g}(h^{'}+k^{'}) \end{aligned}$$

where \(x_{i}^{'}, \,\, y_{j}^{'}\) and \(z_{i,j}^{'}\) are as defined in Sect. 5. Substituting this in the expression for upper bound in (34),

$$\begin{aligned} \mid d_{n,m} \mid&\le \frac{\Big [\sum _{i=1}^{p-1} \sum _{j=1}^{q-1} [4w_{g}(h^{'}+k^{'})]^{2}\Big ]^{1/2}}{T_{h}} \\&= \frac{[[4w_{g}(h^{'}+k^{'})]^{2}]^{1/2} (p-1)^{1/2} (q-1)^{1/2}}{T_{h}} \\&\le w_{g}(h^{'}+k^{'}) T_{H}, \end{aligned}$$

where \(T_{H} = \frac{4(p-1)^{1/2} (q-1)^{1/2}}{T_{h}}.\)

Therefore,

$$\begin{aligned} \Vert g-f\Vert _{\infty }&\le 4 w_{g}(h^{'}+k^{'})+ \frac{w_{g}(h^{'}+k^{'}) T_{H} \Vert h-r\Vert _{\infty }}{1-w_{g}(h^{'}+k^{'}) T_{H}} \\&=w_{g}(h^{'}+k^{'}) \Big [4+\frac{T_{H} \Vert h-r\Vert _{\infty }}{1-w_{g}(h^{'}+k^{'}) T_{H}}\Big ] \end{aligned}$$

As \(h^{'},k^{'} \,\, \rightarrow 0, w_{g}(h^{'}+k^{'}) \rightarrow 0,\) which makes \(\Vert g-f\Vert \le w_{g}(h^{'}+k^{'}) \rightarrow 0.\) Hence, f converges to g uniformly.

Also,

$$\begin{aligned} \mid E_{h^{'}+k^{'}} \mid&\le \mid I \mid \mid J \mid w_{g}(h^{'}+k^{'})\Big [4+\frac{T_{H} \Vert h-r\Vert _{\infty }}{1-w_{g}(h^{'}+k^{'}) T_{H}}\Big ] \\&\rightarrow 0. \end{aligned}$$

Hence, \(\mid E_{h^{'}+k^{'}} \mid \rightarrow 0\) as \(h^{'},k^{'} \rightarrow 0.\) Hence the proof. \(\square\)

7 Examples

Example 1

Consider the following function

$$\begin{aligned} f(x,y)=0.26(x^{2}+y^{2})-0.48xy \end{aligned}$$

defined over the rectangular region \([0,1] \times [0,1].\)

For generating the data set, the length of the rectangle, i.e,  [0, 1] is partitioned into N equal parts, the width of the rectangle [0, 1] is partitioned into M equal parts. Then, the value of the function at the corner points of the subrectangles are calculated. For example, when \(N=2, \,\,\ M=2,\) the data set is:

$$\begin{aligned}{} & {} \Big \{(0,0,0), (0,0.5,0.065), (0,1,0.26), \\{} & {} \quad (0.5,0,0.065),(0.5,0.5,0.01), (0.5,1,0.085), \\{} & {} \quad (1,0,0.26),(1,0.5,0.085), (1,1,0.04)\Big \} \end{aligned}$$

Then, the IFS for this data set consists of 4 maps:

$$\begin{aligned}&W_{n,m}(x,y,z)=\big (L_{n}(x),K_{m}(y), F_{n,m}(x,y,z)\big )\\&\quad =\Big \{\big (L_{1}(x),K_{1}(y), F_{1,1}(x,y,z)\big ),\\&\quad \big (L_{1}(x),K_{2}(y), F_{1,2}(x,y,z)\big ), \\&\quad \big (L_{2}(x),K_{1}(y), F_{2,1}(x,y,z)\big ),\\&\quad \big (L_{2}(x),K_{2}(y), F_{2,2}(x,y,z)\big )\Big \} \end{aligned}$$

where

$$\begin{aligned}&L_{1}(x)=0.5x L_{2}(x)=-0.5x+1 \\&\quad K_{1}(y)=0.5y K_{2}(y)=-0.5y+1\\&\quad F_{1,1}(x,y,z) = 0.25z \\&\quad F_{1,2}(x,y,z)= -0.24x-0.26y+0.24xy+0.26+0.25z \\&\quad F_{2,1}(x,y,z) = -0.26x-0.24y+0.24xy+0.26+0.25z \\&\quad F_{2,2}(x,y,z)= -0.02x-0.02y\\&\qquad \qquad\qquad\quad +0.2084(1.0e-15)xy+0.04+0.25z \end{aligned}$$

The coefficients of these functions are obtained from the endpoint conditions on \(L_{n}, K_{m}\) and \(F_{n,m}.\) The numerical double integral results obtained using the formula (28) is provided in Table 1 along with the actual double integral value. The attractor of this IFS is shown in Fig. 1 along with the original graph of the function mentioned in Example 1.

Table 1 Table comparing the numerical double integral results obtained using the proposed method with the actual double integral value for Example 1
Full size table
Fig. 1
figure 1

Attractor of the IFS for Example 1 (left) and its original graph (right)

Full size image

Example 2

Consider the following function

$$\begin{aligned} f(x,y)=(x+2y-7)^{2}+(2x+y-5)^{2} \end{aligned}$$

defined over the rectangular region \([0,1] \times [0,1].\)

Similar to the above example, the data set is generated by partitioning the rectangle \([0,1] \times [0,1]\) into NM number of subrectangles and taking the corner points of each of the subrectangles along with their function values.

Then, the data set is:

$$\begin{aligned}{} & {} \Big \{(0,0,74), (0,0.5,56.25), (0,1,41), \\{} & {} \quad (0.5,0,58.25),(0.5,0.5,42.5), (0.5,1,29.25), \\{} & {} \quad (1,0,45),(1,0.5,31.25), (1,1,20)\Big \} \end{aligned}$$

The IFS, consisting of the 4 maps are as follows:

$$\begin{aligned} W_{n,m}(x,y,z)& =\big (L_{n}(x),K_{m}(y), F_{n,m}(x,y,z)\big )\\&=\Big \{\big (L_{1}(x),K_{1}(y), F_{1,1}(x,y,z)\big ),\\&\qquad \big (L_{1}(x),K_{2}(y), F_{1,2}(x,y,z)\big ), \\&\qquad \big (L_{2}(x),K_{1}(y), F_{2,1}(x,y,z)\big ),\\&\qquad \big (L_{2}(x),K_{2}(y), F_{2,2}(x,y,z)\big )\Big \} \end{aligned}$$

where

$$\begin{aligned}{} & {} L_{1}(x)=0.5x L_{2}(x)=-0.5x+1 \\{} & {} K_{1}(y)=0.5y K_{2}(y)=-0.5y+1\\{} & {} F_{1,1}(x,y,z) = -8.5x-9.5y+55.5+0.25z \\{} & {} F_{1,2}(x,y,z)= -4.5x+23.5y-4xy+22.5+0.25z \\{} & {} F_{2,1}(x,y,z) = 20.5x-5.5y-4xy+26.5+0.25z \\{} & {} F_{2,2}(x,y,z)= 16.5x+19.5y+10.75+0.25z \end{aligned}$$

The comparison of the numerical double integral value using fractal method is compared with actual value of the double integration in Table 2. Figure 2 shows the attractor of this IFS along with the original graph of the function mentioned in Example 2.

Table 2 Table comparing the numerical double integral results obtained using the proposed method with the actual double integral value for Example 2
Full size table
Fig. 2
figure 2

Attractor of the IFS for Example 2 (left) and its original graph (right)

Full size image

Example 3

Consider the following function

$$\begin{aligned} f(x,y)=(x^{2}+y-11)^{2}+(x+y^{2}-7)^{2} \end{aligned}$$

defined over the rectangular region \([0,1]\times [0,1].\)

Using the procedure explained in Example 1 and 2, the data set generated for this function is given by:

$$\begin{aligned}{} &\quad {} \Big \{(0,0,170), (0,0.5,155.8125), (0,1,136), \\{} & {} \quad (0.5,0,157.8125),(0.5,0.5,144.1250), (0.5,1,125.3125), \\{} & {} \quad (1,0,136),(1,0.5,123.3125), (1,1,106)\Big \}, \\{} & {} \quad \text {considering} \,\, N=2 \,\,\text {and} \,\, M=2. \end{aligned}$$

The IFS is then:

$$\begin{aligned} & W_{n,m}(x,y,z)=\big (L_{n}(x),K_{m}(y), F_{n,m}(x,y,z)\big )\\&\qquad\qquad\qquad =\Big \{\big (L_{1}(x),K_{1}(y), F_{1,1}(x,y,z)\big ),\\&\qquad\qquad\qquad\qquad \big (L_{1}(x),K_{2}(y), F_{1,2}(x,y,z)\big ), \\&\qquad\qquad\qquad\qquad \big (L_{2}(x),K_{1}(y), F_{2,1}(x,y,z)\big ),\\&\qquad\qquad\qquad\qquad \big (L_{2}(x),K_{2}(y), F_{2,2}(x,y,z)\big )\Big \} \\ & L_{1}(x)=0.5x L_{2}(x)=-0.5x+1 \\ & K_{1}(y) =0.5y K_{2}(y)=-0.5y+1\\ & F_{1,1}(x,y,z) = -2.6064x-4.6064y -0.6272xy \\&\qquad\qquad\qquad +122.0943+0.2818z \\ & F_{1,2}(x,y,z) = -2.2993x+28.2007y -1.9868xy\\&\qquad\qquad\qquad +94.0592+0.2467z \\ & F_{2,1}(x,y,z)= 30.2007x-4.2993y -1.9868xy\\&\qquad\qquad\qquad +94.0592+0.2467z \\ & F_{2,2}(x,y,z)= 26.5077x+24.5077y +0.6535xy\\&\qquad\qquad\qquad +70.0241+0.2116z \end{aligned}$$

The numerical double integral value along with the actual integral value is displayed in Table 3. The attractor of this IFS is displayed in Fig. 3 along with the original graph of the function mentioned in Example 3.

Table 3 Table comparing the numerical double integral results obtained using the proposed method with the actual double integral value for example 3
Full size table
Fig. 3
figure 3

Attractor of the IFS for Example 3 (left) and its original graph (right)

Full size image

Example 4

Consider the following function

$$\begin{aligned} f(x,y)=sin(x+y)+(x-y)^{2}-1.5x+2.5y+1 \end{aligned}$$

defined over the rectangular region \([0,1]\times [0,1].\)

With the above mentioned procedures, the data set for this function is:

$$\begin{aligned}{} & {} \Big \{(0,0,1), (0,0.5,2.9794), (0,1,5.3415), \\{} & {} \quad (0.5,0,0.9794),(0.5,0.5,2.3415), (0.5,1,3.9975), \\{} & {} \quad (1,0,1.3415),(1,0.5,1.9975), (1,1,2.9093)\Big \}, \end{aligned}$$

when \(N=2, M=2.\) The IFS, consisting of the 4 maps then will be:

$$\begin{aligned}W_{n,m}(x,y,z) &=\big (L_{n}(x),K_{m}(y), F_{n,m}(x,y,z)\big ) \\& =\Big \{\big (L_{1}(x),K_{1}(y), F_{1,1}(x,y,z)\big ),\\&\qquad \big (L_{1}(x),K_{2}(y), F_{1,2}(x,y,z)\big ), \\&\qquad \big (L_{2}(x),K_{1}(y), F_{2,1}(x,y,z)\big ),\\&\qquad \big (L_{2}(x),K_{2}(y), F_{2,2}(x,y,z)\big )\Big \} \end{aligned}$$

where

$$\begin{aligned}&L_{1}(x)=0.5x L_{2}(x)=-0.5x+1 \\& K_{1}(y)=0.5y K_{2}(y)=-0.5y+1\\& F_{1,1}(x,y,z) = -0.1271x+0.6249y\\&\qquad\qquad\qquad\quad +0.248xy+0.688+0.3120z \\& F_{1,2}(x,y,z)= -1.4258x-3.4028y\\&\qquad\qquad\qquad\quad+1.3709xy+5.1018+0.2397z \\&\quad F_{2,1}(x,y,z) = -0.4439x-0.3847y\\&\qquad\qquad\qquad\quad +1.3709xy+1.1018+0.2397z \\&\quad F_{2,2}(x,y,z)= 1.0170x-1.8173y\\&\qquad\qquad\qquad\quad -0.1657xy+2.7007+0.2086z \end{aligned}$$

The corresponding result of the numerical double integral value is provided in Table 4, with the actual integral value. Figure 4 represents the attractor of this IFS along with the original graph of the function mentioned in Example 4.

Table 4 Table comparing the numerical double integral results obtained using the proposed method with the actual double integral value for Example 4
Full size table
Fig. 4
figure 4

Attractor of the IFS for Example 4 (left) and its original graph (right)

Full size image

On observing the Tables 1, 2, 3 to 4, it can be concluded that:

  • The fractal numerical integral values obtained are much closer to their actual integral values even with fewer number of subrectangles. The error gets reduced even further, as the number of subrectangles in the partition is increased.

  • The close convergence of the fractal numerical integration results also establishes the method of selecting the vertical scaling factor.

Based on the observed results, the research work is significant in the following aspects:

  • Fractal integration techniques are significant in assessing the random characteristics of the signals. For example, in tumor analysis, it is trivial that the degree of vessel branching or the irregularity of a tumor boundary remains consistent over a broad range of magnifications. Since these structures exhibit recurring patterns at various magnifications, the euclidean measurement of specific parameters such as length and area may not give the desired results. However, the evaluation of these parameters can be done quickly and precisely using fractal numerical integration.

  • The obtained formula can be used directly to measure the volume of the region bounded by the fractal interpolation surfaces generated from the given set of values.

Nevertheless, the method of integration proposed here works for the data sets over two-dimensional interpolating domains only. But, the method can be extended to the higher dimensional data sets using multivariate fractal interpolation functions. Further, the method of numerical double integration presented in this work can be used in the volumetric evaluation of physical surfaces. Volume computation is a crucial engineering tool in determining the excavations and the amount of imported soil needed for the construction of tunnels, highways and similar structures. In [29], the authors have compared different interpolation methods in modelling the physical surfaces and in the calculation of volume. The results obtained through these methods can still be enhanced with the proposed method of numerical double integration.

8 Conclusion

Considering the recursive relation satisfied by the bivariate fractal interpolation functions, a novel method of numerical double integration is provided in this paper. The endpoint conditions of the iterated function system are used to derive the coefficients involved in the proposed integration formula. The obtained formula can be used directly to measure the volume of the region bounded by the fractal interpolation surfaces generated from the given set of values. Throughout this study, it has been observed that the bivariate fractal interpolation functions are best approximated with bilinear interpolation functions. It is shown that the recursive relation among these functions can be employed in establishing the closeness of the fractal numerical integration method to the conventional double integration method theoretically. Upon verifying the derived formula of numerical integration through a couple of examples, it is observed that the results obtained using the new formula converge to their actual integral values within a lesser number of iterations. The convergence of the integral results establishes the formulation of the vertical scaling factor provided in this paper. Future scope of this proposed work of numerical double integration involves the volumetric evaluation of irregular physical structures. Further, this research work can be extended to the numerical integration of the data sets over higher dimensional interpolating domains and evaluation of certain higher dimensional parameters.