Iterative numerical method for nonlinear fuzzy Volterra integral equations

1 Introduction

The present paper is devoted to the construction of an iterative numerical method for the nonlinear fuzzy Volterra integral equation x ( t ) = g ( t ) + ( FR ) ∫ 0 t f ( t , s , x ( s ) ) ds , t ∈ [ 0 , T ](1) where f : [ 0 , T ] × [ 0 , T ] × ℝ F → ℝ F and g : [0, T ] → ℝ F are continuous. The study of fuzzy Volterra integral equations is motivated by their applications in control mathematical models (see [13]) and the approach of them starts with the papers of Kaleva (see [22]), Seikkala (see [37]) and Mordeson and Newman (see [27]).

The main problems concerning fuzzy integral equations are related to the study of the existence, existence and uniqueness, and boundedness of the solution. The existence is proved with the Darbo’s fixed point theorem, and the existence and uniqueness of the solution is investigated by using the Banach’s fixed point principle (see [3, 37–39]. Another special tools for the existence of the solution are the Banach-Alaoglu and the stacking theorem (see [1]) and the Adomian decomposition method (see [33]).

Recently, random fuzzy fractional integral equations of Volterra type were investigated for the existence and uniqueness of the solution, in the case of the Caputo’s type fractional integral (see [24]) and of the Riemann-Liouville’s type fractional integral (see [25]). In both these two papers, the tool was the successive approximations method. In [25] two types of equations are approached: with solutions of nondecreasing and nonincreasing diameter, respectively, providing also the boundedness of the solution and its continuous dependence by the data.

The numerical methods for fuzzy integral equations are based on various techniques. The method of successive approximations and other iterative techniques are applied in [7,8,10,14,15, 7,8,10,14,15] and [29]. The analytic-numeric methods like Adomian decomposition, variation iteration, homotopy analysis and homotopy perturbation are used in [2,12,17,16, 2,12,17,16] and [33]. Other techniques used in the construction of the numerical methods for fuzzy Volterra integral equations are: Bernstein polynomials (see [28]), Chebyshev polynomials (see [5]), quadrature rules and Nystrom techniques (see [35]), predictor-corrector techniques (see [36]), fuzzy differential transform methods (see [34]), expansion methods (see [23]) such as: fuzzy collocation and Galerkin methods (see [11]), or Taylor series (see [21]).

In this paper we prove the uniformly boundedness and uniform Lipschitz properties of the sequence of successive approximations associated to the integral Equation (1) and propose an iterative numerical method for approximating the solution of (1). The convergence of the method is proved, and afterthere is tested on two numerical experiments confirming the theoretical results. The paper is organized as follows: in Section 2 we present some preliminary results concerning fuzzy numbers and fuzzy-number-valued functions that will be used in the next sections. Section 3 is devoted to provide some properties of the sequence of successive approximations and to present the iterative algorithm. In Section 4 we prove the convergence of the method, and in Section 5 we present two numerical examples in order to test and to illustrate the accuracy and the convergence of the proposed algorithm.

We remember the notion of fuzzy number and some results involving the structure of the set of fuzzy numbers and the properties of fuzzy-number-valued functions concerning to their integrability.

Definition 1. (see [6] and [18]) A fuzzy subset of the real line u : ℝ → [ 0 , 1 ] is a fuzzy number if it satisfies the following properties:

u is normal, i.e. ∃ x 0 ∈ ℝ such that u (x₀) = 1 ;

The set of all fuzzy numbers is denoted by ℝ F . Since each real number u can be identified with the function u ¯ : ℝ → [ 0 , 1 ] , u ¯ ( x ) = { [ c ] c 1 , x = u 0 , x ≠ u , we infer that each real number is a fuzzy number, and so ℝ ⊂ ℝ F . For a fuzzy number u ∈ ℝ F and for r ∈ (0, 1] are defined the subsets of ℝ, u r = { x ∈ ℝ : u ( x ) ≥ r } being called the r-level sets and the set u 0 = cl { x ∈ ℝ , u ( x ) > 0 } is called the support of u. The 1-level set is called the core of u. If u 1 - = u 1 + then the corresponding fuzzy number is called unimodal.

It is proved that each r-level set u r = [ u r - , u r + ] is a closed interval for any r ∈ [0, 1], and 0 ≤ r₁ ≤ r₂ ≤ 1 implies that u _{r 2}  ⊂ u _{r 1} (see [6]).

From calculus point of view, more convenient is another representation of a fuzzy number, the LU-representation ( u r - , u r + ) (lower-upper representation) introduced in [18]. Therefore, in [18] is established the connection between the properties of a fuzzy number presented in Definition 1 and the LU-representation of a fuzzy number. More precisely, it is proved the following result:

Theorem 1. (see [18]) A subset u of the real line with its r-level sets u r = [ u r - , u r + ] = { x ∈ ℝ : u ( x ) ≥ r } is a fuzzy number if and only if the functions u - , u + : [ 0 , 1 ] → ℝ, defining the end-points of the r-level sets, satisfy the following conditions:

the function u^- given by u - ( r ) = u r - , r ∈ [0, 1] , is bounded, non-decreasing, left-continuous in (0, 1] and it is right-continuous at 0 ;

For u , v ∈ ℝ F , k ∈ ℝ , the addition and the scalar multiplication are defined by [ u + v ] r = [ u ] r + [ v ] r = [ u - r + v - r , u + r + v + r ] , ∀r ∈ [0, 1] [ k · u ] r = k · [ u ] r = { [ ku - r , ku + r ] , if k ≥ 0 [ ku + r , ku - r ] , if k < 0 .

According to [6] we can summarize the following algebraic properties:

u + (v + w) = (u + v) + w and u + v = v + u for any u , v , w ∈ ℝ F ,

The Hausdorff metric (see [6]) is defined onℝ F by v + r | ) } for any u , v ∈ ℝ F .

Lemma 2. (see [6]) This metric has the following properties:

Guided by the property (vi) from Lemma 2, in [16] it is defined a function ‖ · ‖ : ℝ F → ℝ by that has the properties of usual norms:

∀ u ∈ ℝ F and iff u = 0 ˜ ,

Of course, ( ℝ F , + , · , ‖ · ‖ ) is not a normed space because ( ℝ F , + ) it is not a group.

Considering a fuzzy-number-valued function f : I ⊂ ℝ → ℝ F we can define the functions f - r , f + r : I ⊂ ℝ → ℝ, r ∈ [0, 1] by f - r ( t ) = ( f ( t ) ) - r , f + r ( t ) = ( f ( t ) ) + r , ∀t ∈ I, ∀r ∈ [0, 1] . These functions are called the left and right r-level functions of f.

Definition 2. (see [6]) A fuzzy-number-valued function f : [ a , b ] → ℝ F is said to be continuous at t₀ ∈ [a, b] if for each ɛ > 0 there is δ > 0 such that D (f (t) , f (t₀)) < ɛ, whenever

The boundedness of a fuzzy-number-valued function can be expressed in the following form: f : [ a , b ] → ℝ F is bounded iff there is M ≥ 0 such that D ( f ( t ) , 0 ˜ ) ≤ M for all t ∈ [a, b] .

Lemma 3. (see [41]). If f : [ a , b ] → ℝ F is continuous then it is bounded and its supremum sup t ∈ [ a , b ] f ( t ) must exist and is determined by u ∈ ℝ F with u - r = sup t ∈ [ a , b ] f - r ( t ) and u + r = sup t ∈ [ a , b ] f + r ( t ) . A similar conclusion for the infimum is also true.

On the set C ( [ a , b ] , ℝ F ) = { f : [ a , b ] → ℝ F ∣ f continuous } it is defined the metric D ∗ ( f , g ) = sup t ∈ [ a , b ] D ( f ( t ) ,g (t)) , ∀ f , g ∈ C ( [ a , b ] , ℝ F ). We see that (C ( [ a , b ] , ℝ F ) , D ∗ ) is a complete metric space and using the Lemmas 3 and 4 we can derive corresponding properties of the metric D^∗ .

In [40] the notion of Henstock integral for fuzzy-number-valued functions is defined as follows:

Definition 3. (see [40]) Let f : [ a , b ] → ℝ F. For Δ _n  : a = x₀ < x₁ < . . . < x_n-1 < x _n  = b a par-tition of the interval [a, b], we consider thepoints ξ _i  ∈ [x_i-1, x _i ] , i = 1, . . . , n, and the func-tion δ : [ a , b ] → ℝ +. The partition P = {([x_i-1, x _i ] ; ξ _i ) ; i = 1, . . . , n} denoted by P = (Δ _n , ξ) is called δ-fine iff [x_i-1, x _i ] ⊆ (ξ _i  - δ (ξ _i ) , ξ _i  + δ (ξ _i )). For I ∈ ℝ F , the function f is fuzzy Henstock integrable on [a, b] if for any ε > 0 there is a function δ : [ a , b ] → ℝ + such that for any partition δ-fine P, D ( ∑ i = 1 n ( x i - x i - 1 ) · f ( ξ i ) , I ) < ε. The fuzzy number I is named the fuzzy Henstock integral of f and will be denoted by ( FH ) ∫ a b f ( t ) dt.

In the case that the function δ : [ a , b ] → ℝ + is constant, we obtain the Riemann integrability for fuzzy-number-valued functions (see [18]). In this case, I ∈ ℝ F is called the fuzzy Riemann integral of f on [a, b] , being denoted by ( FR ) ∫ a b f ( t ) dt . So, the fuzzy-Riemann integrability is a particular case of the fuzzy-Henstock integrability, and therefore the properties of the integral (FH) will be valid for the integral (FR) , too.

Lemma 4. (see [19]). Let f : [ a , b ] → ℝ F. Then f is (FH) integrable if and only if f - r and f + r are Henstock integrable for any r ∈ [0, 1] . Furthermore, for any r ∈ [0, 1] , [ ( FH ) ∫ a b f ( t ) dt ] r = [ ( H ) ∫ a b f - r ( t )dt , ( H ) ∫ a b f + r ( t ) dt ] .

Remark 1. If f : [ a , b ] → ℝ F is continuous, then f - r and f + r are continuous for any r ∈ [0, 1] and consequently, they are Henstock integrable. According to Lemma we infer that f is (FH) integrable.

Lemma 5. (see [7]). If f and g are fuzzy-Henstock integrable functions and if the function given by D (f (t) , g (t)) is Lebesgue integrable, then

D ( ( FH ) ∫ a b f ( t ) dt , ( FH ) ∫ a b g ( t ) dt ) ≤ ( L ) ∫ a b D ( f ( t ) , g ( t ) ) dt .

Definition 4. (see [7]) For L ≥ 0, a function f : [ a , b ] → ℝ F is L-Lipschitz if D ( f ( x ) , ( f ( x ) ) ≤ L | x − y | for any x, y ∈ [a, b]. A function F : ℝ F → ℝ F is Lipschitz if there exists L′ ≥ 0 such that D (F (u) , F (v)) ≤ L′ · D (u, v), for any u , v ∈ ℝ F .

It is obvious that any Lipschitzian function f : [ a , b ] → ℝ F is continuous.

Lemma 6. (see [7]) Let f : [ a , b ] → ℝ F be a L-Lipschitz function. Then:

D ( ( FR ) ∫ a b f ( t ) dt , [ b - a 2 · f ( a ) + b - a 2 · f ( b ) ] ) ≤ L · ( b - a ) 2 4 .

The trapezoidal quadrature formula, presented above, can be extended for uniform partitions, Δ : a = t 0 < t 1 < . . . < t n = b with t i = a + i · b - a n , ∀ i = 0 , n ¯ , and the corresponding remainder estimate is obtained in [8]:

Lemma 7. (see [8]) For uniform partitions, the following trapezoidal inequality holds:

D ( ( FR ) ∫ a b f ( t ) dt , ∑ i = 0 n - 1 ( t i + 1 - t i ) 2 · [ f ( t i ) + f ( t i + 1 ) ] ) ≤ L ( b - a ) 2 4 n .(2)

In order to obtain a convergent algorithm to approximate the solution of ((1)) we consider the following conditions:

g ∈ C ( [ 0 , T ] , ℝ F ) , f ∈ C ([0, T] × [0, T] × ℝ F , ℝ F ) ;

In [10] we have proved the following result:

Lemma 8. (see [10]) If f ∈ C ( [ 0 , T ] × [ 0 , T ] , ℝ F ) ,g ∈ C ( [ 0 , T ] , ℝ F ) , and a ∈ C ( [ 0 , T ] , ℝ + ) then the functions a · g : [ 0 , T ] → ℝ F and F : [ 0 , T ] → ℝ F given by (a · g) (t) = a (t) · g (t) , ∀t ∈ [0, T] and F ( t ) = ( FH ) ∫ 0 t f ( t , s ) ds are continuous.

Consider the sequence of successive approximations associated to the integral Equation (1):

x m ( t ) = g ( t ) + ( FR ) ∫ 0 t f ( t , s , x m - 1 ( s ) ) ds , t ∈ [ 0 , T ] , m ∈ ℕ ∗(3) and define the sequence of functions ( F m ) m ∈ ℕ ∗, F m : [ 0 , T ] × [ 0 , T ] → ℝ F, given by F _m  (t, s) = f (t, s, x _m  (s)) , s ∈ [0, T] , m ∈ ℕ .

Theorem 9. Under the conditions (i)–(iii) the integral Equation (1) has unique solution in C ( [ 0 , T ] , ℝ F ) , x ∗ ∈ C ( [ 0 , T ] , ℝ F ) and the sequence of successive approximations ( x m ) m ∈ ℕ ⊂ C ( [ 0 , T ] , ℝ F ) , given in ((3)) converges to x^∗ in C ( [ 0 , T ] , ℝ F ) for any choice of x 0 ∈ C ( [ 0 , T ] , ℝ F ) . In addition, the following error estimates hold:

D ( x ∗ ( t ) , x m ( t ) ) ≤ [ α T ] m 1 - α T · D ( x 1 ( t ) , x 0 ( t ) ) , ∀ t ∈ [ 0 , T ] , m ∈ ℕ ∗(4) and

D ( x ∗ ( t ) , x m ( t ) ) ≤ α T 1 - α T · D ( x m ( t ) , x m - 1 ( t ) ) , ∀ t ∈ [ 0 , T ] , m ∈ ℕ ∗ .(5)

Choosing x 0 ∈ C ( [ 0 , T ] , ℝ F ) , x₀ = g, the inequality ((3)) becomes D ( x ∗ ( t ) , x m ( t ) ) ≤ [ α T ] m 1 - α T · M 0 T , ∀ t ∈ [ 0 , T ] , m ∈ ℕ ∗ ,(6) where M₀ ≥ 0 is such that D ( f ( t , s , g ( s ) ) , 0 ˜ ) ≤ M 0 , ∀t, s ∈ [0, T] . Moreover, under the conditions (i)-(v), the sequences ( x m ) m ∈ ℕ ∗ and ( F m ) m ∈ ℕ ∗ are uniformly bounded and uniform Lipschtz.

Proof. Consider the set F ( [ 0 , T ] , R F ) = { f ∣ f : [ 0 , T ] → R F } and define the operator A : C ( [ 0 , T ] , ℝ F ) → F ( [ 0 , T ] , ℝ F ) by

A ( x ) ( t ) = g ( t ) + ( FR ) ∫ 0 t f ( t , s , x ( s ) ) ds , ∀ t ∈ [ 0 , T ] , ∀ x ∈ C ( [ 0 , T ] , ℝ F ) .(7)

Firstly, we intend to show that A (C ([0, T] , ℝ F ) ) ⊂ C ( [ 0 , T ] , ℝ F ) . Therefore, considering arbitrary x ∈ C ( [ 0 , T ] , ℝ F ) , s₀ ∈ [0, T] and ɛ > 0, by the continuity of x, for δ ( ɛ ) = ɛ 2 γ > 0 we infer that D ( x ( s ) , x ( s 0 ) ) < ɛ 2 α for any s ∈ [0, T] with Then, D (f (t, s, x ( s 0 ) ) ≤ γ · ɛ 2 γ + α · ɛ 2 α = ɛ and the function F x : [ 0 , T ] → ℝ F, defined by F _x  (t ; s) = f (t, s, x (s)), is continuous in s₀ for any fixed t ∈ [0, T] . Now, by Lemma 8 we see that the function G x : [ 0 , T ] → ℝ F, defined by G x ( t ) = ∫ 0 t f ( t , s , x ( s ) ) ds is continuous on [0, T] for any x ∈ C ( [ 0 , T ] , ℝ F ) . Since g ∈ C ( [ 0 , T ] , ℝ F ) , we conclude that A (x) is continuous on [0, T] for any x ∈ C ( [ 0 , T ] , ℝ F ) , and consequently, A ( C ( [ 0 , T ] , ℝ F ) ) ⊂ C ( [ 0 , T ] , ℝ F ) . Now, we will apply the Banach’s contraction principle to the operator A : C ( [ 0 , T ] , ℝ F ) → C ( [ 0 , T ] , ℝ F ). For arbitrary u, v ∈ C ([0, T] , ℝ F ), we get D ( A ( u ) ( t ) , A ( v ) ( t ) ) ≤ D ( ( FR ) ∫ 0 t f ( t , s , u ( s ) ) ds , ( FR ) ∫ 0 t f ( t , s , v ( s ) ) ds ) ≤ ∫ 0 t D ( f ( t , s , u ( s ) ) , f ( t , s , v ( s ) ) ) ds ≤ ∫ 0 t α D ( u ( s ) , v ( s ) ) ds ≤ α T · D ∗ ( u , v ) for all t ∈ [0, T] . So, D ∗ ( A ( u ) , A ( v ) ) ≤ α T · D ∗ ( u , v ) , ∀ u , v ∈ C ( [ 0 , T ] , ℝ F ) and according to the condition (iii), A is a contraction mapping. Applying the Banach’s fixed point principle we obtain the existence and uniqueness of the solution, x ∗ ∈ C ( [ 0 , T ] , ℝ F ) , of ((1)) and the uniform convergence of the sequence of successive approximations ((2)) to this solution in ( C ( [ 0 , T ] , ℝ F ) , D ∗ ) , for any choice of the initial term x 0 ∈ C ( [ 0 , T ] , ℝ F ) . The same Banach’s fixed point principle leads to the estimates ((3)) and ((4)). If we choose x₀ = g, then there are M _g , M₀ ≥ 0 such that D ( x 0 ( s ) , 0 ˜ ) = D ( g ( s ) , 0 ˜ ) ≤ M g and D ( F 0 ( t , s ) , 0 ˜ ) = D ( f ( t , s , x 0 ( s ) ) , 0 ˜ ) ≤ M 0 , for all t , s ∈ [ 0 , T ] .(8)

Therefore, D ( x 1 ( t ) , x 0 ( t ) ) ≤ D ( g ( t ) , g ( t ) ) + D ( ( FR ) ∫ 0 t f ( t , s , g ( s ) ) ds , ( FR ) ∫ 0 t 0 ˜ ds ) ≤ ∫ 0 t D ( f ( t , s , g ( s ) ) , 0 ˜ ) ds ≤ ∫ 0 t M 0 ds = M 0 T , ∀ t ∈ [ 0 , T ] , and from ((3)) we obtain now, ((5)). In the purpose to obtain the uniform boundedness of the sequences ( x m ) m ∈ ℕ ∗ and ( F m ) m ∈ ℕ ∗ , by induction, we get D (x _m  (t) , x_m-1 (t))≤∫ 0 t D ( f ( t , s , x m - 1 ( s ) ) , f ( t , s , x m - 2 ( s ) ) ) ds ≤ α ∫ 0 t D ( x m - 1 ( s ) , x m - 2 ( s ) ) ds ≤ α T · D ∗ ( x m - 1 , x m - 2 ) ≤ . . . ≤ ( α T ) m - 1 · D ∗ ( x 1 , x 0 ) , for all t ∈ [0, T] , m ∈ ℕ ∗ . Consequently, D (x _m  (t) , x₀ (t)) ≤ D (x _m  (t) , x_m-1 (t)) + D (x_m-1 (t) , x_m-2 (t)) + . . . + D (x₁ (t) , x₀ (t)) ≤ [(αT) ^m-1 + ( α T ) m - 2 + . . . + 1 ] · D ∗ ( x 1 , x 0 ) ≤ 1 - ( α T ) m 1 - α T · M 0 T , ∀ t ∈ [ 0 , T ] , m ∈ ℕ ∗ and D ( x m ( t ) , 0 ˜ ) ≤ D ( x m ( t ) , x 0 ( t ) ) + D ( x 0 ( t ) , 0 ˜ ) ≤ M 0 T 1 - α T + M g = notation R(9) for any t ∈ [0, T] and for all m ∈ ℕ . So, the sequence ( x m ) m ∈ ℕ is uniformly bounded in C ( [ 0 , T ] , ℝ F ) . On the other hand, D ( F m ( t , s ) , 0 ˜ ) ≤ D ( f ( t , s , x m ( s ) ) , f ( t , s , x 0 ( s ) ) ) + D ( f ( t , s , x 0 ( s ) ) , 0 ˜ ) ≤ α D ( x m ( s ) , x 0 ( s ) ) + M 0 ≤ α M 0 T 1 - α T + M 0 = notation M , ∀ t , s ∈ [ 0 , T ] , m ∈ ℕ ∗ .(10)

Then, the sequence of functions ( F m ) m ∈ ℕ is uniformly bounded in C ( [ 0 , T ] × [ 0 , T ] , ℝ F ), too.

From the condition (iv) we have D (x₀ (t) , x₀ (t′)) for all t, t′ ∈ [0, T] , and for m ∈ ℕ ∗ it follows that

D ( x m ( t ) , x m ( t ′ ) ) ≤ D ( g ( t ) , g ( t ′ ) ) + D ( ( FR ) ∫ 0 t f ( t , s , x m - 1 ( s ) ) ds , ( FR ) ∫ 0 t ′ f ( t ′ , s , x m - 1 ( s ) ) ds )

≤ ( β + μ T + M ) | t − t ′ |(11) for any t, t′ ∈ [0, T] , m ∈ ℕ ∗ , and then the sequence ( x m ) m ∈ ℕ is uniformly Lipschitz with the Lipschitz constant L₀ = β + μT + M . Finally, for arbitrary t, t′, s, s′ ∈ [0, T] we obtain, and D ( F m ( t , s ) , F m ( t ′ , s ′ ) ≤ μ | t − t | + γ | s − s ′ |(12) + α D ( x m s , x m ( s ′ ) ) ≤ μ | t − t ′ | + γ + α L 0 | s − s ′ |(13) for all m ∈ ℕ ∗ . So, the sequence of functions ( F m ) m ∈ ℕ is uniformly Lipschitz. With respect to s the Lipschitz constant is L = γ + αL₀ = γ + α (β + μT + M) .

Now, we consider a uniform partition of [0, T] with the knots t _i  = ih, i = 0 , n ¯ , h = T n , and apply the quadrature rule ((2)) in the computation of the terms of the sequence of successive approximations ((3)), x 0 ( t i ) = g ( t i ) , i = 0 , n ¯ , x m ( t 0 ) = g ( t 0 ) , ∀ m ∈ ℕ ∗ ,

x m ( t i ) = g ( t i ) + ( FR ) ∫ 0 t i F m - 1 ( t i , s ) ds , i = 1 , n ¯ , m ∈ ℕ ∗(14) obtaining x m ( t i ) = g ( t i ) + T 2 n · ∑ j = 1 i [ f ( t i , t j - 1 , x m - 1 ( t j - 1 ) ) + f ( t i , t j , x m - 1 ( t j ) ) ] + R m , i , i = 1 , n ¯ for all m ∈ ℕ ∗ , with the remainder estimate D ( R m , i , 0 ˜ ) ≤ LT 2 4 n , ∀ i = 1 , n ¯ , m ∈ ℕ ∗ .(15)

In this way, it obtains the following iterative algorithm:

Step 1: For given n ∈ ℕ ∗ and for i = 0 , n ¯ compute x 0 ¯ ( t i ) = x 0 ( t i ) = g ( t i ) , and x m ( t 0 ) = g ( t 0 ) , ∀ m ∈ ℕ ∗ ;(16)

Step 2 (the first iterative step): For m = 1 and for all i = 1 , n ¯, compute x 1 ( t i ) ¯ = g ( t i ) + T 2 n · ∑ j = 1 i [ f ( t i , t j - 1 , x 0 ¯ ( t j - 1 ) ) + f ( t i , t j , x 0 ¯ ( t j ) ) ] ;(17)

Steps 3-4 (the generic iterative step): By induction for m ∈ ℕ ∗ , m ≥ 2, we use the values computed at the previous step and for i = 1 , n ¯ , it obtains x m ( t i ) ¯ = g ( t i ) + T 2 n · ∑ j = 1 i [ f ( t i , t j - 1 , x m - 1 ¯ ( t j - 1 ) ) + f ( t i , t j , x m - 1 ¯ ( t j ) ) ] .(18)

This algorithm has the following practical stopping criterion:

Steps 3-4 (a condition of “do-while” type): For previously given ɛ > 0, if D ( x ¯ m ( t i ) , x ¯ m - 1 ( t i ) ) < ɛ , for all i = 1 , n ¯(19) then we stop to this "m" and retain the values x ¯ m ( t i ) , i = 0 , n ¯ , computed at this last iterative step. This condition is active after Step 2.

Step 5: Print “m” and print x ¯ m ( t i ) , i = 0 , n ¯. STOP.

The above presented algorithm will be convergent if we prove that lim m , n → ∞ D ( x ∗ ( t i ) , x ¯ m ( t i ) ) = 0 for all i = 1 , n ¯ . This can be done by providing an apriori error estimate.

Theorem 10. Under the conditions (i)–(v), the solution of the integral Equation ((1)) is approximated on the knots t i = i · T n , i = 0 , n ¯ , by the sequence ( x ¯ m ( t i ) ) m ∈ ℕ , i = 0 , n ¯ , with the following apriori error estimate:

D ( x ∗ ( t i ) , x ¯ m ( t i ) ) ≤ [ α T ] m 1 - α T · M 0 T + LT 2 4 n ( 1 - α T )(20) for all i = 1 , n ¯ , m ∈ ℕ ∗ .

Proof. Because D ( x ∗ ( t i ) , x ¯ m ( t i ) ) ≤ D ( x ∗ ( t i ) , x m ( t i ) ) + D ( x m ( t i ) , x ¯ m ( t i ) ) , and by ((5)) the apriori estimate on the knots t _i , i = 0 , n ¯ , is: D ( x ∗ ( t i ) , x m ( t i ) ) ≤ [ α T ] m 1 - α T · M 0 T , ∀ i = 0 , n ¯ , m ∈ ℕ ∗ .(21)

So, it remains to obtain the estimates for D ( x m ( t i ) , x ¯ m ( t i ) ) , for all i = 0 , n ¯ , m ∈ ℕ ∗ .From ((14)) for m = 1 and by ( (17)) and ( (15)) we get D ( x 1 ( t i ) , x ¯ 1 ( t i ) ) ≤ D ( R 1 , i , 0 ˜ ) ≤ LT 2 4 n , ∀ i = 0 , n ¯ . Now, using ((15)) and ((18)) it follows that: D ( x m ( t i ) , x ¯ m ( t i ) ) ≤ D ( R m , i , 0 ˜ ) + T 2 n · ∑ j = 1 n [ D ( f ( t i , t j - 1 , x m - 1 ( t j - 1 ) ) , f ( t i , t j - 1 , x m - 1 ¯ ( t j - 1 ) ) ) + D ( f ( t i , t j , x m - 1 ( t j ) ) , f ( t i , t j , x m - 1 ¯ ( t j ) ) ) ] ≤ LT 2 4 n + α T 2 n · ∑ j = 1 n [ D ( x m - 1 ( t j - 1 ) , x m - 1 ¯ ( t j - 1 ) ) + D ( x m - 1 ( t j ) , x m - 1 ¯ ( t j ) ) ](22) for all i = 0 , n ¯ and m ∈ ℕ ∗ . Now, from ((22)) for m = 2 we obtain D ( x 2 ( t i ) , x ¯ 2 ( t i ) ) ≤ LT 2 4 n + α T 2 n · ∑ j = 1 n 2 LT 2 4 n = [ 1 + α T ] · LT 2 4 n , ∀ i = 0 , n ¯ . By induction, for m ∈ ℕ ∗ , m ≥ 3, we get D ( x m ( t i ) , x ¯ m ( t i ) ) ≤ [ 1 + α T + . . . + ( α T ) m - 1 ] · LT 2 4 n = 1 - ( α T ) m 1 - α T · LT 2 4 n ≤ LT 2 4 n [ 1 - α T ] , ∀ i = 0 , n ¯ , m ∈ ℕ ∗(23) and together with ((21)) we obtain the inequality ((20)).

In order to confirm the theoretical result (Theorem 10) we test the above presented algorithm on the following examples.

Example 1. Consider the nonlinear fuzzy Volterra integral equation x ( t ) = g ( t ) + ( FR ) ∫ 0 t H ( t , s ) · f ( s , x ( s ) ) ds , t ∈ [ 0 , 0.4 ] where H (t, s) = 2e^-t, f (s, x (s)) = [x (s)] ², (t, s) ∈ [0, 0.4] × [0, 0.4], and g ( t , r ) = [ g - r ( t ) , g + r ( t ) ] , t ∈ [0, 0.4] , r ∈ [0, 1] with g - r ( t ) = e - t - 2 te - t · 0.04 ( 1 - r ) 2 - 0.2 ( 1 - r ) + 0.08 ( 1 - r ) ( 1 - e - t ) , g + r ( t ) = e - t - 2 te - t · 0.04 ( 1 - r ) 2 + 0.2 ( 1 - r ) - 0.08 ( 1 - r ) ( 1 - e - t ) . Here, the fuzzy power with natural exponent is given by [ x ( s ) ] 2 = [ x - r ( s ) · x - r ( s ) , x + r ( s ) · x + r ( s ) ] , the pro-duct being the same as in [9]. The exact solution isx^∗ (t, r) = [(x^∗) ^- (t, r) , (x^∗) ⁺ (t, r)] = [e ^t  - 0.2 (1 - r) , e ^t  + 0.2 (1 - r)] ,  t ∈ [0, 0.4] , r ∈ [0, 1] . For n = 10, with ɛ = 10^-24 in the stop-ping criterion (Steps 3-4), we get m = 16 iterations and the results are presented in Tables 1 and 2, where we are focused on the errors e i - = | ( x ∗ ) - ( t i , r ) - x ¯ m - ( t i , r ) | , e i + = | ( x ∗ ) + ( t i , r ) - x ¯ m + ( t i , r ) | , i = 0 , n ¯, r ∈ {0.25, 0.5, 0.75} . In order to test the convergence, we consider n = 50, and the results e _i  = | (x^∗) ^- (t _i , 1) - x ¯ m - ( t i , 1 ) | = | ( x ∗ ) + ( t i , 1 ) - x ¯ m + ( t i , 1 ) | , for r = 1 and n = 10, n = 50 are presented in Table 3.

Example 2. Now we consider a linear example of fuzzy Volterra integral equation x ( t ) = g ( t ) + ( FR ) ∫ 0 t H ( t , s ) · x ( s ) ds , t ∈ [ 0 , 0.4 ] with H (t, s) = t cos(s - t) , and g ( t , r ) = [ g - r ( t ) , g + r ( t ) ] , t ∈ [0, 0.4] , r ∈ [0, 1], where g - r ( t ) = 2 t ( r 5 + 2 r ) ( 3 - 3 cos t - t 2 ) g + r ( t ) = 6 t ( 2 - r 3 ) ( 3 - 3 cos t - t 2 ) and the exact solution is x^∗ (t, r) = [(x^∗) ^- (t, r) , (x^∗) ⁺ (t, r)] = [t³ (r⁵ +   2r) , t³ (6 - 3r³)] ,  t ∈ [0, 0.4] , r ∈ [0, 1] . For n = 10, with ɛ = 10^-24 in the stopping criterion (Step 3-4), we get m = 7 iterations and the results are presented in Tables 4 and 5, where we are focused on the errors e i - = | ( x ∗ ) - ( t i , r ) - x ¯ m - ( t i , r ) | , e i + = | ( x ∗ ) + ( t i , r ) - x ¯ m + ( t i , r ) | , i = 0 , n ¯, r ∈ {0.25, 0.5, 0.75} . In order to test the convergence, we consider n = 50, and the results e i = | ( x ∗ ) - ( t i , 1 ) - x ¯ m - ( t i , 1 ) | = | ( x ∗ ) + ( t i , 1 ) - x ¯ m + (t _i , 1) |, for r = 1 and n = 10, n = 50 are presented in Table 6.