Numerical methods for impulsive differential equation

Mathematical and Computer Modelling 48 (2008) 46–55 www.elsevier.com/locate/mcm

Numerical methods for impulsive differential equationI X.J. Ran ∗ , M.Z. Liu, Q.Y. Zhu Harbin Institute of Technology, Shenzhen Graduate School, Shenzhen, 518055, PR China Received 4 October 2006; received in revised form 12 September 2007; accepted 19 September 2007

Abstract In this paper, the asymptotical stability of the numerical methods with the constant stepsize for impulsive differential equation x(t) ˙ = αx, 4x = βx,

t 6= k, t > 0, t = k,

x(0 + 0) = x0 , where α 6= 0, β, x0 ∈ R, 1 + β 6= 0, k ∈ N, are investigated. The asymptotical stability conditions of the analytic solution of this equation and the numerical solutions are obtained. Finally, some experiments are given. c 2007 Elsevier Ltd. All rights reserved.

Keywords: Impulsive differential equation; Asymptotical stability; θ-method; Runge–Kutta method

1. Introduction Impulsive differential equations occur in many applications: physics, population dynamics, ecology, biological systems, biotechnology, industrial robotics, pharmacokinetics, optimal control, etc. The theory of the impulsive differential equations has been studied extensively [1,2,4,10,12,14–17]. In 1989, Kulev and Bainov investigated the stability and global stability of systems with impulse by Lyapunov function [5,6]. In 2000, Randelovic gave the algorithm for solving impulsive differential equations [7]. However, in these works the authors did not investigate the stability of the numerical methods for impulsive differential equations. We shall introduce some basic knowledge of the system with impulsive effect at fixed instants of time, considering the following system x˙ = f (t, x),

t 6= τk ,

4x = Ik (x),

t = τk ,

x(t0+ )

(1.1)

= x0 ,

where f : R+ × Ω → Rn , Ik : Ω → Rn , Ω ⊂ Rn is an open subset and τk < τk+1 with limk→∞ τk = 0, as usual 4x(t) = x(t + 0) − x(t − 0), x(t + 0) and x(t − 0) denote the right and left limits of x at t respectively. I This work is supported by the NSF of PR China (No. 10671047). ∗ Corresponding author.

E-mail address: [email protected] (X.J. Ran). c 2007 Elsevier Ltd. All rights reserved. 0895-7177/$ - see front matter doi:10.1016/j.mcm.2007.09.010

47

X.J. Ran et al. / Mathematical and Computer Modelling 48 (2008) 46–55

Definition 1.1 ([1]). The function x(t), t ∈ (t0 , b) is said to be the solution of the system with impulsive effect (1.1), if the following conditions are satisfied: (i) x(t0+ ) = x0 , (t, x(t)) ∈ R+ × Ω for t ∈ (t0 , b), (ii) for t ∈ (t0 , b), t 6= τk , k ∈ N, the function x(t) is differentiable and dx(t) dt = f (t, x(t)), (iii) the function x(t) is left continuous in (t0 , b), if t ∈ (t0 , b) and t = τk , t 6= b, then x(t + 0) = x(t) + Ik (x(t)). Theorem 1.2 ([20]). If f : R+ × Ω → Rn is continuous in (τk , τk+1 ] × Ω , k ∈ N, lim(t,y)→(τ + ,x) f (t, y) is finite k and exists and f is locally Lipschitz continuous with respect to x in R+ × Ω , then the solution x(t) of problem (1.1) is unique. In this paper, we consider the following impulsive differential equation (1.2), construct the numerical schemes for it and discuss the asymptotical stability of these methods. x(t) ˙ = αx, t 6= k, t > 0, 4x = βx, t = k,

(1.2)

