Numerical solution of neutral functional-differential equations with proportional delays

Article History: Received 22 June 2016 Accepted 22 March 2017 Available 15 July 2017 In this paper, homotopy analysis method is improved with optimal determination of auxiliary parameter by use of residual error function for solving neutral functional-differential equations (NFDEs) with proportional delays. Convergence analysis and error estimate of method are given. Some numerical examples are solved and comparisons are made with the existing results. The numerical results show that the homotopy analysis method with residual error function is very effective and simple.


Introduction
Neutral functional-differential equations with proportional delays represent a specific form of delay differential equations.Such equations arise in various fields of science and engineering and play a significant role in the mathematical modeling of real-world phenomena [1].Clearly, most of these equations cannot be solved with well-known exact methods.For this reason, it is necessary to design efficient numerical treatment to approximate solutions.Ishiwata et al. used the rational approximation method and the collocation method to obtain numerical solutions of NFDEs with proportional delays [2,3].Hu et al. applied linear multi step methods to obtain numerical solutions for NFDEs [4].Wang et al. obtained approximate solutions for NFDEs by continuous Runge-Kutta methods and one-leg-θ method [5][6][7].Chen and Wang applied the variational iteration method for solving NFDEs with proportional delays [8].Biazar and Ghanbari applied the homotopy perturbation method to obtain numerical solution of NFDEs with proportional delays [9] and so on [10,11].
Homotopy analysis method (HAM) is proposed by Liao [12,13].This method has been successfully employed to handle a wide variety of scientific and engineering applications.Alomari et al. used modified HAM for solution of delay differential equation in [14].Kumar and Rashidi applied fractional homotopy analysis transform method to obtain approximate analytical solution of nonlinear homogeneous and nonhomogeneous time-fractional gas dynamics equations in [15].Abbasbandy employed HAM to find a family of travelling-wave solutions of the Kawahara equation in [16].Jafari and Seifi used HAM for solution of linear and nonlinear fractional diffusionwave equation in [17].Sakar and Erdogan applied the HAM and Adomian's decomposition method for solving the time-fractional Fornberg-Whitham equation in [18].HAM is different from the perturbation methods it provides the convenient way to control and adjust the convergence region and convergence rate of the series solution.However, for some type of auxiliary operator, in other words some type of base functions, it is generally timeconsuming to get high order approximation, and very large terms appearing in high order approximation [19].The homotopy-series solution is not only dependent upon t but also the convergencecontrol parameter .Convergence-control parameter which supplies us a convenient way to guarantee the convergence of homotopy-series solution.When the approximations contain unknown convergence-control parameters and other physical parameters, it is time-consuming to compute the square residual error at high-order of approximations.To avoid the time-consuming calculation, we will employ averaged square residual error function [19,20].The aim of this paper is to extend the homotopy analysis method with residual error function to obtain the numerical solution of the following neutral functional-differential equations with proportional delays [8,9], with the initial conditions Here, a(t) and b k (t), (k = 0, 1, ..., n − 1) are given analytical functions, and β, p k , c ik , λ i denote given constant with 0 < p k < 1, (k = 0, 1, ..., n).This paper is organized as follows: In Section 2, homotopy analysis method with residual error function is presented.Section 3 is devoted to the convergence analysis of the method.Section 4 contains numerical comparisons between the results obtained by the homotopy analysis method in this work and some existing methods.Finally, concluding remarks are given in the last section.

Homotopy analysis method with residual error function
We consider the following nonlinear differential equation where, N is a nonlinear differential operator, t denotes independent variable, u (t) is an unknown function.By means of generalizing the traditional homotopy method, Liao [12] constructs the socalled zero-order deformation equation here, p ∈ [0, 1] is the embedding parameter, = 0 is a non-zero auxiliary parameter, H (t) = 0 is non-zero auxiliary function, L is an auxiliary linear operator.u 0 (t) is the initial guess of u (t) and φ (t; p) is unknown function, respectively.It is important that one has great freedom to choose auxiliary things such as and L in homotopy analysis method.Obviously, when p = 0 and p = 1 , it holds respectively.Thus, as p increases from 0 to 1, the solution φ (t; p) varies from the initial guesses u 0 (t) to the solution u (x, t).Expanding φ (t; p) in Taylor series with respect to p, we have where If the auxiliary linear operator, the initial guess, the auxiliary operator , and the auxiliary functions are so properly chosen, then the series Eq.( 6) converges at p = 1 and which must be one of the solutions of the original nonlinear equations, as proved by Liao [12].According to Eq.( 7), the governing equations can be deduced from zeroth-order deformation equation Eq.( 4).Define the vector Differentiating Eq.( 4) m-times with respect to the embedding parameter p and then setting p = 0 and finally dividing them by m!, we obtain the mth-order deformation equation, with the assumption (10) where, and The mth-order approximation series solution is given as It is clear from Eq.( 11) that u (t) contains unknown convergence-control parameter which determine the convergence region and rate of the homotopy-series solution.

Selection of optimal value of with residual error function
As given by Liao [19], at the mth-order of approximation, one can define the exact square residual error as, Here, ∆ m contains unknown convergencecontrol parameter.For mth-order approximation, optimal values of the convergence-control parameter is given by the minimum of ∆ m , namely However, it is proven by Liao [19] that the exact residual error ∆ m defined by Eq.( 13) needs too much CPU time to calculate even if the order of approximation is not very high.Thus, to greatly decrease the CPU time, we use here the so-called averaged square residual error √ E m defined by

