H -stability of Runge–Kutta methods with general variable stepsize for pantograph equation

Applied Mathematics and Computation 148 (2004) 881–892 www.elsevier.com/locate/amc

H-stability of Runge–Kutta methods with general variable stepsize for pantograph equation Yang Xu *, MingZhu Liu Department of Mathematics, Harbin Institute of Technology, Harbin 150001, PR China

Abstract This paper deals with H-stability of Runge–Kutta methods with variable stepsize for pantograph equations. It is shown that the Runge–Kutta methods with a regular matrix A are H-stable if and only if the modulus of stability function at infinity is less than 1. Furthermore, the same conclusion can be obtained if the Runge–Kutta methods are stiffly accurate. Ó 2003 Elsevier Inc. All rights reserved. Keywords: Delay differential equations; Infinite lag; Stability; Numerical solution; Runge–Kutta method

1. Introduction We consider the delay differential equation: y 0 ðtÞ ¼ f ðt; yðtÞ; yðaðtÞÞ; yðtÞ ¼ y0 ðtÞ;

t > 0; t 2 ½inf s P 0 faðsÞg; 0 ;

ð1:1Þ

where f ; a and y0 are given functions with aðtÞ 6 t for all t P 0. In recent years, much research has been focused on the numerical solutions of delay differential equations (see, e.g., [1–4,7–22]). These systems can be found in a variety of scientific and engineering fields such as biology, economy,

*

Corresponding author. E-mail address: [email protected] (Y. Xu).

0096-3003/$ - see front matter Ó 2003 Elsevier Inc. All rights reserved. doi:10.1016/S0096-3003(02)00947-5

882

Y. Xu, M.Z. Liu / Appl. Math. Comput. 148 (2004) 881–892

computer aided design, nonlinear dynamical systems and so on, which have a comprehensive list in [10]. It is observed that the delay differential equations play an important role in explaining many different phenomena. In particular, they turn out to be very useful when ordinary differential equations are failed. As far as we know the equations can be classified into two cases according to time lag, one is those with finite time lags, i.e., lim supt!1 ðt aðtÞÞ < 1 and the other is those with infinite time lags, i.e., lim supt!1 ðt aðtÞÞ ¼ 1. There are remarkable differences, both analytically and numerically, between these two classes. Two typical examples are: y 0 ðtÞ ¼ kyðtÞ þ lyðt sÞ; yðtÞ ¼ y0 ðtÞ;

t > 0; t 2 ½ s; 0 ;

ð1:2Þ

and y 0 ðtÞ ¼ kyðtÞ þ lyðqtÞ; yð0Þ ¼ y0 ;

t > 0;

ð1:3Þ

where k; l are complex constants, s 2 ð0; 1Þ and q 2 ð0; 1Þ. While the first class of equations is being used as the test equations for assessing the stability of the numerical methods, especially the Runge–Kutta method, much work has been done in [7–9,21,22]. Theoretical study of the second class of the equations, for instance (1.3) which is called pantograph equation, can be found in [4,10,18]. The numerical methods for the second class have been studied by [3,4], in which the grid is uniform. However, this kind of equation has unbounded time lags, it is usually difficult to investigate numerically the long time dynamical behavior of exact solution due to limited computer memory as shown in [17,19]. There are two kinds of ways to avoid the storage problem. One is to transform Eq. (1.3) into an equation with constant time lag and variable coefficients as shown in [13,16] and apply a numerical method with constant stepsize h ¼ ln q=m to it. As a matter of fact, it seems like applying the numerical method with a grid which is not uniform tn ¼ q n=m and variable stepsizes hn ¼ q n=m ðq 1=m 1Þ. Another way is applying a numerical method with variable stepsizes to Eq. (1.3) directly, which are considered in [1,17]. The nonconstant stepsize strategy considered in [1] is a special case of that introduced by Liu [17]. In this paper, we focus our attention on the H-stability of the Runge–Kutta methods with variable stepsize applied to pantograph equation. In Section 2, we recall some conclusions about asymptotical stability of analytical solution to the system. Furthermore, Runge–Kutta methods with variable stepsize are also given. Finally, Runge–Kutta methods with regular A and stiffly accurate Runge–Kutta methods are applied to pantograph equation, respectively. The same sufficient and necessary condition such that the two kinds of methods are H-stable is presented.

