Solving two-dimensional nonlinear fuzzy Volterra integral equations by homotopy analysis method

1 Introduction

Mathematical models of many physical phenomena and engineering problems are described by means of differential, integral and integro-differential equations. The topic of fuzzy integral equations is one of the important branches of fuzzy analysis theory and plays a major role in numerical analysis. Many mathematical models used in biology, chemistry, physics and engineering are based on integral equations.

One of the first applications of fuzzy integration was given by Wu and Ma [1] who investigated the fuzzy Fredholm integral equation of the second kind. In recent years, many mathematicians have studied a solution to fuzzy integral equations by numerical and analytical methods [2,3, 4,5,6, 7,8].

In 1992, Liao [9] employed the basic idea of the homotopy in topology to propose a general analytic method for nonlinear problems, namely homotopy analysis method (HAM) [10,11,12]. The method is based on the concept of creating function series. If the series converges, its sum is the solution of this system of equations. This method has been successfully applied to solve many types of nonlinear problems [13,14, 15,16]. In [17,18], the HAM is applied to solve two-dimensional linear Volterra fuzzy integral equations and mixed Volterra-Fredholm fuzzy integral equations.

In this paper, we present an application of the HAM for solving the nonlinear fuzzy Volterra integral equations in two variables.

The paper is organized as follows: in Section 2, we give the basic notations of fuzzy numbers, fuzzy functions and fuzzy integrals. In Section 3, we present the basic idea of the HAM. In Section 4, the two-dimensional nonlinear fuzzy Volterra integral equation (2D-NFVIE) and its parametric form are discussed and then we apply HAM to it. In Section 5, we obtain sufficient conditions for convergence of the proposed method and error estimate. In Section 6, we give an example. In Section 7, the conclusion is drawn.

2 Basic concepts

First, we present some notions and results about fuzzy numbers and fuzzy-number-valued functions, and fuzzy integrals.

In [20], the convexity of the fuzzy numbers u ( r x + ( 1 − r ) y ) ≥ min { u ( x ) , u ( y ) } , for any x , y ∈ R , r ∈ [ 0 , 1 ] is shown.

We denote with E 1 the set of all fuzzy numbers and R + = ( 0 , ∞ ) .

The set [ u ] r = { x ∈ R : u ( x ) ≥ r } is called r-level set of the fuzzy number u. Then lower and upper representation of fuzzy number is [ u ] r = [ u ̲ ( r ) , u ¯ ( r ) ] for all r ∈ [ 0 , 1 ] , where u ̲ , u ¯ : [ 0 , 1 ] → R are functions, such that u ̲ is increasing and u ¯ is decreasing.

Let u , v ∈ E 1 , k ∈ R . The addition and the scalar multiplication are defined by

and

With respect to ⊕ in E 1 the neutral element is denoted by 0 ˜ = χ { 0 } .

In [21], the algebraic properties of fuzzy numbers are given. We use the Hausdorff metric to define a distance between fuzzy numbers.

For any fuzzy-number-valued function f : A = [ a , b ] × [ c , d ] → E 1 we define the left and right r -level functions, respectively, f ̲ ( . , . , r ) , f ¯ ( . , . , r ) : A → R for each ( s , t ) ∈ A and r ∈ [ 0 , 1 ] .

Define the metric D ∗ ( f , g ) = max ( s , t ) ∈ A D ( f ( s , t ) , g ( s , t ) ) on the set

We see that ( C ( A , E 1 ) , D ∗ ) is a complete metric space.

In [23], the notion of Henstock integral for fuzzy-number-valued functions is defined as follows.

Let f : A → E 1 , for Δ x n : a = x 0 < x 1 < ⋯ < x n = b , Δ y n : c = y 0 < y 1 < ⋯ < y n = d , be two partitions of the intervals [ a , b ] and [ c , d ] , respectively. Let one consider the intermediate points ξ i ∈ [ x i − 1 , x i ] and η j ∈ [ y j − 1 , y j ] , i = 1 , … , n ; j = 1 , … , n , and δ : [ a , b ] → R + and σ : [ c , d ] → R + . The divisions P x = ( [ x i − 1 , x i ] ; ξ i ) , i = 1 , … , n and P y = ( [ y j − 1 , y j ] ; η j ) , j = 1 , … , n are said to be δ -fine and σ -fine, respectively, if [ x i − 1 , x i ] ⊆ ( ξ i − δ ( ξ i ) , ξ i + δ ( ξ i ) ) and [ y j − 1 , y j ] ⊆ ( η j − σ ( η j ) , η j + σ ( η j ) ) .

Then, I is called the two-dimensional Henstock integral of f and is denoted by

If the above δ and σ are constant functions, then one recaptures the concept of Riemann integral. In this case, I ∈ E 1 will be called two-dimensional integral of f on A and will be denoted by ( F R ) ∫ c d ( F R ) ∫ a b f ( s , t ) d s d t .

Also, if f is continuous, then f ̲ ( . , . , r ) and f ¯ ( . , . , r ) are continuous for any r ∈ [ 0 , 1 ] , and consequently, they are Henstock integrable.

