On a symbolic method for neutral functional-differential equations with proportional delays

Abstract In this paper, we present a new symbolic method for solving neutral functional-differential equations with proportional delays having variable coefficients via homotopy perturbation method. The proposed symbolic method is also applicable to solve the multi-pantograph equations with variable coefficients. Several numerical examples are given to illustrate the proposed method and the results show that the proposed symbolic method is very effective.


Introduction
Functional-differential equations are playing an important role in the mathematical modeling of realworld phenomena and the multi-pantograph equations arise in many applications, such as number theory, probability theory on algebraic structures, electrodynamics, astrophysics, non-linear dynamical systems, quantum mechanics,and cell growth, etc., these equations cannot be solved by one of the known standard methods. Functional-differential equations find use in mathematical models that assume a specified behavior or phenomenon depends on the present as well as the past state of a system. In other words, past events explicitly influence future results. For this reason, Functionaldifferential equations are used to in many applications rather than ordinary differential equations, in which future behavior only implicitly depends on the past. Properties of the analytic solution of these Srinivasarao Thota ABOUT THE AUTHOR Srinivasarao Thota (corresponding author) completed his M.Sc. in Mathematics from Indian Institute of Technology (IIT) Madras, India and Ph.D. in Mathematics from Motilal Nehru National institute of Technology (NIT) Allahabad, India. His research interests are Computer Algebra (Symbolic methods for differential equations), Numerical Analysis (Root finding algorithms) and Mathematical Modelling (Ecology). He has published more than 30 research papers in various international journals and attended and presented research work in several national and international conferences. Also, he had given more than 9 special/invited conferences talks. He has more than 10 years of teaching and research experience. Presently, he is working at Department of Applied Mathematics, School of Applied Natural Sciences, Adama Science and Technology University, Ethiopia.

