  The fundamental aim of this paper is to provide the approximation and
numerical integration of a discrete set of data points with Bernstein fractal
approach. Using Bernstein polynomials in the iterated function system, the
paper initially proposes the numerical integration formula for the data set
corresponding to univariate functions. The proposed formula of integration is
shown to be convergent by examining the data sets of certain weierstrass
functions.
  The paper then extends the Bernstein fractal approximation and numerical
integration technique to two dimensional interpolating regions. Bernstein
polynomials defined over triangular domain has been used for the purpose. The
triangular domain has been partitioned and the newly generated points are
assigned colors in a particular manner to maintain the chromatic number as 3.
Following the above mentioned construction and approximation of bivariate
Bernstein fractal interpolation functions, the paper introduces the numerical
double integration formula using the constructed functions. The convergence of
the double integration formula towards the actual integral value of the data
sets is displayed with the help of some examples including the benchmark
functions. Both the newly introduced iterated function systems are verified for
their hyperbolicity and the resultant fractal interpolation functions are shown
to be continuous.
