Approximate solution of nonlinear fuzzy Fredholm integral equations using bivariate Bernstein polynomials with error estimation

1.   Introduction

Bernstein polynomials have been used recently to solve some linear as well as nonlinear fuzzy integral equations. These polynomials are positive and their sum is unity. The study of fuzzy integral equations was initially considered by Kaleva ^[32] to convert the initial value problem for first-order fuzzy differential equations into the equivalent fuzzy Volterra integral equation. Some authors investigated the existence of a unique solution for fuzzy integral equations using the Banach fixed point theorem (see ^{[9,10,28,36,37]}). Numerical methods for approximating the solution of fuzzy integral equations have been studied by many authors based on quadrature rules and Pcard's approximations (see ^{[11,12,15,16,23,38,43,48,49,53]}). Several numerical iterative approaches have been proposed based on the Picard approximations and other techniques such as the Lagrange interpolation ^[26], divided and finite differences ^[35], hybrid block-pulse functions and Taylor series ^[7], triangular basis functions ^[8] and block-pulse basis functions ^[50,53]. Some other existing numerical and analytical methods in literature are: Adomian decomposition ^[1], Nyström methods ^[2], Bernstein polynomials ^[22], fuzzy Haar wavelets ^[50] and block pulse functions ^[39], homotopy analysis ^[34] and homotopy perturbation ^[5]. Yang and Gong ^[47] introduced the concept of ill-posedness for the fuzzy Fredholm integral equation of the first kind using the Zadehs decomposition theorem of fuzzy number and obtained the approximate solution for this class of integral equations using regularization method based on classic Tikhonov regularization scheme. Early, Ziari et al. ^[55] presented an iterative method based on fuzzy Bernstein polynomials for solving nonlinear fuzzy Volterra integral equations. Recently, Ziari et al. ^[56] proposed an iterative method to solve fuzzy integral equations using generalized quadrature rule. The study of fuzzy integral equations in two dimensions based on iterative technique was started by Sadatrasoul and Ezzati (see ^[41]). Indeed, they developed the iterative method in ^[23] for two dimensional fuzzy integral equations. Also, Ezzati and Ziari in ^[24] proposed a non-iterative numerical method for solving fuzzy Fredholm integral equations in two dimensions based on Bernstein polynomials. Sadatrasoul and Ezzati in ^[42] presented an iterative method to find approximate solution of linear and nonlinear two dimensional Hammerstein fuzzy integral equations by considering optimal quadrature formula for fuzzy-number-valued functions of Lipschitz type. In ^[16], Bica and Popescu constructed the fuzzy trapezoidal cubature rule for fuzzy functions of Lipschitz type to approximate the solution of nonlinear fuzzy Fredholm integral equations in two variables. Recently, Akhavan et al. ^[3] developed an iterative approach and combining trapezoidal and midpoint formulas for such equations. In ^[33], the existence and uniqueness results for fuzzy Fredholm-Volterra integral have been investigated and the authors proposed an iterative numerical approach for these type of equations based on bivariate triangular functions. In this paper, we approximate the integral of two dimensional fuzzy function by the bivariate fuzzy Bernstein polynomials and investigate the error analysis of the numerical method. Also, an iterative procedure is constructed based on two dimensional fuzzy Bernstein polynomials for solving two dimensional nonlinear fuzzy Fredholm integral equations,

The main advantage of proposed method is the reduction of computational cost with respect to existing methods in the literature such as iterative methods based on fuzzy trapezoidal cubature rule in ^[16] and ^[41]. For instance, the number of function evaluations of the proposed method at the knots is $n^2$ for $m = n$ and the same factor for methods in ^[16] and ^[41] equals to $4n^2$ for $m = n$, thus, the proposed approach is reasonable from aspect of reduction of computational cost in comparison with the existing methods. We provide some numerical experiments to confirm the numerical results with theoretical results. Finally, some concluding remarks will be presented.

2.   Preliminaries

Definition 2.1. (See ^[20].) A fuzzy number is a function $w:R\to \left[0, \left. 1\right]\right.$ having the properties:

$(1)$ $w$ is normal, that is $\exists x_{0} \in R$ such that $w(x_{0}) = 1,$

$(2)$ $w$ is fuzzy convex set

$(3)$ $w$ is upper semicontinuous on $R$,