x(0 + 0) = x0 , where α 6= 0, β, x0 ∈ R, 1 + β 6= 0, k ∈ N. Definition 1.3 ([3]). x(t) is said to be the solution of (1.2), if (i) limt→0+ x(t) = x0 = x(0 + 0), (ii) for t ∈ (0, +∞), t 6= k, k ∈ N, x(t) is differentiable and x(t) ˙ = αx(t), (iii) x(t) is left continuous in (0, +∞) and if t = k, then x(k + 0) − x(k) = βx(k), where x(k + 0) = limt→k + x(t). Problem (1.2), in (0, ∞), has a unique solution x(t) = x0 eα{t} ((1 + β)eα )[t] ,

(1.3)

where [t] and {t} denote the greatest-integer function of t and the fractional part of t, respectively. From (1.3), it is easy to obtain the following theorem Theorem 1.4. The solution x ≡ 0 of Eq. (1.2) is asymptotically stable (x(t) → 0 as t → ∞) if and only if |(1 + β)eα | < 1. In this paper, we denote the asymptotical stability region of Eq. (1.2) by H = {(α, β) : |(1 + β)eα | < 1}. We also denote the approximations to the solution x(tkm,l ) and x(k + 0) by xkm,l (1 ≤ l ≤ m) and xkm,0 , respectively. While tkm,l = k + lh, 0 ≤ l ≤ m, k, m ∈ N, h = 1/m. 2. The stability of the θ-method In this section, we investigate the stability of the θ-method for problem (1.2). Applying the θ-method to (1.2), we have xkm,l = xkm,l−1 + h(α(1 − θ )xkm,l−1 + αθ xkm,l ), xkm,0 = (1 + β)x(k−1)m,m ,

l = 1, 2, . . . , m (2.1)

x0,0 = x0 , which can be rewritten as   αh xkm,l = 1 + xkm,l−1 , 1 − θ αh xkm,0 = (1 + β)x(k−1)m,m , x0,0 = x0 .

l = 1, 2, . . . , m (2.2)

48

X.J. Ran et al. / Mathematical and Computer Modelling 48 (2008) 46–55

In fact, if t ∈ (τk , τk+1 ], Eq. (1.2) can be seen as ordinary differential equation. Hence, the θ-method for Eq. (1.2) is convergent of order 1 for θ 6= 12 and order 2 for θ = 21 . In other words, the θ-methods for (1.2) conserve the order of convergence. Definition 2.1. The θ-method for (1.2) (i.e. process (2.1)) is called asymptotically stable at (α, β), if ∃M, for any h = m1 , m > M, such that limk→+∞ X k = 0, where X k = (xkm,0 , xkm,1 , . . . , xkm,m )T . The set Sθ of all pairs (α, β) at which scheme (2.1) is asymptotically stable is called to be the asymptotical stability region. From (2.2) we have m  αh x(k−1)m,l , l = 0, 1, . . . , m. xkm,l = (1 + β) 1 + 1 − θ αh In the following we define M = |α|. Theorem 2.2. The process (2.1) is asymptotically stable at (α, β) if and only if  m (1 + β) 1 + αh < 1, for all m > M, 1 − θ αh where h =

1 m,

[·] is the greatest-integer function.

Lemma 2.3 ([18]). Let ϕ(x) = 1 2 , ϕ(−∞) = 1.

1 x

−

1 ex −1 ,

x ∈ R, then ϕ(x) is a decreasing function and ϕ(+∞) = 0, ϕ(0) =

Theorem 2.4. H ⊆ Sθ , for 0 ≤ θ ≤ 1 if and only if (i) 0 ≤ θ ≤ ϕ(1) for α > 0, (ii) 0 ≤ θ ≤ ϕ(0) for α < 0. Proof. Suppose (α, β) ∈ H. H ⊆ Sθ is equivalent to  m 1 + αh ≤ eα , 1 − θ αh

(2.3)

i.e. 1+ since 1 +

αh ≤ eαh , 1 − θ αh αh 1−θαh

> 0, therefore

