# Bezier Curves Based Numerical Solutions of Delay Systems with Inverse Time

Authors: F. Ghomanjani; M. H. Farahi; A. Kılıçman; A. V. Kamyad; N. Pariz
Journal: Mathematical Problems in Engineering
 (2014)

---

## 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 degreen, 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.


---

## Body

## 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 derivativex˙(t) 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 t=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 att=0 are not equal. Delay systems containing inverse time are an important class of systems:
(1)x˙(t)=x(t-1),x˙(0+)=1,x˙(0-)=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 t=1, as a third degree discontinuity at time t=2, and, more generally, as a discontinuity in the (n+1)st derivative at time t=n.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:(2)x˙(t)=A(t)x(t)+C(t)(x1(t-τ1)⋯xp(t-τp))T+D(t)(x1(tf-t)⋯xp(tf-t))T+G(t)u(t),x(t)=ϕ(t),t∈[-τmax⁡,t0],
where x(t)=(x1(t)⋯xp(t))T∈Rp, u(t)=(u1(t)⋯um(t))T∈Rm are, respectively, state and control functions while ϕ(t)=(ϕ1(t)⋯ϕp(t))T is known vector function and τi’s (i=1,2,…,p) are nonnegative constant time delays, and τmax⁡=max⁡{τi,1≤i≤p}. We assume the matrices A(t)=[aij(t)]p×p, C(t)=[cij(t)]p×p, D(t)=[dij(t)]p×p, and G(t)=[gij(t)]m×m are matrix functions. We need to impose the continuity condition on x(t) 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], Nu¨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–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 Section2, 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[t0,tf] into a set of grid points such that
(3)ti=t0+ih,i=0,1,…,k,
where h=(tf-t0)/k and k is a positive integer. Let Sj=[tj-1,tj] for j=1,2,…,k. Then, for t∈Sj, delay systems containing inverse time (2) can be decomposed to the following problem:
(4)x˙j(t)=A(t)xj(t)+C(t)(x1-k11+j(t-τ1)⋯xp-k1p+j(t-τp))T+D(t)(x1k2-j(tf-t)⋯xpk2-j(tf-t))T+G(t)uj(t),xj(θ)=ϕ(θ),θ∈[-τmax⁡,t0],
where xj(t)=(x1j(t)⋯xpj(t))T and uj(t)=(u1j(t)⋯umj(t))T are, respectively, vectors of x(t) and u(t) which are considered in t∈Sj. We mention that xi-k1i+j(t-τi), 1≤i≤p, is the ith component of (x1-k11+j(t-τ1)⋯xp-k1p+j(t-τp))T where (t-τi)∈[t-k1i+j-1,t-k1i+j], and xik2-j(tf-t), 1≤i≤p, is the ith component of (x1k2-j(tf-t)⋯xpk2-j(tf-t))T where (tf-t)∈[tk2-j-1,tk2-j]. Also
(5)k1i={τih,τih∈N,([τih]+1),τih∉N,1≤i≤p,k2={tfh,tfh∈N,([tfh]+1),tfh∉N,
where [τi/h] and [tf/h] denote the integer part of τi/h and tf/h, respectively.Our strategy is to use Bézier curves to approximate the solutionsxj(t) and uj(t) by vj(t) and wj(t), respectively, where vj(t) and wj(t) are given below. Individual Bézier curves that are defined over the subintervals are joined together to form the Bézier spline curves. For j=1,2,…,k, define the Bézier polynomials of degree n that approximate, respectively, the actions of xj(t) and uj(t) over the interval [tj-1,tj] as follows:
(6)vj(t)=∑r=0narjBr,n(t-tj-1h),wj(t)=∑r=0nbrjBr,n(t-tj-1h),
where
(7)Br,n(t-tj-1h)=(nr)1hn(tj-t)n-r(t-tj-1)r
is the Bernstein polynomial of degree n over the interval [tj-1,tj], arj and brj are, respectively, p and m ordered vectors from the control points (see [22]). By substituting (6) in (4), R1,j(t) for t∈[tj-1,tj] can be defined as follows:
(8)R1,j(t)=v˙j(t)-A(t)vj(t)-C(t)(v1-k11+j(t-τ1)⋯vp-k1p+j(t-τp))T-D(t)(v1k2-j(tf-t)⋯vpk2-j(tf-t))T-G(t)wj(t).Letv(t)=∑j=1k‍χj1(t)vj(t) and w(t)=∑j=1k‍χj2(t)wj(t) where χj1(t) and χj2(t) are, respectively, characteristic function of vj(t) and wj(t) for t∈[tj-1,tj]. Beside the boundary conditions on v(t), at each node, we need to impose the continuity condition on each successive pair of vj(t) to guarantee the smoothness.Since the differential equation is of first order, the continuity ofx (or v) and its first derivative gives
(9)vj(s)(tj)=vj+1(s)(tj),s=0,1,j=1,2,…,k-1,
where vj(s)(tj) is the sth derivative vj(t) with respect to t at t=tj.Thus, the vector of control pointsarj (r=0,1,n-1,n) must satisfy (see the Appendix)
(10)anj(tj-tj-1)n=a0j+1(tj+1-tj)n,(anj-an-1j)(tj-tj-1)n-1=(a1j+1-a0j+1)(tj+1-tj)n-1.
According to the definition of the ti=t0+ih we get that tj-tj-1=h. Therefore,
(11)anj=a0j+1,(anj-an-1j)=(a1j+1-a0j+1).
One may recall that arj is a p ordered vector. This approach is called the subdivision scheme (or h-refinement in the finite element literature). This method is based on the control-point-based method.Remark 1.
By considering theC1 continuity of w, the following constraints will be added to constraints in (10):
(12)bnj(tj-tj-1)n=b0j+1(tj+1-tj)n,(bnj-bn-1j)(tj-tj-1)n-1=(b1j+1-b0j+1)(tj+1-tj)n-1,
where the so-called brj (r=0,1,n-1,n) is an m ordered vector.
Now, the residual function can be defined inSj as follows:
(13)Rj=∫tj-1tj‍∥R1,j(t)∥2dt,
where ∥·∥ is the Euclidean norm (recall that R1,j(t) is a p vector where t∈Sj).Our aim is to solve the following problem overS=⋃j=1k‍Sj:
(14)min⁡∑j=1k‍Rjs.t.anj=a0j+1,(anj-an-1j)=(a1j+1-a0j+1),j=1,2,…,k-1.
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:(15)y˙(t)=R(t)y(t)+Q(t)y(t-α)+S(t)z(t)+a(t),z˙(t)=V(t)y(t)+K(t)z(t+α)+W(t)z(t)+b(t),y(t0)=y0,y(t)=ϕ(t),t∈[-α,t0),z(tf)=z0,z(t)=ψ(t),t∈(tf,tf+α],
where y(t), z(t), a(t), b(t), ϕ(t), and ψ(t) are the vectors of appropriate dimensions. R(t), Q(t), S(t), V(t), K(t), and W(t) are the matrices of appropriate dimensions, and α is nonnegative constant time delay.
Let(16)x(t)=[y(t)Tz(tf-t)T]T,
where T is the transpose; then
(17)x˙(t)=[y˙T(t)-z˙T(tf-t)]T
satisfies that
(18)x˙(t)=A(t)x(t)+C(t)x(t-α)+D(t)x(tf-t)+u(t),t∈[t0,tf],x(t0)=x0=[y0Tz0T]T,
where
(19)A(t)=E11(2)⊗R(t)-E22(2)⊗W(tf-t),C(t)=E11(2)⊗Q(t)-E22(2)⊗K(tf-t),D(t)=E12(2)⊗S(t)-E21(2)⊗V(tf-t),u(t)=[aT(t)-bT(tf-t)]T,
where Eij(f) is the f×f matrix with 1 at its entry (i,j) and zeros elsewhere and ⊗ is Kronecker product (see, e.g., [4, 38, 39]).Remark 3.
Now, the following delay differential equation can be considered:(20)x˙(t)=f(t,x(t),x(t-τ(t,x(t)))),t≥0
with initial condition
(21)x(t)=ϕ(t),t∈[-λ,0],
where λ≡inf⁡{t-τ(t,u):t≥0,u∈R}. In the case when λ is not finite, [-λ,0] denotes the interval (-∞,0].
Furthermore, we assume that(22)τ(t,u)≥0,∀t≥0,u∈R;
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(23)Lx(t)=F(t,x(t),x(τ(t))),t0≤t≤tf
with the following initial conditions:
(24)x(k)(t0)=x0k,k=0,1,…,n-1,x(t)=ϕ(t),t≤t0,
where the differential operator L is defined by L(·)=dn(·)/dtn.


## 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 p=m=1, and the time interval is [0,1]. So, the following problem is considered:
(25)L(x(t),u(t),x(t-τ),x(1-t),dx(t)dt)=dx(t)dt-A(t)x(t)-C(t)x(t-τ)-G(t)u(t)-D(t)x(1-t)=F(t),t∈[0,1],x(t)=x0=a,t≤0,x(1)=xf=b,u(t)=u0=a1,t≤0,
where x(t)∈R, u(t)∈R, and a, b, a1 are given real numbers and A(t), C(t), G(t), D(t), and F(t) are known polynomials for t∈[0,1]. The constant time delay τ is nonnegative.Without loss of generality, we consider the interval[0,1] instead of [t0,tf] since the variable t can be changed with the new variable z by t=(tf-t0)z+t0 where z∈[0,1].Lemma 5.
For a polynomial in Bézier form(26)x(t)=∑i=0n1‍ai,n1Bi,n1(t),
we have
(27)∑i=0n1ai,n12n1+1≥∑i=0n1+1ai,n1+12n1+2≥⋯≥∑i=0n1+m1ai,n1+m12n1+m1+1⟶∫01‍x2(t)dt,m1⟶+∞,
where ai,n1+m1 is the Bézier coefficient of x(t) after degree-elevating to degree n1+m1.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 C1 continuous trajectory solution x¯, C0 continuous control solution u¯, then the approximate solution obtained by the control-point-based method converges to the exact solution (x¯,u¯) 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 Q1,N1(t) of degree N1 and Q2,N2(t) of degree N2 such that ∥diQ1,N1(t)/dti-dix¯(t)/dti∥∞≤ϵ/16, ∥diQ1,N1(t-τ)/dti-dix¯(t-τ)/dti∥∞≤ϵ/16, i=0,1, ∥Q2,N2(t)-u¯(t)∥∞≤ϵ/16, and ∥Q1,N1(1-t)-x¯(1-t)∥∞≤ϵ/16, where ∥·∥∞ stands for the L∞-norm over [0,1]. Especially, we have
(28)∥a-Q1,N1(0)∥∞≤ϵ16,∥b-Q1,N1(1)∥∞≤ϵ16,∥a1-Q2,N2(0)∥∞≤ϵ16.
In general, Q1,N1(t) and Q2,N2(t) do not satisfy the boundary conditions. After a small perturbation with linear and constant polynomials αt+β, γ, respectively, for Q1,N1(t) and Q2,N2(t), we can obtain polynomials P1,N1(t)=Q1,N1(t)+(αt+β) and P2,N2(t)=Q2,N2(t)+γ such that P1,N1(t) satisfies the boundary conditions P1,N1(0)=a, P1,N1(1)=b, and P2,N2(0)=a1. Thus, Q1,N1(0)+β=a, and Q1,N1(1)+α+β=b. By using (28), one has
(29)∥b-Q1,N1(1)∥∞=∥α+β∥∞≤ϵ16,∥a-Q1,N1(0)∥∞=∥β∥∞≤ϵ16.
Since
(30)∥α∥∞-∥β∥∞≤∥α+β∥∞≤ϵ16,
so,
(31)∥α∥∞≤ϵ16+∥β∥∞≤ϵ16+ϵ16=ϵ8.
By the time, froma1=P2,N2(0)=Q2,N2(0)+γ,
(32)∥a1-Q2,N2(0)∥∞=∥γ∥∞≤ϵ16.
Now, we have
(33)∥P1,N1(t)-x¯(t)∥∞=∥Q1,N1(t)+αt+β-x¯(t)∥∞≤∥Q1,N1(t)-x¯(t)∥∞+∥α+β∥∞≤ϵ8<ϵ5,∥dP1,N1(t)dt-dx¯(t)dt∥∞=∥dQ1,N1(t)dt+α-dx¯(t)dt∥∞≤∥dQ1,N1(t)dt-dx¯(t)dt∥∞+∥α∥∞≤3ϵ16<ϵ5,∥P2,N2(t)-u¯(t)∥∞=∥Q2,N2(t)+γ-u¯(t)∥∞≤∥Q2,N2(t)-u¯(t)∥∞+∥γ∥∞≤ϵ8<ϵ5,
so,
(34)∥P1,N1(t-τ)-x¯(t-τ)∥∞<ϵ5,∥dP1,N1(t-τ)dt-dx¯(t-τ)dt∥∞<ϵ5,∥P1,N1(1-t)-x¯(1-t)∥∞<ϵ5.
Now, let LPN(t)=L(P1,N1(t), P2,N2(t), P1,N1(t-τ), P1,N1(1-t), dP1,N1(t)/dt)=dP1,N1(t)/dt-A(t)P1,N1(t)-C(t)P1,N1(t-τ)-G(t)P2,N2(t)-D(t)P1,N1(1-t)=F(t), for every t∈[0,1]. Thus, for N≥max⁡{N1,N2}, an upper bound is found for the following residual:
(35)∥LPN(t)-F(t)∥∞=∥L(dP1,N1(t)dtP1,N1(t),P2,N2(t),P1,N1(t-τ),P1,N1(1-t),dP1,N1(t)dt)-F(t)∥∞≤∥dP1,N1(t)dt-dx¯(t)dt∥∞+∥A(t)∥∞∥P1,N1(t)-x¯(t)∥∞+∥C(t)∥∞∥P1,N1(t-τ)-x¯(t-τ)∥∞+∥G(t)∥∞∥P2,N2(t)-u¯(t)∥∞+∥D(t)∥∞∥P1,N1(1-t)-x¯(1-t)∥∞≤C1(ϵ5+ϵ5+ϵ5+ϵ5+ϵ5)=C1ϵ,
where C1=1+∥A(t)∥∞+∥C(t)∥∞+∥G(t)∥∞+∥D(t)∥∞ is a constant.
Since the residualR(PN):=LPN(t)-F(t) is a polynomial, it can be represented by a Bézier form. Therefore, we have
(36)R(PN)∶=∑i=0m1di,m1Bi,m1(t).
Then, by Lemma 5, there exists an integer M(≥N) such that, when m1>M, we have
(37)|1m1+1∑i=0m1di,m12-∫01(R(PN))2dt|<ϵ,
which gives
(38)1m1+1∑i=0m1di,m12<ϵ+∫01(R(PN))2dt≤ϵ+C12ϵ2.
Suppose x(t) and u(t) are approximated solution of (25) obtained by the control-point-based method of degree m2 (m2≥m1≥M). Let
(39)R(x(t),u(t),x(t-τ),x(1-t),dx(t)dt)=L(x(t),u(t),x(t-τ),x(1-t),dx(t)dt)-F(t)=∑i=0m2‍ci,m2Bi,m2(t),m2≥m1≥M,t∈[0,1].
Define the following norm for difference approximated solution (x(t),u(t)) and exact solution (x¯(t),u¯(t)):
(40)∥(x(t),u(t))-(x¯(t),u¯(t))∥∶=∫01∑j=01|djx(t)dtj-djx¯(t)dtj|2dt+∫01|u(0)-u¯(0)|dt.
By (40), Lemma 5, the boundary conditions x¯(0)=a=P1,N1(0)=x(0), x¯(1)=b=P1,N1(1)=x¯(1), and u¯(0)=a1=P2,N2(0)=u(0), one can show that
(41)∥(x(t),u(t))-(x¯(t),u¯(t))∥≤C((x¯(t),u¯(t),x¯(t-τ),x¯(1-t),dx¯(t)dt))∥22|x(0)-x¯(0)|+|x(1)-x¯(1)|+|u(0)-u¯(0)|+∥R((x(t),u(t),x(t-τ),x(1-t),dx(t)dt)-(dx¯(t)dtx¯(t),u¯(t),x¯(t-τ),x¯(1-t),dx¯(t)dt))∥22)=C∫01‍∑i=0m2‍(ci,m2Bi,m2(t))2dt≤Cm2+1∑i=0m2ci,m22.
The last inequality in (41) is obtained by Lemma 5, where C is a constant positive number. Now
(42)∥(x(t),u(t))-(x¯(t),u¯(t))∥≤Cm2+1∑i=0m2ci,m22≤Cm2+1∑i=0m2di,m22≤Cm1+1∑i=0m1di,m12≤C(ϵ+C12ϵ2)=ϵ1,m1≥M,
where the last inequality in (42) comes from (36). This completes the proof.


### 3.2. Subdivision

Theorem 7.
Let(x,u) 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 (x¯,u¯) and (x¯,u¯) is smooth enough so that the cubic spline T(x¯,u¯) interpolates to (x¯,u¯) and converges to (x¯,u¯) in the order O(hq), (q>2), where h is the maximal width of all subintervals, then (x,u) converges to (x¯,u¯) as h→0.Proof.
We first impose a uniform partition∏d‍=⋃i‍[ti,ti+1] on the interval [0,1] as ti=id, where d=1/(n1+1).
LetId(x¯(t),u¯(t),x¯(t-τ),x¯(1-t),dx¯(t)/dt) be the cubic spline over ∏d‍ which is interpolating to (x¯,u¯). Then, for an arbitrary small positive number ϵ>0, there exists a δ1>0 such that
(43)∥L(x¯(t),u¯(t),x¯(t-τ),x¯(1-t),dx¯(t)dt)-L(Id(x¯(t),u¯(t),x¯(t-τ),x¯(1-t),dx¯(t)dt))∥∞≤ϵ
provided that d<δ1. Let R(Id(x¯(t),u¯(t),x¯(t-τ),x¯(1-t),dx¯(t)/dt))= L(Id(x¯(t),u¯(t),x¯(t-τ),x¯(1-t),dx¯(t)/dt))-F(t) be the residual. For each subinterval [ti,ti+1], R(Id(x¯(t),u¯(t),x¯(t-τ),x¯(1-t),dx¯(t)/dt)) is a polynomial. On each interval [ti,ti+1], we impose another uniform partition ∏i,h‍=⋃j‍[ti,j,ti,j+1] as ti,j=id+jh where h=d/(m1+1), j=0,…,m1. Express R(Id(x¯(t),u¯(t),x¯(t-τ),x¯(1-t),dx¯(t)/dt)) in [ti,j-1,ti,j] as
(44)R(Id(x¯(t),u¯(t),x¯(t-τ),x¯(1-t),dx¯(t)dt))=∑p1=0lrp1i,jBp1,l(t),t∈[ti,j-1,ti,j].
By Lemma 3 in [22], there exists a δ2>0 (δ2≤δ1) such that, when h<δ2, we have
(45)|∑j=1m1(ti,j-ti,j-1)∑p1=0l(rp1i,j)2-(l+1)×∫titi+1‍R2(dx¯(t)dtId(dx¯(t)dtx¯(t),u¯(t),x¯(t-τ),∑p1=0l(rp1i,j)2-(l+1)x¯(1-t),dx¯(t)dt))|dt≤ϵd.
Thus,
(46)|∑i=1n1∑j=1m1(ti,j-ti,j-1)∑p1=0l(rp1i,j)2-(l+1)∫01R2(∑p1=0l(rp1i,j)2-(l+1)Id(∑p1=0l(rp1i,j)2-(l+1)x¯(t),u¯(t),x¯(t-τ),∑p1=0l(rp1i,j)2-(l+1)x¯(1-t),dx¯(t)dt))|dt≤ϵ,
or
(47)∑i=1n1∑j=1m1(ti,j-ti,j-1)∑p1=0l(rp1i,j)2<(l+1)∫01R2(dx¯(t)dtId(dx¯(t)dtx¯(t),u¯(t),x¯(t-τ),x¯(1-t),dx¯(t)dt))dt+ϵ<(l+1)ϵ2+ϵ.
Now combining the partitions ∏d‍ and all ∏i,h‍ gives a denser partition with the length h for each subinterval. Suppose (x(t),u(t)) is the approximate solution by the control-point-based method with respect to this partition, and denote the residual over [ti,j-1,ti,j] by
(48)R(x(t),u(t),x(t-τ),x(1-t),dx(t)dt)=L(x(t),u(t),x(t-τ),x(1-t),dx(t)dt)-F(t)=∑p1=0lcp1i,jBp1,l(t).
Define the following norm for difference approximate solution (x(t),u(t)) and exact solution (x¯(t),u¯(t)):
(49)∥(x(t),u(t))-(x¯(t),u¯(t))∥∶=∑i=1n1∑j=1m1∫ti,j-1ti,j|x(t)-x¯(t)|2dt+∑i=1n1∑j=1m1∫ti,j-1ti,j|dx(t)dt-dx¯(t)dt|2dt+∑i=1n1∑j=1m1∫ti,j-1ti,j|u(0)-u¯(0)|dt.
Then, there is a constant C such that
(50)∥(x(t),u(t))-(x¯(t),u¯(t))∥≤C∥R((x(t),u(t),x(t-τ),x(1-t),dx(t)dt)-(x¯(t),u¯(t),x¯(t-τ),x¯(1-t),dx¯(t)dt))∥2≤Cl+1∑i=1n1∑j=1m1(ti,j-ti,j-1)∑p1=0l(cp1i,j)2;
the last inequality in (50) is obtained by Lemma 5. It can be shown that
(51)Cl+1∑i=1n1∑j=1m1(ti,j-ti,j-1)∑p1=0l(cp1i,j)2≤Cl+1∑i=1n1∑j=1m1(ti,j-ti,j-1)∑p1=0l(rp1i,j)2≤C(ϵ2+ϵl+1)=ϵ2.
By Lemma 3 in [22], we conclude that the approximate solution converges to the exact solution in order o(hq), (q>2). This completes the proof.


## 3.1. Degree Raising

Theorem 6.
If the problem (25) with inverse time in state has a unique C1 continuous trajectory solution x¯, C0 continuous control solution u¯, then the approximate solution obtained by the control-point-based method converges to the exact solution (x¯,u¯) 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 Q1,N1(t) of degree N1 and Q2,N2(t) of degree N2 such that ∥diQ1,N1(t)/dti-dix¯(t)/dti∥∞≤ϵ/16, ∥diQ1,N1(t-τ)/dti-dix¯(t-τ)/dti∥∞≤ϵ/16, i=0,1, ∥Q2,N2(t)-u¯(t)∥∞≤ϵ/16, and ∥Q1,N1(1-t)-x¯(1-t)∥∞≤ϵ/16, where ∥·∥∞ stands for the L∞-norm over [0,1]. Especially, we have
(28)∥a-Q1,N1(0)∥∞≤ϵ16,∥b-Q1,N1(1)∥∞≤ϵ16,∥a1-Q2,N2(0)∥∞≤ϵ16.
In general, Q1,N1(t) and Q2,N2(t) do not satisfy the boundary conditions. After a small perturbation with linear and constant polynomials αt+β, γ, respectively, for Q1,N1(t) and Q2,N2(t), we can obtain polynomials P1,N1(t)=Q1,N1(t)+(αt+β) and P2,N2(t)=Q2,N2(t)+γ such that P1,N1(t) satisfies the boundary conditions P1,N1(0)=a, P1,N1(1)=b, and P2,N2(0)=a1. Thus, Q1,N1(0)+β=a, and Q1,N1(1)+α+β=b. By using (28), one has
(29)∥b-Q1,N1(1)∥∞=∥α+β∥∞≤ϵ16,∥a-Q1,N1(0)∥∞=∥β∥∞≤ϵ16.
Since
(30)∥α∥∞-∥β∥∞≤∥α+β∥∞≤ϵ16,
so,
(31)∥α∥∞≤ϵ16+∥β∥∞≤ϵ16+ϵ16=ϵ8.
By the time, froma1=P2,N2(0)=Q2,N2(0)+γ,
(32)∥a1-Q2,N2(0)∥∞=∥γ∥∞≤ϵ16.
Now, we have
(33)∥P1,N1(t)-x¯(t)∥∞=∥Q1,N1(t)+αt+β-x¯(t)∥∞≤∥Q1,N1(t)-x¯(t)∥∞+∥α+β∥∞≤ϵ8<ϵ5,∥dP1,N1(t)dt-dx¯(t)dt∥∞=∥dQ1,N1(t)dt+α-dx¯(t)dt∥∞≤∥dQ1,N1(t)dt-dx¯(t)dt∥∞+∥α∥∞≤3ϵ16<ϵ5,∥P2,N2(t)-u¯(t)∥∞=∥Q2,N2(t)+γ-u¯(t)∥∞≤∥Q2,N2(t)-u¯(t)∥∞+∥γ∥∞≤ϵ8<ϵ5,
so,
(34)∥P1,N1(t-τ)-x¯(t-τ)∥∞<ϵ5,∥dP1,N1(t-τ)dt-dx¯(t-τ)dt∥∞<ϵ5,∥P1,N1(1-t)-x¯(1-t)∥∞<ϵ5.
Now, let LPN(t)=L(P1,N1(t), P2,N2(t), P1,N1(t-τ), P1,N1(1-t), dP1,N1(t)/dt)=dP1,N1(t)/dt-A(t)P1,N1(t)-C(t)P1,N1(t-τ)-G(t)P2,N2(t)-D(t)P1,N1(1-t)=F(t), for every t∈[0,1]. Thus, for N≥max⁡{N1,N2}, an upper bound is found for the following residual:
(35)∥LPN(t)-F(t)∥∞=∥L(dP1,N1(t)dtP1,N1(t),P2,N2(t),P1,N1(t-τ),P1,N1(1-t),dP1,N1(t)dt)-F(t)∥∞≤∥dP1,N1(t)dt-dx¯(t)dt∥∞+∥A(t)∥∞∥P1,N1(t)-x¯(t)∥∞+∥C(t)∥∞∥P1,N1(t-τ)-x¯(t-τ)∥∞+∥G(t)∥∞∥P2,N2(t)-u¯(t)∥∞+∥D(t)∥∞∥P1,N1(1-t)-x¯(1-t)∥∞≤C1(ϵ5+ϵ5+ϵ5+ϵ5+ϵ5)=C1ϵ,
where C1=1+∥A(t)∥∞+∥C(t)∥∞+∥G(t)∥∞+∥D(t)∥∞ is a constant.
Since the residualR(PN):=LPN(t)-F(t) is a polynomial, it can be represented by a Bézier form. Therefore, we have
(36)R(PN)∶=∑i=0m1di,m1Bi,m1(t).
Then, by Lemma 5, there exists an integer M(≥N) such that, when m1>M, we have
(37)|1m1+1∑i=0m1di,m12-∫01(R(PN))2dt|<ϵ,
which gives
(38)1m1+1∑i=0m1di,m12<ϵ+∫01(R(PN))2dt≤ϵ+C12ϵ2.
Suppose x(t) and u(t) are approximated solution of (25) obtained by the control-point-based method of degree m2 (m2≥m1≥M). Let
(39)R(x(t),u(t),x(t-τ),x(1-t),dx(t)dt)=L(x(t),u(t),x(t-τ),x(1-t),dx(t)dt)-F(t)=∑i=0m2‍ci,m2Bi,m2(t),m2≥m1≥M,t∈[0,1].
Define the following norm for difference approximated solution (x(t),u(t)) and exact solution (x¯(t),u¯(t)):
(40)∥(x(t),u(t))-(x¯(t),u¯(t))∥∶=∫01∑j=01|djx(t)dtj-djx¯(t)dtj|2dt+∫01|u(0)-u¯(0)|dt.
By (40), Lemma 5, the boundary conditions x¯(0)=a=P1,N1(0)=x(0), x¯(1)=b=P1,N1(1)=x¯(1), and u¯(0)=a1=P2,N2(0)=u(0), one can show that
(41)∥(x(t),u(t))-(x¯(t),u¯(t))∥≤C((x¯(t),u¯(t),x¯(t-τ),x¯(1-t),dx¯(t)dt))∥22|x(0)-x¯(0)|+|x(1)-x¯(1)|+|u(0)-u¯(0)|+∥R((x(t),u(t),x(t-τ),x(1-t),dx(t)dt)-(dx¯(t)dtx¯(t),u¯(t),x¯(t-τ),x¯(1-t),dx¯(t)dt))∥22)=C∫01‍∑i=0m2‍(ci,m2Bi,m2(t))2dt≤Cm2+1∑i=0m2ci,m22.
The last inequality in (41) is obtained by Lemma 5, where C is a constant positive number. Now
(42)∥(x(t),u(t))-(x¯(t),u¯(t))∥≤Cm2+1∑i=0m2ci,m22≤Cm2+1∑i=0m2di,m22≤Cm1+1∑i=0m1di,m12≤C(ϵ+C12ϵ2)=ϵ1,m1≥M,
where the last inequality in (42) comes from (36). This completes the proof.