$(4)$ the $\overline{\left\{\left. \, x\in R:w(x) > 0\right\}\right. }$ is compact set.

The set of all fuzzy numbers is denoted by E.

Definition 2.2. (See ^[27].) An arbitrary fuzzy number is represented, in parametric form, by an ordered pair of functions $(\underline{w}(r), \overline{w}(r)), 0 \le r \le 1,$ which satisfy the following requirements:

$(1)$ $\underline{w}(r)$ is a bounded left continuous nondecreasing function over [0, 1],

$(2)$ $\overline{w}(r)$ is a bounded left continuous nonincreasing function over [0, 1],

$(3)$ $\underline{w}(r) \le \overline{w}(r)\; \; \; , \; \; \; 0 \le r \le 1$.

For addition and scalar multiplication of fuzzy numbers in fuzzy set space $E$ we have:

$(1)$ $w_1\oplus w_2 = (\underline{w_1}(r)+\underline{w_2}(r), \overline{w_1}(r)+\overline{w_2}(r)),$

$(2)$ $(\lambda\odot w_1) = \left\{\begin{array}{l} {(\lambda\underline{w_1}(r), \lambda\overline{w_1}(r))}\; \; \; \; \; {\lambda \ge 0, }\\ {(\lambda\overline{w_1}(r), \lambda\underline{w_1}(r))}\; \; \; \; \; {\lambda < 0.} \end{array}\right.$

Definition 2.3. (See ^[31].) For given fuzzy numbers $w_1 = (\underline{w_1}(r), \overline{w_1}(r))$ and $w_2 = (\underline{w_2}(r), \overline{w_2}(r))$, the quantity

is the Hausdorff distance between fuzzy numbers $w_1$ and $w_2$.

Lemma 2.1. (See ^[46].) If $w_1, w_2, w_3, w_4\in E$ and $k\in R$, then we have:

$(1) D(w_1\oplus w_3, w_2\oplus w_3) = D(w_1, w_2),$

$(2) D(k\odot w_1, k\odot w_2) = \left|k\right|D(w_1, w_2),$

$(3) D(w_1\oplus w_2, w_3\oplus w_4)\le D(w_1, w_3)+D(w_2, w_4),$

$(4) D(w_1\oplus w_2, \tilde{0}) \le D(w_1, \tilde{0}) + D(w_2, \tilde{0}).$

In ^[46], it is proved that $(E, D)$ is a complete metric space.

Lemma 2.2. (see ^[4].) For any $k_{1}, k_{2}\in \mathbb{R}$ with $k_{1}\cdot k_{2}\geq 0$ and any $w\in E$ we have $D\left(k_{1}\odot w, k_{2}\odot w\right) = \left\vert k_{1}-k_{2}\right\vert D\left(w, \widetilde{0}\right).$

Remark 2.1. The property (4) in Lemma (2.1) provide the definition of a function $\|.\|:E\rightarrow R^+$ by $\|w_1\| = D(w_1, \tilde{0})$, which has the properties of the usual norms. In ^[5] the properties of this function are exhibited as follows:

(i) $\|w_1\|\ge 0, \; \; \forall\; \; w_1 \in E$, and $\|w_1\| = 0$ iff $w_1 = \tilde{0}$,

(ii) $\|\lambda.w_1\| = |\lambda| \|w_1\|$ and $\|w_1\oplus w_2\|\le \|w_1\|+\|w_2\|, \; \; \forall \; \; w_1, w_2 \in E, \; \; \forall\; \; \lambda \in R$,

(iii) $|\|w_1\|-\|w_2\||\le D(w_1, w_2)$ and $D(w_1, w_2)\le\|w_1\|+\|w_2\| \; \; \forall \; \; w_1, w_2 \in E$.

Definition 2.4. (See ^[25].) A fuzzy-real-number-valued function $f:[a, b]\to E$ is said to be continuous in $t_{0} \in [a, b]$, if for given $\varepsilon > 0$ there exists a $\delta > 0$ such that $\left|x-x_{0} \right| < \delta.$ implies that $D(f(x), f (x_{0})) < \varepsilon$, for all $x\in [a, b]$ and We say that $f$ is fuzzy continuous on $[a, b]$ if $f$ is continuous at each $t_{0} \in [a, b]$, and we will use the notation $C([a, b], E)$ for the space of all continuous fuzzy functions.

Definition 2.5. (See ^[29].) Let $f:[a, b]\to E$ be given. The function $f$ is fuzzy-Riemann integrable to $I(f)\in E$ if for any $\varepsilon > 0$, there exists $\delta > 0$ such that for any division $P = \left\{\left[u, v\right]; \xi \right\}$ of $[a, b]$ with the norms $\Delta (p) < \delta$, we have

where $\sum$ denotes the fuzzy summation. In this case it is denoted by $I(f) = (FR)\int _a^b f(x)dx.$

Lemma 2.3. (See ^[29].) If $f, \; g:[a, b]\subseteq{R}\to{E}$ are fuzzy continuous functions, then the function $F:[a, b]\to R_+$ by $F(x) = D(f(x), g(x))$ is continuous on $[a, b]$ and

Theorem 2.1. (See ^[30].) If $f, g : [a, b]\rightarrow E$ are (FR) integrable fuzzy functions, and $\alpha$, $\beta$ are real numbers, then

Definition 2.6. (See ^[11].) For $L\geq 0,$ a function $f:\left[a, b \right] \rightarrow E$ is $L$-Lipschitz if

for any $x, y\in \left[a, b\right]$. A function $F:E \rightarrow E$ is Lipschitz if there exists $L^{\prime }\geq 0$ such that $D\left(F\left(u\right), F\left(v\right) \right) \leq L^{\prime }\cdot D\left(u, v\right)$, for any $u, v\in E.$

According to ^[11], any Lipschitz function is continuous.

For $f\in C_F [a, b]$, let us consider the Bernstein-type fuzzy polynomials

where $P_{n, k} (x)\; (k = 0, ..., n)$ are Bernstein polynomials of degree n defined on $[a, b]$ as:

It is obvious that $P_{n, k} (x)\ge 0, \; \; \; \forall \; \; x\in [a, b]$ and $P_{n, 0} (x), P_{n, 1} (x), ..., P_{n, n} (x)$ are linearly independent algebraic polynomials of degree $n$ and $\sum _{k = 0}^n p_{n, k} (x) = 1$.

The following definitions are related to fuzzy-number-valued functions in two variables.

Definition 2.9.(See ^[3].) Let $f:[a, b]\times[c, d]\rightarrow E$ be given$.$ The function $f$ is fuzzy-Riemann double integrable to $I(f)\in E$ if for any $\varepsilon > 0$, there exists $\delta > 0$ such that for any partition $\Delta = (\Delta_x, \Delta_y)$ of the region $[a, b]\times[c, d],$ with $\Delta_x:$ $a = x_{0} < x_{1} < \ldots x_{m-1} < x_{m} = b$ and $\Delta_y : c = y_{0} < y_{1} < \ldots < y_{n-1} < y_{n} = d$ and any set of intermediate points $\left\{(\xi_i, \eta_j):i = 1, 2, ..., m\; , j = 1, 2, ..., n\right\}$, with norm $\upsilon (\Delta) = \max\{\left\vert x_{i}-x_{i-1}\right\vert+\left\vert y_{j}-y_{j-1}\right\vert :i = 1, \ldots, m, \; j = 1, \ldots, n\} < \delta$, we have

In this case, the fuzzy number $I(f)$ will be denoted by

Definition 2.10. (See ^[40].) A function $f:\left[a, b\right] \times \lbrack c, d]\rightarrow E$ is called:

(i) continuous in $\left(x_{0}, y_{0}\right) \in \left[a, b\right] \times \lbrack c, d]$ if for any $\varepsilon > 0$ there exists $\delta > 0$ such that for any $\left(x, y\right) \in \left[a, b\right] \times \lbrack c, d]$ with $\left\vert x-x_{0}\right\vert < \delta,$ $\left\vert y-y_{0}\right\vert < \delta$ we have

The function $f$ is continuous on $\left[a, b\right] \times \lbrack c, d]$ if it is continuous in each $\left(x, y\right) \in \left[a, b\right] \times \lbrack c, d].$

(ii) bounded if there exists $M\geq 0$ such that

The set of all continuous functions $f:\left[a, b\right] \times\lbrack c, d]\rightarrow E$ is denoted by $C\left([a, b]\times\lbrack c, d], E\right).$

Lemma 2.4. (See ^[40].) If $f\in C\left([a, b]\times \lbrack c, d], E\right)$ then

exists and

Lemma 2.5. If $f, g:\left[a, b\right] \times \lbrack c, d]\rightarrow E$ are continuous fuzzy functions then the function $\varphi:\left[a, b\right] \times \lbrack c, d]\rightarrow \mathbb{R}_{+}$ defined by $\varphi\left(s, t\right) = D\left(f\left(x, y\right), g\left(x, y\right) \right)$ is continuous on $\left[a, b\right] \times \lbrack c, d]$ and