3 Basic idea of HAM

According to [6], we will give a brief overview of the main used method. HAM transforms the considered equation into the corresponding deformation equation. Using this method we solve the operator equation

where N is the nonlinear operator, u is the unknown function and Ω is the domain of z . Define the homotopy operator ℋ in the following way:

where p ∈ [ 0 , 1 ] is an embedding parameter, h ≠ 0 is the convergence control parameter [24,25] and u 0 is the initial approximation of the solution of equation (1). The linear operator L is the auxiliary with property L ( 0 ) = 0 . This operator can be arbitrarily selected. The practice is to choose L so that the equations, obtained in the next stages of the procedure, would be as simple to solve as possible.

From the equation ℋ ( Φ ( z ; p ) ) = 0 , we get the so-called zero-order deformation equation:

For p = 0 we have Φ ( z ; 0 ) = u 0 ( z ) and for p = 1 we get the searched solution of equation (1). By expanding function Φ ( z ; p ) : Ω × [ 0 , 1 ] → R into the Maclaurin series with respect to parameter p we have

where u 0 ( z ) = Φ ( z ; 0 ) and

If series (4) converges for p = 1 , then we get the sought solution

In order to determine the function u m we differentiate m -times, with respect to parameter p , the left and right-hand side of formula (3), then the obtained result is divided by m ! and substituted with p = 0 , which gives the so-called m th-order deformation equation ( m ≥ 1 )

where u ˜ m − 1 ( z ) = { u 0 ( z ) , u 1 ( z ) , … , u m − 1 ( z ) } .

Applying L − 1 on both sides of equation (7), we get

where

and

If we are not able to determine the sum of series in (6), then as the approximate solution of considered equation we accept the partial sum of this series

By appropriate selection of the convergence control parameter h , we can influence the convergence region of series and the rate of this convergence [15,26]. One of the methods for selecting the value of convergence control parameter is the so-called h -curve. To obtain this curve we need to investigate the behavior of a certain quantity of the exact solution as a function of parameter h [27]. Another method is the so-called “optimization method” proposed in the study by Liao [15]. In this method, we define the squared residual of the governing equation:

The optimum value of the convergence control parameter is obtained by finding the minimum of this squared residual, whereas the effective region of the convergence control parameter h is defined as

To speed up the calculations, Liao [15] suggested to replace the integral in formula 11 by its approximate value obtained by applying the quadrature rules. By choosing a different value of the control parameter h than the optimal one, but still belonging to the effective region, we also obtain the convergent series, only the rate of convergence is lower. A version of the method with the above described selection of optimal value of the convergence control parameter is called the basic optimal HAM [15].

4 Using of HAM

We consider 2D-NFVIE of the following form:

where g , u : A = [ a , b ] × [ c , d ] → E 1 are continuous fuzzy-number valued functions, k : A × A → R + is continuous function and G : E 1 → E 1 is continuous function on E 1 .

According to [2], we introduce the parametric form of the integral equation (13). Let u ( s , t , r ) = ( u ̲ ( s , t , r ) , u ¯ ( s , t , r ) ) , g ( s , t , r ) = ( g ̲ ( s , t , r ) , g ¯ ( s , t , r ) ) , 0 ≤ r ≤ 1 and ( s , t ) ∈ A be the parametric form of the functions u ( s , t ) and g ( s , t ) .

So, the parametric form of equation (13) is as follows:

Let for ( x , y ) ∈ A , we have

Then,

for a ≤ x ≤ s ≤ b , c ≤ y ≤ t ≤ d and 0 ≤ r ≤ 1 .

Let the function G ( β ) be increasing for β ∈ [ u ̲ ( x , y , r ) , u ¯ ( x , y , r ) ] and k ( s , t , x , y ) ≥ 0 for all a ≤ x ≤ s ≤ b , c ≤ y ≤ t ≤ d and 0 ≤ r ≤ 1 . Then the parametric form of equation (13) is as follows:

We consider the operators L and N can be defined in the following way:

Then we get the following formula for the function u ̲ m ( s , t , r )

where χ m is defined by (8).

For the operator R m , m ≥ 1 from (9) we obtain

By using (16) and the definitions of appropriate operators for m ≥ 1 we obtain

where u ̲ 0 ∈ C ( A × [ 0 , 1 ] , R ) and for m ≥ 2

5 Existence and convergence

In this section, we prove convergence of the presented method and obtain the error estimation between the exact and the approximate solution.

We introduce the following conditions:

1. g ∈ C ( A , E 1 ) , k ∈ C ( A × A , R + ) and G ∈ C ( E 1 , E 1 ) ;

2. there exists L G ≥ 0 such that D ( G ( u ) , G ( v ) ) ≤ L G D ( u , v ) for all u , v ∈ E 1 ;

3. α = M k L G Δ < 1 , where ∣ k ( s , t , x , y ) ∣ ≤ M k for all ( s , t ) , ( x , y ) ∈ A , according to the continuity of k and Δ = ( b − a ) ( d − c ) .