Y. Xu, M.Z. Liu / Appl. Math. Comput. 148 (2004) 881–892

883

2. Runge–Kutta method with variable stepsize In this paper, we consider the pantograph equation: y 0 ðtÞ ¼ kyðtÞ þ lyðqtÞ;

t > 0;

ð2:1Þ

where k; l 2 C, q 2 ð0; 1Þ with the initial condition such that yð0Þ ¼ y0 2 C. The existence and uniqueness of solutions to Eq. (2.1) has been studied by Kuang and Feldstein [14], Iserles [10] and Liu [18], etc. It is demonstrated in Liu [18] that the solutions tend to zero (algebraically) if and only if Rk < 0;

jlj < jkj;

ð2:2Þ

where Rk is the real part of k. Here we give out a definition similar to the concepts in [20]. Definition 2.1. Eq. (2.1) is called asymptotically stable if, for any q 2 ð0; 1Þ and any initial value, the solution yðtÞ of this system satisfies condition lim yðtÞ ¼ 0:

t!1

Thus the analytical stability region can be defined as follows, S ¼ fðk; lÞ 2 C2 jðk; lÞ satisfy condition ð2:2Þg: Then we formulate the definition of numerical asymptotical stability. Definition 2.2. Let q 2 ð0; 1Þ and H ¼ ft0 ; t1 ; . . . ; tn ; . . .g be an assigned mesh. The asymptotical stability region of a numerical method for Eq. (2.1) is the set SðHÞ of ðk; lÞ, such that the numerical solution yn satisfies lim yn ¼ 0:

n!1

Definition 2.3. The numerical method applied to Eq. (2.1) is called H-stable, if SðHÞ  S; for any q 2 ð0; 1Þ and H ¼ ft0 ; t1 ; . . . ; tn ; . . .g. In Sections 3 and 4, we focus our attention on H-stability of Runge–Kutta method. Now we recall the Runge–Kutta method and introduce variable stepsize schemes. For the general pantograph equation: y 0 ðtÞ ¼ f ðt; yðtÞ; yðqtÞÞ;

884

Y. Xu, M.Z. Liu / Appl. Math. Comput. 148 (2004) 881–892