Proof. For simplicity we consider only $m = n$, where the $m$ and $n$ are the number of points in the partition of $[a, b]$ and $[c, d]$, respectively, the case of $m\neq n$ follows similarly. Firstly, we prove that the function $\varphi\left(s, t\right) = D(f(s, t), g(s, t))$ be continuous in every point $(s_0, t_0)\in[a, b]\times[c, d]$, for this purpose, we let $\{(s_n, t_n)\}_{n\ge1}, \; (s_n, t_n)\in [a, b]\times[c, d],$ such that $\lim_{n\rightarrow \infty} (s_n, t_n) = (s_0, t_0)$. In this case, we have:

and on the other hand, we have:

Taking to the limit when $n\rightarrow \infty$, and according to the continuity of $f\text{ and } g$ we obtain:

that is $\varphi$ is continuous at each $(s_0, t_0)\in[a, b]\times[c, d]$. Now, let $P_{n} = \{([s_{n-1}, s_{n}]; \xi _{n})\}, \; n\in N$ and $Q_{n} = \{([t_{n-1}, t_{n}]; \eta_{n})\}, \; n\in N$ be two sequences of partitions of $[a, b]$ and $[c, d]$ with $\Delta(P_n)\rightarrow 0,$ when $n\rightarrow \infty$ and $\Delta(Q_n)\rightarrow 0,$ when $n\rightarrow \infty$, respectively.

So, we have:

Taking to the limit when $n\rightarrow \infty$, we get:

Thus, the proof is complete.

Moreover, it can be proved that any continuous function $f:\left[a, b\right] \times \lbrack c, d]\rightarrow E$ is bounded.

Consider a fuzzy-number-valued function $f:[a, b]\times \lbrack c, d]\rightarrow E$ having the following Lipschitz property: there exist $L_{1}, L_{2}\geq 0$ such that

for all $x_{1}, x_{2}\in \lbrack a, b],$ $y_{1}, y_{2}\in \lbrack c, d].$

Similar to the one dimensional case, for $f\in C_{F} \left([a, b]\times[c, d] \right)$ we express the fuzzy two- dimensional Bernstein operators as follows:

Definition 2.10. (See ^[4].) Two dimensional fuzzy Bernstein operator on the region $[a, b]\times[c, d]$ is defined as follows:

where $x_i = a+\frac{i(b-a)}{m}, \; (i = 0, ..., m), \; y_j = c+\frac{j(d-c)}{n}, \; (j = 0, ...n)$ and $P_{m, i} (x), \; (i = 0, ..., m), \; P_{n, j} (x), \; (j = 0, ..., n)$ are Bernstein polynomials of degrees $m$ and $n$ respectively, moreover the polynomials $P_{m, i}(x) \text{ and } P_{n, j}\; \; \forall i, j,$ are defined on intervals $[a, b] \text{ and } [c, d],$ respectively.

Let $\varphi_{i, j}(x, y) = P_{m, i} (x).P_{n, j} (y)$, it is obvious that $\phi_{i, j}(x, y)\ge 0, \; \; \; \forall \; \; (x, y)\in [a, b]\times[c, d], \; \forall i, j, \text{ and }\varphi_{0, 0}(x, y), \varphi_{0, 1}(x, y), ..., \varphi_{m, n}(x, y)$ are linearly independent, and

3.   Approximation of integral of fuzzy function

The fuzzy function $f\in C_{F} \left([a, b]\times[c, d] \right)$ can be approximated using fuzzy Two dimensional Bernstein-type polynomials as

