Abstract
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.
1. 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 . 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 .
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–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–10]).
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 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 [23], Nrnberger 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–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–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.
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 [22]). 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 [40].
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 [32]).
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 [41]), 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
Since
so,
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.
3.2. 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 , (), 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 [22], 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 [22], we conclude that the approximate solution converges to the exact solution in order , (). This completes the proof.
4. 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 1. Consider the delay system containing inverse time described by (see [4]) where 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 1 and 2.
Example 2. Consider the boundary value problem described by (see [4]) From (18), we have (see [4]) 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.
Example 3. Consider the time-varying delay system described by (see [42])
The exact solutions are [42]
Here, this problem is solved by choosing and . the following approximate solutions and are found. In Tables 1 and 2, exact, numerical results of this method and obtained results in [42] are shown, respectively:
Example 4. Consider the following system described by (see [40]) Analytic solution of the initial value problem (IVP) is . By choosing and (degree raising), we obtain the following solution: In Table 3, exact and presented methods are shown, respectively.
Example 5. Consider the following system described by (see [40]) where the exact solution is on and on . By choosing and (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 6. Consider the following system described by (see [40]) The solution of this problem is By choosing and (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 7. 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 and . The graph of error is shown in Figure 5, and the following approximate solution is found:
Example 8. Consider the second-order linear decay differential equation: The exact solution of this problem is . Here, this problem is solved by choosing and . The following approximate solution is found. In Table 6, exact, numerical results of this method, error of presented method, and error of the method in [43] are shown, respectively.
5. 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.
Appendix
In this Appendix, we specify the derivative of Bézier curve.
By (6), we have where .
Now, we have (see [44]) where , and By using (A.2), the first derivative is shown as Now, we specify the procedure of derivation of (10) from (9).
By (6), we have by substituting into (A.5), one has To preserve the continuity of the Bézier curves at the nodes, one needs to impose the condition ; so from (A.6), we have From (A.4), the first derivatives of and are, respectively, By substituting into (A.8), we have and to preserve the continuity of the first derivative of Bézier curves at nodes, by equalizing (A.9), we have where it shows the equality (10).
Conflict of Interests
The authors declare that they have no conflict of interests regarding publication of this paper.
Acknowledgment
The authors would like to thank the anonymous reviewers for their careful reading, constructive comments, and nice suggestions which have improved the paper very much.