Convergence analysis and error estimate
In this section we present convergence analysis and error estimate for our method.
The proofs of Theorem 1. and Theorem 2. can be found in [12].
Theorem 4. Assume that the series solution ∞ n=0 u n (t) defined in Eq. (11), is convergent to the solution u(t).If the truncated series m n=0 u n (t) is used as an approximation to the solution u(t) of Eq.( 3), then the maximum absolute truncated error is estimated as, Proof.From Theorem 3. and Eq.( 16), we have for n ≥ m.Now, as n → ∞ then S n → u (t).So, Since 0 < γ < 1, we have 1 − γ n−m < 1. Herewith the above inequality becomes, This completes the proof.

Numerical examples
Now, we apply the homotopy analysis method with residual error function which presented in Section 2-3 to some NFDEs with proportional delay.
We select auxiliary linear operator, with property in which c 1 is an integral constant to be determined by initial condition Eq. (20).
Furthermore, Eq.( 19) suggest that we define a nonlinear operator as From mth-order deformation equation we can obtain following components, Then the solution is, where n is the number of terms.
We define the following residual error function, for obtaining optimal value of .Figure 1(b) shows averaged square residual error function for the 7th-order approximation, i.e., with respect to for t ∈ [0,1].
To determine the region of validity of the convergence-control parameter , we plot the values of u ′ (0.2), and u ′′ (0.2) in Figure 1(a).It appears that should at least lie within the interval [−2, −1].For the best possible value within this region, the averaged square residual error at the 7th-order approximation was evaluated from Eq.( 26) which gives rise to the optimal value of = -1.18613,resulting in a averaged square residual error 4.27 × 10 −4 .For 10th-order approximation we found = −1.12058and averaged square residual error 5.20 × 10 −5 .The accuracy is improved by optimal choice of .In Table 1, we compare the absolute errors of the homotopy analysis method (n = 7 and n = 10) with those of the Runge-Kutta method (R-K) of [1,8], variational iteration method (VIM) [8] with n i = 7 and the one-leg θ method [5,6] with θ = 0.8, where h = 0.01 and homotopy perturbation method (HPM) [9] with n = 7.  Example 2. Consider the second-order NFDE with proportional delay [8,9]: with initial conditions, which has the exact solution u (t) = t 2 [8].
We select auxiliary linear operator, with property in which c 1 and c 2 are integral constants to be determined by initial condition Eq.( 28).Furthermore, Eq.( 27) suggest that we define a nonlinear operator as From mth-order deformation equation, we can obtain following components, . . .We define the following residual function, for obtaining optimal value of .Figure 2(b) shows averaged square residual error function for the 6th-order approximation, i.e., with respect to for t ∈ [0,1].
To determine the region of validity of the convergence-control parameter , we plot the values of u ′ (0.1), and u ′′ (0.1), in Figure 2(a).It appears that should at least lie within the interval [−2, −1].For the best possible value within this region, the averaged square residual error at the 6th-order approximation was evaluated from Eq.( 33) which gives rise to the optimal value of = -1.49346,resulting in a averaged square residual error 8.16 × 10 −4 .For 10th-order approximation we found = −1.45885and averaged square residual error 5.93 × 10 −6 .In Table 2, we compare the absolute errors of the homotopy analysis method (n = 6 and n = 10) with those of the Runge-Kutta method of [1,8], variational iteration method [8] with n i = 6 and the one-leg θ method [5,6] with θ = 0.8, where h = 0.01 and homotopy perturbation method [9] with n = 6.
Example 3. Consider the third-order neutral functional differential equation with proportional delay [8,9]: with initial conditions, which has the exact solution u (t) = t 4 [8].
We select auxiliary linear operator, with property for obtaining optimal value of .To determine the region of validity of the convergence-control parameter , we plot the values of u ′ (0.2), and u ′′ (0.2) in Figure 3(a).It appears that should at least lie within the interval [−1.5, −0.5].For the best possible value within this region, the averaged square residual error at the 4th-order approximation was evaluated from Eq.( 40) which gives rise to the optimal value of = -1.0932155,resulting in a averaged square residual error 1.77 × 10 −5 .For 7th-order approximation we found = −1.08382and averaged square residual error 2.05×10 −8 .In Table 3, we compare the absolute errors of the homotopy analysis method (n = 4 and n = 7) with those of the Runge-Kutta method of [1,8], variational iteration method [8] with n i = 4 and the one-leg θ method [5,6] with θ = 0.8, where h = 0.01 and homotopy perturbation method [9] with n = 5.

Conclusion
In this paper, we have demonstrated the suitability of the homotopy analysis method with residual error function for solving neutral functionaldifferential equations with proportional delays.We obtain high-accuracy approximate solutions after only a few iterations.The numerical results also show that the HAM with residual error function is more effective than Runge-Kutta method, HPM, VIM and other methods for solving NFDEs with proportional delays.

Table 1 .
Comparison of absolute errors for Example 1.

Table 2 .
Comparison of absolute errors for Example 2.

Table 3 .
Comparison of absolute errors for Example 3.