## 3.2. Subdivision

Theorem 7.
Let(x,u) 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 (x¯,u¯) and (x¯,u¯) is smooth enough so that the cubic spline T(x¯,u¯) interpolates to (x¯,u¯) and converges to (x¯,u¯) in the order O(hq), (q>2), where h is the maximal width of all subintervals, then (x,u) converges to (x¯,u¯) as h→0.Proof.
We first impose a uniform partition∏d‍=⋃i‍[ti,ti+1] on the interval [0,1] as ti=id, where d=1/(n1+1).
LetId(x¯(t),u¯(t),x¯(t-τ),x¯(1-t),dx¯(t)/dt) be the cubic spline over ∏d‍ which is interpolating to (x¯,u¯). Then, for an arbitrary small positive number ϵ>0, there exists a δ1>0 such that
(43)∥L(x¯(t),u¯(t),x¯(t-τ),x¯(1-t),dx¯(t)dt)-L(Id(x¯(t),u¯(t),x¯(t-τ),x¯(1-t),dx¯(t)dt))∥∞≤ϵ
provided that d<δ1. Let R(Id(x¯(t),u¯(t),x¯(t-τ),x¯(1-t),dx¯(t)/dt))= L(Id(x¯(t),u¯(t),x¯(t-τ),x¯(1-t),dx¯(t)/dt))-F(t) be the residual. For each subinterval [ti,ti+1], R(Id(x¯(t),u¯(t),x¯(t-τ),x¯(1-t),dx¯(t)/dt)) is a polynomial. On each interval [ti,ti+1], we impose another uniform partition ∏i,h‍=⋃j‍[ti,j,ti,j+1] as ti,j=id+jh where h=d/(m1+1), j=0,…,m1. Express R(Id(x¯(t),u¯(t),x¯(t-τ),x¯(1-t),dx¯(t)/dt)) in [ti,j-1,ti,j] as
(44)R(Id(x¯(t),u¯(t),x¯(t-τ),x¯(1-t),dx¯(t)dt))=∑p1=0lrp1i,jBp1,l(t),t∈[ti,j-1,ti,j].
By Lemma 3 in [22], there exists a δ2>0 (δ2≤δ1) such that, when h<δ2, we have
(45)|∑j=1m1(ti,j-ti,j-1)∑p1=0l(rp1i,j)2-(l+1)×∫titi+1‍R2(dx¯(t)dtId(dx¯(t)dtx¯(t),u¯(t),x¯(t-τ),∑p1=0l(rp1i,j)2-(l+1)x¯(1-t),dx¯(t)dt))|dt≤ϵd.
Thus,
(46)|∑i=1n1∑j=1m1(ti,j-ti,j-1)∑p1=0l(rp1i,j)2-(l+1)∫01R2(∑p1=0l(rp1i,j)2-(l+1)Id(∑p1=0l(rp1i,j)2-(l+1)x¯(t),u¯(t),x¯(t-τ),∑p1=0l(rp1i,j)2-(l+1)x¯(1-t),dx¯(t)dt))|dt≤ϵ,
or
(47)∑i=1n1∑j=1m1(ti,j-ti,j-1)∑p1=0l(rp1i,j)2<(l+1)∫01R2(dx¯(t)dtId(dx¯(t)dtx¯(t),u¯(t),x¯(t-τ),x¯(1-t),dx¯(t)dt))dt+ϵ<(l+1)ϵ2+ϵ.
Now combining the partitions ∏d‍ and all ∏i,h‍ gives a denser partition with the length h for each subinterval. Suppose (x(t),u(t)) is the approximate solution by the control-point-based method with respect to this partition, and denote the residual over [ti,j-1,ti,j] by
(48)R(x(t),u(t),x(t-τ),x(1-t),dx(t)dt)=L(x(t),u(t),x(t-τ),x(1-t),dx(t)dt)-F(t)=∑p1=0lcp1i,jBp1,l(t).
Define the following norm for difference approximate solution (x(t),u(t)) and exact solution (x¯(t),u¯(t)):
(49)∥(x(t),u(t))-(x¯(t),u¯(t))∥∶=∑i=1n1∑j=1m1∫ti,j-1ti,j|x(t)-x¯(t)|2dt+∑i=1n1∑j=1m1∫ti,j-1ti,j|dx(t)dt-dx¯(t)dt|2dt+∑i=1n1∑j=1m1∫ti,j-1ti,j|u(0)-u¯(0)|dt.
Then, there is a constant C such that
(50)∥(x(t),u(t))-(x¯(t),u¯(t))∥≤C∥R((x(t),u(t),x(t-τ),x(1-t),dx(t)dt)-(x¯(t),u¯(t),x¯(t-τ),x¯(1-t),dx¯(t)dt))∥2≤Cl+1∑i=1n1∑j=1m1(ti,j-ti,j-1)∑p1=0l(cp1i,j)2;
the last inequality in (50) is obtained by Lemma 5. It can be shown that
(51)Cl+1∑i=1n1∑j=1m1(ti,j-ti,j-1)∑p1=0l(cp1i,j)2≤Cl+1∑i=1n1∑j=1m1(ti,j-ti,j-1)∑p1=0l(rp1i,j)2≤C(ϵ2+ϵl+1)=ϵ2.
By Lemma 3 in [22], we conclude that the approximate solution converges to the exact solution in order o(hq), (q>2). This completes the proof.