1 1 − αh = ϕ(αh). (2.4) αh e −1 In view of Lemma 2.3, (2.3) is equivalent to 0 ≤ θ ≤ ϕ(1) for α > 0 and 0 ≤ θ ≤ ϕ(0) for α < 0, respectively. The proof is complete.  0≤θ ≤

3. Runge–Kutta method We shall discuss s-stage Runge–Kutta method c

A bT

where   c1 c2    c =  . ,  ..  cs

 a11 a21  A= .  ..

a12 a22 .. .

as1

as2

 · · · a1s · · · a2s   ..  , .. . .  · · · ass

  b1 b2    b =  . .  ..  bs

49

X.J. Ran et al. / Mathematical and Computer Modelling 48 (2008) 46–55

Using the Runge–Kutta method to (1.2), we obtain xkm,l = xkm,l−1 + αh

s X

bi Yi ,

l = 1, 2, . . . , m

i=1

Yi = xkm,l−1 + αh

s X

ai j Y j ,

i = 1, 2, . . . , s

(3.1)

j=1

xkm,0 = (1 + β)x(k−1)m,m , x0,0 = x0 . Let Y = (Y1 , Y2 , . . . , Ys )T , h¯ = αh, e = (1, 1, . . . , 1)T . It follows that l = 1, 2, . . . , m,

xkm,l = xkm,l−1 + h¯ bT Y, Y = xkm,l−1 e + h¯ AY,

(3.2)

xkm,0 = (1 + β)x(k−1)m,m , x0,0 = x0 . For small h, matrix I − αh A is invertible, then (3.2) has a simple form xkm,l = (1 + h¯ bT (I − h¯ A)−1 e)xkm,l−1 ,

l = 1, 2, . . . , m,

xkm,0 = (1 + β)x(k−1)m,m , x0,0 = x0 .

(3.3)

If the Runge–Kutta method is of order p with the stability function R(h¯ ) = 1 + h¯ bT (I − h¯ A)−1 e, then there exists constant C > 0, such that |eh¯ − R(h¯ )| ≤ Ch p+1 . For scheme (3.1), on the assumption that xkm,l = x(tkm,l ), 1 ≤ l ≤ m, we have if l = m, |x(t(k+1)m,1 ) − x(k+1)m,1 | = |eαh x(t(k+1)m,0 ) − (1 + h¯ bT (I − h¯ A)−1 e)x(k+1)m,0 | = eαh |1 + β||eαh x(tkm,m ) − ((1 + h¯ bT (I − h¯ A)−1 e)xkm,m )| = eαh |1 + β||eh¯ − R(h¯ )||xkm,m | < e|1 + β||xkm,m |Ch p+1 = C0 h p+1 ,

(3.4)

if 1 ≤ l ≤ m − 1, |x(tkm,l+1 ) − xkm,l+1 | = |eαh x(tkm,l ) − ((1 + h¯ bT (I − h¯ A)−1 e)xkm,l )| = |eh¯ − R(h¯ )||xkm,l | ≤ |xkm,l |Ch p+1 = C1 h p+1 ,

(3.5)

where C0 = e|1 + β||xkm,m |C, C1 = |xkm,l |C. Both (3.4) and (3.5) imply that for Eq. (1.2) the Runge–Kutta method is also convergent of p order. Lemma 3.1 ([8,9,11,19]). There exists a unique (m, n)-pad´e approximation Rmn (z) = N × N. Furthermore Z (−1)n m+n+1 1 n ez Pn (z) − Q m (z) = z u (1 − u)m euz du, (m + n)! 0

Q m (z) Pn (z)

to ez , for (m, n) ∈

50

X.J. Ran et al. / Mathematical and Computer Modelling 48 (2008) 46–55

where Q m (z) = Pn (z) =

m X m! (m + n − j)! j z , (m + n)! j=0 j!(m − j)! n X (m + n − j)! n! (−z) j . (m + n)! j=0 j!(n − j)!

