PECULIAR SPLINE COLLOCATION METHOD FOR SOLVING ROUGH AND STIFF DELAY DIFFERENTIAL PROBLEMS

. As well known, solutions of delay diﬀerential equations (DDEs) are characterized by low regularity. In particular solutions of neutral delay diﬀeren-tial equations (NDDEs) frequently exhibit discontinuities in the ﬁrst derivative so that the diﬀerential problems become rough. The aim of this paper is to approximate the solutions of such rough delay diﬀerential problems by means of a peculiar deﬁcient spline collocation method. Signiﬁcant numerical examples are provided to enlighten the features of the proposed method.


INTRODUCTION
We consider the following neutral delay differential equation (NDDE): it can be proved that the problem (1) has a unique solution The proof is analogous to that reported in [2] and we omit it here.
It is well known that discontinuities in low-order derivatives are almost always present in the solutions of delay differential problems: generally there is at least a discontinuity in the first derivative of solution y at the initial point τ 0 =a, even though f , ϕ are analytical vectorial functions in their arguments.Discontinuities for DDEs are serious because they propagate.If the solution has a discontinuity in a derivative somewhere, then discontinuities appear in the rest of the interval at a spacing given by the delay, that is at points {τ j = a + jτ }, (j = 1, 2, ..., M ) with τ M = b.However, under general conditions, if the solution is smooth enough, the propagated discontinuities are smoothed and they appear in derivatives of increasing order; once the orders are high enough, the discontinuities do not interfere with the numerical method any longer.
On the contrary, for NDDEs discontinuity in the first derivative at the initial point a very often remains even when j increases.In this case the conditions of existence and uniqueness of the solution still hold in subintervals I j := (τ j , τ j+1 ) (j = 0, 1, ..., M − 1).Even discontinuities in the history function ϕ propagate in the same way and also in this case the conditions of existence and uniqueness of the solution still hold in convenient subintervals.In particular, NDDEs often present discontinuities in the first derivative characterized by variation of sign with a very high module or a sudden increasing of module.
We remark that in order to preserve the accuracy of the numerical approximation, the ends of those convenient subintervals have to coincide with some of the discretization mesh points.Obviously the points where the discontinuities in derivatives happen are known in advance, even if the kind of irregularity can not be evident.
The aim of this work is to present a method for the numerical approximation of solutions of NDDEs so that the order of regularity of the analytical solution is preserved when the problem is rough.
We remark that sometime very stiff DDEs problems can be viewed numerically as rough because of the discretization step, even though it is very small.We will show that our method works efficiently even for these cases as well as for the extreme case when the solution itself is discontinuous.
In Section 2 we present and discuss the proposed numerical method; in Section 3 we present the algorithm; in Section 4 we report some numerical examples showing significant applications of the algorithm.