the Runge–Kutta method, presented by [5], gives out the recurrence relation: ( P ðnþ1Þ h ynþ1 ¼ yn þ hnþ1 si¼1 bi f ðtn þ ci hnþ1 ; yi ; y ðtn þ ci hnþ1 ÞÞ; P ð2:3Þ ðnþ1Þ ðnþ1Þ h yi ¼ yn þ hnþ1 sj¼1 aij f ðtn þ cj hnþ1 ; yj ; y ðtn þ cj hnþ1 ÞÞ; where ðA; b; cÞ denotes Runge–Kutta method, with matrix A ¼ ðaij Þss , vectors ðnþ1Þ T T b ¼ ðb1 ; b2 ; . . . ; bs Þ , c ¼ ðc1 ; c2 ; . . . ; cs Þ . And yi , y h ðtn þ ci hnþ1 Þ can be interpreted as an approximation to yðtn þ ci hnþ1 Þ and yðqðtn þ ci hnþ1 ÞÞ, for i ¼ 1; . . . ; s, respectively. Furthermore, we suppose to have the numerical solution available till some point t0 > 0. Here, variable stepsize is introduced as follows. For simplicity, without loss of generality, we assume that t0 ¼ 1. Let Tl ¼

1 1 t0 ¼ l ; ql q

Tlþ1 ¼

1 1 t0 ¼ lþ1 ; qlþ1 q

where l ¼ 0; 1; 2; . . . If we are interested in the values of yðtÞ at points tð1Þ ; tð2Þ ; . . . ; tðm 1Þ with ð1Þ t < tð2Þ <    < tðm 1Þ , then there exist k1 ; k2 ; . . . ; km 1 such that qki tðiÞ 2 ½T0 ; T1 Þ;

i ¼ 1; 2; . . . ; m 1:

We define the grid points and variable stepsizes as follows: t0 ¼ T0 ;

t m ¼ T1 ;

ti ¼ qki tðiÞ

and tkmþi ¼ q k ti ;

hnþ1 ¼ tnþ1 tn ;

for k ¼ 1; 2; 3; . . . and i ¼ 1; 2; . . . ; m 1. It is easy to see that for any n, we have qtn ¼ tn m ;

qhn ¼ hn m

and

lim hn ¼ 1

n!1

ð2:4Þ

and the stepsizes increase geometrically. If h1 ¼ h2 ¼    ¼ hm 1 , then the above definition about grid points and variable stepsizes turns out to be that described in [1]. If ti ¼ q i=m ; i ¼ 0; 1; 2; . . . ; m 1, then the above definition is shown in [17] (Section 5) and appears to be the constant step in [13], by a change of the variable.

Y. Xu, M.Z. Liu / Appl. Math. Comput. 148 (2004) 881–892

885

3. H-stability of Runge–Kutta method with regular A This section deals with the Runge–Kutta method with nonsingular A. Using the results (2.4) of variable stepsize, we find that for pantograph equation (2.1), the Runge–Kutta method (2.3) reads (

Ps ðnþ1Þ ðn mþ1Þ ynþ1 ¼ yn þ hnþ1 i¼1 bi ðkyi þ lyi Þ; P ðnþ1Þ ðnþ1Þ ðn mþ1Þ s yi ¼ yn þ hnþ1 j¼1 aij ðkyj þ lyj Þ:

ð3:1Þ

Let ðnÞ

ðnÞ

T

Yn ¼ ðyn ; y1 ; y2 ; . . . ; ysðnÞ Þ ; then the difference equation (3.1) can be written as ðnÞ

ðnÞ

ðnÞ

D0 Ynþ1 ¼ D1 Yn þ D2 Yn mþ1 ; where knþ1 ¼ khnþ1 , lnþ1 ¼ lhnþ1 and ðnÞ D0

 ¼

1 0

 knþ1 bT ; I knþ1 A

ðnÞ D1

 ¼

 0 ; 0

1 e

ðnÞ D2

with e ¼ ð1; 1; . . . ; 1ÞT 2 Rs . After some simple calculations, we can obtain that ðnÞ

ðnÞ

Ynþ1 ¼ G1 Yn þ G2 Yn mþ1 ; where 1

ðnÞ G1

¼

1 þ knþ1 bT ðI knþ1 AÞ e 0 1

ðI knþ1 AÞ e 1

ðnÞ G2

¼

0 lnþ1 bT ðI knþ1 AÞ 0

lnþ1 ðI knþ1 AÞ 1 A

0

!

! ;

:

Let T

T T ; . . . ; Yn mþ1 Þ ; Zn ¼ ðYnT ; Yn 1

then difference equation (3.1) can be rewritten as Znþ1 ¼ Qn Zn ;

 ¼

0 0

 lnþ1 bT ; lnþ1 A

886

Y. Xu, M.Z. Liu / Appl. Math. Comput. 148 (2004) 881–892

where

0

ðnÞ

0 0

G1 B I B Qn ¼ B B ... @ 0 0

0 0 0 0

... ... ... ... ...

0 0

I 0 0 I

ðnÞ

G2 0 ... 0 0

1 C C C C A

; mðsþ1Þ

with 0 and I denote zero matrix and identity matrix with appropriate dimension, respectively. Let ðnÞ

G1 ¼ lim G1 ; n!1

ðnÞ

G2 ¼ lim G2 ;

Q ¼ lim Qn ;

n!1

ð3:2Þ

n!1

then  G1 ¼ 0

T

1

0



1 b A e 0 ; 0 0

G1 B B I B Q¼B B... B @ 0 0

B G2 ¼ @

0

...

0

0

0

... ...

0

0

0 0

... ...

I 0

G2

1

0 0

1 l bT A 1 C k A; l Is k

C 0 C C ...C C: C 0 0 A I 0

It is easily seen that det½fImðsþ1Þ Q ¼ det½fm Isþ1 fm 1 G1 G2 2 ðf ð1 bT A 1 eÞÞfm 1 6 ¼ det 4 0

l T 1 3 b A k 7 :

l 5 m f þ Is k

ð3:3Þ

Theorem 3.1. Assume the matrix A is regular, then Runge–Kutta method ðA; b; cÞ, applied to any Eq. (2.1) with condition (2.2) for any initial value, is Hstable if and only if jr1 j < 1; ð3:4Þ where r1 ¼ 1 bT A 1 e is the limit of stability function of Runge–Kutta method. Proof. Assume that condition (3.4) is satisfied, then we can see from (3.3) with (2.2), qðQÞ < 1, where qðÞ denotes the spectral radius. According to Lemma 5.6.10 in [6], there exists a norm k  k such that kQk < 1. In view of (3.2), for any fixed e with 1 kQk > e > 0, there exists N1 > 0 such that for n P N1 ,

Y. Xu, M.Z. Liu / Appl. Math. Comput. 148 (2004) 881–892

887

kQn Qk < 1 kQk e and kZnþ1 k 6 ðkQk þ kQn QkÞ  kZn k < ð1 eÞkZn k; which implies limn!1 yn ¼ 0. Conversely, assume that jr1 j P 1, then similar to the analysis of Theorem 3.2 in [1], we choose l ¼ 0 and Rk < 0, which satisfies condition (2.2) and Eq. (3.1) reduces to ynþ1 ¼ rðknþ1 Þyn ; 1

where rðknþ1 Þ ¼ 1 þ knþ1 bT ðI knþ1 AÞ e and limn!1 rðknþ1 Þ ¼ r1 . If jr1 j > 1, then for any fixed e with jr1 j 1 > e > 0, there exists N2 such that for n P N2 , jynþ1 j ¼ jðr1 þ rðknþ1 Þ r1 Þyn j P ðjr1 j jrðknþ1 Þ r1 jÞ  jyn j P ðjr1 j eÞ  jyn j; which implies yn does not vanish for jr1 j e > 1. If jr1 j ¼ 1, it is obviously seen that yn is bounded, but does not vanish. The proof is completed.  Corollary 3.1. The Gauss–Legendre methods are not H-stable. Corollary 3.2. The RadauIA methods, RadauIIA methods and LobattoIIIC methods are H-stable. Corollary 3.3. The one-leg h-method is H-stable if and only if 1=2 < h 6 1:

4. H-stability of stiffly accurate method This part is concerned with the Runge–Kutta method with singular A. Here we just consider the stiffly accurate Runge–Kutta method, i.e., asj ¼ bj ;

1 6 j 6 s:

ð4:1Þ

Since A is singular, without loss of generality, we assume that a1j ¼ 0, 1 6 j 6 s and write   0 0 A¼ ;  a A where  a ¼ ða21 ; a31 ; . . . ; as1 ÞT , A ¼ ðaij Þsi;j¼2 . In fact, LobattoIIIA methods and linear h-method are typical examples of such methods.

888

Y. Xu, M.Z. Liu / Appl. Math. Comput. 148 (2004) 881–892

Now applying this special method to Eq. (2.1), we can also arrive at difference equation (3.1). Furthermore, in [13], some properties of stiffly accurate method are given. Lemma 4.1 [13]. If the Runge–Kutta method is stiffly accurate and A is regular, then 1

lim rðzÞ ¼ eTs 1 A  a;

ð4:2Þ

z!1

1

T

with es 1 ¼ ð0; 0; . . . ; 0; 1Þ 2 Rs 1 and rðzÞ ¼ 1 þ zbT ðI zAÞ e be the stability function of Runge–Kutta method. Using (4.1), we find that, for recurrence relation (3.1), ynþ1 ¼ ysðnþ1Þ . Hence, the difference equation (3.1) can be written as (

Ps ðnþ1Þ ðn mþ1Þ ¼ yn þ hnþ1 j¼1 aij ðkyj þ lyj Þ; P ðnþ1Þ ðn mþ1Þ ¼ yn þ hnþ1 sj¼1 asj ðkyj þ lyj Þ;

ðnþ1Þ

yi

ynþ1

ð4:3Þ

for i ¼ 1; . . . ; s 1. Let ðnÞ

ðnÞ

ðnÞ

T

Un ¼ ðy1 ; y2 ; . . . ; ys 1 ; yn Þ ; then Eq. (4.3) reads ðnÞ

ðnÞ

Unþ1 ¼ M1 Un þ M2 Un mþ1 ; where ðnÞ M1 ðnÞ

M2



 1 ¼ ; 1 1 knþ1 ðI knþ1 AÞ  a þ ðI knþ1 AÞ e   0 0 ¼ : 1  lnþ1 ðI knþ1 AÞ 1 A lnþ1 ðI knþ1 AÞ a 0 0

Putting T

T T ; . . . ; Un mþ1 Þ ; Wn ¼ ðUnT ; Un 1

then the difference equation (4.3) can be written in the form Wnþ1 ¼ Sn Wn ;

ð4:4Þ

Y. Xu, M.Z. Liu / Appl. Math. Comput. 148 (2004) 881–892

889

where 0

ðnÞ

0 0

M1 B I B Sn ¼ B B ... @ 0 0

0 0

ðnÞ

0 0 0 0

... ... ... ... ...

M2 0 ... 0 0

I 0 0 I

1 C C C: C A

Let ðnÞ

ðnÞ

M1 ¼ lim M1 ; n!1

M2 ¼ lim M2 ;

S ¼ lim Sn ;

 ;

0

n!1

ð4:5Þ

n!1

where  M1 ¼ 0

0

1

0

A  a

M1 B B I B S¼B B... B @ 0 0

1

M2 ¼

0

...

0

0

0

... ...

0

0

0 0

... ...

I 0

0 I

M2

0 l I k

l 1 A a 1k

! ;

C 0 C C ...C C: C 0 A 0

Theorem 4.1. Assume the matrix A is regular, then the stiffly accurate Runge– Kutta method ðA; b; cÞ, applied to any Eq. (2.1) with condition (2.2) for any initial value, is H-stable if and only if jr1 j < 1;

ð4:6Þ 1

where r1 ¼ eTs 1 A  a. Proof. Assume that condition (4.6) is satisfied, detðfI SÞ ¼ det½fm I fm 1 M1 M2 ;

ð4:7Þ 1

T we denote a ¼ ð 1; r1 ; . . . ; rs 2 Þ with ri ¼ eTi A a and T

ei ¼ ð0; . . . ; 0 ; 1; 0; . . . ; 0Þ 2 Rs 1 ; |fflfflfflffl{zfflfflfflffl} i 1

then by condition (4.6) and Lemma 4.1, for jfj P 1, we have   a fIs 1 det½fIs M1 ¼ det 6¼ 0; 0 f þ rs 1 which means that the matrix ðfIs M1 Þ is invertible whenever jfj P 1.

890

Y. Xu, M.Z. Liu / Appl. Math. Comput. 148 (2004) 881–892

Moreover, ðfI M1 Þ 1 M2

0

rs 1 B B r1 f

l B 1 B ... ¼ k ½fðf þ rs 1 Þ B B @ rs 2 f rs 1 f

0 f þ rs 1

... ...

0 0

... 0 0

... ...

f þ rs 1 0

1 1 C r1 C C ... C C; C rs 2 A f

which means that the characteristic polynomial is  s 1 l1 1 det½zI ðfI M1 Þ M2 ¼ z z þ : kf Under condition (2.2), we can arrive at     l 1 1  < 1;  sup q½ðfI M1 Þ M2 ¼ max 0;  k f whenever jfj ¼ 1 and qðÞ denotes the spectral radius. According to Corollary 1.2 in [9] (see also Lemma 2.1 in [12]), the above two results imply that the zeros f of characteristic polynomial (4.7) satisfy jfj < 1. By Lemma 5.6.10 in [6], there exists a norm k  k such that kSk < 1. In view of (4.5), for any fixed e with 1 kSk > e > 0, there exists N 1 > 0 such that for n P N1 kSn Sk < 1 kSk e and kWnþ1 k 6 ðkSk þ kSn SkÞ  kWn k < ð1 eÞkWn k; which implies limn!1 yn ¼ 0. Conversely, assume that jr1 j P 1, then we choose l ¼ 0 and Rk < 0, which satisfies condition (2.2) and Eq. (4.3) reduces to ynþ1 ¼ rðknþ1 Þyn ; where limn!1 rðknþ1 Þ ¼ r1 . If jr1 j > 1, then for any fixed e with jr1 j 1 > e > 0, there exists N 2 such that for n P N 2 , jynþ1 j ¼ jðr1 þ rðknþ1 Þ r1 Þyn j P ðjr1 j jrðknþ1 Þ r1 jÞ  jyn j P ðjr1 j eÞ  jyn j; which implies yn does not vanish. If jr1 j ¼ 1, it is obviously seen yn is bounded, but does not vanish. The proof is completed. 

Y. Xu, M.Z. Liu / Appl. Math. Comput. 148 (2004) 881–892

891

Corollary 4.1. The LobattoIIIA methods are H-stable. Corollary 4.2. The linear h-method is H-stable if and only if 1=2 < h 6 1.

Acknowledgements This paper is supported by the National Natural Science Foundation of China (10271036).

References [1] A. Bellen, N. Guglielmi, L. Torelli, Asymptotic stability properties of h-methods for the pantograph equation, Appl. Numer. Math. 24 (1997) 279–293. [2] A. Bellen, N. Guglielmi, M. Zennaro, On the contractivity and asymptotic stability of systems of delay differential equations of neutral type, BIT 39 (1999) 1–24. [3] M.D. Buhmann, A. Iserles, On the dynamics of a discretized neutral equation, IMA J. Numer. Anal. 12 (1992) 339–364. [4] M.D. Buhmann, A. Iserles, Stability of the discretized pantograph differential equation, Math. Comput. 60 (1993) 575–589. [5] K. Dekker, J.G. Verwer, Stability of Runge–Kutta Methods for Stiff Nonlinear Differential Equations, North-Holland, Amsterdam, 1984. [6] R.A. Horn, C.R. Johnson, Matrix Analysis, Cambridge University Press, Cambridge, 1985. [7] K.J. inÕt Hout, A new interpolation procedure for adapting Runge–Kutta methods to delay differential equations, BIT 32 (1992) 634–649. [8] K.J. inÕt Hout, Stability analysis of Runge–Kutta methods for systems of delay differential equations, IMA J. Numer. Anal. 17 (1997) 17–27. [9] K.J. inÕt Hout, M.N. Spijker, Stability analysis of numerical methods for delay differential equations, Numer. Math. 59 (1991) 807–814. [10] A. Iserles, On the generalized pantograph functional-differential equation, European J. Appl. Math. 4 (1993) 1–38. [11] A. Iserles, Numerical analysis of delay differential equations with variable delays, Ann. Numer. Math. 1 (1994) 133–152. [12] T. Koto, NP-stability of Runge–Kutta methods based on classical quadrature, BIT 37 (1997) 870–884. [13] T. Koto, Stability of Runge–Kutta methods for the generalized pantograph equation, Numer. Math. 84 (1999) 233–247. [14] Y. Kuang, A. Feldstein, Monotonic and oscillatory solutions of a linear neutral delay equation with infinite lag, SIAM J. Math. Anal. 21 (1990) 1633–1641. [15] M.Z. Liu, M.N. Spijker, The stability of h-method in the numerical solution of delay differential equations, IMA J. Numer. Anal. 10 (1990) 31–48. [16] Y.K. Liu, Stability analysis of h-methods for neutral functional-differential equations, Numer. Math. 70 (1995) 473–485. [17] Y.K. Liu, On the h-method for delay differential equations with infinite lag, J. Comput. Appl. Math. 71 (1996) 177–190. [18] Y.K. Liu, Asymptotic behaviour of functional-differential equations with proportional time delays, European J. Appl. Math. 7 (1996) 11–30.

892

Y. Xu, M.Z. Liu / Appl. Math. Comput. 148 (2004) 881–892

[19] Y.K. Liu, Numerical investigation of the pantograph equation, Appl. Numer. Math. 24 (1997) 309–317. [20] L.H. Lu, Numerical stability of the h-methods for systems of differential equations with several delay terms, J. Comput. Appl. Math. 34 (1991) 291–304. [21] G. Van den Heuvel, New stability results for Runge–Kutta methods adapted to delay differential equations, Appl. Numer. Math. 34 (2000) 293–308. [22] M. Zennaro, P -stability properties of Runge–Kutta methods for delay differential equations, Numer. Math. 49 (1986) 305–318.