In recent years, there has been a growing interest in the nonlinear Volterra–Fredholm integrodifferential equations, NVFID, in short, which are a combination of differential and Fredholm-Volterra integral equations. These equations occur frequently as a model in mathematical biology and the physical sciences. Volterra–Fredholm integrodifferential equations is sometimes difficult to solve analytically. Therefore, finding either the analytical approximation or numerical solution of such equations is of great interest. The interested reader can refer to [1–11] (and the references therein) for more research works.

We aim in this work to study the convergence of the homotopy analysis method, for solving NVFID equations, of the form: (1)u′(x)+q(x)u(x)+λ1∫01k1(x,t)(u(t))pdt   +λ2∫0xk2(x,t)(u(t))qdt=f(x), p,q∈N,u(0)=u0, x∈[0,1], where μ is a real number, the kernel k ₁(x, t), k ₂(x, t) are known continuous functions over [a, b] × [a, b], and f(x), q(x) are given continuous function defined over [a, b].

The homotopy analysis method (HAM) [12] has been proved to be one of the useful techniques to solve numerous linear and nonlinear functional equations [13–17].

As mentioned in [13, 14], unlike all previous analytic techniques, the homotopy analysis method provides great freedom to express solutions of a given nonlinear problem by means of different base functions. Also this method provides a way to adjust and control the convergence region and the rate of convergence of solution series, by introducing the auxiliary parameter ℏ.

By properly choosing the base functions, initial approximations, auxiliary linear operators, and auxiliary parameter ℏ, HAM gives rapidly convergent successive approximations of the exact solution. A systematic description of this analytic technique, for nonlinear problems, can be found in [13].

This paper is organized as follows. In Section 2, a short description of the basic ideas of the homotopy analysis method will be stated. In Section 3, the homotopy analysis method is applied to construct approximate solution of (1). Section 4 is devoted to the convergence analysis of the method. In Section 5, our numerical findings are reported and demonstrate the accuracy of the proposed scheme, by considering three numerical examples. Finally, conclusions are stated in the last section.

Let us assume the following nonlinear differential: (2)N[u(τ)]=0, where N is a nonlinear operator, τ is an independent variable, and u(τ) is the solution of the equation.

By utilizing the concept of homotopy in topology, Liao [13] constructs the so-called zero-order deformation equation (3)(1−q)L[Φ(x;q)−u0(x)]−ℏqN[Φ(x;q)]=0, where q ∈ [0,1] is the embedding parameter, ℏ ≠ 0 is a nonzero auxiliary parameter, L is an auxiliary linear operator, u ₀(x) is an initial guess of u(x), and Φ(x; q) is an unknown function, respectively.

In the view of HAM the solution of original equation is assumed to be expanded in terms of embedding parameter q as (4)Φ(x;q)=u0(x)+∑i=0∞ui(x)qi, where u _i(x)'s are obtained as follows.