THE NUMERICAL METHOD
Basically the proposed method was already presented in [4], [5], [6], but here we present, in particular, some grid modifications in order to tackle efficiently with rough problems and interesting and significant applications.
2.1.Basic theoretical results.The main features of the proposed numerical method are the following.
As a class of approximating functions, we choose the class of deficient splines of degree m ≥ 2 and deficiency 2, denoted by s : [a, b] → R n , that is s belongs to the vectorial space of splines of degree m with s ∈ C m−2 .
Without loss of generality, we refer to the generic interval I j := [τ j , τ j+1 ] (j = 0, 1, ..., M −1) and suppose that unique discontinuity point is t * j such that t * j ∈ I j ; then in a subinterval [t k , t k+1 ] ⊂ [τ j , τ j+1 ] , (k = 0, 1, ..., 2N − 1), where , the spline function s is defined as: We choose to determine the vectors of coefficients a k , b k by the following system of collocation conditions: taking into account that: Thus our numerical model (step by step) is reduced to compute the solution of system (2) of 2 vectorial equations with 2 vectors unknowns, which can be linear or non-linear.Some theoretical results referring to existence and uniqueness of numerical solution, consistency and convergence are analogous to those in [5], where we studied the scalar case only.
In particular existence and uniqueness can be easily proved as in Theorem 1 in [5], which is an extension to deficient spline case of results in [1].
Moreover, experimentally, we have that for m < 4, consistency and convergence are guaranteed, at contrary for m ≥ 4 zero-stability is not assured; therefore divergence happens, as it was already pointed out in [10], for the classic splines.
Such an upper bound on m is not a disadvantage for our numerical method, because our aim is to approximate solutions characterized by low regularity, therefore our choice is m = 2, so that s ∈ C 0 , in such a way that the numerical approximation belongs to the same class of regularity of solution of considered rough problems.
We remind that for m = 2 our method guarantees a second order convergence [5].
2.2.Rough problems.As discretization over discontinuities in low-order derivatives can make the numerical approximation less accurate with a consequent loss of the order of convergence, we need to use convenient mesh points.This is the main point of the proposed method: here we give details about the choice of the mesh according to different types of discontinuities.
1. Discontinuities appear at τ j only, j = 0, 1, 2, ..., and then they are propagated and smoothed; in particular at τ 1 there is a discontinuity in y (k) with k ≥ 2. As these problems are not rough, our numerical method tackles with this kind of discontinuities using a convenient subgrid only in the first integration interval [a, a + h] where h is the chosen integration step.In general this is the case of DDEs. 2. Discontinuities appear at τ j only, j = 0, 1, 2, ..., but the propagated discontinuities are not smoothed; this is the case of NDDEs.Then the method: a) detects the discontinuities at least in the first derivative by means of a control of the first derivative of the approximating numerical solution (always available together with the solution itself); b) uses a subgrid in the first integration interval after τ j , if the discontinuity is severe.3. Discontinuities at least in the first derivative appear in the integration interval when history itself exhibits discontinuities in low-order derivatives.The history is the solution prior the initial point and its discontinuities must also be taken into account as they propagate into the interval of integration because of the delays.This kind of discontinuities are handled like discontinuities at known points during the integration, as they happen at each successive interval given by the delays.Therefore they can be numerically treated as in the previous case, with the caution of making the discontinuity points coincide with convenient mesh points.4. Discontinuities in the solution itself complicate matters; also in this case discontinuities happen at points known in advance, therefore the numerical treatment is analogous to the previous case.5. When the solution behavior is stiff, the solution can be viewed numerically discontinuous in the first derivative, even though it is continuous, because of the integration step when it is not small enough; in this case our collocation method provides very good results as it can use deficient spline functions of class C 0 which approximate efficiently the stiff behavior of solution.Indeed the main feature of our numerical model is that it builds a piecewise polynomial function approximating the analytical solution, which exhibits similarly a piecewise behavior.
Moreover we point out that our method presents some main advantages: -It provides an accurate global approximation to the solution of (1), which allows to take into account the historical behavior of the approximating spline.This is really relevant to the "delayed" nature of the considered equation.
-The step-size can be changed at any step, if necessary, without added complications.
-The use of deficient splines allows to relax the regularity of the function at the connection points without any loss of order of convergence.On the other hand, the deficient splines functional class takes into account and follows the low regular functions which are solutions of the problem.
-The numerical treatment of rough problems by means of spline functions belonging to class C 0 , performs better than other methods.
-As at each step, our method provides not only the global solution in the considered interval, but also its global first derivative, by which we can adjust the mesh grid to our problem and the choice is easy and straightforward.

ALGORITHM
Our purpose is just to show that a simple, easy to implement, very fast algorithm in MATLAB can provide efficiently a numerical approximation of the solution of DDEs when the solution exhibits very low regularity (in particular for NDDEs).Actually what is a difficulty for other methods, is an advantage for our method.
We do not propose any numerical method for DDEs in competition with the many general purpose software packages, which have been developed in the recent years.We remind here that up-to-date links to those programs can be found in the web site: www.cs.kuleuven.ac.be/˜koen/delay/software.shtml.
Obviously k − 1 refers to the initial conditions, when k = 0.
Here we report an outline of the used algorithm.
Algorithm 1. Input: discontinuity points τ j , t * j , τ j+1 , N , number of subgrid intervals d Output: global piecewise numerical solution, given in particular at each discretization point.