Lemma 3.2. Denote g j =

(m+n− j)! j!(n− j)!

( j ≤ n), m, n, j ∈ N, then {g j } is monotonically decreasing.

Proof. It is clear that g j > 0,  gj (m + n − j)! (m + n − ( j + 1))! = g j+1 j!(n − j)! ( j + 1)!(n − ( j + 1))! (m + n − j)( j + 1) = n− j   m = 1+ ( j + 1) > 1.  n− j Lemma 3.3. If |z| ≤ 1, z ∈ R, then Pn (z) > 0 and Rmn (z) > 0. Proof. By Lemma 3.1, we know that  n!   [(g0 − g1 z) + z 2 (g2 − g3 z) + · · · + z n−1 (gn−1 − gn z)] n is odd (m + n)! Pn (z) = n!   [(g0 − g1 z) + z 2 (g2 − g3 z) + · · · + z n−2 (gn−2 − gn−1 z) + gn z n ] n is even. (m + n)! Since g j > 0 and g j > g j+1 (Lemma 3.2), if −1 ≤ z ≤ 1, then for all even k we have gk z k − gk+1 z k+1 = z k (gk − gk+1 z) ≥ z k (gk − gk+1 ) > 0, which implies Pn (z) > 0. Similarly, Q m (z) is greater than zero, for |z| ≤ 1. Hence Rmn = Q m (z)/Pn (z) > 0.



Definition 3.4. The Runge–Kutta method for Eq. (1.2) (scheme (3.1)) is called asymptotically stable at (α, β) if and only if ∃M > 0, for all m > M, h = m1 (i) (I − αh A) is invertible, (ii) for any given xkm,l (0 ≤ l ≤ m) by relation (3.3), such that limk→∞ X k = 0 where X k = (xkm,0 , xkm,1 , . . . , xkm,m ). The set S of all pairs (α, β) at which scheme (3.1) is asymptotically stable is called to be the asymptotical stability region. From (3.3), we know that xkm,l = (1 + β)(1 + h¯ bT (I − h¯ A)−1 e)m x(k−1)m,l , 0 ≤ l ≤ m. In the following we choose M = |α|. Theorem 3.5. Scheme (3.1) is said to be asymptotically stable at (α, β) if and only if |(1 + β)R m (h¯ )| < 1,

for all m > M

where R(h¯ ) = 1 + h¯ bT (I − h¯ A)−1 e is the stability function of the Runge–Kutta method. The following theorem will be as an immediate consequence. Theorem 3.6. H ⊆ S if and only if eh¯ ≥ |R(h¯ )|.

51

X.J. Ran et al. / Mathematical and Computer Modelling 48 (2008) 46–55 Table 1 The higher-order Runge–Kutta methods

(m, n) α > 0, H ⊆ S α < 0, H ⊆ S

Gauss

Lobatto IIIC

Lobatto IIIA, IIIB

Radau IA, IIA

(s, s) s is even s is odd

(s − 2, s) s is even s is odd

(s − 1, s − 1) s is odd s is even

(s − 1, s) s is even s is even

Lemma 3.7. Suppose Rmn (z) = Q m (z)/Pn (z) is the (m, n)-pad´e approximation to ez (i) if 1 > z > 0, then ez > Rmn (z) > 0 if and only if n is even, (ii) if 0 > z > −1, then ez > Rmn (z)mn > 0 if and only if m is odd. Proof. It is known from Lemma 3.3 Pn (z) > 0, Rmn (z) > 0 for −1 < z < 1, therefore ez − Rmn (z) has the same sign as Pn (z)ez − Q m (z), by Lemma 3.1, Z (−1)n m+n+1 1 n z u (1 − u)m euz du, (3.6) ez Pn (z) − Q m (z) = (m + n)! 0 the lemma is proved from (3.6).