Appling recently proposed “mth-order homotopy-derivative operator” [14]: (5)Dm(ϕ)=1m!∂mϕ∂qm|q=0 to both sides of (3), one reads (6)L[um(x)−χmum−1(x)]=hRm−1(x), subject to initial condition (7)um(0)=0, where (8)Rm−1(x)=Dm−1(N[u(x;q)])=1(m−1)!∂m−1N[u(x;q)]∂qm−1|q=0,χm={1,m>1,0,m≤1. In this way u _m's, for m ≥ 1, are dependent upon not only x but also the auxiliary parameter ℏ.

As proved in [13], whenever (4) is a convergent series at q = 1, its limit must satisfy the original equation N[u(τ)] = 0. To find a proper convergence-control parameter ℏ, to get a convergent series solution or to accelerate its convergence rate, there is a classic way of plotting the so-called ‘‘ℏ-curves” or ‘‘curves for convergence-control parameter”. Such a region can be found, although approximately, by plotting the curves of these unknown quantities versus ℏ.

However, it is a pity that curves for convergence-control parameter (i.e., ℏ-curves) give us only a graphical region and cannot tell us which value of ℏ ₀ ∈ R _ℏ gives the fastest convergent series. Moreover, recently in [18] a misinterpreted usage of ℏ-curves has been reported.

To investigate the exact solution of (1), by means of the homotopy analysis method, let us consider the nonlinear operator (9)N[Φ(x;q)]=∂Φ(x;q)∂x+q(x)Φ(x;q)+λ1∫01k1(x,t)(Φ(x;q))pdt  +λ2∫0xk2(x,t)(Φ(x;q))qdt−f(x), and linear operator (10)Lϕ=dϕdx, with the following property: (11)L[c]=0, where c is constant.

The initial guess u ₀(x) is chosen such that it satisfies the initial condition of problem; that is, (12)u0(0)=u(0)=u0. from (8) and by using (9), we have (13)Rm[ϕ→m−1(x)]=∂ϕm−1∂x+q(x)ϕm−1−(1−χm)f(x) +λ1∫01k1(x,t)[∑r1=0m−1ϕm−1−r1(x)×∑r2=0r1ϕr1−r2(x)∑r3=0r2ϕr2−r3(x)⋯∑rp−2=0rp−3ϕrp−3−rp−2(x)×∑rp−1=0rp−2ϕrp−2−rp−1(x)ϕrp−1(x)]dt +λ2∫0xk2(x,t)[∑r1=0m−1ϕm−1−r1(x)∑r2=0r1ϕr1−r2(x)×∑r3=0r2ϕr2−r3(x)⋯∑rp−2=0rp−3ϕrq−3−rp−2(x)×∑rp−1=0rp−2ϕrq−2−rp−1(x)ϕrq−1(x)]dt. From (13), it should be noted that the right-hand side of (6) is only dependent upon ϕ→m-1. Thus, ϕ ₀(r, t), ϕ ₁(r, t), ϕ ₂(r, t),…, will be computed recursively, by means of solving the linear high-order deformation equation (6) subject to (7). So, the mth-order approximation of u(x) is given by (14)um(x)=∑k=0mϕk(x).

In this section, some theorems and conditions in the framework of convergence of the homotopy analysis method are stated.

Whereas the series (15)ϕ0(x)+∑m=1+∞ϕm(x) converges, where ϕ _m(x)'s are resulted from (6), (7), and (13), the limit of the series is an exact solution of (1).

Since, by hypothesis, the series is convergent, it holds (16)s(x)=∑m=0+∞ϕm(x). So, the necessary condition for the convergence of the series is valid; that is, (17)lim⁡m→∞ϕm(x)=0. Using (6) and (17), we have (18)ℏH(x)∑m=1+∞Rm[ϕ→m−1(x)] =lim⁡n→∞∑m=1nL[ϕm(x)−χmϕm−1(x)] =L{lim⁡n→∞∑m=1n[ϕm(x)−χmϕm−1(x)]} =L{lim⁡n→∞ϕn(x)} =0. Since ℏ ≠ 0, we must have (19)∑m=1+∞Rm[ϕ→m−1(x)]=0. On the other hand, we have (20)∑m=1+∞Rm[ϕ→m−1(x)]   =∑m=1n[∂ϕm−1∂x+q(x)ϕm−1−(1−χm)f(x)  +λ1  ∫01k1(x,t)[∑r1=0m−1ϕm−1−r1(x)∑r2=0r1ϕr1−r2(x)×∑r3=0r2ϕr2−r3(x)…∑rp−2=0rp−3ϕrp−3−rp−2(x)×∑rp−1=0rp−2ϕrp−2−rp−1(x)×ϕrp−1(x)]dt+λ2∫0xk2(x,t)[∑r1=0m−1ϕm−1−r1(x)∑r2=0r1ϕr1−r2(x)×∑r3=0r2ϕr2−r3(x)…∑rp−2=0rp−3ϕrq−3−rp−2(x)×∑rp−1=0rp−2ϕrq−2−rp−1(x)×ϕrq−1(x)]dt] =∑m=0∞ϕm′(x)+q(x)∑m=0∞ϕm(x)−f(x)  +λ1∑m=1∞{∫01k1(x,t)×[∑r1=0m−1ϕm−1−r1(x)∑r2=0r1ϕr1−r2(x)×∑r3=0r2ϕr2−r3(x)…∑rp−2=0rp−3ϕrp−3−rp−2(x)×∑rp−1=0rp−2ϕrp−2−rp−1(x)×ϕrp−1(x)]dt}  +λ2∑m=1∞{∫0xk2(x,t)[∑r1=0m−1ϕm−1−r1(x)∑r2=0r1ϕr1−r2(x)×∑r3=0r2ϕr2−r3(x)…∑rq−2=0rq−3ϕrq−3−rp−2(x)×∑rq−1=0rq−2ϕrq−2−rq−1×(x)ϕrq−1(x)]dt} =∑m=0∞ϕm′(x)+q(x)∑m=0∞ϕm(x)−f(x)  +λ1∫01k1(x,t)[∑rp−1=0∞ϕm−1−r1(x)∑rp−2=rp−1∞ϕr1−r2(x)×∑rp−2=rp−1∞ϕr2−r3(x)…∑r2=r3∞ϕrp−3−rp−2(x)×∑r1=r2∞ϕrp−2−rp−1(x)×∑m=r1∞ϕm−r1(x)]dt  +λ2∫0xk2(x,t)[∑rq−1=0∞ϕm−1−r1(x)∑rq−2=rq−1∞ϕr1−r2(x)×∑rq−2=rq−1∞ϕr2−r3(x)…∑r2=r3∞ϕrp−3−rq−2(x)×∑r1=r2∞ϕrq−2−rq−1(x)∑m=r1∞ϕm−r1(x)]dt =∑m=0∞ϕm′(x)+q(x)∑m=0∞ϕm(x)−f(x)    +λ1∫  0  1k1(x,t)[∑i1=0∞ϕi1(x)∑i2=0∞ϕi2(x)∑i3=0∞ϕi3(x)…∑ip−2=0∞ϕip−2(x)×∑ip−1=0∞ϕip−1(x)∑ip=0∞ϕip(x)]dt  +λ2  ∫  0xk2(x,t)[∑j1=0∞ϕj1(x)∑j2=0∞ϕj2(x)×∑j3=0∞ϕj3(x)…∑jq−2=0∞ϕjq−2(x)×∑jq−1=0∞ϕjq−1(x)∑jq=0∞ϕjq(x)]dt =s′(x)+q(x)s(x)−f(x)  +  λ1∫  abk1(x,t)(s(t))pdt    +λ2∫0xk2(x,t)(s(t))qdt.  

So, from (19), we obtain (21)s′(x)+q(x)s(x)+λ1∫01k1(x,t)(s(x))pdt  +λ2∫0xk2(x,t)(s(x))qdt=f(x). Also, from initial conditions (7) and (12), the following holds: (22)s(0)=∑i=0∞ϕi(0)=ϕ0(0)=u0(0)=u0.

Since s(x) satisfies (21) and (22), we conclude that it is an exact solution of (1). This completes the proof.

Suppose that ℑ ⊂ ℝ is a Banach space with a suitable norm ||·||, say ||·||_∞, over which the sequence ϕ _m(x) is defined for a prescribed value of ℏ. Assume also that the initial approximation ϕ ₀(x) remains inside the domain of the solution u(x). Taking r ∈ R as a constant, the following statements holds.

If there exists some r ∈ [0,1], such that for all k ∈ ℕ we have ||ϕ _k+1(x)|| ≥ r||ϕ _k(x)||, then the series solution u(x) = ∑_k=0 ^∞ ϕ _k(x)q ^k converges absolutely at q = 1, over the domain of x.

Indeed, this is a special case of Banach' fix point theorem, which can be found in standard texts on real analysis such as in [19].

Theorems 1 and 2 state that the homotopy series solution ∑_k=0 ^∞ ϕ _k(x), of nonlinear problem (1), converges to an exact solution u(x), under the condition that ∃γ, 0 < γ < 1 such that ∀k ≥ k ₀, ||ϕ _k+1(x)|| ≤ γ||ϕ _k(x)||, for some k ₀ ∈ ℕ. In other words, if, for every i ≥ k ₀, we consider the parameters β _i's, as were defined in [20], (23)βi+1={||ϕi+1||||ϕi||,||ϕi||≠0,0,||ϕi||=0, for i ∈ ℕ ∪ {0}, then the series solution∑_k=0 ^∞ ϕ _k(x) of problem (1) converges to an exact solution u(x), when ∀i ≥ k ₀, 0 ≤ β _i < 1.

To illustrate the convergence study of HAM, two examples are presented.

Let us first consider the following linear Fredholm integral equation: (24)u′(x)+u(x)−14∫01t(u(t))3dt  +12∫0xx(u(t))2dt=f(x), where f(x) = 2x + x ² + (1/10)x ⁶ − (1/32), with the initial condition u(0) = 0. The exact solution is u(x) = x ² [6].

To solve (24) by means of HAM and according to the initial condition, we can choose the initial guess of u(x) as follows: (25)ϕ0=x The linear operator and the zero-order deformation equations are as defined by (10) and (3), respectively.

Therefore, the mth-order deformation equations (6) read as (26)ddx[um(x)−χmum−1(x)]=hRm−1(x), subject to initial condition (27)um(0)=0, where (28)Rm−1(x)=∂ϕm−1∂x+ϕm−1−(1−χm)f(x)−14∫01t[∑r1=0m−1ϕm−1−r1(x)∑r2=0r1ϕr1−r2(x)ϕr2(x)]dt+12∫0xx[∑r1=0m−1ϕm−1−r1(x)ϕr1(x)]dt. Now, we can successively obtain the solution of the mth-order deformation equation (26) and (27), for m ≥ 1.

Table 1 shows some values of β _i's, defined as in (23), which were obtained from the HAM series solution (14) by using different values of ℏ. Form Table 1, since β _i < 1 for ℏ = −0.74, ℏ = −0.9, and ℏ = −1, we can conclude that the HAM approach converges to the exact solution for ℏ = −0.74, ℏ = −0.9, and ℏ = −1. Also, one can observe that β _i's are not less than one for ℏ = 0.4. So, the HAM approach may be divergent for ℏ = 0.4.

In Table 2 the relative errors δ _n of the n terms approximation of HAM, which is defined as (29)δ(xj)=|uexact(xj)−uappr(xj)uexact(xj)|, for different values of ℏ, are presented. It is evident that the auxiliary parameter ℏ can also be effectively implemented to adjust and control the rate of convergence of the series solutions by HAM.

Consider now the following equation [7, 11]: (30)u′(x)−∫0xcos⁡(x−t)(u(t))2dt, =−2sin(x)−13cos⁡(x)−23cos⁡(2x), with initial condition u(0) = 1. The exact solution is u(x) = cos⁡(x) − sin(x).

To obtain the approximate solution of this example, if we choose ϕ ₀(x) = 1, then the solutions of mth-order deformation equations can be resulted in (6) and (7), where (31)Rm−1(x)=∂ϕm−1∂x+ϕm−1−(1−χm)×[−2sin(x)−13cos⁡(x)−23cos⁡(2x)]−∫01cos⁡(x−t)[∑r1=0m−1ϕm−1−r1(x)ϕr1(x)]dt. Since HAM provides a family of solution expression in the auxiliary parameter h, the convergence region of solution series depends upon the value of ℏ.

Table 3 shows the computed values of β _i's for different values of ℏ, in this example. Also, for different values of ℏ, the relative errors of approximate analytic solutions by HAM solution are reported in Table 4.

Clearly, one can observe that the approximate solution for ℏ = −1, is more accurate than the approximate solutions obtained for ℏ = −1.1, ℏ = −0.9, and ℏ = −0.8. Therefore, it can be claimed that the auxiliary parameter ℏ plays an important role to adjust and control the convergence of the series solution. It can be seen from Tables 3 and 4 that, for smaller values of β _i, our approximate solutions are in very excellent agreement with the exact solutions.

In this paper, we presented the application of the homotopy analysis method (HAM) for solving a special form of nonlinear Volterra-Fredholm integro-differential equation.

The sufficient condition for the convergence of the method is illustrated and then verified for two examples.

As we can see in Tables 2 and 4, HAM solutions have a good agreement with the numerical results provided that appropriate values for the convergence control parameter ℏ are chosen.

The ability of the HAM is mainly due to the fact that the method provides a way to ensure the convergence of series solution.

Convergence conditions in applying HAM for other equations, and systems of ordinary and partial differential equations, integral equations, and integro-differential equations are also under study in our research team.