The considered interval is [τ j , τ j+1 ]
Step 0: (1) , h * * = h (2)  Step 1: collocation in intervals [t k , t k+1 ], (k = 0, ..., 2N − 1) Step 2: ∆ = (left first derivative at t * j ) -(right first derivative at t In Sect. 2 different types of discontinuities were taken under consideration; here we present some numerical examples referring to those discontinuities.We always used deficient spline functions s of degree m = 2, s ∈ C 0 in order to follow the low degree of regularity of the analytical solutions.Our algorithm was implemented in MATLAB and is available upon request from the authors About the type 1., several examples were reported in [5] and here we do not provide any further examples.
About discontinuities of type 2. we report the following two examples; they are both artificially built in order to show how the numerical method works.
Here a convenient subgrid was used in the first integration subinterval in both [0, 1] and [1,2].
Table 1 reports some used integration steps h and the corresponding absolute errors of the computed solution in t = 2, indicated by Err A .We remark that in this case we have to use small integration steps, because in the first lag we are approximating an increasing curve of exponential type, which is to be considered an unstable solution (i.e.[1, p. 25]).For h = 0.001, in t = 1, we have numerically that the computed first derivative is s (1) = 11.575 using the last left interval, whereas s (1 + h/d) = −6.2915(d = 25), using the first right interval.The concordance with analytical left and right first derivatives in t = 1 is encouraging: indeed their rounded values are respectively y (1) − = 11.575 and y (1) + = −6.2899.Fig. 1 shows the numerical solution of Eq. (A), computed with integration step h = 0.001, which from a graphical point of view, coincides with the exact solution.is to be solved on [0, 2] with history y(t) = e −t for t ≤ 0. The analytical solution is Therefore the solution is again continuous with discontinuity in the first derivative at t = 1.
As in the previous Example the algorithm uses a convenient subgrid in the first integration subinterval in both [0, 1] and [1,2].Comparing our computed solution with the exact one in t = 2, we obtained an absolute error equal to 9E − 15 when the integration step is h = 0.5; obviously is not worth reducing the integration step any more.Figure 2 shows the exact solution (solid line) together with the integration intervals (squares) and the used subgrids (circles).
We emphasize that these two examples present different stiffness features and Eq.(B) is to be considered the worse case from this point of view.On the other hand, in the first lag Eq. (A) exhibits an unstable solution and this compels small integration steps, whereas in the first lag Eq. (B) exhibits an asymptotically stable solution, so very large integration steps can be used with very satisfactory results.In Fig. 2 the excellent overlapping of the numerical solution of Eq. (B) with the analytical one is shown.From these examples we conclude that the stability properties of the solution in the first lag look relevant for the order of convergence; indeed our method does not comply with its own (second) order for Eq.(A), whereas it does for Eq.(B) as the stability properties of the two solutions are very different.
Next we consider the following example relating to discontinuities of type 3.
which is to be solved on [0, 1] with history y 1 (t) = t + 1 2 , y 2 (t) = 1 for t ≤ 0. This is again an example built artificially, so no comparisons with results in literature can be done.
The history function exhibits a discontinuity in t = − 1 2 and this propagates in the integration interval; so we have that the analytical solution is y 1 (t) = −0.5t 2 +1.5t+0.5, y 2 (t) = 1 for t ∈ [0, 1  2 ], whereas it is y . The solution belongs to class C 0 , as the first derivative is discontinuous in t = 1  2 .In this case the algorithm uses a convenient subgrid only in the first subinterval [0, 1  2 ]; actually differences between left first derivatives and right first derivatives in t = 1  2 are very small.So the subgrid in the first interval in [ 1 2 , 1] is not implemented by the algorithm.Nevertheless the solution computed with integration steps h = 0.5 and h = 0.25 (very large integration steps) exhibits an error O(10 −16 ), which is the machine precision.Here the reason of such good results is that the analytical solution is given by polynomials of degree less or equal to 2 and we approximate them by means of splines of degree m = 2.
About discontinuities of type 4., we report the following nonlinear example We consider this problem in order to get a comparison with applied models.Eq. (D) is presented by Shampine and Thompson as Example 6, but originally it was proposed as a model for life cycle of population of lemmings ( [12] and references therein).This problem is characterized by discontinuity in the solution itself at the initial point and by discontinuity in the first derivative at t = 0.74.We notice that the equation has a constant solution y(t) = 19.01 in the first lag; then the solution moves away from that constant value because of the discontinuity at the initial point and presents an increasing oscillating behavior.Shampine and Thompson use y(0) = 19.00001,whereas here we use y(0) = 19.01,so as to see the cyclic behavior sooner.
The exact (rounded) solution in t = 1.48 is y(1.48)= 18.9841040225841; the following Table reports the absolute error of our numerical solution referring to that exact value, together with the corresponding used integration steps.From Table 2 we see that results are really satisfactory.Therefore our method looks to deal well with the discontinuity in the solution itself at the initial point.However we remark that when the total integration interval becomes large, the increasing oscillating behavior of the solution is not well approximated by our piecewise polynomials of degree 2.
The last two examples are known as stiff problems.The linear DDE equation is to be solved in [0, 10] with history y(t) = e −t for t ≤ 0. The reference solution [14] is y(10) = 0.10954547858196.In t = 1 the solution is decreasing and continuous with its first derivative, but when the used integration step is not small enough, the numerical solution can be viewed as rough.Table 3 reports the absolute errors of our results (indicated by Err E ) compared with those by Zennaro [14] (indicated by Err Z ).Comparison shows that our results are very good for large integration steps, and sometime even better than those reported in literature.Figure 3 shows the exact solution (solid line) together with our numerical integration intervals (squares) and the used subgrids (circles), referring just to the first two lags.
At last for Eq.(E), Table 5 reports a comparison of computing time between our algorithm and dde23 by Shampine and Thompson [12], using MATLAB 6 and a Personal Computer equipped with a AMD Athlon Processor (750 MHz).The first column reports the order of magnitude of both the absolute errors at t = 10; T 1 is the time (in seconds) taken by our algorithm and analogously T 23 the time by dde23.The saving of computing time achieved by our algorithm is evident.The last example we present is the following equation which is to be solved in [0, 2] with history y(t) = cos t for t ≤ 0. This is a very stiff equation already studied by Guglielmi and Hairer [9].The analytical solution is where K = −4.5538739E+15.The coefficients are rounded to 8 digits.The behavior of the solution is similar to that one of the previous equation, moreover the problem looks ill-conditioned with respect to the number of significant digits considered in the coefficients.In t = 1 the solution is increasing and continuous with its first derivative; however it has a relative maximum in t = 1.00005, so for any integration step greater than 0.00005, the problem looks numerically rough.Results are good and analogous to those reported for Eq.(E); their absolute errors are reported in Table 6, referring to the analytical solution in t = 2. Figure 4 shows the exact solution (solid line) together with the numerical integration intervals (squares) and the used subgrids (circles).We remark that even in the last two examples, results are very satisfactory in particular for large integration step size, when subgrids follow the exact solution in a very nice way.
Summarizing our numerical results, we notice that if the class of low regularity of analytical solutions is preserved by the numerical method, then results can be very good even when the used numerical method looks easy and simplified, as in the case of the numerical method we propose.

Table 4
[3]orts the absolute errors of our results (indicated by Err E ) compared with those by Brugnano and Trigiante[3](indicated by Err BT )