Bezier Curves Based Numerical Solutions of Delay Systems with Inverse Time

This paper applied, for the first time, the Bernstein’s approximation on delay differential equations and delay systems with inverse delay that models these problems. The direct algorithm is given for solving this problem. The delay function and inverse time function are expanded by the Bézier curves. The Bézier curves are chosen as piecewise polynomials of degree n, and the Bézier curves are determined on any subinterval by n + 1 control points. The approximated solution of delay systems containing inverse time is derived. To validate accuracy of the present algorithm, some examples are solved.


Introduction
Delay differential equations (DDEs) differ from ODEs in that the derivative at any time depends on the solution at prior times (and in the case of neutral equations on the derivative at prior times).
DDEs often arise when traditional pointwise modeling assumptions are replaced by more realistic distributed assumptions, for example, when the birth rate of predators is affected by prior levels of predators or prey rather than by only the current levels in a predator-prey model.
Because the derivative ẋ () depends on the solution at previous time(s), it is necessary to provide an initial history function to specify the value of the solution before time  = 0.In many common models the history is a constant; but nonconstant history functions are encountered routinely.
For most problems there is a jump derivative discontinuity at the initial time.In most models, the DDE and the initial function are incompatible: for some derivative order, usually the first, the left and right derivatives at  = 0 are not equal.Delay systems containing inverse time are an important class of systems: ẋ () =  ( − 1) , ẋ (0 A fascinating property is how such derivative discontinuities are propagated in time.For the equation and history just described, for example, the initial first discontinuity is propagated as a second degree discontinuity at time  = 1, as a third degree discontinuity at time  = 2, and, more generally, as a discontinuity in the ( + 1)st derivative at time  = .Delay differential equations are type of differential equations where the time derivatives at the current time depend on the solution and possibly its derivatives at previous times (see [1][2][3][4]).
The basic theory concerning the stable factors, for example, existence and uniqueness of solutions, was presented in [1,3].Since then, DDEs have been extensively studied in recent decades and a great number of monographs have been published including significant works on dynamics of DDEs by Hale and Lunel [5] and on stability by Niculescu [6].The interest in study of DDEs is caused by the fact that many processes have time delays and have been models for better representations by systems of DDEs in science, engineering, economics, and so forth.Such systems, however, are still not feasible to actively analyze and control precisely; thus the study of systems of DDEs has actively been conducted over the recent decades (see [7][8][9][10]).

Mathematical Problems in Engineering
Wu et al. [11] developed a computational method for solving an optimal control problem which is governed by a switched dynamical system with time delay.Kharatishivili [12] has approached this problem by extending Pontryagin's maximum principle to time delay systems.The actual solution involves a two-point boundary-value problem in which advances and delays are presented.In addition, this solution does not yield a feedback controller.Optimal-time control of delay systems has been considered by Oguztoreli [13] who obtained several results concerning bang-bang controls which are parallel to those of LaSalle [14] for nondelay systems.For a time-invariant system with an infinite upper limit in the performance measure, Krasovskii [15] has developed the forms of the controller and the performance measure.Ross [16] has obtained a set of differential equations for the unknowns in the forms of Krasovskii.However, Ross's results are not applicable to time-varying systems with a finite limit in the performance measure.
Basin and Perez [17] presented an optimal regulator for a linear system with multiple state and input delays and a quadratic criterion.The optimal regulator equations were obtained by reducing the original problem to the linearquadratic regulator design for a system without delays (see [17,18]).
This paper aims at solving delay systems containing inverse time of the following form: +  () u () , where Piecewise polynomial functions are often used to represent the approximate solution in the numerical solution of differential equations (see [19][20][21][22]).B-splines, due to numerical stability and arbitrary order of accuracy, have become popular tools for solving differential equations (where Bézier form is a special case of B-splines).There are many papers and books dealing with the Bézier curves or surface techniques.
Harada and Nakamae [23], N ü rnberger and Zeilfelder [24] used the Bézier control points in approximated data and functions.Zheng et al. [22] proposed the use of control points of the Bernstein-Bézier form for solving differential equations numerically and also Evrenosoglu and Somali [25] used this approach for solving singular perturbed two-point boundary-value problems.The Bézier curves are used in solving partial differential equations; as well, Wave and Heat equations are solved in Bézier form (see [26][27][28][29]), the Bézier curves are used for solving dynamical systems (see [30]), and also the Bézier control points method is used for solving delay differential equation (see [31,32]).
The Bézier curves method was presented which was stated for solving the optimal control systems with pantograph delays (see [33]).The method was computationally attractive and also reduced the CPU time and the computer memory and at the same time keeps the accuracy of the solution.The algorithm had been successfully applied to the pantograph equations.Comparing with other methods, the results of numerical examples demonstrated that this method was more accurate than some existing methods (see [33]).
Using Bezier curve, Ghomanjani et al. [34] had used least square method for numerical solutions of time-varying linear optimal control problems with time delays in state and control.
Some other applications of the Bézier functions and control points are found in [35][36][37] that are used in computer aided geometric design and image compression.
The use of the Bézier curves is a novel idea for solving delay systems containing inverse time.The approach used in this paper reduces the CPU time and the computer memory comparing with existing methods (see the numerical results).Although the method is very easy to use and straightforward, the obtained results are satisfactory (see the numerical results).We suggest a technique similar to the one used in [22,25] for solving delay systems containing inverse time.The current paper is organized as follows.
In Section 2, Function approximation will be introduced.Convergence analysis will be stated in Section 3. In Section 4, some numerical examples are solved which show the efficiency and reliability of the method.Finally, Section 5 will give a conclusion briefly.

Function Approximation
Divide the interval [ 0 ,   ] into a set of grid points such that where ℎ = (  −  0 )/ and  is a positive integer.Let   = [ −1 ,   ] for  = 1, 2, . . ., .Then, for  ∈   , delay systems containing inverse time (2) can be decomposed to the following problem: where [  /ℎ] and [  /ℎ] denote the integer part of   /ℎ and   /ℎ, respectively.Our strategy is to use Bézier curves to approximate the solutions x  () and u  () by k  () and w  (), respectively, where k  () and w  () are given below.Individual Bézier curves that are defined over the subintervals are joined together to form the Bézier spline curves.For  = 1, 2, . . ., , define the Bézier polynomials of degree  that approximate, respectively, the actions of x  () and u  () over the interval [ −1 ,   ] as follows: where is the Bernstein polynomial of degree  over the interval [ −1 ,   ], a   and b   are, respectively,  and  ordered vectors from the control points (see [22]).By substituting (6) in (4),  1, () for  ∈ [ −1 ,   ] can be defined as follows: Since the differential equation is of first order, the continuity of x (or k) and its first derivative gives where k ()  (  ) is the th derivative k  () with respect to  at  =   .
Thus, the vector of control points a   ( = 0, 1,  − 1, ) must satisfy (see the Appendix) According to the definition of the   =  0 + ℎ we get that   −  −1 = ℎ.Therefore, One may recall that a   is a  ordered vector.This approach is called the subdivision scheme (or ℎ-refinement in the finite element literature).This method is based on the controlpoint-based method.
Remark 1.By considering the  1 continuity of w, the following constraints will be added to constraints in (10): where the so-called b   ( = 0, 1,  − 1, ) is an  ordered vector.
Now, the residual function can be defined in   as follows: where ‖ ⋅ ‖ is the Euclidean norm (recall that  1, () is a  vector where  ∈   ).
Our aim is to solve the following problem over  = ⋃  =1   : The mathematical programming problem ( 14) can be solved by many subroutine algorithms.Here, we used Maple 12 to solve this optimization problem.

Convergence Analysis
In this section, without loss of generality, we analyze the convergence of the control-point-based method applied to the problem ( 2

Lemma 5. For a polynomial in Bézier form
we have where  , 1 + 1 is the Bézier coefficient of () after degreeelevating to degree  1 +  1 .
The convergence of the approximate solution could be done in two ways: (1) degree raising the Bezier polynomial approximation, (2) subdivision of the time interval.
In the following, the convergence in each case can be proven, although in numerical examples, we used only subdivision case (see [32]).

Degree Raising
Theorem 6.If the problem (25) with inverse time in state has a unique  1 continuous trajectory solution ,  0 continuous control solution , then the approximate solution obtained by the control-point-based method converges to the exact solution (, ) as the degree of the approximate solution tends to infinity.
Proof.Given an arbitrary small positive number  > 0, by the Weierstrass theorem (see [41]), one can easily find polynomials Now, let where Since the residual (  ) :=   () − () is a polynomial, it can be represented by a Bézier form.Therefore, we have Then, by Lemma 5, there exists an integer (≥ ) such that, when  1 > , we have which gives Suppose () and () are approximated solution of (25) obtained by the control-point-based method of degree  2 The last inequality in ( 41) is obtained by Lemma 5, where  is a constant positive number.Now ‖( () ,  ()) − ( () ,  ())‖ ≤   2 + 1 where the last inequality in (42) comes from (36).This completes the proof.

Subdivision
Theorem 7. Let (, ) be the approximate solution of the problem (25) with inverse time obtained by the subdivision scheme of the control-point-based method.If (25) has a unique solution (, ) and (, ) is smooth enough so that the cubic spline (, ) interpolates to (, ) and converges to (, ) in the order (ℎ  ), ( > 2), where ℎ is the maximal width of all subintervals, then (, ) converges to (, ) as ℎ → 0.

Numerical Examples
Applying the presented method, in Examples 1, 2, and 3, the Bézier curves are chosen as piecewise polynomials of degree 3.
Example 8. Consider the delay system containing inverse time described by (see [4]) where we have the following exact solution: Let () = 1.Then, by ( 14) and choosing  = 3,  = 6 we have the approximate solution 1.000000001 + 0.00001180813339 + 0.9996447663 2 + 0.0023695 The graphs of approximate trajectories are shown in Figures 1 and 2.
By choosing  = 1 and  = 7 (degree raising), we obtain the following solution: In Table 4, exact, numerical results of this method, method in [40], error of presented method, and error of the method in [40] are shown, respectively.
Example 13.Consider the following system described by (see [40]) (66) The solution of this problem is (67) By choosing  = 1 and  = 7 (degree raising), we obtain the following solution: In Table 5, exact, numerical results of this method, method in [40], error of presented method, and error of the method in [40] are shown, respectively.
Example 14.Consider the following LDDE described by with the initial conditions where the exact solution of this example is () =  − .Here, this problem is solved by choosing  = 10 and  = 3.
The graph of error is shown in Figure 5, and the following approximate solution () is found: The exact solution of this problem is () = In Table 6, exact, numerical results of this method, error of presented method, and error of the method in [43] are shown, respectively.

Conclusions
Using the Bézier curves, the general algorithm is provided for the delay systems containing inverse time.Numerical examples show that the proposed method is efficient and very easy to use.Now, we specify the procedure of derivation of (10) from (9).By (6), we have by substituting  =   into (A.5), one has   where it shows the equality (10).

(𝑡) and u(𝑡) which are considered
() and w  () for  ∈ [ −1 ,   ].Beside the boundary conditions on k(), at each node, we need to impose the continuity condition on each successive pair of k  () to guarantee the smoothness.
2. Here, this problem is solved by choosing  = 1 and  = 7.The following approximate solution () is found.

Table 1 :
Exact and estimated values of  1 () for Example 3.

Table 3 :
Exact and estimated values of () for Example 4.

Table 4 :
Exact and estimated values of () for Example 5.

Table 5 :
Exact and estimated values of () for Example 6.

Table 6 :
Exact and estimated values of () for Example 8.