The following theorem is obviously true by Lemma 3.7. Theorem 3.8. If the stability function of Runge–Kutta method is Rmn (h¯ ) = Q m (h¯ )/Pn (h¯ ) which is given by the (m, n)-pad´e approximation to exponential function eh¯ , h¯ = αh, then (i) if α > 0, H ⊆ S if and only if n is even, (ii) if α < 0, H ⊆ S if and only if m is odd. As a consequence of Theorem 3.8, we formulate the following corollary. Corollary 3.9. For the higher-order Runge–Kutta methods, it is easy to see from Theorem 3.8 (see Table 1): (i) for the s-stage Gauss and Lobatto IIIC methods, H ⊆ S if and only if s is even for α > 0, s is odd for α < 0, (ii) for the s-stage Lobatto IIIA and Lobatto IIIB methods, H ⊆ S if and only if s is odd for α > 0, s is even for α < 0, (iii) for the s-stage Radau IA and Radau IIA methods, H ⊆ S if and only if s is even. Proof. At the beginning of this section, we have discussed that the Runge–Kutta methods preserve the order of convergence. Hence Gauss scheme, Lobatto IIIC, Lobatto IIIA, IIIB and Radau IA, IIA, which belong to Runge–Kutta methods, should conserve order. As we know the stability functions of them are the pad´e approximation to exponential function ez , denoted by Rmn (z) = Q m (z)/Pn (z). Suppose all the above method are s-stage, it is true that Gauss scheme is of order m + n = 2s, where m = s, n = s, from Theorem 3.8, we obtain H ⊆ S if and only if s is even for α > 0 and s is odd for α < 0, respectively. Similarly, the results for Lobatto IIIC, Lobatto IIIA, IIIB and Radau IA, IIA can be gained which are displayed in Table 1.  As we know that all the s-stage explicit Runge–Kutta methods of order p = s = 1, 2, 3, 4 possess stability 3 p 2 function (see [8,11,13]) R(h¯ ) = 1 + h¯ + h¯2! + h¯3! + · · · + h¯p! which are (s, 0) approximations to exponential function eh¯ . Therefore, s-stage explicit Runge–Kutta methods preserve the asymptotical stability of the exact solution if and only if s is even for α < 0 and s = 1, 2, 3, 4 for α > 0. 4. Numerical experiments We consider the following example  ˙ = −5x, t 6= k, t > 0, x(t) 4x = −80x, t = k,  x(0 + 0) = 10,

(4.1)

52

X.J. Ran et al. / Mathematical and Computer Modelling 48 (2008) 46–55

Fig. 1. The numerical solution and the analytic solution of (4.1), the dashed curve is the numerical solution, the solid is the analytical solution.

which has the analytic solution in the following form x(t) = 10e−5{t} ((1 + (−80))e−5 )[t] ,

t ∈ (0, +∞).

(4.2)

We know the exact solution of (4.1) is asymptotically stable. Generally speaking, if θ-method with θ = 1 is applied to ODE, the numerical solution should be convergent. However, for impulsive differential equation (4.1), let θ = 1 and h = 1/m(m = 15), the numerical solution of (4.1) is divergent, see Fig. 1. This example indicates that the θ-method keeps the order of convergence but some properties have changed for system (1.2). The following examples will be helpful to illustrate the main results of this paper. Example 4.1.  t 6= k, t > 0  y˙ (t) = 1.2y 4y = −0.9y t = k  y(0 + 0) = 1. Example 4.2.  ˙ = −3.9x t 6= k, t > 0 x(t) 4x = 10x t =k  x(0 + 0) = 3.

(4.3)

(4.4)

