A Galerkin-like method for solving linear functional differential equations under initial conditions

In this paper, we present a weighted residual Galerkin method to solve linear functional differential equations. We consider the problem with variable coefficients under initial conditions. Assuming the exact solution of the problem has a Taylor series expansion convergent in the relevant domain, we seek a solution of the given problem in the form of a polynomial having degree N of our choice. Substituting this polynomial with unknown coefficients in the given equation yields an expression linear in these coefficients. We then proceed as in the weighted residual method and take inner product of this expression with monomials up to degree N , resulting in N +1 linear algebraic equations. Appropriately incorporating the initial conditions and solving the resulting linear system, we obtain the approximate solution to the given problem. Additionally, we present a way of estimating the absolute error of the obtained approximation, which is then used to improve the original approximation through a method called residual correction. We also show that the upper bound for the error of the proposed method depends on the Taylor truncation error of the exact solution. The proposed scheme and the residual correction technique are illustrated in several example problems.


Introduction
The first order delay differential equation (DDE) with constant coefficients and proportional delay u ′ (x) = au(x) + bu(qx), x ∈ I = [0, T ], 0 < q < 1 (1.1) is called the pantograph equation and first appeared in the mathematical modelling of the wave motion in the current line between an electric locomotive and its overhead catenary wire [1,2]. from C m+1 (I) can be found which also satisfies the same initial condition [3]. interested reader is referred to [4,5] for analytical and approximate solutions of the pantograph equation and for DDEs in general.
The pantograph equation plays a role in a wide range of applications. This has brought about an increasing interest on the numerical solution of DDEs of pantograph type and researchers have made use of numerous methods. To name a few, Adomian decomposition method was used for this purpose by Dehghan et al. in [6], where the convergence of the approach was also established. Another popular method, homotopy * Correspondence: syuzbasi@akdeniz.edu.tr perturbation method, was used by Feng [7], where the multipantograph equation was considered with variable coefficients. The authors of this paper used a weighted residual Galerkin approach [8] to obtain approximate solutions of the same problem. Ishiwata and Muroya used two kinds of collocation methods based on special mesh distributions in [9]. Other collocation methods were presented in conjunction with Bernstein [10] and Bessel [11] polynomials for the same purpose. In [12], pantograph type DDEs were solved by shifted Chebyshev polynomials and an applicable error analysis was presented. Another collocation method that is based on exponential functions rather than polynomials can be found in [13].
In this paper, we consider the linear nonhomogeneous m-th order generalized pantograph equation with variable coefficients, given by under the initial conditions Here, the given functions P jk (x) and f (x) in equation ( Bernoulli polynomials [14], a Taylor polynomial approach [15], a collocation method involving Lucas polynomials [16], variational iteration method [17], and an operational matrix method using polynomials in the standard basis [18]. In addition to these, pantograph equations have been solved using shifted orthonormal Bernstein polynomials [19]. There are also studies interested in solving a nonlinear version of equation (1.2); the reader can see [20], for example, where Lucas polynomials are employed in conjunction with collocation points to solve nonlinear delay differential equations of this type. In addition, residual correction technique described in Section 3.1 is applied to two example problems and comparisons with other methods are made whenever possible. Finally, conclusions of the paper are summarized in Section 5.

Solution method
In this section, we outline the Galerkin-like procedure we will use in order to solve equation (1.2). The method, which relies on taking inner product of an expression with monomials, also called the method of moments, was employed to obtain approximate solutions of high-order Fredholm integro-differential equations [21] and high-order integro-differential equations with weakly singular kernel [22]. In this study, we develop the method to compute the approximate solutions of generalized pantograph-type functional differential equations.
To begin with, we assume that the unique solution y(x) of equation (1.2) can be expressed as a power series of the form We then truncate this power series after the (N + 1) st term so that Here, the coefficients a i will be obtained as the output of the procedure. The derivative y holds for any nonnegative integer k , where y (0) n denotes the function y n itself. In order to express the delayed terms in equation (1.2), equation (2.2) can be utilized in a straightforward manner. For this purpose, given two constants λ and µ , we replace the variable x in X(x) with λx + µ and write equation (2.2) as It is worth noting that the vector X(λx + µ) can be expressed in a form that is more convenient for computer implementations. Namely, if we define the After that, substituting the matrix expressions for y taking the inner product of each element of Φ with both sides of equation (2.4) results in a linear system WA = F of N + 1 equations in the unknowns a 0 , a 1 , . . . , a N . Explicitly, the coefficient matrix W and the column vector F on the right-hand side of this system is given by where the inner product is defined by Here, f and g are functions from the Hilbert Since the initial conditions in equation (1.3) should also be satisfied, m rows of W and the corresponding m entries of F should be modified accordingly. For the sake of being deterministic, let us fix the modified entries to be the last m rows of W and the last m entries of F . Determining the matrix equivalents of the initial conditions is easy in view of equation (2.2). Namely, ( i + 1 )-st initial condition is given by . Thus, feeding the m initial conditions into the system for the last m rows of W and the lengthm vector [ α 0 α 1 . . . α m−1 ] T for the last m entries of F . After performing this step, we are left with the modified systemWA =F explicitly given bỹ from which we obtain the matrix of unknown coefficients A =W −1F and the approximate solution

Error analysis and residual correction
In this section, we first present a way of improving the accuracy of an already obtained approximate solution. Then, we give an upper bound for the absolute error in terms of the Taylor truncation error of the exact solution.

Error estimation and residual correction
In situations where it is not possible to measure the accuracy of an approximate solution, the residual function is useful since it gives an idea about the efficiency of this approximation. Furthermore, when the equation is linear in the unknown function, it is possible to exploit the residual function of the approximate solution in order to obtain a more accurate solution. In our case, since the generalized pantograph equation given by equations (1.2) is linear in the unknown function y(x) , the method of residual correction can be applied in a rather straightforward manner. This section is about how the residual function gives rise to an error estimation of an approximate solution to equations (1.2) obtained as the output of our method. We also exhibit how this estimation is used to obtain better approximate solutions.
To this end, let us consider the residual function of equation (1.2). We now replace y(x) by the approximate solution y N (x) to obtain as the residual function of y N (x) . Subtracting equation (

Error bound for the solution
In this part, we relate the error bound for the approximate solution y N (x) to the truncation error of the Taylor polynomial corresponding to the exact solution.

Theorem 3.1 Let y N (x) and y(x) denote the approximate and the exact solutions of problem (1.2), respectively.
If where y T N (x) denotes the N th degree Taylor polynomial of y around the point x = q ∈ [0, b] and R T N represents its remainder term.
Proof Since y is (N + 1) -times continuously differentiable, it can be represented by its Taylor series as is the reminder term of the Taylor expansion of y . Thus, y(x) − y T N (x) = R T N (x). By using this and the triangle inequality, we obtain 2 As a result, we have found an upper bound of the absolute error in terms of the Taylor truncation error of the exact solution. Note that this is not an a priori error bound; it only serves as a means to compare the actual error to this Taylor truncation error.

Illustrative examples
In this section, we apply the method explained in Section 2 to several generalized pantograph equations and compare the resulting approximate solutions with some other methods present in the literature. We also apply residual correction to obtain better approximate solutions using the existing ones. All the calculations have been performed using MATLAB. [23] and then in [14,24,25]:

Example 4.1 Let us first consider the third-order generalized pantograph equation studied in
The exact solution of this problem is y(t) = e −t . Table 1 compares the absolute errors of the solutions obtained by the present method, the Taylor series method [23], the Hermite method [24], the Chebyshev-Gauss grid method [25], and the collocation method based on Bernoulli operational matrix [14] for  Table 2 makes it clear that the accuracy of the improved solutions is a direct consequence of the accuracy of error estimations, which is in parallel with the remark made at the end of Section 3.1. It is also seen that the Chebyshev method CMBOM method Hermite method Present method Hermite method with N=8 Taylor method with N=5 Chebyshev method with N=8 CMBOM method with N=8 Present method with N=5 Present method with N=8 estimates corresponding to M = 13 are more accurate than those corresponding to M = 10 , which is because the proposed method is more accurate when implemented using higher-degree polynomials. These remarks can also be confirmed from Figure 2, where the data in Table 2 are depicted graphically.
and its absolute truncation error is given by As for the second part of the error bound (3.6), we first compute the maximum distance between the fifth degree Taylor polynomial of e −x around x = 0 and the approximate solution  gives |y(x) − y 5 (x)| < 0.0026 . This procedure has also been carried out for the values N = 8 and N = 11 . The obtained bounds are shown in Table 3 together with the maximum actual errors |e N (x)| corresponding to these N values. It is understood that the error bounds for this problem are rather loose since the actual maximum error values are much smaller in reality. It should be stressed once more that the error bound (3.6) is not an a priori error bound; it only shows that the actual error is bounded in part by the Taylor truncation error of the actual solution.

Example 4.2
Our second example is the following third order generalized pantograph equation from [13], whose exact solution is y(x) = e −x .
We applied residual correction to y 12 (x) for the values M = 4 and M = 7 . Table 4 and Figure 3 illustrate the absolute error values of the original approximation y 12 (x) and of the improvements y 12,4 (x) and y 12,7 (x) for several values of x . It can be commented looking at the table that residual correction greatly reduces the absolute error also for this problem. When we examine the two graphs in the figure, just like as in the first example problem, we see that the estimated absolute errors |e 12,4 (x)| and |e 12,7 (x)| are very close to the actual error |e 12 (x)|, which explains the significant decrease in the actual corrected errors |E 12,4 (x)| and |E 12,7 (x)|.    Implementing the present method with N = 2 yields the solution y 2 (x) = x 2 , which is the exact solution. In fact, for any choice of N > 2 we get a k = 0 for 2 < k ≤ N as a result of the algorithm, which means any N ≥ 2 yields the exact solution y(x) = x 2 . This is not a surprise since the scheme described in Section 2 makes it clear that the unknown coefficients of the approximate solution y N (x) obtained as a result are equal to the actual coefficients of the exact solution. Therefore, this example application shows that the present method yields the exact solution in case this solution is a polynomial.

Conclusions
In this paper, we have presented a Galerkin-like approach for the approximate solution of linear functional differential equations. The method includes taking inner product of an expression derived from equation (1.2) with a set of monomials. The results of this approach are compared with some popular methods present in the literature and it is revealed that the present method performs slightly better than the others in terms of absolute error. Residual error correction technique to improve the accuracy of approximate solutions has also been discussed. Simulation results show that significant improvements in the approximate solutions can be achieved as a result of employing this technique. It is stressed that this fact can be attributed to the accuracy of our error estimation related to the approximate solution. Lastly, the presented method implemented with parameter N gives rise to the exact solution in case the exact solution is a polynomial of degree at most N .