Recent Advances on Methods and Applications of Nonlinear Differential EquationsView this Special Issue
Research Article | Open Access
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 , and the Bézier curves are determined on any subinterval by control points. The approximated solution of delay systems containing inverse time is derived. To validate accuracy of the present algorithm, some examples are solved.
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 . 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 are not equal. Delay systems containing inverse time are an important class of systems: 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 , as a third degree discontinuity at time , and, more generally, as a discontinuity in the st derivative at time .
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  and on stability by Niculescu . 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–10]).
Wu et al.  developed a computational method for solving an optimal control problem which is governed by a switched dynamical system with time delay. Kharatishivili  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  who obtained several results concerning bang-bang controls which are parallel to those of LaSalle  for nondelay systems. For a time-invariant system with an infinite upper limit in the performance measure, Krasovskii  has developed the forms of the controller and the performance measure. Ross  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  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 linear-quadratic regulator design for a system without delays (see [17, 18]).
This paper aims at solving delay systems containing inverse time of the following form: where , are, respectively, state and control functions while is known vector function and ’s () are nonnegative constant time delays, and . We assume the matrices , , , and are matrix functions. We need to impose the continuity condition on and its first derivative where these constraints appeared in Section 2.
Piecewise polynomial functions are often used to represent the approximate solution in the numerical solution of differential equations (see [19–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 , Nrnberger and Zeilfelder  used the Bézier control points in approximated data and functions. Zheng et al.  proposed the use of control points of the Bernstein-Bézier form for solving differential equations numerically and also Evrenosoglu and Somali  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–29]), the Bézier curves are used for solving dynamical systems (see ), 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 ). 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 ).
Using Bezier curve, Ghomanjani et al.  had used least square method for numerical solutions of time-varying linear optimal control problems with time delays in state and control.
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.
2. Function Approximation
Divide the interval into a set of grid points such that where and is a positive integer. Let for . Then, for , delay systems containing inverse time (2) can be decomposed to the following problem: where and are, respectively, vectors of and which are considered in . We mention that , , is the th component of where , and , , is the th component of where . Also where and denote the integer part of and , respectively.
Our strategy is to use Bézier curves to approximate the solutions and by and , respectively, where and are given below. Individual Bézier curves that are defined over the subintervals are joined together to form the Bézier spline curves. For , define the Bézier polynomials of degree that approximate, respectively, the actions of and over the interval as follows: where is the Bernstein polynomial of degree over the interval , and are, respectively, and ordered vectors from the control points (see ). By substituting (6) in (4), for can be defined as follows:
Let and where and are, respectively, characteristic function of and for . Beside the boundary conditions on , at each node, we need to impose the continuity condition on each successive pair of to guarantee the smoothness.
Since the differential equation is of first order, the continuity of (or ) and its first derivative gives where is the th derivative with respect to at .
Thus, the vector of control points () must satisfy (see the Appendix) According to the definition of the we get that . Therefore, One may recall that is a ordered vector. This approach is called the subdivision scheme (or -refinement in the finite element literature). This method is based on the control-point-based method.
Remark 1. By considering the continuity of , the following constraints will be added to constraints in (10):
where the so-called () is an ordered vector.
Now, the residual function can be defined in as follows: where is the Euclidean norm (recall that is a vector where ).
Our aim is to solve the following problem over : The mathematical programming problem (14) can be solved by many subroutine algorithms. Here, we used Maple 12 to solve this optimization problem.
Remark 2. Consider the following boundary value problem:
where , , , , , and are the vectors of appropriate dimensions. , , , , , and are the matrices of appropriate dimensions, and is nonnegative constant time delay.
Let where is the transpose; then satisfies that where where is the matrix with 1 at its entry and zeros elsewhere and is Kronecker product (see, e.g., [4, 38, 39]).
Remark 3. Now, the following delay differential equation can be considered:
with initial condition
where . In the case when is not finite, denotes the interval .
Furthermore, we assume that that is, (20) is a delay differential equation. The existence and uniqueness of the solution of initial value problem (20)-(21) was stated in .
Equation (20) is converted into a nonlinear programming problem (NLP) by applying Bézier control points method, whereas the MATLAB optimization routine FMINCON is used for solving resulting NLP. Numerical example shows that the proposed method is efficient and very easy to use.
Remark 4. Now, we limit ourselves to consider the following nonlinear delay differential equation in the type with the following initial conditions: where the differential operator is defined by .
3. Convergence Analysis
In this section, without loss of generality, we analyze the convergence of the control-point-based method applied to the problem (2) with time delay in state when , and the time interval is . So, the following problem is considered: where , , and , , are given real numbers and , , , , and are known polynomials for . The constant time delay is nonnegative.
Without loss of generality, we consider the interval instead of since the variable can be changed with the new variable by where .
Lemma 5. For a polynomial in Bézier form we have where is the Bézier coefficient of after degree-elevating to degree .
Proof. See [22, page 245].
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 ).
3.1. Degree Raising
Theorem 6. If the problem (25) with inverse time in state has a unique continuous trajectory solution , 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 , by the Weierstrass theorem (see ), one can easily find polynomials of degree and of degree such that , , , , and , where stands for the -norm over . Especially, we have
In general, and do not satisfy the boundary conditions. After a small perturbation with linear and constant polynomials , , respectively, for and , we can obtain polynomials and such that satisfies the boundary conditions , , and . Thus, , and . By using (28), one has
By the time, from , Now, we have so, Now, let , , , , , for every . Thus, for , an upper bound is found for the following residual: where is a constant.
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 , we have which gives Suppose and are approximated solution of (25) obtained by the control-point-based method of degree (). Let Define the following norm for difference approximated solution and exact solution : By (40), Lemma 5, the boundary conditions , , and , one can show that The last inequality in (41) is obtained by Lemma 5, where is a constant positive number. Now where the last inequality in (42) comes from (36). This completes the proof.
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 , (), where is the maximal width of all subintervals, then converges to as .
Proof. We first impose a uniform partition on the interval as , where .
Let be the cubic spline over which is interpolating to . Then, for an arbitrary small positive number , there exists a such that provided that . Let be the residual. For each subinterval , is a polynomial. On each interval , we impose another uniform partition as where , . Express in as By Lemma 3 in , there exists a () such that, when , we have Thus, or Now combining the partitions and all gives a denser partition with the length for each subinterval. Suppose is the approximate solution by the control-point-based method with respect to this partition, and denote the residual over by Define the following norm for difference approximate solution and exact solution : Then, there is a constant such that the last inequality in (50) is obtained by Lemma 5. It can be shown that By Lemma 3 in , we conclude that the approximate solution converges to the exact solution in order , (). This completes the proof.
4. Numerical Examples
Example 1. Consider the delay system containing inverse time described by (see ) where we have the following exact solution:
Let . Then, by (14) and choosing , we have the approximate solution
Example 2. Consider the boundary value problem described by (see ) From (18), we have (see ) where , and we have the following exact solution: Let . Then, by (14) and choosing , we have the approximate solution : The graphs of approximate trajectories are shown in Figures 3 and 4.