The coefficient of problem (4.3) is α = 1.2, which is greater than zero. According to Table 1, we shall choose θ-methods, 3-Lobatto IIIB and 2-Radau IA methods to obtain the numerical solution of (4.3) with stepsize h = 1/m. In Eq. (4.4), we observed that α = −3.9 < 0, so that 2-Lobatto IIIB and 3-Gauss methods are applied. In Tables 2–5 the relative errors (R E) and absolute errors (AE) are listed. The true solutions of (4.3) at t = 8 and (4.4) at t = 4 are y ≈ 1.215104e − 001, x ≈ 6.703558e − 004, respectively. In order to obtain convergent solution of these two problems, by Theorem 2.4, we know that for (4.3), θ must satisfy 0 ≤ θ ≤ ϕ(1) ≈ 0.4180233, for (4.4), 0 ≤ θ ≤ ϕ(0) = 0.5. The ratios give a sound support to show that the methods preserve their order of convergence. There are four pictures in Fig. 2, in which the dashed curves are numerical solution, the solid curves are obtained by analytical expressions. The upper left picture is for Example 4.1, using θ-method with θ = 0; the upper right one is for Example 4.2, using θ-method with θ = 0.5; the lower left is for Example 4.1 using 2-Radau IA method (the absolute errors are so small that two curves are almost overlapped); the lower right is for Example 4.2 using 2-Lobatto IIIB method.

53

X.J. Ran et al. / Mathematical and Computer Modelling 48 (2008) 46–55

Fig. 2. The numerical solutions and analytic solution.

Table 2 Errors of problem (4.3) θ =0 AE

RE

θ = 0.4 AE

RE

2 4 10 20 100 200

7.856074e−002 5.496876e−002 2.845945e−002 1.571442e−002 3.422626e−003 1.730344e−003

6.465351e−001 4.523790e−001 2.342141e−001 1.293257e−001 2.816734e−002 1.424029e−002

1.636205e−002 1.226849e−002 6.081804e−003 3.262087e−003 6.901199e−004 3.474957e−004

1.346555e−001 1.009666e−001 5.005171e−002 2.684615e−002 5.679512e−003 2.859802e−003

Ratio

1.9780

1.9780

1.9860

1.9860

θ = 0.2 AE

RE

θ = 0.5 AE

RE

8 16 100 200 250 500

6.428826e−004 4.928432e−004 1.141830e−004 5.910432e−005 7.376674e−009 3.989841e−009

9.590170e−001 7.351964e−001 1.703320e−001 8.816857e−002 7.102379e−002 3.600510e−002

1.838000e−004 5.025963e−005 1.324487e−006 3.313107e−007 2.120533e−007 5.301817e−008

2.741827e−001 7.497455e−002 1.975797e−003 4.942310e−004 3.163295e−004 7.908959e−005

Ratio

1.9726

1.9726

3.9996

3.9996

m

Table 3 Errors of problem (4.4) m

54

X.J. Ran et al. / Mathematical and Computer Modelling 48 (2008) 46–55

Table 4 Errors of problem (4.3) 3-Lobatto IIIB AE

RE

2-Radau IA AE

RE

2 4 5 10 20 40

1.071861e−004 6.596533e−006 2.696801e−006 1.681199e−007 1.050075e−008 6.561920e−010

8.821141e−004 5.428780e−005 2.219399e−005 1.383584e−006 8.641851e−008 5.400295e−009

6.253585e−003 7.146806e−004 3.594719e−004 4.340316e−005 5.335311e−006 2.204902e−007

1.715514e−002 1.960547e−003 9.861209e−004 1.190657e−004 1.463609e−005 1.814579e−006

Ratio

16.0026

16.0026

8.0658

8.0658

2-Lobatto IIIB AE

RE

3-Gauss AE

RE

5 10 20 40 50 100

3.900707e−004 1.228127e−004 3.251276e−005 8.245050e−006 5.285872e−006 1.324487e−006

5.818861e−001 1.832053e−001 4.850075e−002 1.229951e−002 7.885173e−003 1.975797e−003

2.392048e−008 3.672177e−010 5.712409e−012 8.916153e−014 2.336651e−014 3.713392e−016

3.568325e−005 5.477952e−007 8.521458e−009 1.330063e−010 3.485687e−011 5.539435e−013

