Applied Mathematics and Computation 187 (2007) 725–740 www.elsevier.com/locate/amc
Segmented Tau approximation for test neutral functional differential equations Luis F. Cordero a, Rene´ Escalante
b,*,1
a
b
Departamento de Matema´tica y Estadı´stica, Facultad de Ciencias Econo´micas y Sociales, Universidad de Carabobo, Venezuela Departamento de Co´mputo Cientı´fico y Estadı´stica, Divisio´n de Ciencias Fı´sicas y Matema´ticas, Universidad Simo´n Bolı´var, Ap. 89000, Caracas 1080-A, Venezuela
Abstract We use the segmented formulation of the Tau method to approximate the solutions of the neutral delay differential equation y 0 ðtÞ ¼ ayðtÞ þ byðt sÞ þ cy 0 ðt sÞ þ f ðtÞ; yðtÞ ¼ WðtÞ; t 6 0;
t P 0;
which represents, for different values of a, b, c and s, a family of functional differential equations that some authors have considered as test equations in different numerical experimentations. The Tau method introduced by Lanczos is an important example of how to get approximations of functions defined by a differential equation. In the formulation of a step by step Tau version is expected that the error is minimized at the matching points of successive steps. Through the study of recent papers it seems to be demonstrated that the step by step Tau method is a natural and promising strategy for the numerical solution of functional differential equations. In preliminary experimentation significant improvements have been obtained when compared with the numerical results obtained elsewhere. 2006 Elsevier Inc. All rights reserved. Keywords: Functional differential equations; Neutral delay differential equations; Polynomial approximation; Step by step Tau method approximation
1. Introduction and preliminaries In this paper the segmented Lanczos–Tau method is used to find numerical solutions of the linear nonhomogeneous functional differential problem of neutral type
*
Corresponding author. Address: CCS: 91646, P.O. Box 025323, Miami, FL 33102 5323, USA. E-mail addresses:
[email protected] (L.F. Cordero),
[email protected] (R. Escalante). 1 This author was partially supported by the Centro de Estadı´stica y Software Matema´tico (CESMa) and the Decanato de Investigacio´n y Desarrollo (DID) at USB, project DI-CB-020-05. 0096-3003/$ - see front matter 2006 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2006.08.085
726
L.F. Cordero, R. Escalante / Applied Mathematics and Computation 187 (2007) 725–740
y 0 ðtÞ ¼ ayðtÞ þ byðt sÞ þ cy 0 ðt sÞ þ f ðtÞ; yðtÞ ¼ WðtÞ; t 6 0;
t P 0;
ð1Þ
where s > 0, W is a continuous function in [s, 0] and f is a continuous function in [0, 1). This equation has as particular cases: the test neutral (c 5 0 case) and test delay (c = 0 case) functional differential equations. For example, if in (1) we do • a = b = 1, c = 2, s = 1, W(t) = t and f(t) = 0, we obtain test equation used in [1] with an unstable difference operator (see Section 3, Experiment 2), • a = c = 1, b = 0, s = 1, W(t) = 1 and f(t) = 0, we obtain the test equation used in [2] (see Section 3, Experiment 3), • a = 1, b = 1, c = 1/4, s = 1, W(t) = t and f(t) = sin(t), we obtain the test equation considered in [1]. Neglecting the forcing term, the resulting problem is unstable (see Section 3, Experiment 5), • a = p e3pp/2, b = 1, c = 0, s = 3p/2, W(t) = ept + sin(t) and f(t) = (p e3pp/2)sin(t), we obtain the test equation considered in [3]. When p takes large negative values the stiffness of Eq. (1) is increased substantially (see Section 3, Experiment 4), which all have available analytical solutions [4] (we hope that this fact will help us to show the ‘kindness’ of the method here proposed). Although the class of equations considered here is relatively simple, the method can be extended to more general problems (e.g. nonlinear systems). See [5–7] where some extensions to nonlinear case are considered using the step by step Tau method. In [3] the particular case c = 0 was studied applying a spectral method type Tau [8] (our numerical results are consistent with the theoretical and practical results reported there, see Section 3, Experiment 4). The case c = 0 is the so-called linear stability DDE test equation, and for a particular choice of parameters we can produce a stiff delay differential equation (e.g., see Section 3, Experiment 4) or an unstable problem (for example, if a = 5, b = 1, W(t) = 1 and f(t) = 0 [9,10]). The theory of delay differential equations has undergone rapid development in the last fifty years. We suggest to see [11,12], where the readers will be able to find a wealth of reference materials on the subject. In more general models, the derivative y 0 (t) may depend on y and y 0 itself at some past value t s. These type of equations are the so-called neutral differential equations. Recently these equations have received considerable attention by many authors (see some references in [12]). For a review of some basic concepts concerning neutral and delay differential equations we refer to [13,14]. The Tau method (or the s-method), first introduced by Lanczos [15–17] and later used by Luke [18] to obtain rational and polynomial approximations mostly for hypergeometric functions, is an important example of how to get approximations of functions defined by a differential equation. A convergence analysis and error bounds for Tau method was considered by Lanczos [15,19], Luke [18], and Ortiz and Pham [20,21]. The recursive form of the Tau method, formalized by Ortiz in [22], was extended to the case of systems of ordinary differential equations in [23], and also an error analysis was given there. The basic philosophy of the method was extended to the numerical solution of linear and nonlinear initial value, boundary value, and mixed problems for ordinary differential equations (see [24,21,25]), consequently applied to the eigenvalue problems [26–28], to ‘‘stiff’’ problems [24], and to the partial differential equations [29], among others. The Tau method has also been used as an analytic tool in the discussion of equivalence results across numerical methods (see [30,31]). It is an important feature of the Tau method that no trial solutions, approximate quadratures or large matrix inversions are required [24]. In the formulation of a step by step Tau version it is allowed to construct piecewise polynomial approximations of a given function which can be used to start a refining process (see [32] for details). The Ortiz’ Step by Step Tau method (or SST method to abbreviate) was later applied efficiently to the solution of linear and nonlinear boundary value problems [24]. More recently, computational strategies for a parallel implementation of the SST method were proposed in [33]. In recent papers [5–7] the step by step Tau method was applied to find polynomial approximations to the solution of particular nonlinear delay differential equations. These papers seem to show that the segmented Tau method is a natural and promising strategy in the numerical solution of functional differential equations.
L.F. Cordero, R. Escalante / Applied Mathematics and Computation 187 (2007) 725–740
727
Next we briefly sketch the Ortiz’s recursive formulation of the Tau method and its segmented version (the SST method). The notation used here is, with slight modifications, that of [22]. Let Pj be the class of polynomials of degree less than or equal to j. Let us consider the equation defined by the differential operator L m X LyðxÞ pi ðxÞy ðiÞ ðxÞ ¼ f ðxÞ; ð2Þ i¼0
where pi(x) and f(x) are polynomials of finite degree (on the contrary, they are polynomial approximations of given functions), and y(i) represents the ith derivative of y(x). We assume also that the solution y(x) of (2) satisfies conditions ðfk ; yÞ ¼ sk ;
k ¼ 1; . . . ; m;
x 2 ½a; b;
ð3Þ
where [a, b] is a compact interval, fk are linear functionals acting on y(x), and represents the supplementary conditions (i.e., the initial, boundary or mixed conditions) to be satisfied by the required solution y(x). In order to express the Tau method approximations as a weighted arithmetic mean of successive partial sums of the power series solutions and obtaining approximations of high precision, Lanczos proposed in 1956 [17] the concept of the canonical polynomials Qn(x), n 2 N0 N [ f0g, associated with a linear differential operator. Afterwards, Ortiz introduced the more workable definition of Qn(x) [22], LQn(x) = xn + Rn, where Rn(x) is a polynomial generated by {xi}, i 2 S (set of indices for which canonical polynomials remain undefined), and is called the residual polynomial of Qn(x). Note: RS spani2S{xi}. Another important concept is that of generating polynomials [22], which are obtained from applying L to xn, where n 2 N0 . From them we can find a recursive relation for the canonical polynomials. Each differential operator L (that belongs to the class characterized by (2)) is uniquely associated with a sequence {Qn}, n 2 N0 S, of canonical polynomials. For further details about how to generate this sequence and precise definitions of Ortiz’ algebraic theory of the Tau method we refer to [22]. We now follow the notation used in [24]. Let v = {vi(x)} = Vx be a polynomial basis defined by a lower ^ ¼ fQ ^ n ðxÞg, n 2 N S such that triangular matrix V = (vij), with i; j 2PN, acting on x = (1, x, x2, . . .)t. Let Q n ^ ^ ^ LQn ðxÞ ¼ vn ðxÞ þ Rn ðxÞ, and Qn ðxÞ ¼ j¼0 vn Qj ðxÞ, where j 2 S. We consider the perturbed equation ð4Þ Ly n ðxÞ ¼ f ðxÞ þ H n ðxÞ; Pm ðnÞ where H n ðxÞ ¼ i¼0 si vni ðxÞ 2 P n ; which is called the perturbation term, is expressed in terms of the basis v. Pr ðnÞ The si ’s are parameters that we wish to find. Further we assume that f ðxÞ ¼ i¼0 ai vi ðxÞ; where r 6 n. Hence, m r X X ðnÞ ^ ^ i ðxÞ; si Q ai Q ð5Þ y n ðxÞ ¼ ni ðxÞ þ i¼0
i¼0
where i 62 S. Eq. (5) is called the Tau approximant of order n of y(x), and satisfies exactly (4) and (3). m Tau parameters are chosen such that yn satisfies exactly conditions (3). Additional Tau parameters (say s) are chosen such that the residuals of Lyn(x) match the components of f(x) + Hn(x) belonging to RS. If there exist t exact polynomial solutions of (2), then s + m t 1 = m (see [24,22] for details). We observe that the generation of the canonical polynomials depend neither on the conditions (3) nor on the interval in which the solution is required. These features allow us to introduce the concept of Ortiz’s segmented approximations (i.e. SST method). Let P be a partition of the interval [a, b] into subintervals Ij = [xj1, xj], j = 1, . . . , E, where x0 = a < x1 < x2 < < xE = b. We consider, for j = 1, . . . , E, the Tau approximants y nj ðxÞ, with x 2 Ij, which define a piecewise Tau approximant yn(x) of order n of the solution y(x) of problem (2) and (3), if each of the y nj ’s satisfies (2) with a polynomial perturbation term H jn ðxÞ, with x 2 Ij, and, for k = 1, . . . , m, ðfk ; y nr Þ ¼ sk ðiÞ (r = 1, E), and, for j = 2, . . . , E, y ðiÞ nj1 ðxj1 Þ ¼ y nj ðxj1 Þ, i = 0, . . . , m 1. The construction of a piecewise Tau approximation yn(x) of the solution y(x) of problem (2) and (3) depends only on one matrix V and one canonical sequence Q [24, Theorem 2]. There is a range of possibilities in the choice of a basis v. For instance, if v ¼ fxn gn2N , the Tau method carries out the power series expansion method, where a high accuracy is to be expected near the point of expansion. On the other hand, Lanczos in [17,19] suggested the choice of Chebyshev and Legendre polynomials to
728
L.F. Cordero, R. Escalante / Applied Mathematics and Computation 187 (2007) 725–740
obtain a better distribution of errors over the interval in which the original problem is defined and the approximate solution is required. In particular, Lanczos’ remark concerning the approximations obtained using a Legendre polynomial perturbation term, allows us to obtain accurate estimations at the end point of that interval [19], which is the key fact to construct the step by step formulation of the Tau method in which the error is minimized at the matching point of successive steps [32]. This paper is organized as follows. In Section 2, we find the piecewise polynomial approximation of the involved neutral problem using the Lanczos–Tau method, and we solve the outlined test neutral differential problem using the recursive formulation of the Tau method. In Section 3, we present preliminary numerical experimentation to illustrate the advantages of the proposed strategy. Finally, in Section 4, we present some concluding remarks. 2. Solving test neutral differential equation using the segmented formulation We will use the following notation. The symbol I will denote the unit interval [0, 1] and, for each integer k, k P 1, the interval [ks, (k + 1)s] will be denoted by Ik. Motivated by some recent papers [6,7] we consider the application of the segmented Tau method approach to problem (1). In particular, we apply the Tau method to each finite interval Ik of the domain and we solve the given problem in each interval separately. Starting with the continuous function W given on [s, 0], which will be a polynomial function (otherwise, we will consider an accurate polynomial approximation of W), the method generates a polynomial approximation on the next interval [0, s]. In the same way it proceeds into subsequent intervals, each one with the same length s. The initial condition is satisfied by the approximate solution in the first interval, and the initial condition in each subsequent interval is taken to be the endpoint value of the previous solution. Thus, Eq. (1) is defined over [0, 1) in a piecewise manner. Here, we notice that our segmented strategy will force to shift the intervals Iks to the unit interval I. It is clear that if t 2 [s, 1), then t 2 Ik for some k P 1. Hence, t = s(x + k), with x 2 I. Then, we can define y k ðxÞ ¼ yðsðx þ kÞÞ ¼ yðtÞ;
x 2 I; k P 0;
y 1 ðxÞ ¼ wðsðx 1ÞÞ; x 2 I; F k ðxÞ ¼ fk ðsðx þ kÞÞ; x 2 I; where w is equal to W if the latter is a polynomial (on the contrary, it is polynomial approximation in [s, 0]). Moreover, if f is a polynomial, then for each k P 0, fk is equal to f when t is restricted to lie in Ik (on the contrary, it is a polynomial approximation of f in Ik). We will suppose that rk r X X ð1Þ ðkÞ y 1 ðxÞ ¼ aj xj and F k ðxÞ ¼ cj xj ; k P 0; j¼0
j¼0
where r and rk are nonnegative integers. We will distinguish two cases: a 5 0 and a = 0. Case a 5 0: Eq. (1) leads to the sequence of differential equations y 0k ðxÞ asy k ðxÞ ¼ bsy k1 ðxÞ þ cy 0k1 ðxÞ þ sF k ðxÞ; r X ð1Þ y 1 ðxÞ ¼ aj xj ; x 2 I;
x 2 I; k P 0; ð6Þ
j¼0
subject to conditions y k ð0Þ ¼ y k1 ð1Þ;
k P 0:
ð7Þ
Let n be a fixed integer greater than or equal to max{r, rk} (where r = degree(y1(x)) and rk = degree(Fk(x)), for each k P 0). Let us note that if there exists a countable set of fk’s, then it is required that the set of all its degrees is bounded above. From Eq. (4), for each k P 0, we will consider the following perturbed form of Eq. (6):
L.F. Cordero, R. Escalante / Applied Mathematics and Computation 187 (2007) 725–740
Y 0k ðxÞ asY k ðxÞ ¼ bsY k1 ðxÞ þ cY 0k1 ðxÞ þ sF k ðxÞ þ H ðkÞ n ðxÞ;
x 2 I; k P 0;
729
ð8Þ
where Y1(x) = y1(x) (x 2 I), Yk1(x), for each k P 1, is the polynomial solution from the previous interval, and the perturbation term H ðkÞ n ðxÞ is defined in terms of a shifted Legendre polynomial as n X C j xj ; H ðkÞ n ðxÞ ¼ sk j¼0
where the Cjs are the coefficients of the shifted Legendre polynomial P n ðxÞ of degree n in I. The choice of Legendre over Chebyshev polynomials is motivated by their superior endpoint accuracy (see comments at the end of Section 1, and [19]). The polynomials we seek are assumed to be of the form n X ðkÞ aj xj ; k P 1; Y k ðxÞ ¼ j¼0
with the conditions Y k ð0Þ ¼ Y k1 ð1Þ;
ð9Þ
k P 0:
Next, let us define the linear operator L to be Lð:Þ ¼
d ðÞ asðÞ: dx
So, the perturbed Eq. (8) becomes LY k ðxÞ ¼ bsY k1 ðxÞ þ cY 0k1 ðxÞ þ
rk X
ðkÞ
scj xj þ sk
n X
j¼0
C j xj :
ð10Þ
j¼0
Now, since bsY k1 ðxÞ þ cY 0k1 ðxÞ ¼ bsaðk1Þ xn þ n
n1 X
ðk1Þ
bsaj
ðk1Þ
þ cðj þ 1Þajþ1
xj
j¼0
and denoting, for rk < n, 8 ðk1Þ ðk1Þ ðkÞ > þ cðj þ 1Þajþ1 þ scj ; 0 6 j 6 rk ; > < bsaj ðk1Þ bðk1Þ ¼ bsaðk1Þ þ cðj þ 1Þajþ1 ; rk þ 1 6 j 6 n 1; j j > > : ; j ¼ n; bsaðk1Þ n and for rk = n, ( ¼ bðk1Þ j
ðk1Þ
bsaj
bsanðk1Þ
ðk1Þ
ðkÞ
þ cðj þ 1Þajþ1 þ scj ; þ
scðkÞ n ;
0 6 j 6 n 1; j ¼ n;
Eq. (10) becomes LY k ðxÞ ¼
n X ðbðk1Þ þ sk C j Þxj : j
ð11Þ
j¼0
On the other hand, in view of the linearity of L, the canonical polynomials can be obtained from the generating polynomials (see Section 1) Lxm ¼ mxm1 asxm ¼ mLQm1 ðxÞ asLQm ðxÞ ¼ LðmQm1 ðxÞ asQm ðxÞÞ; from which xm ¼ mQm1 ðxÞ asQm ðxÞ:
730
L.F. Cordero, R. Escalante / Applied Mathematics and Computation 187 (2007) 725–740
Thus, we obtain a recursive definition of the canonical polynomials 1 Qm ðxÞ ¼ ðxm þ mQm1 ðxÞÞ; m P 0 as and by induction on m, the above relation reduces to m m! X 1 xi ; m P 0: Qm ðxÞ ¼ as i¼0 ðasÞmi i!
ð12Þ
Note that there are no undefined canonical polynomials. Therefore, the exact polynomial solution of the perturbed Eq. (11) is set to be n X ðbðk1Þ þ sk C j ÞQj ðxÞ Y k ðxÞ ¼ j j¼0
(note that Eq. (11) can be retrieved by applying the operator L to both sides of this expression). Now, if we substitute (12) in the above expression we have 0 1 j n bðk1Þ þ sk C j j! X j 1 X @ Y k ðxÞ ¼ x i A: ji as j¼0 i¼0 ðasÞ i! Pn Pn Pn Pj i j Moreover, since j¼0 i¼0 Aij x ¼ j¼0 i¼j Aji x (where Aij denotes the ijth element in an array), the exact polynomial solution of the perturbed equation (11) is reduced to 1 0 ðk1Þ n n b þ s C i! X X k i i 1 @ Ax j : Y k ðxÞ ¼ ij as j¼0 i¼j ðasÞ j!
ð13Þ
Now, for each k P 0, since the parameter sk was introduced in order to account for the initial condition imposed at x = 0, it is calculated using (9). Therefore, by setting x = 0 in (13) we solve for each k, the equation Yk(0) = Yk1(1), to obtain !1 ! n n X X C i i! biðk1Þ i! sk ¼ asY k1 ð1Þ þ : i i ðasÞ i¼0 ðasÞ i¼0 Finally, from (13) we have the desired polynomial approximations to the solution of the functional differential equation (1) t yðtÞ Y k ðxÞ ¼ Y k k ; t 2 I k ; k P 0; s or also, yðtÞ Y ½ t s
t s
h t i ; s
t P 0;
where [w] denotes the integer part of w. Case a = 0: Eq. (1) leads to the sequence of differential equations: y 0k ðxÞ ¼ bsy k1 ðxÞ þ cy 0k1 ðxÞ þ sF k ðxÞ; r X ð1Þ y 1 ðxÞ ¼ aj xj ; x 2 I;
x 2 I; k P 0; ð14Þ
j¼0
subject to conditions y k ð0Þ ¼ y k1 ð1Þ;
k P 0:
ð15Þ
L.F. Cordero, R. Escalante / Applied Mathematics and Computation 187 (2007) 725–740
731
We can summarize this case in the following theorem. P k1 ðk1Þ j Theorem 2.1. For all k P 0, if y k1 ¼ nj¼0 aj x , then problem (14) and (15) has an exact polynomial solution given by ! ðk1Þ ðkÞ nk1 rk X X bsanðk1Þ bsaj1 scj jþ1 ðk1Þ k1 þ caj xnk1 þ1 þ x ; x 2 I: ð16Þ y k ðxÞ ¼ y k1 ð1Þ þ xj þ j nk1 þ 1 jþ1 j¼1 j¼0 Proof. It follows by induction on k.
h
Finally, with the yk’s given by (16), we have the polynomial approximations to the solution of the functional differential Eq. (1) yðtÞ y k ðxÞ ¼ y k
t s
k ;
t 2 I k ; k P 0;
which we can also express as yðtÞ y ½st
t s
h t i s
;
t P 0:
.Remarks • In the case a = 0 we do not need to apply the Tau method directly (see Section 3, Experiment 7). • In case that W and all the fk’s are polynomials, exact solutions are obtained instead of approximate solutions. 3. Numerical experiments In this section we compare the performance of the piecewise polynomial approximation obtained in Section 2 to solve problem (1) with the numerical results reported by other authors. The numerical examples illustrate that feasibility and higher-order accuracy can be achieved. Moreover, the proposed method requires less computing time and can be used widely. As we observed at the end of Section 1, Lanczos demonstrated that when Chebyshev polynomials are replaced by Legendre polynomials as a basis a significant increase in accuracy at the endpoint of the approximation range of a Tau approximant is attained [19]. It follows from this result that more accurate information is provided to each interval by its predecessor by way of supplying it with a more accurate initial value. For this reason we preferred to use in (8) Legendre polynomials. The results in this section were obtained using MATLAB 7.0 package and were run on an Intel Pentium 4 Processor (with 16 significant digits). Here, we will take into account some computational parameters obtained as the processing time and the accuracy of the approximations. It already had been reported in [34] that the original Tau method is, in many cases, comparable to the accuracy of best uniform approximations or near optimal polynomial approximations of the same degree. In each case we compare the numerical results obtained by us with the best values attained from the direct evaluation of the analytical solution (i.e., see YEXAC in tables below). In the tables of numerical results we also show the nodes ‘t’ in which we evaluate the approximate solutions represented by YSST (i.e., the approximate solutions achieved from to apply the SST method). Likewise, the tables show the corresponding absolute errors (i.e., ERROR) and computing time CPUTIME (except Tables 1, 3 and 7). We use the symbol ‘ < ’ in some tables to indicate that a higher precision than 1016 was reached in any node by a particular numerical solution (although we cannot know exactly how precise it is). In the reported experiments we also compare our numerical results with the those obtained by other authors elsewhere. Next we show some of the carried out numerical experiments. For our experimentation we select, for different values of the parameters, some particular equations of the type (1) taken from [35,4]. Let us note that
732
L.F. Cordero, R. Escalante / Applied Mathematics and Computation 187 (2007) 725–740
Table 1 Results of Experiment 1 showing the errors involved by the application of the SST method and spline strategies reported in [1] t
0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0
YEXACT
0.2553506895400424. . . 0.5229561744103176. . . 0.8055297000976271. . . 1.1063852321231171. . . 1.4295704571147614. . . 1.7025852818153557. . . 2.0904677160858514. . . 2.6208949716308472. . . 3.3281691659926915. . . 4.2547941531425408. . .
n=3
n=7
ERROR(YSST)
ERROR(YSST)
ERROR(d16 S )
ERROR(d16 H)
1.42 · 1004 1.39 · 1004 1.78 · 1004 9.36 · 1005 7.01 · 1006 1.62 · 1003 1.51 · 1003 1.96 · 1003 1.06 · 1003 1.16 · 1004
5.60 · 1011 1.39 · 1012 6.36 · 1011 9.44 · 1012 2.22 · 1016 8.08 · 1010 4.10 · 1012 9.20 · 1010 1.16 · 1010 3.55 · 1015
1.59 · 103 3.89 · 103 5.09 · 103 1.60 · 103 4.05 · 103 2.57 · 103 5.49 · 104 2.90 · 103 2.28 · 103 2.29 · 104
3.10 · 104 7.22 · 104 9.25 · 104 1.33 · 103 4.04 · 103 1.10 · 103 2.47 · 103 2.44 · 103 3.64 · 103 4.64 · 103
Experiments 4 and 5 involve nonhomogeneous problems, and Experiments 1, 2, 3, 5 consider functional differential equations of neutral type (i.e., c 5 0). MATLAB function dde23 is a solver based on the explicit Runge–Kutta (2, 3) pair used by the ODE solver ode23 [36]. We test the case c = 0 through the Experiments 4 and 6 and we demonstrate the reliability, precision and speed of our method as compared to dde23. We always ran DDE solver dde23 using the optional integration argument options to set the scalar relative error tolerance RelTol in 1e8 and the vector of absolute error tolerances AbsTol in 1e16 (in all components). Although MATLAB does not currently possess a built-in routine to solve neutral differential equation problems directly, we made use of dde23 in our programming software to treat numerically the neutral case (i.e. c 5 0). This was made in two ways. In the first one, we use in the current step the derivative information, provided by the same dde23, in the immediately previous interval Ik (see ‘‘approxder’’ at top of the neutral case tables). In the second way, we will suppose that the analytical solution of the given problem is previously known, so that we could calculate the derivative corresponding to the current interval Ik directly, and use this information to compute the numerical solution in the following interval (see ‘‘exactder’’ at top of the neutral case tables). Since in this paper we treat with a simple equation (i.e., Eq. (1)), these computing strategies do not present bigger difficulty (nevertheless, we were monitoring any theoretical or practical difficulty that could exist [37, Section 4.5]). Finally, in Experiment 6, we consider the particular case a = 0, which was treated in the last part of Section 2 (see Theorem 2.1). Experiment 1. We consider Eq. (1) with a = b = 1, c = 1/4, s = 1, W(t) = t and f(t) = 0 (see [1]). It has in [0, 2] the analytical solution ( 1 þ t þ 14 et ; 0 6 t 6 1; yðtÞ ¼ 1 4 1 t 17 t1 3 t1 t þ 4 e þ 16 e þ 16 te ; 1 6 t 6 2; 2 with a first-order discontinuity at t = k for k P 0. Table 1 shows the numerical results obtained. The last two columns of Table 1 show the errors obtained by Kappel and Kunisch in [1]. They correspond to the best approximations attained for this experiment using 16 cubic spline (d16 S ) and cubic Hermite spline (dH ) approximations. It is interesting to observe that it is sufficient to use polynomial approximations of degree 3 to attain similar errors to those obtained in [1], which are also shown in Table 1. Unfortunately, we cannot make an appropriate comparison with the results obtained by these authors since the hardware and software platforms used by them were very different to the one used by us. Table 2 presents numerical results when the MATLAB function dde23 is used. Here we only show segmented Tau approximations of degree 7. It is clear that the results obtained were favorable to our approach.
L.F. Cordero, R. Escalante / Applied Mathematics and Computation 187 (2007) 725–740
733
Table 2 Results for Experiment 1 t
0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0
YEXACT
n=7
Approxder
Exactder
ERROR(YSST)
ERROR(dde23)
ERROR(dde23)
0.2553506895400424. . . 0.5229561744103176. . . 0.8055297000976271. . . 1.1063852321231171. . . 1.4295704571147614. . . 1.7025852818153557. . . 2.0904677160858514. . . 2.6208949716308472. . . 3.3281691659926915. . . 4.2547941531425408. . .
5.60 · 1011 1.39 · 1012 6.36 · 1011 9.44 · 1012 2.22 · 1016 8.08 · 1010 4.10 · 1012 9.20 · 1010 1.16 · 1010 3.55 · 1015
2.73 · 1010 1.19 · 1009 2.93 · 1009 5.72 · 1009 9.80 · 1009 2.23 · 1008 2.33 · 1008 2.36 · 1008 2.25 · 1008 1.90 · 1008
2.73 · 1010 1.19 · 1009 2.93 · 1009 5.72 · 1009 9.80 · 1009 2.22 · 1008 2.30 · 1008 2.27 · 1008 2.06 · 1008 1.56 · 1008
CPUTIME (s)
0.320461
1.552232
0.680979
Experiment 2. Now we consider Eq. (1) with a = b = 1, c = 2, s = 1, W(t) = t and f(t) = 0 (see [1]). It has in [0, 2] the analytical solution 2 þ t þ 2et 0 6 t 6 1; yðtÞ ¼ t t1 4 t þ 2e 2ðt þ 1Þe 1 6 t 6 2; with a first-order discontinuity at t = k for k P 0. Table 3 shows the numerical results obtained. The observations made in the previous experiment can also be applied here. However, in this case we have an unstable difference operator. For that reason in [1] it is required to use the weighting function g(s) = 1 8s, s 2 [1, 0] (i.e., for an unstable difference operator they cannot use g 1 instead of the appropriate weighting function). In our case it is also enough to use polynomial approximations of degree 3 to attain similar errors to those obtained in [1]. Table 3 shows the errors corresponding to segmented Tau approximations of degrees 7 and 16. Numerical results obtained using the DDE solver dde23 are shown in Table 4. Here we also report the step-by-step Tau approximations of degree 7. It is clear the superior accuracy obtained from our piecewise polynomial approximations. In Fig. 1, we plot the piecewise polynomial approximations obtained. Experiment 3. We consider Eq. (1) with a = c = 1, b = 0, s = 1, W(t) = 1 and f(t) = 0 (see [2]). It has in [0, 4] the analytical solution 8 t e; 0 6 t 6 1; > > > < ðt 1Þet1 þ et ; 1 6 t 6 2; yðtÞ ¼ 1 2 t2 t1 t > ðt 2tÞe þ ðt 1Þe þ e ; 2 6 t 6 3; > > : 21 3 2 t3 1 2 t2 t1 t ðt 3t 3t þ 9Þe þ ðt 2tÞe þ ðt 1Þe þ e ; 3 6 t 6 4; 6 2 Table 3 Results for Experiment 2 showing the errors obtained by the application of the SST method and spline strategies reported in [1] t
YEXACT
n=7
n = 16
ERROR(YSST) 0.25 0.50 0.75 1.00 1.25 1.50 1.75 2.00
0.8180508333754827. . . 1.7974425414002564. . . 2.9840000332253496. . . 4.4365636569180911. . . 3.9525715398288463. . . 3.2197717871754872. . . 2.1157052606417484. . . 0.4684212271070258. . .
09
2.41 · 10 3.74 · 1009 2.66 · 1009
