Application of Vieta–Lucas Series to Solve a Class of Multi-Pantograph Delay Differential Equations with Singularity

: The main focus of this paper was to ﬁnd the approximate solution of a class of second-order multi-pantograph delay differential equations with singularity. We used the shifted version of Vieta– Lucas polynomials with some symmetries as the main base to develop a collocation approach for solving the aforementioned differential equations. Moreover, an error bound of the present approach by using the maximum norm was computed and an error estimation technique based on the residual function is presented. Finally, the validity and applicability of the presented collocation scheme are shown via four numerical test examples.


Introduction
Delay differential equations (DDEs) [1][2][3] play a pivotal role in mathematical modeling when the rate of change of some time-dependent phenomena does not only depend on its current state but also depends on a certain past state. In the current state, it might possible that we are subject to the poor performance of our model but when we include the performance of some former state, we might have better results.
In other words, delay differential equations arise whenever a phenomenon has a delayed effect in a differential equation. Since this happens quite frequently, DDEs are ubiquitous in nature. Their most common applications include models of population dynamics [4], control theory [5], dynamics of neural networks [6], models on the spread of diseases [7] and applications to feedback problems in electrical and electronics engineering [8]. Among other applications are modeling problems in quantum dot lasers [9] and topics in climate science [10]. Interested readers are referred to [11] for an extensive survey of the uses of DDEs in mathematical biology, where they can also find a monograph of the numerical solutions of DDEs [12].
We will now have a brief introduction of pantograph equations. A pantograph equation is actually the most prominent delay differential equation with proportional delay which frequently appears in many scientific and engineering models, and plays an important role in diverse applications such as number theory, astrophysics, electrodynamics, probability theory on algebraic structures, quantum mechanics, nonlinear dynamical systems, and cell growth [13][14][15][16][17]. The name pantograph originated in 1851 from the work of Ockendon and Tayler on the collection of current by the pantograph head of an electric locomotive [14].
The main subject of this research paper was to consider a polynomial-based technique to obtain an approximate solution to the following second-order singular multi-pantograph delay differential equation (SMDDE) in the form [18][19][20]: associated to the initial conditions χp0q " χ 0 , χ 1 p0q " χ 1 , where the pantograph parameters are 0 ă τ s ă 1 for s " 1, 2, . . . , m and χ 0 , χ 1 are the known real numbers. Here, the functions p s pηq, s " 1, 2, . . . , m and qpηq are assumed to be continuous. In recent years, many studies have been carried out on the approximate and numerical solutions of differential equations of pantograph type. We reviewed some of the previous studies in this section. Recently, the generalized pantograph type differential equations have been approximately solved by using the shifted Legendre polynomial solutions, the special matrix approach, the Taylor operational matrix method, and an approach based on the combination of the variational iteration method and Laplace transform scheme [21][22][23][24]. Furthermore, singular-type differential equations with delay have been solved by the Bernoulli, Bessel, reproducing kernel and machine learning approaches, as well as by neural-swarming heuristic methods in [2,3,[18][19][20]. In addition, the collocation approach based on the Jacobi rational-Gauss collocation scheme, the stable Runge-Kutta methods, the shifted orthonormal Bernstein polynomials method and the least squares-Epsilon-Ritz algorithm were introduced for solving generalized pantograph equations [25][26][27][28].
In recent decades, collocation techniques based on (orthogonal) polynomials have been successfully applied to various areas of physical science and engineering. They have been proven to be efficient, robust, and provide exponential convergence in many application areas, as can be seen in [29][30][31][32][33][34][35]. The main goal of this study was to give an efficient collocation-based polynomial approximation for the solution to the SMDDE (1). In this respect, a new class of orthogonal functions called Vieta-Lucas polynomials is considered. This technique is based on taking the truncated VL expansion of the function in the SMDDE. Hence, the resulting matrix equation can be solved, and thus the unknown VL coefficients can be approximately determined.
The rest of this research paper is structured as given below. In Section 2, we give a detailed overview of Vieta-Lucas (VL) functions having definitions of VL polynomials and shifted VL polynomials. Some of their important and symmetry properties were also reviewed. Section 3 is devoted to presenting a collocation approach called a VL-collocation scheme to solve a class of singular multi-pantograph delay differential equations. By exploiting the shifted VL basis functions together with a set of discrete collocation points, this technique reduces the underlying model (1) into an algebraic linear fundamental matrix equation. Moreover, an error bound in the maximum norm for the VL-collocation scheme applied to the SMDDE was obtained and an error estimation method was introduced in Section 4. In the numerical Section 5, three numerical test examples were simulated to illustrate the accuracy and efficiency of the proposed VL-collocation method. We finally summarize the study in Section 6.
We clearly have VL 0 pξq " 2. The explicit forms of the remaining polynomials are obtained as follows: where r¨s denotes the ceiling function. A simple calculation using j " 1 in (4) shows that VL j pξq " ξ. Using VL 0 pξq " 2 and VL 1 pξq " ξ, we can also generate the VL polynomials by the following recursive relation: On account of the former recursion, the next five VL polynomials are obtained as Furthermore, the class of VL polynomials forms an orthogonal system on Ω. The orthogonality condition is given by where ωpξq " p4´ξ 2 q´1 {2 is the related weighting function. Lemma 1. The VL polynomial upξq :" VL j pξq satisfies the following differential equation: Proof. If we differentiate upξq with respect to ξ, we have: Similarly, for the second derivative, we obtain: We then insert the expressions of derivatives u 1 , u 2 , u " 2 cospjϕq and ξ " 2 cos ϕ into the given differential equation. After simplifying, all terms are canceled out. Thus, the proof is completed.

Remark 1.
We note that the VL differential Equation (6) is symmetric with respect to the transformation ξ Ñ´ξ.

The Shifted VL Polynomials
In order to utilize the family of VL polynomials for our model problem, we are required to transform them on the interval r0, 1s. In this respect, we apply the change of variable ξ " 4η´2 on the differential Equation (6). It follows that: where VL ‹ j pηq " VL j p4η´2q is the shifted VL polynomials defined on Ω ‹ :" r0, 1s. In view of (4), we can express the explicit form of shifted VL polynomials as the (unique) solution of (7). After some manipulations, we obtain: Accordingly, a few shifted VL polynomials on r0, 1s are constructed as follows: Accordingly, the weight function ωpξq for the shifted VL polynomials on Ω ‹ is expressed as ω ‹ pηq " pη´η 2 q´1 {2 . Hence, the orthogonal condition (5) for the shifted VL polynomials is written as A function gpηq P L 2 ω ‹ pΩ ‹ q may be expanded in terms of VL functions as gpηq " where the unknown coefficients e r can be computed with the aid of orthogonal property (9). In practice, we express the solution to the SMDDE (1) as a combination of the finite number of pR`1q VL functions as In order to devise the VL matrix scheme, we construct the vector of VL basis functions as Here, G R pηq is the vector of the monomial basis function defined by Additionally, the matrix H R in (11) is an upper-triangular structure of size pR`1qˆpR`1q: For instance, by setting R " 5, we arrive at: Our ultimate goal is to determine the R`1 unknown coefficients e r , r " 0, 1, . . . , R in (10) by utilizing the shifted VL-collocation matrix technique. To this end, let us define the vector: E R " re 0 e 1 . . . e R s T .
Thus, we can rewrite the approximate solution χ R pηq compactly as

The VL-Collocation Approach
Supposedly, the considered SMDDE (1) has a solution in terms of the truncated series expansion of shifted VL functions. Our aim was to systematically estimate the unknown coefficients of the assumed solution. Thus, we are required to express χ R pηq and all its derivatives in the matrix representation forms. To begin, by using (11) and (12), the truncated series solution χ R pηq in (10) can be rewritten as We now associate a set of collocation points tη r u R r"0 on Ω ‹ to our approximation technique. In this respect, the following Chebyshev points on Ω ‹ are taken as where ν r " cos´2 rπ´π 2R`2¯a re the Chebyshev nodes on p´1, 1q. Now, the collocation points (14) are inserted into (13). The resulting system of equations is: By introducing the vectors G R and X, we can represent the former equations in a compact form as We then consider the evaluation of the derivatives χ R pηq in (1). Their matrix forms at the collocation points are illustrated in the subsequent lemma. Lemma 2. The matrix expressions of χ 1 R pηq and χ 2 R pηq evaluated at the collocation points are given by Proof. Differentiating once from (13), we obtain: A straightforward calculation shows that the derivatives of G R pηq can be presented as a product of itself and the powers of a (sparse) matrix M R " pm i,j q R i,j"0 as Therefore, we combine (18) and (19) to express the first derivative of the unknown solution in the form: In an analogue way, we can write the second-order derivative of χ R pηq by the following form: Consequently, the Chebyshev nodes (14) are substituted into the relations (20) and (21). Finally, we define two vectors X r1s and X r2s as Finally, we consider the delayed-like terms χ 1 R pτ s ηq appeared in (1) to be written in the matrix forms. For this purpose, we first utilize (20) and then replace η by τ s η to obtain: To treat the first term G R pτ s ηq on the right-hand side, we use the fact that it can be represented as G R pτ s ηq " G R pηq N s τ , s " 1, 2, . . . , m, where the delayed matrices N τ s depend on the parameters τ s which have diagonal representations such as Thus, we can update the expression for d dη χ R pτ s ηq in (22) as Let us now put the collocation points (14) into (22). It follows that the matrix representations for d dη χ R pτ s ηq, s " 1, 2, . . . , m can be expressed as To date, all involved terms in the model problem (1) are written in matrix form. To proceed, we replace the Chebyshev nodes (14) as the collocation points into (1). This gives us the system: To continue, we introduce the following matrices with diagonal structures as These allow us to rewrite the Equations (25) in the following matrix notation: Lemma 3. Suppose that the approximate solution to the SMDDE (1) can be written as (13). To determine the unknown E R , the following fundamental matrix of the equation is obtained: where: Proof. The proof is straightforward if we put the relations (15), (17), and (24) into the matrix Equation (26).
The above Lemma indicates that by solving the linear system of simultaneous Equation (27), we obtain the unknown Vieta-Lucas coefficients e r for r " 0, 1, . . . , R. However, the initial conditions (2) have to date been left unspecified. To do so, let us first pay attention to χp0q " χ 0 . By placing η to 0 in (13), we have the matrix expression for the first initial condition as For the second initial condition χ 1 p0q " χ 1 , we utilize the relation (20). Similarly, we write η to 0 in (20) to have: In the final step, the replacement of two (arbitrary) rows of the fundamental matrix Equation (27) was carried out. In other words, two row matrices r s Z 0,R ; χ 0 s and r s Z 1,R ; χ 1 s enter the augmented matrix rZ R ; F R s in (27). Thus, we obtain the modified fundamental matrix equation r s Z R ; s F R s. After solving this equation, the unknown VL coefficients in (13) will be found.

L 8 -Error Bound and Error Analysis
Here, we analyze the error for the presented collocation technique. Firstly, an upper boundary of errors is given and then a theorem to estimate the error function is presented.

Theorem 1 (Upper Boundary of Errors).
We assume that χpηq and χ R pηq " V R pηqE R denote the exact and the VL polynomial solutions with R-th degree for model problem (1) in the interval 0 ď η ď L, L ą 0, respectively. Let χ Mac R pηq " G R pηq p E R show the Maclaurin expansion series by R-th degree of χpηq. Then, the error of the VL polynomial solution χ R pηq satisfies: where G R pηq " " 1 η η 2 . . . η R ‰ , V R pηq and the matrix H R are defined in (11), and we also have k R " maxtL R , 1u.
Proof. To begin, we add and subtract the Maclaurin expansion χ Mac R pηq with the R-th degree in the error }χpηq´χ R pηq} 8 . Utilizing the triangle inequality give us: We know from relation (13) that the VL polynomial solution χ R pηq " V R pηqE R can be expressed as χ R pηq " G R pηqH R E R . In addition, according to the hypothesis of the Theorem, we can express the Maclaurin series of χpηq as χ Mac R pηq " G R pηq p E R . From these data, we can obtain the following inequality: The norm }G R pηq} on r0, Ls is bounded by Hence, we can arrange Equation (30) as follows: From the calculus, we know that the remainder error of the Maclaurin expansion series χ Mac R pηq by R-th degree is given by From here, we can write the inequality: Finally, we insert the relations (31) and (32) into (29). Thus, we obtain the inequality: Hence, the proof is finished.
Theorem 2 (Error Estimation). Let χpηq and χ R pηq be the exact solution as well as the VL polynomial solution with the degree-R of the initial value problems (1) and (2), respectively. Then, the following error problem to estimate the error function Φ R pηq " χpηq´χ R pηq is gained: where F R pηq represents the residual function related to the initial value problems (1) and (2) for the VL solution χ R pηq.
Proof. Let us substitute the VL polynomial solution with the R-th degree in (1). Then, we have: By arranging Equation (35), we can write: Since the VL polynomial solution (10) satisfies the initial conditions (2), we can write the expression: Now, we subtract the relations (36) and (37) from problems (1) and (2) and thus we have the error problem (34) as required.

Corollary 1.
By solving the error problem in Theorem 2 by using the proposed VL-collocation technique in Section 3, the error function can be estimated. This kind of error estimation, which is known as the residual correction technique, can be utilized for an improvement of the accuracy of the approximate solutions, as can be seen in [38].

Numerical Test Examples
We checked the accuracy of the present Vieta-Lucas collocation scheme by solving four numerical examples. A comparison between the obtained computational results and the true exact solutions was performed. For a validation study, we compared our numerical results with those reported in the existing published works. MATLAB R2021a software was used for simulations.
By computing the approximated solution χ R pηq for this example with R " 5, we obtain: . . .
We highlight the difference between the above exact and approximate solutions in Figure 1. In addition to R " 5, 10, the graph of absolute error using R " 15 is further visualized in Figure 1. We then reported the accuracy of the proposed VL-collocation approach by calculating the values of absolute errors via: E R pηq :" |χ R pηq´χpηq|.
(38) Table 1 shows the numerical outcomes obtained at some points η P Ω ‹ utilizing R " 5 and R " 10. Additionally, the absolute errors are presented in Table 1. For validation, the results of an existing numerical procedure are shown. For this purpose, the outcomes of the called ANN-PSO-AS approach are based on the stochastic method via artificial neural networks [19]; the statistical values performed after sixty executions utilizing the median (MED), minimum (Min), semi-interquartile range (SI-R), and mean are reported. One can readily observe that our computational results with R " 5 are comparable with those obtained via ANN-PSO-AS method, and for R " 10, are superior and significantly more accurate. Test Example 4. In the second test problem, we consider the following second-order SMDDE given by [19] The initial conditions are χp0q " 1 and χ 1 p0q " 0. The exact solution of this problem is given by χpηq " cospηq.
For the second example, we set R " 5 and R " 10. The approximated solutions χ R pηq obtained via VL-collocation procedure are given by . . , which can be used for comparison with our approximate solutions obtained above. The approximate solution obtained with R " 5 along with the exact one are plotted in Figure 2.
The graphs of the achieved absolute errors using R " 5, 10 and R " 15 are further shown in Figure 2. It can be clearly seen that there is excellent agreement between our approximate solutions and the related exact solution. Some comparisons in terms of achieved (absolute) errors were also made in Table 2 for the Test Example 4. Similar to Table 1, the method ANN-PSO-AS was utilized for comparisons. Numerical solutions for R " 5, 10 evaluated at points η j " j{20 for j " 0, 1, . . . , 20 are further shown in Table 2. Clearly, the results of the proposed VL-collocation procedure are more accurate than those obtained via ANN-PSO-AS. This shows the advantage of our approach and also due to its ease in implementation.

Test Example 5.
In the third test example, we pay our attention to the following second-order SMDDE as [19] The initial conditions are χp0q " 0 and χ 1 p0q " 1. The analytical solution of the above problem is χpηq " sinhpηq.
Let us first consider R " 5. The approximated solution χ 5 pηq over Ω ‹ is gained as Similarly, by taking R " 10, the approximate solution obtained via VL-collocation method is as follows χ 10 pηq " 0.0000001429027479 η 10`0 .00000241264259 η 9`0 .0000004596522829 η 8`0 .0001980345147 η 7 0.0000001967391815 η 6`0 .008333269637 η 5`0 .00000001189264844 η 4`0 .1666666656 η 3 The graphs of absolute errors E R pηq for diverse values of R are visualized in Figure 3. The exponential convergence of the proposed scheme can be clearly seen in Figure 3. Furthermore, we present the benefits of the proposed method and validate our results. To this end, some comparisons with respect to the absolute errors are performed in Table 3 with R " 5, 10. Our results are compared with the outcomes of ANN-PSO-AS [19]. Clearly, the performance of the VL-collocation approach outperforms, ANN-PSO-AS [19] particularly when R is becomes large.
Utilizing diverse values of R " 2 i , i " 0, 1, 2, 3, 4, the computed L 8 error norms together with the related NOC are shown in Table 4. The graph of E R,8 errors are also depicted in Figure 4. It can be clearly seen that the numerical order of the convergence rate has an exponential behavior. This shows the advantage of the present method and its ability to achieve an arbitrary order of convergence when R is increased.  Test Example 6. The last test case is devoted to a slightly different second-order multi-pantograph differential equation as [23,24] where: The initial conditions are χp0q " 1 and χ 1 p0q "´1. The analytical solution of this model has a rational structure and is given by χpηq " 1 η 2`η`1 .
In Figure 5, we visualize the graphs of absolute errors E R pηq for different values of R " 8, 12, 16, 20, 24, R " 28, which again shows the exponential convergence of the VL-collocation procedure.  For comparison, let us consider the approximate analytical solution obtained via a combination of the Laplace transform and variational iteration method (LVIM) [24] and the operational matrix method based on the Taylor polynomials (OMMTP) [23]. The numerical results using R " 10 and R " 15 along with the achieved absolute errors are reported in Table 5. As shown in Figure 5, more accurate results were obtained by increasing R.

Conclusions
A collocation approach based on the novel family of orthogonal functions was developed for the numerical treatment of a class of second-order singular multi-pantograph delay differential equations (SMDDEs). The Vieta-Lucas (VL) functions we reviewed as the bases from the differential perspective and their error in the maximum norm was analyzed. The effectiveness of the developed VL-collocation procedure was demonstrated with the aid of four test examples. The experimental results reported by the developed collocation technique revealed that the VL matrix approach was very accurate, robust and convenient for the SMDDE and related model problems in comparison with some existing computational and analytical approaches.