## 4. Numerical Examples

Applying the presented method, in Examples1, 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])
(52)x˙(t)=[t2+1-t20-9]x(t)+[1-190]x(t-13)+[10-11]x(1-t)+[4t+38t+15]u(t),ϕ(t)=[t2-1t2+1],t∈[-13,0],
where we have the following exact solution:
(53)x(t)=[x1(t)x2(t)]T=[t2-1t2+1]T.Letu(t)=1. Then, by (14) and choosing n=3, k=6 we have the approximate solution x(t)=[x1(t)x2(t)]T(54)x1(t)={-1.000000001+8.333333337×10-9t+0.9999999669t2+10-7t3,0≤t≤16,-0.9999999988+8.13333333×10-9t+0.9999999829t2,16≤t≤13,-0.9999999997+2.00×10-10t+t2,13≤t≤12,-0.9999999927-2.202222223×10-8t+1.000000017t2,12≤t≤23,-0.9999999902-1.504444443×10-8t+0.9999999963t2+10-8t3,23≤t≤56,-1.000000032+1.120666667×10-7t+0.9999998702t2+5×10-8t3,56≤t≤1,x2(t)={1.000000001+0.000011825t+0.9996447669t2+0.0023693t3,0≤t≤16,1.000000001+0.00001180813339t+0.9996447663t2+0.0023695t3,16≤t≤13,0.9999999645+0.00001211131104t+0.9996439669t2+0.0023702t3,13≤t≤12,1.000000063+0.00001151408882t+0.9996452169t2+0.0023693t3,12≤t≤23,0.9581187057+0.1594325022t+0.8040813829t2+0.0783674t3,23≤t≤56,0.9581181451+0.1594344559t+0.8040791002t2+0.0783683t3,56≤t≤1.The graphs of approximate trajectories are shown in Figures1 and 2.Figure 1
The graph of approximated trajectoryx1(t) for Example 1.Figure 2
The graph of approximated trajectoryx2(t) for Example 1.Example 2.
Consider the boundary value problem described by (see [4])
(2)y˙(t)=16ty(t-14)-16z(t)+8t2+17t+16,z˙(t)=64ty(t)-64z(t+14)+51t2+76t+65,y(t)=t2-1,-14≤t≤0,z(t)=t3+1,1≤t≤54.
From (18), we have (see [4])
(56)x˙(t)=[16t0064]x(t-14)+[0-1664t-640]x(1-t)+[8t2+17t+16-51t2+178t-62],ϕ(t)=[t2-1-t3+3t2-3t+1],t∈[-14,0],
where x(t)=[x1(t)x2(t)]T=[y(t)z(1-t)]T, and we have the following exact solution:
(57)x(t)=[x1(t)x2(t)]T=[t2-1t3+1]T.
Let u(t)=1. Then, by (14) and choosing n=3, k=4 we have the approximate solution x(t)=[x1(t)x2(t)]T:
(58)x1(t)={t2-1,0≤t≤14,t2-1,14≤t≤12,t2-1,12≤t≤34,-1.000000006+2.0625×10-8t+0.999999975t2+10-8t3,34≤t≤1,x2(t)=t3+1.
The graphs of approximate trajectories are shown in Figures 3 and 4.Figure 3
The graph of approximated trajectoryx1(t) for Example 2.Figure 4
The graph of approximated trajectoryx2(t) for Example 2.Example 3.
Consider the time-varying delay system described by (see [42])
(59)[x˙1(t)x˙2(t)]=[01-25-5t][x1(t-14)x2(t-14)]+[01],[x1(t)x2(t)]=[00],t∈[-14,0].
The exact solutions are [42](60)x1(t)={0,t∈[0,14],132-14t+12t2,t∈[14,12],132-1996t+316t2+58t3-512t4,t∈[12,34],-964132768+3739124576t-31831024t2+70652304t3-135384t4-8596t5+518t6,t∈[34,1],x2(t)={t,t∈[0,14],-5384+t+58t2-53t3,t∈[14,12],7751536-178t+1295192t2-11524t3-7532t4+53t5,t∈[12,34],87997132120-10511024t-9575549152t2+215151536t3-553253072t4+33596t5+2125576t6-2521t7,t∈[34,1].
Here, this problem is solved by choosingk=8 and n=3. the following approximate solutions x1(t) and x2(t) are found. In Tables 1 and 2, exact, numerical results of this method and obtained results in [42] are shown, respectively:(61)x1(t)={-0.001524977445t+0.04981148910t2-0.3456171465t3,t∈[0,18],-0.002668294207+0.06251408351t-0.4625009986t2+1.020549487t3,t∈[18,14],0.006613889339-0.04887212012t-0.01695618114t2+0.4264897281t3,t∈[14,38],0.01307452454-0.1005572014t+0.1208707015t2+0.303976944t3,t∈[38,12],0.1271590458-0.7850643303t+1.489884961t2-0.608699230t3,t∈[12,58],0.06579667219-0.4905249419t+1.018621948t2-0.357358960t3,t∈[58,34],0.3247255416-1.526240419t+2.399575918t2-0.9711162800t3,t∈[34,78],0.6384881122-2.601997790t+3.629012898t2-1.439473220t3,t∈[78,1],,x2(t)={1.003041110t-0.09123330000t2+0.6082219700t3,t∈[0,18],0.003925049727+0.9088399072t+0.6623763650t2-1.401403820t3,t∈[18,14],0.003925049727+0.9088399072t+0.6623763650t2-1.401403820t3,t∈[14,38],-0.02462216250+1.075221794t+0.4666746125t2-1.558091100t3,t∈[38,12],0.3991598156-1.467470069t+5.552058325t2-4.948346900t3,t∈[12,58],0.00006281562500+0.4481955219t+2.486993388t2-3.313645600t3,t∈[58,34],-1.159405308+5.086068009t-3.696836582t2-0.5652767300t3,t∈[34,78],-5.634050302+20.42770799t-21.23013942t2+6.114076730t3,t∈[78,1].Table 1
Exact and estimated values ofx1(t) for Example 3.
t
Exactx1(t)
Presentx1(t)
x
1
(
t
) in [42]
0.00
0.000000
0.0000000000000000
-
0.000088
0.05
0.000000
0.0000050777072400
-
0.000046
0.10
0.000000
0.0000000000000000
0.000021
0.15
0.000000
-
0.000253099630375
0.000083
0.20
0.000000
-
0.000501121553000
-
0.000128
0.25
0.000000
-
1.06875
×
1
0
-
11
-
0.000024
0.30
0.001250
0.0019414196591000
0.001400
0.35
0.005000
0.0057172621996375
0.004987
0.40
0.011250
0.0116454806360000
0.011157
0.45
0.020000
0.0199999999857500
0.019968
0.50
0.031250
0.0310107172150000
0.031304
0.55
0.044971
0.04479153044625000
0.045021
0.60
0.061000
0.06099999990000000
0.060991
0.65
0.079086
0.0791835285950000
0.079044
0.70
0.098917
0.0989798441000000
0.098901
0.75
0.120117
0.1201170006000000
0.120143
0.80
0.142244
0.1422502585600000
0.142266
0.85
0.164728
0.1647280007500000
0.164710
0.90
0.186819
0.1868145712000000
0.186803
0.95
0.207606
0.2076060001475000
0.207623
1.00
0.226030
0.2260300002000000
0.226030Table 2
Exact and estimated values ofx2(t) for Example 3.
t
Exactx2(t)
Presentx2(t)
x
2
(
t
) in [42]
0.00
0.000000
0.000000000000000
0.001169
0.05
0.050000
0.050000000000000
0.049923
0.10
0.100000
0.099999999970000
0.100294
0.15
0.150000
0.150424766127000
0.149740
0.20
0.200000
0.200976855207000
0.199902
0.25
0.250000
0.250636614672500
0.250170
0.30
0.298229
0.298229000000000
0.298294
0.35
0.345083
0.342083000067500
0.342098
0.40
0.380313
0.380416662700000
0.380186
0.45
0.411667
0.411748202343750
0.411593
0.50
0.434896
0.434896000125000
0.435025
0.55
0.448306
0.4482677054750000
0.448326
0.60
0.448532
0.4485758408000000
0.448483
0.65
0.432078
0.432134688390000
0.432080
0.70
0.395846
0.395846000275000
0.395868
0.75
0.337199
0.337199000906250
0.337171
0.80
0.254052
0.254052000960000
0.254038
0.85
0.145303
0.145637497343750
0.145354
0.90
0.011316
0.011635894970000
0.011295
0.95
-
0.145872
-
0.14587200166625
-
0.145924
1.00
-
0.322405
-
0.32240500200000
-
0.322386Example 4.
Consider the following system described by (see [40])
(62)x˙(t)=8t+1x(t-(t2+12)),t≥0,x(t)=(t+1)2,t∈[-12,0].
Analytic solution of the initial value problem (IVP) is x(t)=(t+1)2. By choosing k=1 and n=16 (degree raising), we obtain the following solution:
(63)x(t)=1+0.2018032795×10-4t12-0.01572515756t7-0.008572702573t5+0.01741959010t6-0.0001540665901t11-0.1834453040×10-5t13+1.101285958×10-7t14+0.008669328894t8+1.999552507t+6.306939519×10-11t16-3.928281389×10-9t15-0.003213347229t9+0.9993525856t2+0.0008342736689t10+0.004438985657t3-0.002620448442t4.
In Table 3, exact and presented methods are shown, respectively.Table 3
Exact and estimated values ofx(t) for Example 4.
t
Exact
Presented method
0.5
2.25
2.24991525903163
1.0
4
4.00000000000000
1.5
6.25
6.24995700258759
2
9
9.00000000128046Example 5.
Consider the following system described by (see [40])
(64)x˙(t)=x(t-1-1t+1),t≥0,x(t)={23(t+2),-2≤t≤-0.5,1,-0.5≤t≤0,
where the exact solution is x(t)=1+(2/3)t+t3/3-(2/3)log⁡(t+1) on [0,1] and x(t)=1-(2/3)log⁡2+t on [1,2]. By choosing k=1 and n=7 (degree raising), we obtain the following solution:
(65)x(t)=1+5.424427795t5-1.611981446t6-2.552250886t2+7.963903747t3+0.3574277875t-9.236517482t4+0.1928923646t7.
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.Table 4
Exact and estimated values ofx(t) for Example 5.
t
Exact
Presented method
Method in [40]
Error of presented method
Error of the method in [40]
0.5
1.10468992792789
1.10468992817860
1.1451
2.50709000000000
×
1
0
-
10
1.2232
×
1
0
-
3
1.0
1.53790187962670
1.53790188062000
1.5361
9.93297
×
1
0
-
10
1.7685
×
1
0
-
3
1.4
1.93790187962670
1.93768171138582
1.9361
0.220168240883
×
1
0
-
3
1.7685
×
1
0
-
3
1.5
2.03790187962670
2.03790188078453
2.0362
1.157827
×
1
0
-
9
1.6125
×
1
0
-
3
2.0
2.53790187962670
2.53790188032000
2.5870
6.93297
×
1
0
-
10
4.9096
×
1
0
-
2Example 6.
Consider the following system described by (see [40])
(66)x˙(t)=-x(t-τ(t)),t∈[0,2],x(0)=1,τ(t)≡{t-2+4-2t,0≤t≤2,0,t>2.
The solution of this problem is
(67)x(t)={(t-2)24,0≤t≤2,0,t>2.
By choosing k=1 and n=7 (degree raising), we obtain the following solution:
(68)x(t)=1-1.000000002t+3.207267830×10-9t6+0.2500000112×t2-3.416339151×10-10t7-1.204800000×10-8t5-2.304000000×10-8t3+2.296000000×10-8t4.
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.Table 5
Exact and estimated values ofx(t) for Example 6.
t
Exact
Presented method
Method in [40]
Error of presented method
Error of the method in [40]
1.0
0.25
0.250000000017634
0.250013
1.7634
×
1
0
-
11
1.28346
×
1
0
-
5
2.0
0.0
0.0
5.26486
×
1
0
-
7
0.0
5.26486
×
1
0
-
7Example 7.
Consider the following LDDE described by(69)d3x(t)dt3=-x(t)-x(t-0.3)+e-t+0.3,0≤t≤1,
with the initial conditions
(70)x(0)=1,dx(0)dt=-1,d2x(0)dt2=1,x(t)=e-t,hhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhht≤0,
where the exact solution of this example is x(t)=e-t. Here, this problem is solved by choosing k=10 and n=3. The graph of error is shown in Figure 5, and the following approximate solution x(t) is found:(71)x(t)={1-t+0.5t2-0.172928t3,t∈[0,0.1],0.9999767558-0.999302674t+0.49302674t2-0.1496838t3,t∈[0.1,0.2],0.9998081522-0.996773576t+0.48038104t2-0.1286073t3,t∈[0.2,0.3],0.9992953871-0.991645877t+0.46328853t2-0.1096154t3,t∈[0.3,0.4],0.9982244623-0.983613861t+0.44320829t2-0.0928817t3,t∈[0.4,0.5],0.9964070488-0.972709334t+0.42139914t2-0.0783422t3,t∈[0.5,0.6],0.9937493164-0.9594206t+0.39925114t2-0.0660377t3,t∈[0.6,0.7],0.9903114379-0.944686777t+0.37820273t2-0.0560146t3,t∈[0.7,0.8],0.9863822587-0.929952279t+0.35978451t2-0.0483403t3,t∈[0.8,0.9],0.9825547252-0.917193744t+0.34560826t2-0.0430898t3,t∈[0.9,1].Figure 5
The graph of error for Example7.Example 8.
Consider the second-order linear decay differential equation:(72)x¨(t)=34x(t)+x(t2)-t2+2,0≤t≤1,x(0)=0,x˙(0)=0.
The exact solution of this problem is x(t)=t2. Here, this problem is solved by choosing k=1 and n=7. The following approximate solution x(t) is found. 							(73)x(t)=1.882848000t2-5.072623999t3+15.56400000t4-28.14240000t5+30.84000000t6-t9+7t8-0.06182400000t-20.01000000t7.
In Table 6, exact, numerical results of this method, error of presented method, and error of the method in [43] are shown, respectively.Table 6
Exact and estimated values ofx(t) for Example 8.
t
Exact
Presented method
Error Of presented method
Error of the method in [43]
0.2
0.04
0.0400000000049152
4.9152
×
1
0
-
12
1.73
×
1
0
-
6
0.4
0.16
0.1600000000193540
1.9354
×
1
0
-
11
1.10
×
1
0
-
5
0.6
0.36
0.3600000000221180
2.2118
×
1
0
-
11
1.26
×
1
0
-
4
0.8
0.64
0.6400000000073730
7.373
×
1
0
-
12
7.07
×
1
0
-
4


## 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.



---

*Source: 602641-2014-02-27.xml*