Applied Mathematics and Computation 190 (2007) 866–881 www.elsevier.com/locate/amc
Segmented Tau approximation for a parametric nonlinear neutral differential equation Luis F. Cordero a
a,*
, Rene´ Escalante
b,1
Depto. de Matema´tica y Estadı´stica, Facultad de Ciencias Econo´micas y Sociales, Universidad de Carabobo, Maracay, Venezuela Depto. 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
b
Abstract The segmented formulation of the Tau method is used to approximate the solutions of the parametric nonlinear neutral differential equation y 0 ðtÞ ¼ ryðtÞða þ byðt sÞ þ cy 0 ðt sÞÞ; yðtÞ ¼ WðtÞ;
t P 0;
t 6 0;
which represents, for different values of the parameters r, a, b, c and s, a family of functional differential equations with some of its members arising in areas as different as the number theory, mathematical biology, and population dynamics. For this equation no closed form of analytical solution is available. The numerical results obtained are consistent with the theoretical and practical results reported elsewhere. 2007 Elsevier Inc. All rights reserved. Keywords: Delay differential equations; Functional differential equations; Neutral differential equations; Polynomial approximations; Step by step Tau method approximation
1. Introduction In this paper the segmented Lanczos-Tau method is used to find numerical solutions of the nonlinear functional differential problem of neutral type y 0 ðtÞ ¼ ryðtÞða þ byðt sÞ þ cy 0 ðt sÞÞ; yðtÞ ¼ WðtÞ;
t P 0;
t 6 0;
*
ð1Þ
Corresponding author. 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 2007 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2007.01.081
L.F. Cordero, R. Escalante / Applied Mathematics and Computation 190 (2007) 866–881
867
where r 6¼ 0; b 6¼ 0; s > 0 and W is a continuous function in ½s; 0. Eq. (1) has as particular cases some known neutral and delay differential equations. For example, if in this equation we take • r > 0; b ¼ 1, and replace c by c ðc > 0Þ, we obtain the neutral delay logistic equation y 0 ðtÞ ¼ ryðtÞða yðt sÞ cy 0 ðt sÞÞ;
t P 0;
which arises in a food-limited population model [1,2]. • r ¼ 1; c ¼ 0, and replace b by b ðb > 0Þ, we obtain the delay logistic equation y 0 ðtÞ ¼ ða byðt sÞÞyðtÞ;
t P 0;
which arises in population dynamics [2,3]. Let us notice that if in this last equation we take b ¼ 1; s ¼ 1 and do xðtÞ ¼ a1 yðtÞ 1, we obtain the delay differential equation x0 ðtÞ ¼ axðt 1Þ½1 þ xðtÞ;
tP0
that arises in the study of population dynamics [3] and the prime number distributions [4]. Since for Eq. (1) no closed form of analytical solution is available (and in particular for none of the equations of the examples above mentioned), we obtain polynomial approximations to its solutions by applying a spectral technique, the step by step Tau method with the adopted approach in the papers [3,5,6]. These seem to show that the segmented Tau method is a natural and promising strategy in the numerical solution of neutral and delay differential equations. Lanczos’ work [7] provides the background of the Tau method, likewise papers by Ortiz [8,9] provide the background of the recursive form of the Tau method and the step by step Tau method (or SST method to abbreviate). In [5] a brief sketch of the recursive formulation of the Tau method and the SST method is given. For a review of some basic concepts concerning neutral and delay differential equations we refer to [10,11]. In the Tau method, a perturbator polynomial is introduced into the differential equation and from the perturbed equation an exact polynomial solution is obtained (Tau solution). This solution is a approximation to the solution of the original differential equation. In the segmented version of the Tau method (the SST method), the interval under consideration is divided into subintervals and the Tau method is applied separately in each subinterval. The Tau solution obtained in one interval is used as an input in the next interval. In [5,6] the differential equation in each one of the subintervals is shifted to a corresponding equation in the interval [0, 1] (we apply this strategy here). This way, a sequence of differential equations defined in [0, 1] is established with the Tau solution for each equation providing information to its successor. Lanczos in [12] established that the use of a perturbator polynomial defined in terms of the Legendre polynomials, leads to obtain accurate estimations at the end points of the interval in which the original problem is given. This 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 the successive steps [9]. In light of this fact, in this paper we will define the perturbator polynomial in terms of the Legendre polynomials. 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 parametric nonlinear 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. To conclude, in Section 4, we present some final remarks. 2. Solving the parametric neutral 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 I k . Motivated by recent work [3,5,6] we consider the application of the segmented Tau method approach to problem (1). In particular, we apply the strategy used in [5,6] according to which the differential in each Ik is shifted to a corresponding equation in I and each new problem obtained this way is solved separately, but taking into account the information from the previous equation.
868
L.F. Cordero, R. Escalante / Applied Mathematics and Computation 190 (2007) 866–881
Starting with the shift of the continuous function W given in ½s; 0 to I, which is a polynomial function (otherwise, we will consider a polynomial approximation of W), the method generates a polynomial defined in I that approximates the solution of (1) in I0. With this last polynomial a polynomial defined in I is generated that approximates the solution of (1) in I1. We proceed in similar fashion into subsequent intervals I k ; k P 2. To match the polynomial approximations obtained, there are conditions imposed at endpoints of their domain. Thus, Eq. (1) is defined over ½0; 1Þ in a piecewise manner. It is clear that if t 2 ½s; 1Þ, then t 2 I k , 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; where w is the same function W given in (1) when this is a polynomial, or a polynomial approximation in ½s; 0 of W if it is not a polynomial. Thus, (1) leads to the sequence of differential equations y 0k ðxÞ rðas þ bsy k1 ðxÞ þ cy 0k1 ðxÞÞy k ðxÞ ¼ 0; y 1 ðxÞ ¼ wðsðx 1ÞÞ; x 2 I
x 2 I; k P 0;
ð2Þ
and to match the solutions, at the endpoints of the intervals, we impose the conditions y k ð0Þ ¼ y k1 ð1Þ;
ð3Þ
k P 0:
Let n be a fixed integer greater or equal than the degree of w (n is the desired degree for the approximating polynomials that we seek). For each k P 0, let us consider the following perturbed form of Eq. (2): Y 0k ðxÞ rðas þ bsY k1 ðxÞ þ cY 0k1 ðxÞÞY k ðxÞ ¼ H ðkÞ n ðxÞ;
x 2 I;
ð4Þ
where Y k1 ðxÞ, k P 1, is the polynomial solution of the previous differential equation. Here, Y 1 ðxÞ ¼ wðsðx 1ÞÞ for all x 2 I, and the perturbation term H ðkÞ n ðxÞ is defined in terms of the shifted Legendre polynomial in I, P n ðxÞ, as ! gk X ðkÞ i ðkÞ H n ðxÞ ¼ si x P n ðxÞ i¼0
with gk ð0 6 gk 6 nÞ denoting the degree of Y k1 ðxÞ (gk will have the same meaning in the rest of the paper). We note here that one of the motives for the choice of this particular form of H ðkÞ n ðxÞ is to obtain (if possible) the polynomials Y k ðxÞ ðk P 0Þ to be of degree n. Remark 1. Eq. (4) is considered only while Y k1 ðxÞ does not have the constant value ab1 in I. When Y k1 ðxÞ ¼ ab1 , for all x 2 I, Eq. (4) is not solved anymore. This is because in this case, if we do y k1 ðxÞ ¼ ab1 in (2), then it is easy to establish that the exact polynomial solution of (2) is y k ðxÞ ¼ ab1 , for all x 2 I. If we apply the same procedure successively, it is followed that for all j P k 1, the exact polynomial solution of (2) is y j ðxÞ ¼ ab1 ; x 2 I: ðnÞ
ðnÞ
We denote by C 0 ; C 1 ; . . . ; C nðnÞ the coefficients of the shifted Legendre polynomial P n ðxÞ. Using properties of polynomials multiplication to the expression ! ! gk n X X ðkÞ i ðnÞ v ðkÞ si x Cv x H n ðxÞ ¼ i¼0
v¼0
we have H ðkÞ n ðxÞ
¼
gk i X X i¼0
v¼0
! ðnÞ sðkÞ v C iv
i
x þ
ng Xk
gk X
i¼1
v¼0
! ðnÞ sðkÞ v C gk þiv
x
gk þi
Next, for each k P 0, we define the linear differential operator DðkÞ ðÞ :¼
d ðÞ rðas þ bsY k1 ðxÞ þ cY 0k1 ðxÞÞðÞ dx
þ
gk X
ni X
i¼1
v¼ngk
! ðnÞ sðkÞ nv C iþv
xnþi :
ð5Þ
L.F. Cordero, R. Escalante / Applied Mathematics and Computation 190 (2007) 866–881
869
thus, with this definition and from (5), for each k P 0: the perturbed equation (4) becomes DðkÞ Y k ðxÞ ¼
gk i X X i¼0
! ðnÞ sðkÞ v C iv
xi þ
v¼0
ng Xk
gk X
i¼1
v¼0
! ðnÞ sðkÞ v C gk þiv
xgk þi þ
gk X
ni X
i¼1
v¼ngk
! ðnÞ sðkÞ nv C iþv
xnþi ;
x 2 I:
ð6Þ
In the following, we will suppose that the polynomials that we seek are of the form: Y k ðxÞ ¼
n X
ðkÞ
aj x j ;
k P 1
j¼0
and to be in concordance with (3), they are subject to the conditions Y k ð0Þ ¼ Y k1 ð1Þ;
ð7Þ
k P 0:
We now follow the recursive formulation presented by Ortiz [8]. We define a set of canonical polynomials ðkÞ m fQðkÞ by means of the relation DðkÞ QðkÞ m ðxÞg associated with D m ðxÞ ¼ x ; m P 0. Due to Remark 1, we assume that Y k1 ðxÞ does not have the constant value ab1 in I. By applying the operator DðkÞ to xm, we have d m ðx Þ r as þ bsY k1 ðxÞ þ cY 0k1 ðxÞ ðxm Þ dx ! gk gX k 1 X ðk1Þ j ðk1Þ j m1 r as þ bs aj x þ c ðj þ 1Þajþ1 x xm : ¼ mx
DðkÞ xm ¼
j¼0
ð8Þ
j¼0
Here we will distinguish two cases (depending of the degree of Y k1 ðxÞ): gk ¼ 0 and gk 6¼ 0: • For gk ¼ 0, from (8) and the linearity of DðkÞ , it is followed that: h i ðk1Þ m DðkÞ xm ¼ mxm1 rs a þ ba0 x h i ðkÞ ðk1Þ ¼ mDðkÞ Qm1 ðxÞ rs a þ ba0 DðkÞ QðkÞ m ðxÞ h i ðkÞ ðk1Þ ðxÞ ¼ DðkÞ mQm1 ðxÞ rs a þ ba0 QðkÞ m thus, h i ðkÞ ðk1Þ xm ¼ mQm1 ðxÞ rs a þ ba0 QðkÞ m ðxÞ ðk1Þ
and since a0
¼ Y k1 ðxÞ 6¼ ab1 , we obtain the following recursive definition of the canonical polynomials:
QðkÞ m ðxÞ ¼
1 rs½a þ
ðkÞ
ðk1Þ ba0
ðxm mQm1 ðxÞÞ;
m P 0:
By induction on m this last relation reduces to QðkÞ m ðxÞ ¼
m X
m! ðk1Þ
rs½a þ ba0
1 ðk1Þ
j¼0
j!ðrs½a þ ba0
Þmj
xj ;
m P 0:
ð9Þ
870
L.F. Cordero, R. Escalante / Applied Mathematics and Computation 190 (2007) 866–881
• For gk 6¼ 0, from (8) and the linearity of DðkÞ , it is followed that: h i ðk1Þ ðk1Þ m þ ca1 DðkÞ xm ¼ mxm1 r as þ bsa0 x
gX k 1
h i ðk1Þ ðk1Þ r bsaj þ cðj þ 1Þajþ1 xjþm rbsaðk1Þ xgk þm gk
j¼1 ðkÞ
ðk1Þ
ðk1Þ
¼ mDðkÞ Qm1 ðxÞ r½as þ bsa0 þ ca1 DðkÞ QðkÞ m ðxÞ gX 1 h i k ðk1Þ ðk1Þ ðkÞ r bsaj þ cðj þ 1Þajþ1 DðkÞ Qjþm ðxÞ rbsaðk1Þ DðkÞ QðkÞ gk þm ðxÞ gk j¼1 ðkÞ
ðk1Þ
¼ DðkÞ mQm1 ðxÞ r½as þ bsa0
gX k 1
ðk1Þ
þ ca1
QðkÞ m ðxÞ !
ðk1Þ r½bsaj
ðk1Þ ðkÞ 1Þajþ1 Qjþm ðxÞ
þ cðj þ
rbsaðk1Þ QðkÞ gk þm ðxÞ gk
j¼1
thus, ðkÞ
ðk1Þ
xm ¼ mQm1 ðxÞ r½as þ bsa0
gX k 1
ðk1Þ
r½bsaj
ðk1Þ
þ ca1 ðk1Þ
QðkÞ m ðxÞ
ðkÞ
þ cðj þ 1Þajþ1 Qjþm ðxÞ rbsaðk1Þ QðkÞ gk þm ðxÞ: gk
j¼1
From this last expression we obtain the following recursive definition of the canonical polynomials: QðkÞ gk þm ðxÞ ¼ þ
1
ðkÞ
ðk1Þ
rbsagk gX k 1
ðk1Þ
xm mQm1 ðxÞ þ r½as þ bsa0
ðk1Þ
þ ca1
QðkÞ m ðxÞ
!
ðk1Þ r½bsaj
þ cðj þ
ðk1Þ ðkÞ 1Þajþ1 Qjþm ðxÞ
;
ð10Þ
m P 0:
j¼1 ðkÞ
ðkÞ
ðkÞ
Note that the first gk canonical polynomials Q0 ðxÞ; Q1 ðxÞ; . . . ; Qgk 1 ðxÞ remain undefined. If we denote 1ðk1Þ by cðkÞ and if for all s P 0 rbsag k 8 ðkÞ ; j ¼ s 1; sc > < ðk1Þ ðk1Þ ðkÞ þ ca1 ; j ¼ s; aj;s ¼ rcðkÞ ½as þ bsa0 > : ðkÞ ðk1Þ ðk1Þ rc ½bsajs þ cðj s þ 1Þajsþ1 ; j ¼ s þ 1 ð1Þ s þ gk 1; then the recursive expression (10) becomes ðkÞ m QðkÞ gk þm ðxÞ ¼ c x þ
gkX þm1
ðkÞ
ðkÞ
aj;m Qj ðxÞ;
ð11Þ
m P 0:
j¼m1 ðkÞ For each m P 0, we define bðkÞ m;m ¼ 1 and bm;s , for s ¼ m 1; m 2; . . . ; 0, in terms of the following decreasing recursive forms: w1 X ðkÞ bðkÞ agk þsþj;wþs bðkÞ s ¼ m 1 ð1Þ 0; m;s ¼ m;wþsj ; j¼0
where w ¼ minfm s; gk þ 1g. With this, we can establish the following proposition. ðkÞ
Proposition 1. For each m P 0, if Qgk þm ðxÞ is given by (11), then ! min fm2;g gX jþ1 m k 1 X X k 1g X ðkÞ ðkÞ ðkÞ ðkÞ s ðkÞ ðkÞ Qgk þm ðxÞ ¼ c bm;s x þ aj;s bm;s Qj ðxÞ þ s¼0
j¼0
s¼0
j¼m1
m X s¼0
! ðkÞ aj;s bðkÞ m;s
ðkÞ
Qj ðxÞ; m P 0:
ð12Þ
L.F. Cordero, R. Escalante / Applied Mathematics and Computation 190 (2007) 866–881
Proof. It follows by induction on m.
871
h ðkÞ
We point out that the practical utility of (12) lies in that it defines any canonical polynomial Qj ðxÞ, j P gk , only in terms of the first gk canonical polynomials. The exact polynomial solution of the perturbed equation (6), for each k P 0, is set to be ! ! ! gk ng gk gk i ni X X X Xk X X ðkÞ ðkÞ ðnÞ ðkÞ ðkÞ ðnÞ ðkÞ ðnÞ ðkÞ Y k ðxÞ ¼ sv C iv Qi ðxÞ þ sv C gk þiv Qgk þi ðxÞ þ snv C iþv Qnþi ðxÞ; i¼0
i¼1
v¼0
v¼0
i¼1
v¼ngk
x 2 I:
ð13Þ
Note that applying the operator DðkÞ to both sides of (13), Eq. (6) is retrieved. For any k P 0, two possibilities exist (i) gk ¼ 0. (ii) gk P 1. ðkÞ
Case (i). The unknown s0 was introduced in order to account for the initial condition imposed at x ¼ 0 given in (7). Thus, we have Theorem 2. If gk ¼ 0, the exact polynomial solution of the perturbed equation (6), comes given by Y k ðxÞ ¼
ðk1Þ
rs½a þ ba0
!
ðnÞ
n n X X
ðkÞ
s0
C j j! ðk1Þ
j¼i
i¼0
i!ðrs½a þ ba0
xi ;
Þji
x 2 I;
ð14Þ
where ðkÞ s0
¼ rs½a þ
ðk1Þ ba0 Y k1 ð1Þ
!1
ðnÞ
n X
C j j!
j¼0
ðrs½a þ ba0
ðk1Þ
Þj
:
Proof. By taking gk ¼ 0 in (13) and after using (9) in the expression that is obtained, we have n X ðkÞ
Y k ðxÞ ¼ s0
i¼0
0 1 ðnÞ n i X X C i i! B C ðnÞ ðkÞ i C i Qi ðxÞ ¼ h @ h iij xj A ðk1Þ ðk1Þ rs a þ ba0 i¼0 j¼0 j! rs a þ ba 0 ðkÞ s0
Pn Pi Pn Pn and since i¼0 ð j¼0 Aij xj Þ ¼ i¼0 ð j¼i Aji Þxi (where Aij denotes the ijth element in an array), (14) follows. ðkÞ Finally, due to (7), we obtain s0 by setting x ¼ 0 in (14) and solving for this unknown the equation Y k ð0Þ ¼ Y k1 ð1Þ: h ðkÞ
Case (ii). The gk þ 1 unknowns sj were introduced in order to account for the initial condition imposed at ðkÞ ðkÞ x ¼ 0 given in (7) in addition to the gk undefined canonical polynomials. Due to this, s0 ; s1 ; . . . ; sðkÞ gk are calculated, solving the linear system with gk þ 1 equations Y k ð0Þ ¼ Y k1 ð1Þ coefficient of Q0 ðxÞ ¼ 0 coefficient of Q1 ðxÞ ¼ 0 .. . coefficient of Qgk 1 ðxÞ ¼ 0: Thus, we can establish the following theorem.
872
L.F. Cordero, R. Escalante / Applied Mathematics and Computation 190 (2007) 866–881
Theorem 3. If gk P 1, the exact polynomial solution of the perturbed equation (6), comes given by " ng !# ng gk ng v Xk ðkÞ ðnÞ Xk ðkÞ ðnÞ Xk X X ðkÞ ðnÞ ðkÞ Y k ðxÞ ¼ cðkÞ s0 bj;i C gk þj þ sðkÞ bj;i C gk þjv þ bngk þj;i C nþjv xi v j¼i
i¼0
þ
"
n X
cðkÞ
gk X
v X
sðkÞ v
v¼iþgk n
i¼ngk þ1
j¼i
v¼1
j¼1
!# ðnÞ
bðkÞ ngk þj;i C nþjv
xi ;
x 2 I;
ð15Þ
j¼iþgk n
ðkÞ
where the sj ’s are obtained according to the following: • When, 2gk 6 n 1, the linear system of equations is solved ng Xk
! ðkÞ ðnÞ bj;0 C gk þj
ðkÞ s0
þ
gk ng Xk X
j¼0
v¼1
ðkÞ ðnÞ bj;0 C gk þjv
v X
þ
j¼0
! ðkÞ ðnÞ bngk þj;0 C nþjv
sðkÞ v ¼
j¼1
Y k1 ð1Þ cðkÞ
For each i ¼ 0 ð1Þ gk 1 : " ðnÞ Ci
þ
þ
ðkÞ C gðnÞ ai;0 k
i X
j iþ1 X X
þ
ðnÞ C iv
þ
ðkÞ C ðnÞ gk v ai;0
v¼1 ng Xk
iþ1 X
j¼iþ2
s¼0
þ
þ
j iþ1 X X
þ
v iþ1 X X j¼1
! ðkÞ ai;s bðkÞ j;s
þ
sðkÞ v þ
ðnÞ
C gk þjv þ ai;s bðkÞ j;s
ðnÞ C nþjv
ðkÞ C gðnÞ a k v i;0
ng Xk
iþ1 X
j¼iþ2
s¼0
ðkÞ
s0
þ
j iþ1 X X
þ
v iþ1 X X j¼1
! ðkÞ
! ðkÞ ai;s bðkÞ j;s
s¼0
j¼1
"
gk X
# ðnÞ C gk þj
s¼0
v¼iþ1
!
! ðkÞ ai;s bðkÞ j;s
!
ðkÞ ai;s bðkÞ ngk þj;s
s¼0
#
ðnÞ C gk þjv
ðkÞ
iþ1 X
ng Xk j¼iþ2
s¼0
j¼1
ðnÞ C gk þj
s¼0
j¼1
"
! ðkÞ ai;s bðkÞ j;s
ðnÞ
C gk þjv
!
ðkÞ ai;s bðkÞ ngk þj;s
ðnÞ
C nþjv
s¼0
#
ðnÞ
C gk þjv sðkÞ ai;s bðkÞ j;s v ¼ 0
• When, 2gk > n 1, the linear system of equations is solved ng Xk
! ðkÞ ðnÞ bj;0 C gk þj
ðkÞ s0
þ
gk ng Xk X
j¼0
v¼1
ðkÞ ðnÞ bj;0 C gk þjv
v X
þ
j¼0
! ðkÞ ðnÞ bngk þj;0 C nþjv
sðkÞ v ¼
j¼1
Y k1 ð1Þ cðkÞ
For each i ¼ 0 ð1Þ n gk 2 : " ðnÞ Ci
þ
þ
ðkÞ C gðnÞ ai;0 k
i X
j iþ1 X X
þ
j¼1
" ðnÞ C iv
þ
þ
þ
iþ1 X
j¼iþ2
s¼0
j iþ1 X X
ðkÞ C ðnÞ gk v ai;0
j¼1
s¼0
þ
þ
ðkÞ ai;s bðkÞ j;s
ng Xk j¼iþ2
v iþ1 X X j¼1
!
ðnÞ C gk þjv
ðkÞ ai;s bðkÞ ngk þj;s
s¼0
#
sðkÞ v
þ
gk X
ðnÞ
ai;s bðkÞ C gk þjv þ j;s
ng Xk
iþ1 X
j¼iþ2
s¼0
iþ1 X
! ðkÞ ai;s bðkÞ j;s
# ðnÞ C gk þj
ðnÞ C nþjv
þ
j iþ1 X X j¼1
ðkÞ C gðnÞ a k v i;0
þ
v iþ1 X X j¼1
! ðkÞ
ðkÞ
s0
s¼0
!
"
v¼iþ1
! ðkÞ
ðnÞ C gk þj
s¼0
v¼1 ng Xk
! ðkÞ ai;s bðkÞ j;s
ðnÞ
#
! ðkÞ ai;s bðkÞ j;s
s¼0 ðkÞ ai;s bðkÞ ngk þj;s
s¼0
ai;s bðkÞ C gk þjv sðkÞ j;s v ¼ 0
ðnÞ
C gk þjv
! ðnÞ
C nþjv
L.F. Cordero, R. Escalante / Applied Mathematics and Computation 190 (2007) 866–881
For i ¼ n gk 1: " ðnÞ Ci
þ
ðkÞ C gðnÞ ai;0 k
þ
ngk j X X j¼1
þ þ
ng Xk
j X
j¼1
s¼0
ng Xk
j X
j¼1
s¼0
!
ðkÞ
! ðkÞ ai;s bðkÞ j;s
s¼0
#
ðnÞ C gk þj
#
ðnÞ
ai;s bðkÞ C gk þjv sðkÞ j;s v þ ! ðkÞ
gk X
ðkÞ s0
þ
i X
" ðnÞ C iv
þ
v¼1
"
ng v Xk X j¼1
ðkÞ
C gðnÞ a þ k v i;0
v¼iþ1
#
þ
ðkÞ C gðnÞ a k v i;0
v X
ng Xk
j¼1
s¼0
!
ðkÞ
873
! ðkÞ ai;s bðkÞ ngk þj;s
ðnÞ
C nþjv
s¼0 ðnÞ
ai;s bðkÞ ngk þj;s C nþjv
ðnÞ
C gk þjv sðkÞ ai;s bðkÞ j;s v ¼ 0
For each i ¼ n gk ð1Þ gk 1: " ! # " ! ngk iþg ng j j k n X X X Xk X ðnÞ ðkÞ ðkÞ ðnÞ ðkÞ ðnÞ ðkÞ ðkÞ ðkÞ ðnÞ ðnÞ ðkÞ ðnÞ C i þ C gk ai;0 þ ai;s bj;s C gk þj s0 þ C iv þ C gk v ai;0 þ ai;s bj;s C gk þjv j¼1
þ
þ
þ þ
v X
ng k þj X
j¼1
s¼0
s¼0
ðkÞ ai;s bðkÞ ngk þj;s
iþgX k nþ1
ng k þj X
j¼1
s¼0
iþgX k nþ1
ng k þj X
j¼1
s¼0
v X
iþ1 X
j¼iþgk nþ2
s¼0
!
# ðnÞ C nþjv
sðkÞ v
þ
! ðkÞ ai;s bðkÞ ngk þj;s
v¼1
"
ðnÞ ðkÞ C ni1 ai;0
ðkÞ siþgk nþ1
gk X
þ
v¼iþgk nþ2 ðnÞ
ðkÞ
a þ C nþjv þ C gðnÞ k v i;0 !
ðnÞ C ij
þ
j¼iþgk nþ1
!
ðkÞ
þ
#
ðnÞ C 2nþjgk i1
ðkÞ ai;s bðkÞ ngk þj;s
j¼1 min fi;iþg Xk nþ1g
#
min fi;vg X
j X
j¼1
s¼0
"ng j Xk X j¼1
s¼0
ng Xk
ðkÞ ai;s bðkÞ j;s
! ðkÞ ai;s bðkÞ j;s
ðnÞ
C nþji1
! ðnÞ
C gk þjv
s¼0
ðnÞ
C ij
j¼v
ðnÞ
ðkÞ ai;s bðkÞ ngk þj;s C nþjv sv ¼ 0
Proof. See Appendix. h Finally, we have the polynomial approximations to the solution of the nonlinear neutral differential equation (1) t yðtÞ Y k ðxÞ ¼ Y k k ; t 2 I k ; k P 0; ð16Þ s where for each k P 0, Y k ðxÞ is found according to the following: (i) If the degree of Y k1 ðxÞ is equal to 0 and Y k1 ðxÞ does not have the constant value ab1 , then formula (14) is used. (ii) If the degree of Y k1 ðxÞ is not equal to 0, then formula (15) is used. (iii) If Y k1 ðxÞ has the constant value ab1 , then we do not seek Y k ðxÞ, instead we use Remark 1 to conclude that yðtÞ ¼ ab1 ;
t 2 ½ðk 1Þs; 1Þ:
We note that (16) can also be expressed as t h t i yðtÞ Y ½ t ; t P 0; s s s where ½w denotes the integer part of w. 3. Numerical results In this section we compare the performance of the piecewise polynomial approximation obtained in Section 2 to solve problem (1) with the theoretical and numerical results reported by other authors. Unfortunately, as
874
L.F. Cordero, R. Escalante / Applied Mathematics and Computation 190 (2007) 866–881
far as we know, there does not exist in the literature extensive well-known numerical experimentation with which we can carry out appropriate numerical comparisons. Also we do not have the analytical solutions of (1). 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. The results of the experiments in this section were obtained using MATLAB 7.2 package and were run on an Intel Pentium D processor (with 16 significant digits). In the tables of numerical results we 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). The tables also show the computing time CPUTIME. Due to the ðkÞ ill conditioned matrices involved in the sj ’s parameters computation, we have applied to some of our experiments some form of scaling strategy. MATLAB function dde23 is a solver of ordinary and delay differential equations based on the explicit Runge–Kutta (2, 3) pair used by the ODE solver ode23 [13]. We test the case c ¼ 0 through the Experiments 2, 3 and 4, and we demonstrate 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 tolerance ‘RelTol’ in 1e8 and the vector of absolute error tolerances ‘AbsTol‘ in 1e-16 (in all components). Next we show some of the carried out numerical experiment. For our experimentation we selected some particular cases of equations of the type (1) taken from [2]. pffiffi Experiment 1: We consider Eq. (1) with r ¼ ppffiffi3 þ 201 ; a ¼ 1; b ¼ 1; c ¼ 2p3 þ 251 ; s ¼ 1 and WðtÞ ¼ t þ 2. This equation is a particular case of the neutral delay logistic equation and models a food-limited population (see Section 1). It has a first-order discontinuity at t ¼ k for k P 0 integer. Table 1 shows some of the numerical results obtained by using a segmented Tau approximation of degree 7. In Fig. 1, we plot the piecewise polynomial approximations obtained. A reference solution for t ¼ 40 reported in [2] is 0:804413836191349. Thus, from the penultimate row of Table 1 it follows that for this node, the error obtained by the application of the SST method is 4:2016 107 : Experiment 2: We consider Eq. (1) with r ¼ 1; a ¼ 75 ; b ¼ 1; c ¼ 0; s ¼ 1 and WðtÞ ¼ 0:01. This equation is a particular case of the delay logistic equation (see Section 1). It has a ðk þ 1Þst order discontinuity at t ¼ k for k P 0 integer. Table 2 shows some of the numerical results obtained by using a segmented Tau approximation of degree 8 and using the DDE solver dde23. In Fig. 2, we plot the piecewise polynomial approximations obtained. Table 1 Results for Experiment 1 t
YSST
2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40
0.592758483441024. . . 0.847357153802180. . . 1.425085420368763. . . 0.499981021082671. . . 1.326667587733823. . . 1.051372017509054. . . 0.565384239897487. . . 1.526702949318229. . . 0.792765049327441. . . 0.758029610899032. . . 1.465747828783425. . . 0.639692985784929. . . 1.043299902552723. . . 1.256381746215828. . . 0.586527510055331. . . 1.319806734683968. . . 1.012381941623081. . . 0.636247886193021. . . 1.462791192028606. . . 0.804413416033268. . .
CPUTIME (s)
1.609375
L.F. Cordero, R. Escalante / Applied Mathematics and Computation 190 (2007) 866–881
875
2 1.8 1.6 1.4
y(t)
1.2 1 0.8 0.6 0.4 0.2
0
5
10
15
20 time t
25
30
35
40
Fig. 1. Segmented Tau approximation in [0,40] corresponding to Experiment 1. Table 2 Results for Experiment 2 showing also results with dde23 t
YSST
dde23
0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 6.5 7.0 7.5 8.0 8.5 9.0 9.5 10.0
0.020037090740435. . . 0.040148500529942. . . 0.080267448421081. . . 0.159316935661258. . . 0.311664854810063. . . 0.592425328370563. . . 1.064775041578119. . . 1.722001455064376. . . 2.311988759617619. . . 2.331477414785061. . . 1.690949940365590. . . 1.032428729872882. . . 0.747922040383385. . . 0.772118890293693. . . 1.011955555500156. . . 1.408125191926395. . . 1.829572427863560. . . 2.021479460402627. . . 1.805543693653312. . . 1.367208017731018. . .
0.020037090597705. . . 0.040148499963922. . . 0.080267446727054. . . 0.159316931118312. . . 0.311664843672819. . . 0.592425302542852. . . 1.064774963019523. . . 1.722001343405133. . . 2.311988678713552. . . 2.331477427578766. . . 1.690950076249193. . . 1.032428814972972. . . 0.747922062244889. . . 0.772118888595608. . . 1.011955512384420. . . 1.408125133064433. . . 1.829572371038602. . . 2.021479441647275. . . 1.805543720932559. . . 1.367208073909903. . .
CPUTIME (s)
1.03125
1.15625
A reference solution for t ¼ 10 reported in [2] is 1:367208017754907. Thus, from the penultimate row of Table 2 it follows that for this node, the errors obtained by the application of the SST method and the DDE solver dde23 are 2:3889 1011 and 5:6155 108 , respectively. Experiment 3: We consider Eq. (1) with r ¼ 1; a ¼ 32 ; b ¼ 1; c ¼ 0; s ¼ 1 and WðtÞ ¼ 32 ðt þ 1Þ. It has a ðk þ 1Þst order discontinuity at t ¼ k for k P 0 integer. Table 3 shows some of the numerical results obtained by using a segmented Tau approximation of degree 9 and using the DDE solver dde23. In Fig. 3, we plot the piecewise polynomial approximations obtained. By doing xðtÞ ¼ 23 yðtÞ 1, we obtain the problem 3 x0 ðtÞ ¼ xðt 1Þ½1 þ xðtÞ; 2 xðtÞ ¼ t; t 6 0:
t P 0;
876
L.F. Cordero, R. Escalante / Applied Mathematics and Computation 190 (2007) 866–881 2.5
2
y(t)
1.5
1
0.5
0
0
2
4
6
8
10
time t
Fig. 2. Segmented Tau approximation in [0,10] corresponding to Experiment 2.
Table 3 Results for Experiment 3 showing also results with dde23 t
YSST
dde23
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
3.175500024919042. . . 1.130860716492482. . . 0.505044661027307. . . 1.185123684686417. . . 2.478027169600124. . . 1.735171972405921. . . 0.751459530503598. . . 1.115837612191472. . . 2.136709208686827. . . 1.891342638734044. . . 0.940251538980253. . . 1.108343745704871. . . 1.935473980922958. . . 1.922526411135078. . . 1.084566539315653. . . 1.123070114330643. . . 1.803100792429180. . . 1.909460222635392. . . 1.195585192273148. . . 1.147223062932988. . .
3.175500026681722. . . 1.130860709405948. . . 0.505044651532585. . . 1.185123592919753. . . 2.478027064149757. . . 1.735172107469037. . . 0.751459699802538. . . 1.115837511866635. . . 2.136708978170111. . . 1.891342802084340. . . 0.940251673675954. . . 1.108343691444265. . . 1.935473812084545. . . 1.922526486671191. . . 1.084566662333494. . . 1.123070095788933. . . 1.803100666138099. . . 1.909460241797122. . . 1.195585296315836. . . 1.147223075455078. . .
CPUTIME (s)
1.359375
2.28125
A reference solution reported in [2] for this last problem is xð20Þ ¼ 0:235184625529838. Hence we get yð20Þ ¼ 32 ðxð20Þ þ 1Þ ¼ 1:147223061705243. From the penultimate row of Table 3 it follows that for t ¼ 20, the errors obtained by the application of the SST method and the DDE solver dde23 are 1:2277 109 and 1:3750 108 , respectively. Experiment 4: We consider Eq. (1) with r ¼ 72 ; a ¼ 1; b ¼ 191 ; c ¼ 0; s ¼ 37 and WðtÞ ¼ 19:001. Thus we 50 obtain a particular case of the equation that models the Lemmings cycle of about 4 years [2]. It has a ðk þ 1Þst order discontinuity at t ¼ ks for k P 0 integer. Table 4 shows some of the numerical results obtained by using a segmented Tau approximation and using the DDE solver dde23. In this case, motivated by the ill conditioned matrices involved, we have considered the following strategy of approximation: we initially used poly-
L.F. Cordero, R. Escalante / Applied Mathematics and Computation 190 (2007) 866–881 3.5
3
2.5
y(t)
2
1.5
1
0.5
0
0
5
10 time t
15
20
Fig. 3. Segmented Tau approximation in [0,20] corresponding to Experiment 3.
Table 4 Results for Experiment 4 showing also results with dde23 t
YSST
dde23
4.44 8.88 13.32 17.76 22.20 26.64 31.08 35.52 40.00
19.005181701590246. . . 18.905491861931360. . . 19.199451310179281. . . 24.952770201313484. . . 0.495733815547412. . . 1.162713290109195. . . 3.202036613854786. . . 8.677124972490198. . . 24.765723535481836. . .
19.005199064580463. . . 18.904572222588627. . . 19.203903271535101. . . 24.990741853727513. . . 0.474942233776180. . . 1.118567314523729. . . 3.079277123134402. . . 8.359540990541770. . . 24.765984275794626. . .
CPUTIME (s)
2.078125
26.921875
100 90 80 70
y(t)
60 50 40 30 20 10 0
0
5
10
15
20 time t
25
30
35
40
Fig. 4. Segmented Tau approximation in [0,40] corresponding to Experiment 4.
877
878
L.F. Cordero, R. Escalante / Applied Mathematics and Computation 190 (2007) 866–881
nomials of degree 4 in ½0; 10s (this corresponds to the first 10 segmented approximations), and starting from 10s on, we use polynomials of degree 9. In Fig. 4, we plot the piecewise polynomial approximations obtained. A reference solution for t ¼ 40 reported in [2] is 24:7645486398692. Thus, from the penultimate row of Table 4 it follows that for this node, the errors obtained by the application of the SST method and the DDE solver dde23 are 1:1749 103 and 1:4356 103 respectively. Also from Table 4 we observe that although the precisions reached for both approaches are similar, the CPU times are very different (in order to reach with dde23 the same precision given by our approach, the CPU time should be approximately 1200% longer than the one reported by YSST).
4. Final remarks This paper together with some recent papers [3,5,6] demonstrates that the step by step Tau method is a natural approach, easy to apply, numerically efficient, and promising strategy in the numerical solution of functional differential equations. Moreover, the approximate solutions are consistent with the results reported elsewhere.
Appendix Proof of Theorem 3. Taking in (12) m ¼ 0, the expression obtained is used in (13) and then operating on the subindexes of some canonical polynomials, we have
Y k ðxÞ ¼
gX k 1
i X
i¼0
v¼0
þ
! ðnÞ sðkÞ v C iv
ðkÞ Qi ðxÞ
gk X
þ
!" ðnÞ sðkÞ v C gk v
c
ðkÞ
þ
gk X
i¼1
v¼0
! ðnÞ sðkÞ v C gk þiv
ðkÞ Qgk þi ðxÞ
þ
# ðkÞ ðkÞ ai;0 Qi ðxÞ
i¼0
v¼0
ng Xk
gX k 1
!
gk X
ni X
i¼1
v¼ngk
ðnÞ C iþv
sðkÞ nv
ðkÞ
Qnk þðnnk þiÞ ðxÞ
and now using (12) with m ¼ i and m ¼ n gk þ i, we have Y k ðxÞ ¼
gX k 1
i X
i¼0
v¼0
þc
ðkÞ
! ðnÞ sðkÞ v C iv
þ
þ
þc
ðkÞ
ng Xk
gk i X X
ðnÞ sðkÞ v C gk v
! ðnÞ sðkÞ v C gk þiv
s bðkÞ i;s x
þ
þ
gk X
i¼1 j¼i1
v¼0
! ðnÞ sðkÞ v C gk þiv
min fnnX k þi2;gk 1g
ni X
i¼1
j¼0
v¼ngk
i¼1 j¼nnk þi1
gk X
i¼0
v¼0
ni X
gk X
i¼1
j¼0
v¼0
! ðkÞ aj;s bðkÞ i;s
ðkÞ Qj ðxÞ
þc
ðkÞ
v¼ngk
!
jþ1 X
!
nn k þi X s¼0
! ðkÞ aj;s bðkÞ nnk þi;s
ðkÞ
Qj ðxÞ
jþ1 X
! ðkÞ aj;s bðkÞ i;s
ðkÞ
Qj ðxÞ
s¼0
s¼0
ðkÞ
Qj ðxÞ
s¼0
ðkÞ
!
! ðkÞ aj;s bðkÞ nnk þi;s
ðkÞ
ai;0 Qi ðxÞ
ðnÞ sðkÞ v C gk þiv
gk nn k þi X X i¼1
ðnÞ sðkÞ nv C iþv
ðnÞ sðkÞ nv C iþv
! ðnÞ sðkÞ v C gk v
min fX i2;gk 1g
s¼0
gk X
gX k 1
i X
gX k 1
ngk X
v¼0
s¼0
ng k 1 Xk gX
gk X
gk X v¼0
i¼1
þ
ðkÞ Qi ðxÞ
ni X v¼ngk
! ðnÞ sðkÞ nv C iþv
s bðkÞ nnk þi;s x
L.F. Cordero, R. Escalante / Applied Mathematics and Computation 190 (2007) 866–881
879
By rearranging and simplifying ! gk ng gk i X Xk X X ðkÞ ðkÞ ðnÞ ðkÞ ðkÞ ðnÞ ðkÞ sv C gk v þ c bi;s sv C gk þiv xs Y k ðxÞ ¼ c v¼0
þc
ðkÞ
i¼1
gk nn k þi X X i¼1
þ þ
þ
i¼0
v¼0
v¼0 gk jþ1 X X
gk X
þ
j¼0 i X
gk X
i¼1 j¼i1
s¼0
v¼0
ðkÞ ðnÞ sðkÞ v C gk v ai;0
ðkÞ
Qi ðxÞ !
ðkÞ ðkÞ ðnÞ aj;s bðkÞ i;s sv C gk þiv
v¼0
s¼0
ng k 1 Xk gX
xs
!
ng i2;gk 1g Xk min fX
ðnÞ sðkÞ v C iv
ðkÞ
Qj ðxÞ
!
ðkÞ ðkÞ ðnÞ aj;s bðkÞ i;s sv C gk þiv
min fnnX k þi2;gk 1g
i¼1
þ
v¼ngk
s¼0
ðkÞ
Qj ðxÞ !
jþ1 X ni X
ðkÞ ðnÞ ðkÞ aj;s bðkÞ nnk þi;s snv C iþv
s¼0 v¼ngk
j¼0
gk X
!
ðnÞ ðkÞ bðkÞ nnk þi;s snv C iþv
i X
gk X
v¼0
s¼0
gk 1 X
i¼1
þ
ni X
gX k 1
nn k þi X
ni X
s¼0
v¼ngk
i¼1 j¼nnk þi1
ðkÞ
Qj ðxÞ
!
ðkÞ ðnÞ ðkÞ aj;s bðkÞ nnk þi;s snv C iþv
ðkÞ
Qj ðxÞ
To continue, we use the results ! ng ng ng i Xk Xk X Xk Xk ng Ai;s xs ¼ As;0 þ As;i xi ; i¼1
s¼1
s¼0
gk ng k þi X X i¼1 ng Xk
Ai;s x ¼
s¼0 min fi2;g X k 1g
i¼1 ng Xk i¼1 gk X
s
gk X
i¼0
s¼1
! ðkÞ
Ai;j Qj ðxÞ
!
! ¼
j¼i1
gk X
gX k 1
i¼1
j¼ngk þi1
ðkÞ
s¼iþgk n ng Xk
i¼0
j¼iþ2
!
! ðkÞ
Aj;i Qi ðxÞ;
ðkÞ
Aj;i Qi ðxÞ;
¼
min fng k 1;gk 1g X
gk X
i¼0
j¼1
gX k 1
iþgX k nþ1
i¼ngk
j¼1
¼
As;i xi ;
j¼1
!
! Ai;j Qj ðxÞ
i¼ngk þ1 min fng k 2;gk 1g X
ðkÞ Ai;j Qj ðxÞ
j¼0
i¼1
gk X
min fiþ1;ng X kg
i¼0
min fngX k þi2;gk 1g
!
n X
As;i x þ
gX k 1
¼
s¼i i
ðkÞ Ai;j Qj ðxÞ
j¼0 gX k 1
i¼1
ng Xk
!
! Aj;i
ðkÞ Qi ðxÞ
!
gX k 1
gk X
i¼ngk
j¼iþgk nþ2
þ
ðkÞ
Aj;i Qi ðxÞ;
ðkÞ
Aj;i Qi ðxÞ;
(where Au;w denotes the uwth element in an array). Thus, Y k ðxÞ ¼ cðkÞ
gk X
ðnÞ ðkÞ sðkÞ v C gk v þ c
v¼0
þc
ðkÞ
ng Xk i¼0
þ
gk 1 X
i X
i¼0
v¼0
nk ns X X
ng Xk
gk X
s¼1
v¼0
ðkÞ ðnÞ bs;0 sðkÞ v C gk þsv
þ cðkÞ
i¼1
! ðnÞ ðkÞ bðkÞ nnk þs;i snv C sþv
i
x þc
ðkÞ
ðnÞ
sðkÞ v C iv þ
v¼0
! ðkÞ
n X i¼ngk þ1
s¼1 v¼ngk gk X
ngk ng gk X Xk X
ðkÞ
ðnÞ sðkÞ v C gk v ai;0 Qi ðxÞ
s¼i
! ðkÞ ðnÞ i bðkÞ s;i sv C gk þsv x
v¼0
gk X
ns X
s¼iþgk n v¼ngk
! ðnÞ ðkÞ bðkÞ nnk þs;i snv C sþv
xi
880
L.F. Cordero, R. Escalante / Applied Mathematics and Computation 190 (2007) 866–881
þ
ng Xk
min fng k 2;gk 1g X
j¼iþ2 s¼0
i¼0
þ þ
gX k 1
min fX iþ1;ngk g
i¼0
j¼1
s¼0
þ
nk X
þ
ðkÞ
Qi ðxÞ
v¼0
!
nj X
ðkÞ ðnÞ ðkÞ ai;s bðkÞ nnk þj;s snv C jþv
nj iþ1 X X
gX k 1
iþnX k þj k nþ1 nn X
i¼nnk
j¼1
nj X
ðkÞ
Qi ðxÞ
!
ðkÞ ðnÞ ðkÞ ai;s bðkÞ nnk þj;s snv C jþv
j¼iþgk nþ2 s¼0 v¼ngk
i¼ngk
ðkÞ
Qi ðxÞ
!
ðkÞ ðkÞ ðnÞ ai;s bðkÞ j;s sv C gk þjv
s¼0 v¼ngk
j¼1
i¼0 gX k 1
ðkÞ ðkÞ ðnÞ ai;s bðkÞ j;s sv C gk þjv
v¼0
gk j X X
nk X iþ1 X
min fng k 1;gk 1g X
!
gk iþ1 X X
ðkÞ
Qi ðxÞ !
ðkÞ ðnÞ ðkÞ ai;s bðkÞ nnk þj;s snv C jþv
ðkÞ
Qi ðxÞ:
v¼ngk
s¼0
By rearranging and simplifying " g ! # ng nk ns k X X Xk ðkÞ X ðkÞ ðnÞ ðnÞ ðnÞ ðkÞ ðkÞ ðkÞ Y k ðxÞ ¼ c sv C gk v þ bs;0 C gk þsv þ bnnk þs;0 snv C sþv v¼0
þ cðkÞ þc
ng Xk
gk X
i¼1
s¼i
v¼0
n X
ðkÞ
þ þ
gX k 1
i X
i¼0
v¼0
ðkÞ ðnÞ bðkÞ s;i sv C gk þsv þ
ðnÞ ðkÞ i bðkÞ nnk þs;i snv C sþv x
!
ns X
þ
gk X
ðnÞ ðkÞ bðkÞ nnk þs;i snv C sþv
ðkÞ ðnÞ sðkÞ v C gk v ai;0
xi
þ
min fX iþ1;ngk g
gk j X X
j¼1
s¼0 v¼0
v¼0 ng iþ1 Xk X
min fng k 2;gk 1g X
gk X
i¼0
j¼iþ2 s¼0 v¼0
min fng k 1;gk 1g X
nj nk X iþ1 X X
ðkÞ
!
ðnÞ
nk X
ðkÞ
! ðkÞ ðnÞ ðkÞ ai;s bðkÞ nnk þj;s snv C jþv
ðkÞ
Qi ðxÞ
nj iþ1 X X
ðkÞ
ðnÞ
ðkÞ ai;s bðkÞ nnk þj;s snv C jþv þ
iþnX k þj k nþ1 nn X j¼1
s¼0
!
nj X
ðkÞ
• If gk 1 6 n gk 2: " g ! ng nk ns k X X Xk ðkÞ X ðkÞ ðnÞ ðnÞ ðnÞ ðkÞ ðkÞ sv C gk v þ bs;0 C gk þsv þ bnnk þs;0 sðkÞ Y k ðxÞ ¼ c nv C sþv v¼0
þ
ng Xk
gk X
i¼1
s¼i
v¼0
gX k 1 i¼0
þ
ng Xk
ðkÞ ðnÞ bðkÞ s;i sv C gk þsv
i X
!
ðnÞ ðkÞ bðkÞ nnk þs;i snv C sþv
xi
ns X
! #
ðnÞ ðkÞ bðkÞ nnk þs;i snv C sþv
xi
s¼iþgk n v¼ngk ðnÞ
sðkÞ v C iv þ
gk X
ðkÞ
ðnÞ sðkÞ v C gk v ai;0 þ
v¼0
v¼0 gk iþ1 X X
j¼iþ2 s¼0
þ
nk ns X X s¼1 v¼ngk
gk X
n X i¼ngk þ1
þ
s¼1 v¼ngk
s¼1
ng Xk
v¼0
ðkÞ ðkÞ ðnÞ ai;s bðkÞ j;s sv C gk þjv
gk j X iþ1 X X j¼1
þ
s¼0
nj nk X iþ1 X X j¼1
s¼0 v¼ngk
ðkÞ
ðnÞ
ðkÞ ai;s bðkÞ j;s sv C gk þjv
v¼0 ðkÞ ðnÞ ðkÞ ai;s bðkÞ nnk þj;s snv C jþv
ðnÞ
ðkÞ
ðkÞ ai;s bðkÞ nnk þj;s snv C jþv Qi ðxÞ:
v¼ngk
Two possibilities exist: either gk 1 6 n gk 2 or gk 1 > n gk 2.
þ
ðkÞ
Qi ðxÞ
ðkÞ ai;s bðkÞ j;s sv C gk þjv Qi ðxÞ
j¼iþgk nþ2 s¼0 v¼ngk
i¼ngk
! ðkÞ ðkÞ ðnÞ ai;s bðkÞ j;s sv C gk þjv
j¼1 s¼0 v¼ngk
i¼0
þ
!
s¼iþgk n v¼ngk
ðnÞ sðkÞ v C iv
gk 1 X
nk ns X X s¼1 v¼ngk
gk X
i¼ngk þ1
þ
s¼1 v¼ngk
s¼1
ng Xk
! ðkÞ
Qi ðxÞ
L.F. Cordero, R. Escalante / Applied Mathematics and Computation 190 (2007) 866–881
881
• If gk 1 > n gk 2 : " g ! ng nk ns k X X Xk ðkÞ X ðkÞ ðnÞ ðnÞ Y k ðxÞ ¼ cðkÞ sðkÞ C gðnÞ þ b C bnnk þs;0 sðkÞ þ s;0 g v nv C sþv v þsv k k v¼0
þ
s¼i
i¼1
þ
n X
þ
i X
i¼0 ng Xk
v¼0
ðnÞ
þ þ
gk X
xi gk j X iþ1 X X
ðkÞ
ðnÞ sðkÞ v C gk v ai;0 þ
j¼1
ðkÞ ðkÞ ðnÞ ai;s bðkÞ j;s sv C gk þjv
v¼0
þ
ðnÞ
sðkÞ v C ngk 1v þ
gk X
ðkÞ
ðnÞ
ðkÞ ai;s bðkÞ j;s sv C gk þjv
v¼0
!
ðkÞ ðnÞ ðkÞ ai;s bðkÞ nnk þj;s snv C jþv
ðkÞ
Qi ðxÞ
s¼0 v¼ngk ðkÞ
ðnÞ sðkÞ v C gk v angk 1;0 þ
v¼0 ðkÞ
s¼0
nj nk X iþ1 X X j¼1
nj nk ng Xk X X j¼1
ðnÞ ðkÞ bðkÞ nnk þs;i snv C sþv
v¼0
v¼0
þ
ðnÞ ðkÞ i bðkÞ nnk þs;i snv C sþv x
! #
ns X
sðkÞ v C iv þ
gk iþ1 X X
j¼iþ2 s¼0
þ
!
s¼1 v¼ngk
v¼0 gk X
ng k 2 X
ng k 1 X
ðkÞ ðnÞ bðkÞ s;i sv C gk þsv þ
nk ns X X
s¼iþgk n v¼ngk
i¼ngk þ1
þ
s¼1 v¼ngk
s¼1
ngk ng gk Xk X X
! ðnÞ
ng Xk
gk j X X
j¼1
s¼0
ðkÞ
ðnÞ
ðkÞ angk 1;s bðkÞ j;s sv C gk þjv
v¼0
ðkÞ
ðkÞ angk 1;s bðkÞ nnk þj;s snv C jþv Qngk 1 ðxÞ
s¼0 v¼ngk
gk 1 X
i X
i¼ngk
v¼0
nk X
ðnÞ sðkÞ v C iv
þ
nj iþ1 X X
j¼iþgk nþ2 s¼0 v¼ngk
gk X
ðkÞ ðnÞ sðkÞ v C gk v ai;0
þ
v¼0 ðkÞ ðnÞ ðkÞ ai;s bðkÞ nnk þj;s snv C jþv
ng Xk
gk j X X
j¼1
s¼0
þ
ðkÞ
v¼0
iþnX k þj k nþ1 nn X j¼1
ðnÞ
ðkÞ ai;s bðkÞ j;s sv C gk þjv
s¼0
nj X
! ðkÞ ðnÞ ðkÞ ai;s bðkÞ nnk þj;s snv C jþv
ðkÞ
Qi ðxÞ
v¼ngk
Finally, in each one of these two last expressions for Y k ðxÞ, we equal to zero the coefficients of the gk undefined canonical polynomials. Next, after we adjust some terms, we obtain (15) and the required corresponding linear systems of equations. h References [1] Y. Kuang, A. Feldstein, Boundedness of solutions of a nonlinear nonautonomous neutral delay equation, J. Math. Anal. Appl. 156 (1991) 293–304. [2] C.A. Paul, A test set of functional differential equations, Numerical Analysis Report 243, Mathematics Department, University of Manchester, 1994. [3] H.G. Khajah, Tau method treatment of a delayed negative feedback equation, Comput. Math. Appl. 49 (2005) 1767–1772. [4] E.M. Wright, A non-linear difference-differential equation, J. Reine Angew. Math. 194 (1955) 66–87. [5] L.F. Cordero, R. Escalante, Segmented Tau approximation for test neutral functional differential equations, Appl. Math. Comput. (2006), doi:10.1016/j.amc.2006.08.085. [6] E.L. Ortiz, H.G. Khajah, On a differential-delay equation arising in number theory, Appl. Numer. Math. 21 (1996) 431–437. [7] C. Lanczos, Applied Analysis, Pretince-Hall, Englewood Cliffs, NJ, 1956. [8] E.L. Ortiz, The Tau method, SIAM J. Numer. Anal. 6 (1969) 480–492. [9] E.L. Ortiz, Step by step Tau method I, Comput. Math. Appl. 1 (1975) 381–392. [10] A. Bellen, M. Zennaro, Numerical Methods for Delay differential Equations, Oxford University Press, New York, 2003. [11] O. Diekmann, S. van Gil, S.V. Lunel, H.O. Walther, Delay Equations, Springer, New York, 1995. [12] C. Lanzos, Legendre versus Chebyshev polynomials, in: J.J.H. Miller (Ed.), Studies in Numerical Analysis, Academic Press, London, 1973, pp. 191–201. [13] L. Shampine, S. Thompson, Solving DDEs in MATLAB, Appl. Numer. Math. 37 (2001) 441–458.