where $s_i = a+\frac{i(b-a)}{m}, \; (i = 0, ..., m), \; t_j = c+\frac{j(d-c)}{n}, \; (j = 0, ...n).$ Hence, the approximate value of the integral of fuzzy function can be obtained as follows:

In ^[21], the derivative of Bernstein polynomials on interval $[a, b]$ for fixed $m$ is expressed as follows:

By integrating from the above equality, we obtain:

namely, for fixed $m$ all Bernstein basis functions have the same definite integral over the interval [a, b]. Since $\sum _{i = 0}^m p_{m, i} (s) = 1$ we have:

Similarity for two dimensional Bernstein polynomials, we obtain:

So, we have the following approximation for double integral of two dimensional fuzzy function:

where $s_i = a+\frac{i(b-a)}{m}, \; (i = 0, ..., m), \; t_j = c+\frac{j(d-c)}{n}, \; (j = 0, ...n).$

In the following theorem, we obtain the error estimate of the above approximation.

Theorem 3.1. Let $f\in C_{F} \left([a, b]\times[c, d] \right)$ be a Lipschitzian function. Then we have:

Proof. According to (2.5), we have:

where $\varphi_{i, j}(s, t) = P_{m, i} (s).P_{n, j} (t)$. Then, applying the parts of 2 and 3 of Lemma 2.1, we obtain:

As regards, $f$ satisfies in Lipschitz condition namely, inequality (2.3) holds for $f$, we obtain:

It is well known following inequality which is expressed in ^[29]:

According to $\varphi_{i, j}(s, t) = P_{m, i} (s).P_{n, j} (t)$ and considering the $\sum _{i = 0}^mP_{m, i} (s) = 1, \; \sum _{j = 0}^nP_{n, j} (t) = 1,$ and inequality (3.3), we obtain:

4.   Two dimensional fuzzy integral equations

Here, we consider the two dimensional nonlinear fuzzy Fredholm integral Eq. (1.1), where H(x, y, s, t) is a crisp kernel function over $([a, b]\times[c, d])^2$, $f, \; F$ are continuous fuzzy-number-valued functions and $G:E\rightarrow E$ is a continuous fuzzy function. We assume that $H$ is continuous and therefore it is uniformly continuous with respect to $(s, t)$ and there exists $M_H > 0$, such that $M_H = \max\limits_ {a\le x, s\le b, \; c\le y, t\le d} |H(x, y, s, t)|$.

Let $\Omega = \{f : [a, b]\times[c, d]\rightarrow E; f$ is continuous$\}$ be the space of the two-dimensional fuzzy continuous functions with the metric $D^{\ast}(f, g) = \sup \limits_{a\leq s\leq b, \; c\leq t\leq d}D(f(s, t), g(s, t))$, for $f, g\in \Omega$. In the following theorem, sufficient conditions for the existence of an unique solution of Eq (1.1) are given.

Theorem 4.1. (See ^[41].) Let $H(x, y, s, t)$ be continuous and positive for $a\le x, s\le b, \; c\le y, t\le d$ and $f : [a, b]\times[c, d]\rightarrow E$ be continuous on $[a, b]\times[c, d]$. Moreover assume that there exists $L > 0$, such that

If $B = M_HL(b -a)(d-c) < 1$, then the iterative procedure for $k\in\bf{N}$

converges to the solution of $F$ of (1.1). In addition, the following error bound holds:

where $\Omega$ is the set of fuzzy continuous functions on $I = [a, b]\times[c, d]$.

Remark 4.1. If $F_0 = f$ the error estimate (4.2), becomes:

where $M_{0} = \sup\limits_{(x, y)\in I}\|G(\tilde 0)\|_F$.

Proof. It is observed that

and hence, it obtains the inequality (4.3).

Now, we introduce the numerical method to find the approximate solution of the two dimensional nonlinear fuzzy Fredholm integral equation (1.1). In this way, we consider the following uniform partitions of the region $[a, b]\times[c, d]$:

with $s_{i} = a+\frac{i(b-a)}{m}\; \; , \; \; t_{j} = c+\frac{j(d-c)}{n}\; \; , \; \; 0\le i\le m, \; 0\le j\le n$. Then according to (3.2), the following iterative procedure gives the approximate solution of Eq (1.1) in the point $(x, y)\in [a, b]\times[a, b]$ using two dimensional fuzzy Bernstein-type polynomials:

5.   Error analysis

In this section, we investigate the convergence of the iterative proposed method to the solution of Eq (1.1) under the following conditions:

(i) $f : [a, b]\times[c, d]\rightarrow E$ is fuzzy continuous;

(ii) $H : ([a, b]\times[c, d])^2\rightarrow R^{+}$ is continuous;

(iii) There exist $\alpha, \beta \geq 0$ such that

for any $s^\prime, s^{\prime\prime }\in [a, b], t^\prime, t^{\prime\prime }\in [c, d];$

(iv) There exists $L > 0$ such that

where $F_{1}, F_{2} : [a, b]\times[c, d]\rightarrow E$;

(v) $B = M_{H}L(b -a)(d-c) < 1,$ where $L$ is as given in the above item and $M_{H}\geq 0$ is such that $M_{H} = \max\limits_{a\le x, s\le b, \; c\le y, t\le d} H(x, y, s, t)$;

(vi) There exist $\mu, \lambda \geq 0$ such that

for any $x, s^{\prime }, s^{\prime\prime }\in \lbrack a, b], \text{ and } y, t^{\prime }, t^{\prime\prime } \in \lbrack c, d];$

(vii) There exist $\gamma, \eta \geq 0$ such that

for any $s, s^{\prime }, s^{\prime\prime }\in \lbrack a, b], t, t^{\prime }, t^{\prime\prime } \in \lbrack c, d];$

Firstly, we prove that functions $z_k\left(s, t\right)$ for all $k$ are Lipschitzian.

Lemma 5.1. Consider the iterative procedure (4.5). Under the conditions (i)–(vii) then functions $z_k$ for all $k\ge 0$ are Lipschitzian.

Proof. Using the properties of the distance between fuzzy numbers stated in Lemma 2.1 we obtain

Consequently, using conditions (iii) and (vii), we obtain:

By applying the properties of the norm function $\Vert.\Vert_{\digamma}$ in Remark 1, we obtain

By using again the properties of the norm function $\Vert.\Vert_{\digamma}$ in Remark 1, we have

By applying the inequality (5.2) into the above inequality, we have:

Taking supremum from the above inequality for $(s, t)\in [a, b]\times[c, d]$, it follows that

By successive substitutions on the above inequality, we obtain

Since, $\frac{1-B^{k}}{1-B}\leq\frac{1}{1-B},$ for all $k\in\bf{N}$, we obtain:

Taking supremum from the inequality (5.2), it follows that

Now, by taking into account the inequality (5.4), from the above inequality we get:

Finally, by considering the above inequality, from (5.1) we obtain:

By supposing $L_1 = \alpha+\gamma(b-a)(d-c)\frac{L\Vert f\Vert_{\digamma}+M_{0}}{1-B}$ and $L_2 = \beta+\eta(b-a)(d-c)\frac{L\Vert f\Vert_{\digamma}+M_{0}}{1-B}$, we have:

Thus, the functions $z_k(s, t)$ for all $k$ are Lipschitzian.

Now, we prove an interesting result about the satisfying functions $H\left(x, y, s, t\right)\odot G\left(z_k\left(s, t\right)\right)$ in Lipschitz condition which is used in the proof of the main result.

Lemma 5.2. Consider the iterative procedure (4.5). Under the conditions (i)–(vii) and supposing functions $z_k$ are Lipschitzian then the functions $\varphi_k(s, t) = H\left(x, y, s, t\right)\odot G\left(z_k\left(s, t\right)\right)$ are Lipschitzian.

Proof. Using fuzzy Distance, we have:

Using part (2) of Lemma 2.1, Lemma 2.2 and condition (iii), we obtain:

According to Lemma 5.1 and using condition (vi), we get:

By applying the inequality (5.5) from the proof of Lemma 5.1, we obtain

Now, by taking into account the above inequality, from (5.6) we get:

By assuming $L_1^\prime = \left(M_{H}L L_1+\frac{\mu(L\Vert f\Vert_{\digamma}+M_{0})}{1-B}\right)$ and $L_2^\prime = \left(M_{H}L L_2+\frac{\lambda(L\Vert f\Vert_{\digamma}+M_{0})}{1-B}\right).$ we have:

Hence, the functions $\varphi_k(s, t)$ for all $k$ are Lipschitzian.

Theorem 5.1. Under the conditions (i)–(vii) the iterative procedure Eq (4.5) converges to the unique solution of Eq (1.1), F, and its error estimate is as follows:

where $M_{0} = \sup\limits_{0\le x, y \le 1}\|F(\tilde 0)\|_F$.

Proof. Let $\varphi_k(s, t) = H(x, y, s, t)\odot G(F_k(s, t))$ and $\psi_k(s, t) = H(x, y, s, t)\odot G(z_k(s, t))$ for all $k\in \mathbb{N^\ast}$. Since

we have

Regarding to part (2) of Lemma 2.1, Theorem 3.1 and condition (iii), we have:

Applying Lemma 5.1 and Theorem 3.1 for second part of above expression, we obtain:

Taking supremum from the above inequality for $(x, y)\in [a, b]\times[c, d]$, and according to $D\left(F_{0} (s, t)), z_{0} (s, t)\right) = 0$, we deduce:

Now, since

we conclude that

Using again part (2) of Lemma 2.1, Theorem 3.1, condition (iii) and Lemma 5.1, we obtain:

Then, we have:

By induction, for $k\geq3$, we obtain:

Therefore, we have:

Since $B < 1$, we conclude that

Using Eqs (4.3) and (5.7), we obtain:

Remark 5.1. Since $B < 1$, it is easy to show that

this result shows that the proposed method is convergent.

6.   Numerical stability analysis

For the iterative numerical method, it is more suitable investigating the stability of the obtained numerical solution with respect to the choice of the first iteration. So, in order to study the numerical stability of the iterative method (4.5) with respect to small changes in the starting approximation, we consider an another initial iteration term $\hat{F}_0\in \Omega$ such that there exists $\epsilon > 0$ for which $D(F_0(x, y), \hat{F}_0(x, y) < \epsilon$, for all $(x, y)\in [a, b]\times[c, d].$ So, we have the new sequence of successive approximations as:

Using the same iterative method, the terms of produced sequence are:

Definition 6.1. We say that the method of successive approximations applied to solve Eq. (1.1) is numerically stable with respect to the choice of the first iteration iff there exist positive numbers $p, \; q$ and constants $k_1, k_2, k_3 > 0$ which are independent by step-sizes $h_1 = \frac{b-a}{m}$ and $h_2 = \frac{d-c}{n}$ respectively, such that

Theorem 6.1. Under the conditions of Theorem 5.1, the iterative method (6.2) is numerically stable with respect to the choice of the first iteration.

Proof. Clearly, we observe that:

Using Eq (5.7), we obtain:

Also, by similar reasoning we have:

In order to obtain the upper bound of $D(F_k(x, y), \hat{F}_k(x, y))$, we observe that

and thus

Therefor, we obtain:

By induction for $k \ge 2$, we get:

According to $B < 1$, this inequality becomes:

Then, from eq-stab1, we get:

By comparing inequalities (6.3) and (6.7), we deduce that

So, the proof is complete.

7.   Numerical experiments

In this section, we present numerical results obtained by applying the proposed algorithm on two examples in order to demonstrate the validity and accuracy of the method. The approximate solution is calculated for different values of $k$, $m$ and $n$.

Example 1. Consider

where

The above fuzzy equation has the exact solution

The results of the proposed iterative algorithm in given points are shown in terms of errors in Tables 1 and 2.

This method well performs for linear two dimensional fuzzy integral equation.

Example 2. Consider

where

The above fuzzy equation has the exact solution

The results of the proposed iterative algorithm in given points are shown in terms of errors in Tables 3 and 4.

8.   Conclusions

In this paper, we have developed an numerical iterative method for the solution of two dimensional fuzzy Fredholm integral equations of the second kind. The method presented in this work was applied Picard iterations and numerical approximation to evaluate the iterated integral based on bivariate Bernstein polynomials. Moreover, we have considered the conditions to prove the sequence of approximations of the solution at given point converge to the exact solution. The concept of numerical stability is defined based on the choice of the first iteration and then the numerical stability of the proposed iterative algorithm is proven. Finally, some numerical examples are included to examine the accuracy and the convergence of the proposed method.

Conflict of interest

The authors declare that there is no conflict of interest.