PUBLIC INTEREST STATEMENT
Functional-differential equations (FDEs) and multi-pantograph equations (MPEs) are playing an important role in the mathematical modeling of real-world phenomena. FDEs find use in mathematical models that assume a specified behavior or phenomenon depends on the present as well as the past state of a system. In other words, past events explicitly influence future results. For this reason, FDEs are used to in many applications rather than ordinary differential equations, in which future behavior only implicitly depends on the past. FDEs and the MPEs arise in many applications, such as number theory, probability theory on algebraic structures, electrodynamics, astrophysics, non-linear dynamical systems, quantum mechanics, cell growth, etc.
In this paper, we focus on a new symbolic method to obtain an analytic or approximate solution of a given neutral functional-differential equations in easy and efficient manner. The main idea is that, we follow, similar to the symbolic methods for differential equations discussed in (Thota, 2019;Thota & Kumar, 2016, 2015Thota & Kumar, 2016a, 2016b, 2015, 2017Thota et al., 2012;Thota, 2017Thota, , 2019aThota, , 2018aThota, , 2018bThota, , 2018cThota, , 2018dThota, , 2019b to present the proposed method symbolically. Indeed, in this paper, we propose a new symbolic method for finding an analytic or approximate solution of neutral functional-differential equations as well as the multi-pantograph equations with variable coefficients by using the concept of homotopy perturbation method (HPM). The approach to analytic solution of this method is easy and computationally very attractive.
The HPM was introduced first by Dr. He in 1998(He, 1999He, 2000). This HPM is one of the a powerful and efficient techniques to find the solutions of the given non-linear equations without using a linearization process. This method is based on a combination of the perturbation and homotopy methods. HPM can take the advantages of the conventional perturbation method while eliminating its restrictions. Many engineers, scientists, and researchers have been applied successfully this method to solve several kinds of non-linear and linear equations in science and engineering (Biazar & Ghazvini, 2008;He, 1999He, , 2019He, , 2004He, 2000;He & Ain, 2020;He & Jin, 2020;Shang, 2012). For example, in (He, 2019), the authors discussed the simpler, the better: Analytical methods for nonlinear oscillators and fractional oscillators; in (Shang, 2012), authors presented a Lie algebra approach to susceptible-infected-susceptible epidemics, Lie algebraic discussion for affinity-based information diffusion in social networks, analytical solution for an in-host viral infection model with time-inhomogeneous rates; in (Biazar & Ghazvini, 2008), new promises and future challenges of fractal calculus: From two-scale thermodynamics to fractal variational principle are presented; and in (He & Jin, 2020) the authors focused a short review on analytical methods for the capillary oscillator in a nanoscale deformable tube.
The main objective of this paper is to present a new easy and efficient symbolic method to obtain an analytic or approximate solution of a given neutral functional-differential equation. In detail, we propose a symbolic method to find an analytic or approximate solution of a given neutral functional-differential equations as well as the multi-pantograph equations with variable coefficients by using the concept of HPM. In this paper, we concern with the neutral functionaldifferential equation with proportional delays (Jafari & Aminataei, 2012) of the following form where λ; c 0 ; c 1 ; . . . ; c nÀ 1 2 C; μ j ðxÞ and f ðxÞ are analytic functions; 0<q 1 < � � � <q nÀ 1 � 1, and the superscripts indicate differentiation. In order to propose a new method, we rewrite equation (1) as y ðnÞ ðxÞ ¼ λyðxÞ À ðβðxÞyðq n xÞÞ ðnÞ þ ∑ nÀ 1 j¼0 μ j ðxÞy ðjÞ ðq j xÞ þ f ðxÞ; x � 0; y ðiÞ ð0Þ ¼ c i ; i ¼ 0; 1; . . . ; n À 1: (2) The multi-pantograph equation is a kind of delay differential equations of the form Since we are interested to find an analytical or approximate solution of the given neutral functionaldifferential equations with proportional delays as well as the multi-pantograph equations with variable coefficients in symbolic manner, we use a powerful and efficient method namely HPM. This method is coupling of the homotopy method, a basic concept of topology, and the classical perturbation technique, provides us with a suitable way to obtain approximate or analytic solution for wide variety of problems arising in different fields. One of the advantages of this HPM is, good approximate solutions can be obtained with only a few iterations and HPM yields a very rapid convergence of the solution series in most cases, and another advantage of HPM is that the perturbation equation can be freely constructed in many ways by homotopy in topology and the initial approximation can also be freely selected. In Section 2 we describe the symbolic formulation of the proposed method. Several examples are given in Section 3 to illustrate the efficiency of the symbolic method, and comparisons are made to confirm the reliability of the symbolic method.

Basic ideas of homotopy perturbation method
In this section, we recall the basic ideas of the HPM, see, for example (He, 1999;He, 2000;Shakeri & Dehghan, 2008) for more details. Consider the following non-linear differential equation: with the initial conditions where L is a differential operator, I c is a initial condition operator, Ω is the set of initial conditions of the domain S, and f ðxÞ is a forcing analytical function. One can divide the differential operator L into two parts L l and L n , where L l is a linear operator and L n ia a non-linear operator. Now the equation (4) can be written as Using HPM, we can construct a homotopy vðx; pÞ : S � ½0; 1� ! R which satisfies, for p 2 ½0; 1�; x 2 S, Hðv; pÞ ¼ ð1 À pÞ½L l ðvÞ À L l ðy 0 Þ� þ p½LðvÞ À f ðxÞ� ¼ 0; ¼ L l ðvÞ À L l ðy 0 Þ þ pL l ðy 0 Þ þ p½L n ðvÞ À f ðxÞ� ¼ 0; where y 0 is an initial approximation and p is an embedding parameter. One can easily observe that Hðv; 0Þ ¼ L l ðvÞ À L l ðy 0 Þ ¼ 0 and Hðv; 1Þ ¼ LðvÞ À f ðxÞ ¼ 0. According to the HPM, we can take the parameter p as a small, and assume that the solution of (7) and (8) can be written as a power series in p v ¼ ∑ 1 n¼0 p n v n : Setting p ¼ 1, the approximate solution of non-linear differential equation (4) is In the following section, we present the proposed method.

Symbolic formulation
In this section, we illustrate the basic idea of the proposed method to find the solution of the neutral functional-differential equations corresponding to the special choice of S ¼ C 1 ½0; 1� for simplicity. Recall equation (2) as follows: LðyðxÞ þ βðxÞyðq n xÞÞ ¼ λyðxÞ þ ∑ nÀ 1 j¼0 μ j ðxÞy ðjÞ ðq j xÞ þ f ðxÞ; where L : S ! S is a linear differential operator, i.e., L ¼ D n ¼ d n dx n , y 2 S, and c 0 ; . . . ; c nÀ 1 are constants of C. Analytically speaking, the differential operator L acts on the Banach space (C½0; 1�; k �k 1 ) with dense domain of definition C n ½0; 1� and algebraically speaking, the domain of L is the complex vector space C 1 ½0; 1� without any prescribed topology.

OAMA A 1813961
From replacement lemma (Collins, 2006), we have Thus, equation (12) can be written in terms of I as follows I 2 f ðxÞ ¼ xIf ðxÞ À Ixf ðxÞ; and as an operator we have I 2 ¼ xI À Ix. We can easily check that D 2 I 2 ¼ 1 and also D 2 ðxI À IxÞ ¼ 1. The operator xI À Ix is called the normal form of I 2 . Now replacement lemma can be restated in operator-based notations as follows: Lemma 2.1. If f ðxÞ is integrable, then I 2 f ðxÞ ¼ xIf ðxÞ À Ixf ðxÞ for all x 2 ½0; 1�.
Proof. Result easily follows from the integration by parts as Ixf ðxÞ ¼ xIf ðxÞ À I 2 f ðxÞ.
Since D n I n ¼ 1, I n is the right inverse of L ¼ D n , and by Lemma 2.1, one can easily compute the normal form of I n by induction on n. Let the normal form of right inverse, I n , be denoted by L y and then LL y ¼ 1. We can also compute the normal form of the right inverse of the differential operator L by using the variation of parameters formula (see (Coddington & Levinson, 1955)) as follows: Lemma 2.2. For a given differential operator L ¼ D n with the fundamental system 1; x; x 2 ; . . . ; x nÀ 1 , the normal form of the right inverse of L is given by where I is the integral operator, i.e., If ¼ ð x 0 fdx, W i is the Wronskian matrix associated with 1; x; . . . ; x nÀ 1 whose i-th column is the n-th unit vector.
The solution of the semi homogeneous FDEs is computed as follows: It is clear that the solution of the semi non-homogeneous FDEs (see, e.g. (Thota, 2019;Thota & Kumar, 2016, 2015Thota & Kumar, 2016a, 2016b, 2015, 2017Thota et al., 2012;Thota, 2017Thota, , 2019a) depends on the interpolation technique. Indeed, to compute the solution of the IVP (11), it is enough to compute the interpolating function which is satisfying the initial data. Therefore the solution is Now the symbolic formulation of the solution of (9) in terms of normal form is given by, where ϕ ¼ À βðxÞyðq n xÞ þ ∑ nÀ 1 i¼0 c i x i i! ; L y is the right inverse of the differential operator L as given in the Lemma 2.2; and gðx; y; y 0 ; . . . ; y ðnÀ 1Þ Þ ¼ λyðxÞ þ ∑ nÀ 1 j¼0 μ j ðxÞy ðjÞ ðq j xÞ þ f ðxÞ.
As mentioned, in introduction, that one of the advantages of HPM is that the perturbation equation can be freely constructed in many ways by homotopy in topology and the initial approximation can also be freely selected. Here, we are going to discuss two ways to construct the perturbation equation. One way to construct the perturbation equation is as follows.
According to HPM (He, 1999;He, 2000), we can write the solution (14) as a power series in p, Setting p ¼ 1, the approximate solution of equation (15) is To solve the given equation (1) using HPM, we have the convex homotopy from equation (7) as follows: Hðv; pÞ ¼ ð1 À pÞ½v ðnÞ ðxÞ À y . . . Identifying the components v i 's, the nth approximation of the exact solution can be obtained as Another way to construct the perturbation equation is as follows.
Using the HPM, we can construct the following homotopy for the equation (1). The functions v n , for n ¼ 0; 1; . . . , are determined as follows (see, e.g. (He, 1999), (He, 2000)): The functions v n are obtained by equating the corresponding coefficient of power of p on both sides, for n ¼ 1; 2; 3; . . . , Here, we are mainly concerned about the perturbation equation that can be constructed easily in many ways by homotopy, and the initial guess that can also be selected freely by HPM. Once one chooses these parts, the homotopy equation is completely determined, because the remaining part is actually the original equation and we have less freedom to change it.
Remark: This theory can also be easily extended to fractal calculus.

Error involving due to approximation
If y k ðxÞ is an approximate solution at level n ¼ k, then the error involved due to the approximation is calculated as follows and it must be approximately zero, for i ¼ 0; 1; . . . , In the following section, we present several numerical examples to illustrate the simplicity and efficiency of the proposed method.

Numerical examples
EXAMPLE 3.1. Consider the following multi-pantograph equation where λ ¼ a 2 and μ ¼ a 2 e ax 2 . In symbolic notations, we have LyðxÞ ¼ f ðx; yÞ where L ¼ D, f ðx; yÞ ¼ a 2 e ax 2 yð x 2 Þ þ a 2 yðxÞ and c 0 ¼ 1. The right inverse of L, computed as discussed in Section 2, is L y ¼ I.
By HPM, we have the following convex homotopy Comparing the corresponding coefficient of power of p on both sides, for n ¼ 1; 2; 3; . . . , and so on, in this manner the rest of components of the homotopy perturbation solution can be obtained. Now we can approximate the analytical solution yðxÞ by truncated series as y n ! e ax asn ! 1: Figure 1 shows the graphical comparison between the exact solution and the approximate solution obtained by the proposed method for a ¼ 1. It can be seen that the solution obtained by the present method nearly identical to the exact solution. The simplicity and accuracy of the proposed method is illustrated by computing the absolute error E 7 ¼ jy exact ðxÞ À y appox ðxÞj. Figure 2 shows the absolute error between the exact and approximate solution at level n ¼ 7 which is significantly small; thus, it indicates the convergence of series solution very rapidly. The accuracy of the results can be improved by introducing more terms of the approximate solutions. In Table 1, the solution is compared with the well-known exact solution of the given equation (20). From the numerical results in Table 1, it is clear that the approximate solution is in high agreement with the exact solution. The table also shows the absolute error between the exact solution and approximate solution. The convergence rate of the solution series is very fast. EXAMPLE 3.2. Consider a second-order neutral functional-differential equation with proposed delay y 00 ðxÞ ¼ y 0 ð x 2 Þ À x 2 y 00 ð x 2 Þ þ 2; x 2 ½0; 1� yð0Þ ¼ 1; y 0 ð0Þ ¼ 0: In symbolic notations, we have The symbolic representation of (22) is  and so on, in this manner the rest of components of the homotopy perturbation solution can be obtained. Now we can approximate the analytical solution yðxÞ by truncated series as y n ! x 2 asn ! 1: In Figure 3, we show the graphical comparison between the exact solution and the approximate solution obtained by the proposed method at the level n ¼ 5 and the absolute errors E 4 ; E 5 , and E 6 are given in Table 2. EXAMPLE 3.4. Consider a following NFDE with proportional delay (Biazar & Ghanbari, 2012) yð0Þ ¼ 1: Following the proposed method, similar to Example 3.1, the seventh order of HPM approximate solution is x 2 À 82677 524288 x 3 þ 1240155 33554432 x 4 À 1736217 268435456 x 5 þ 1736217 2147483648 x 6 À 248031 4294967296 x 7 ; and the exact solution of the given NFDE is yðxÞ ¼ e À x . The graphical comparison of the approximate solution at seventh iteration with the exact solution is shown in Figure 4. From Figure 4, it is clear that the approximate solution almost same as the exact solution.