Ratio

3.9909

3.9909

64.0681

64.0681

m

Table 5 Errors of problem (4.4) m

Acknowledgements We would like to thank the National Natural Science Foundation of China No. 10671047. Meanwhile, this work is supported by the NSF of China under Grant 10572154. References [1] D.D. Bainov, P.S. Simenov, Systems with Impulse Effect Stability Theory and Applications, Ellis Horwood Limited, Chichester, 1989. [2] D.D. Bainov, S.I. Kostadinov, N. van Minh, P.P. Zabreiko, A topological classification of differential equations with impulse effect, Tamkang J. Math. 25 (1994) 15–27. [3] D.D. Bainov, P.S. Simenov, Systems with Impulse Effect, Stability Theory and Applications, vol. 40, Ellis Horwood Limited, Chichester, 1989. [4] S.D. Borisenko, V.I. Kosolapov, Yu.A. Obolenskii, Stability of Processes Under Continuations and Discrete Disturbances, Naukova Dumka, Kiev, 1988 (in Russian). [5] G. Kulev, D.D. Bainov, On the stability of systems with impulsive by sirect method of Lyapunov, J. Math. Anal. Appl. 140 (1989) 324–340. [6] D.D. Bainov, G. Kulev, Application of Lyapunov’s functions to the investigation of global stability of solutions of system with impulses, Appl. Anal. 26 (1) (1988) 255–270. [7] B.M. Randelovic, L.V. Stefanovic, B.M. Dankovic, Numerical solution of impulsive differential equations, Facta Univ. Ser. Math. Inform. 15 (2000) 101–111. [8] J.C. Butcher, The Numerical Analysis of Ordinary Differential Equations: Runge–Kutta Methods, Wiley, New York, 1987. [9] K. Dekker, J.G. Verwer, Stability of Runge–Kutta Methods for Stiff Nonlinear Differential Equations, North-Holland, Amsterdam, 1984. [10] Shen Jianhua, New maximum principles for first-order impulsive boundary value problems, Appl. Math. Lett. 16 (2003) 105–112. [11] E. Hairer, S.P. Nøsett, G. Wanner, Solving Ordinary Differential Equations II, Stiff and Differential Algebraic Problems, Springer, New York, 1993. [12] M.U. Akhmetov, A. Zafer, Successive approximation method for quasilinear impulsive differential equations with control, Appl. Math. Lett. 13 (2000) 99–105. [13] J.D. Lambert, Computational Methods in Ordinary Differential Equations, John Wiley, London, 1972. [14] V. Lakshmikantham, D.D. Bainov, P.S. Simenov, Theory of Impulsive Differential Equations, World Scientific, Singapore, New Jersey, London, Hong Kong, 1989, p. 273. [15] Alberto dˆaaˆ aˆ Onofrio, On pulse vaccination strategy in the SIR epidemic model with vertical transmission, Appl. Math. Lett. 18 (2005) 729–732.

X.J. Ran et al. / Mathematical and Computer Modelling 48 (2008) 46–55

55

[16] J.J. Nieto, Impulsive resonance periodic problems of first order, Appl. Math. Lett. 15 (2002) 489–493. [17] A.M. Samoilenko, N.A. Perestyuk, Differential Equations with Impulse Effect, Visca Skola, Kiev, 1987, p. 286 (in Russian). [18] M.H. Song, Z.W. Yang, M.Z. Liu, Stability of θ-methods for advanced differential equations with piecewise continuous arguments, J. Comput. Math. Appl. 49 (2005) 1295–1301. [19] X.Y. Xu, J.K. Li, G.L. Xu, Introduction to Pad`e Approximation, Shanghai Scientific & Technical Publishers, 1990, pp. 24–33 (in Chinese). [20] X.L. Fu, B.Q. Yan, Introduction to the Impulsive Differential System, Scientific Publishers, 2005, pp. 1–33 (in Chinese).