Variational iteration method for singular perturbation initial value problems

Computer Physics Communications 181 (2010) 947–956

Contents lists available at ScienceDirect

Computer Physics Communications www.elsevier.com/locate/cpc

Variational iteration method for singular perturbation initial value problems ✩ Yongxiang Zhao, Aiguo Xiao ∗ School of Mathematics and Computational Science, Hunan Key Laboratory for Computation and Simulation in Science and Engineering, Xiangtan University, Xiangtan, Hunan 411105, PR China

a r t i c l e

i n f o

a b s t r a c t

Article history: Received 25 October 2009 Received in revised form 2 January 2010 Accepted 30 January 2010 Available online 2 February 2010

In this paper, the variational iteration method (VIM) is applied to solve singular perturbation initial value problems (SPIVPs). The obtained sequence of iterates is based on the use of Lagrange multipliers. Some convergence results of VIM for solving SPIVPs are given. Moreover, the illustrative examples show the efficiency of the method. © 2010 Elsevier B.V. All rights reserved.

Keywords: Variational iteration method Initial value problems Singular perturbation Convergence

1. Introduction

Singular perturbation initial value problems often arise in many scientific and engineering fields such as automatic control, biology, medical science, economics, etc. (cf. [4,19]). These problems are characterized by a small parameter  multiplying the derivatives, and there exist initial boundary layers where the solutions change rapidly. The variational iteration method was first proposed by He [9–11], and has been extensively worked out over a number of years by numerous authors. This method solves the problems without any need to discrete the variables. Therefore, there is no need to compute the round off errors and one is not faced with necessity of large computer memory and time. Applications of the method have been enlarged due to its flexibility, convenience and efficiency. The VIM was first applied to autonomous ordinary differential systems, delay differential equations, and fractional differential equations by He et al. [7–9]. Salkuyeh [3] studied the convergence of VIM for solving linear system of ODEs with constant coefficients. Wazwaz [1] used the VIM for analytic treatment of linear and nonlinear ODEs. Saadatmandi, Dehghan [2] and Yu [21] applied the VIM to solve pantograph equations. Xu [14], Saadati, Dehghan, Vaezpour and Saravi [20] considered the convergence of VIM for solving integral equations. Though Darvishi, Khani and Soliman [16] applied the VIM to some stiff ODEs, but these stiff problems do not contain the singular perturbation problems. Tatari and Dehghan [17], Mamode [18] considered the VIM for solving second order initial value problems. Lu [13] applied VIM to solve two point boundary value problems. Rafei, Ganji, Daniali and Pashaei [15], Marinca, Herisanu and Bota [22] applied the VIM to Oscillations. For more comprehensive survey on this method and its applications, the reader is refer to the review articles [5,6] and the references therein. In this paper, we apply the VIM to SPIVPs to obtain the analytical or approximate analytical solutions. The convergence results of VIM for solving SPIVPs are obtained. Some illustrative examples confirm the theoretical results. In the rest parts of the text, we denote x(t ) = x(t ,  ), y (t ) = y (t ,  ) for simplicity, here,  is the singular perturbation parameter. The vectors (a1 , a2 , . . . , an ) T = a  b = (b1 , b2 , . . . , bn ) T means each component ai  b i (i = 1, 2, . . . , n).  ·  denotes the standard Euclidean norm of a vector.

✩ This work is supported by projects NSF of China (10971175), Specialized Research Fund for the Doctoral Program of Higher Education of China (20094301110001), NSF of Hunan Province (09JJ3002) and Hunan Provincial Innovation Foundation for Postgraduate (S2008yjscx02). Corresponding author. E-mail addresses: [email protected] (Y. Zhao), [email protected] (A. Xiao).

*

0010-4655/$ – see front matter doi:10.1016/j.cpc.2010.01.007

©

2010 Elsevier B.V. All rights reserved.

948

Y. Zhao, A. Xiao / Computer Physics Communications 181 (2010) 947–956

2. Convergence 2.1. Case 1 Consider the following singular perturbation initial value problem

⎧    ⎪ ⎨ x (t ) = f t , x(t ), y (t ) , 0  t  T ,    y  (t ) = g t , x(t ), y (t ) , 0    1, ⎪ ⎩ x(0) = x0 , y (0) = y 0 ,

(2.1)

where t ∈ [0, T ] is the ‘time’, x ∈ R n1 and y ∈ R n2 are the state variables,  is the singular perturbation parameter, 0 > 0 is a sufficiently small constant. f : [0, T ] × R n1 × R n2 → R n1 , g : [0, T ] × R n1 × R n2 → R n2 are given continuous mappings which satisfy the following Lipschitz conditions

   f (t , x1 , y 1 ) − f (t , x2 , y 2 )  l1 (t )x1 − x2  + l2 (t ) y 1 − y 2 ,    g (t , x1 , y 1 ) − g (t , x2 , y 2 )  k1 (t )x1 − x2  + k2 (t ) y 1 − y 2 ,

(2.2a) (2.2b)

where l1 (t ), l2 (t ), k1 (t ), k2 (t ) are continuous bounded functions. According to VIM we can construct the correction functionals as follows

t

   λ1 (s, t ) xn (s) − ˜f s, xn (s), yn (s) ds,

xn+1 (t ) = xn (t ) +

(2.3a)

0

  1   λ2 (s, t ) yn (s) − g˜ s, xn (s), yn (s) ds,

t yn+1 (t ) = yn (t ) +



(2.3b)

0

where λ1 (s, t ), λ2 (s, t ) are general Lagrange multipliers, which can be defined optimally via variational theory, and ˜f , g˜ denote the restrictive variation, i.e., δ ˜f = δ g˜ = 0. Thus, we have

t δ xn+1 (t ) = δ xn (t ) +

  λ1 (s, t ) δ xn (s) ds,

0

t δ yn+1 (t ) = δ yn (t ) +

  λ2 (s, t ) δ yn (s) ds,

0

and the stationary conditions are obtained as



∂λ1 (s, t ) = 0, ∂s ∂λ2 (s, t ) = 0. ∂s

1 + λ1 (s, t ) s=t = 0,



1 + λ2 (s, t ) s=t = 0,

Moreover, the general Lagrange multiplier can be readily identified by

λ1 (s, t ) = λ2 (s, t ) = −1. Therefore, the variational iteration formulas can be written as

t xn+1 (t ) = xn (t ) −







xn (s) − f s, xn (s), yn (s)

ds,

(2.4a)

0

t

yn+1 (t ) = yn (t ) −

yn (s) −



 1  g s, xn (s), yn (s) ds.



(2.4b)

0

Now, we show that the sequences {xn (t )}n∞=1 , { yn (t )}n∞=1 defined by (2.4) with x0 (0) = x0 , y 0 (0) = y 0 converge to the solution of (2.1). Theorem 1. Let x(t ), xi (t ) ∈ (C 1 [0, T ])n1 , y (t ), y i (t ) ∈ (C 1 [0, T ])n2 , i = 0, 1, . . . . The sequences defined by (2.4) with x0 (0) = x0 , y 0 (0) = y 0 converge to the solution of (2.1). Proof. Obviously from system (2.1) we have

Y. Zhao, A. Xiao / Computer Physics Communications 181 (2010) 947–956

t



x(t ) = x(t ) −





x (s) − f s, x(s), y (s)

ds,

949

(2.5a)

0

t

y (t ) = y (t ) −



 1  y (s) − g s, x(s), y (s) ds. 

(2.5b)



0

Introduce E i x(t ) = xi (t ) − x(t ), E i y (t ) = y i (t ) − y (t ), i = 0, 1, . . . , where E j x(0) = E j y (0) = 0, j = 0, 1, . . . . Now from (2.4)–(2.5) we obtain

t



E n+1 x(t ) = E n x(t ) −

 







E n x (s) − f s, xn (s), yn (s) − f s, x(s), y (s)

ds,

0

t

E n+1 y (t ) = E n y (t ) −

E n y  (s) −



   1  g s, xn (s), yn (s) − g s, x(s), y (s) ds.



0

Moreover, we can derive

t E n+1 x(t ) =

 







f s, xn (s), yn (s) − f s, x(s), y (s)

ds,

(2.6a)

0

E n+1 y (t ) =

1

t

 







g s, xn (s), yn (s) − g s, x(s), y (s)



ds.

(2.6b)

0

From the Lipschitz conditions (2.2), we have



 E n+1 x(t )  E n+1 y (t )





 t   E n x(s) ds 0 ,

k2 t  E n y (s) ds  0

l1



l2

k1



(2.7)

where li = max0s T li (s), ki = max0s T ki (s), i = 1, 2. Therefore



 E 1 x(t )  E 1 y (t )







 t   E 0 x(s) ds

0t k2   0  E 0 y (s) ds  l2 max0s T  E 0 x(s)t . k2 max0s T  E 0 y (s)t 

l1

l2

k1



l1 k1





Moreover, we can derive



 E n x(t )  E n y (t )





Tn



n! T



l2

k1

k2



n

( 0 ) n!

l1



n



max0s T  E 0 x(s)



max0s T  E 0 y (s)

O(0 ) O(0 ) O(1) O(1)



max0s T  E 0 x(s) max0s T  E 0 y (s)

 .

0 , T , max0sT  E 0 x(s), max0sT  E 0 y (s), li , ki , i = 1, 2, are constants. By using the Stirling’s formula, we have  



( T en/0 )n  E n x(t ) O(0 ) O(0 ) max0s T  E 0 x(s) √ , max0s T  E 0 y (s)  E n y (t ) 2π n(1 + O (1/n)) O (1) O (1)

(2.8)

Noting that



thus, ( E n x(t ),  E n y (t )) T → 0 as n → ∞.

(2.9)

2

2.2. Case 2 Consider the special case of (2.1)

⎧      ⎪ ⎨ x (t ) = Ax(t ) + F t , x(t ), y (t ) = f t , x(t ), y (t ) , 0  t  T ,      y  (t ) = B y (t ) + G t , x(t ), y (t ) = g t , x(t ), y (t ) , 0    1, ⎪ ⎩ x(0) = x0 , y (0) = y 0 ,

(2.10)

where F : [0, T ] × R n1 × R n2 → R n1 , G : [0, T ] × R n1 × R n2 → R n2 are given continuous mapping which satisfy the Lipschitz conditions (2.2), the matrices A = (ai j ) ∈ R n1 ×n1 , B = (b i j ) ∈ R n2 ×n2 can be decomposed into A = A 0 + A 1 , B = B 0 + B 1 respectively, where A 0 = diag(a11 , a22 , . . . , an1 n1 ) and B 0 = diag(b11 , b22 , . . . , bn2 n2 ). It’s easy to show that the right-hand sides of (2.10) also satisfy the Lipschitz conditions. If the right-hand sides of (2.10) are consider as nonlinear terms, then, we can also use the correction functionals constructed in Case 1, and get the similar results as Theorem 1. Now, we construct the following correction functionals

950

Y. Zhao, A. Xiao / Computer Physics Communications 181 (2010) 947–956

t

   Λ1 (s, t ) xn (s) − A 0 xn (s) − A 1 x˜ n (s) − F˜ s, xn (s), yn (s) ds,

xn+1 (t ) = xn (t ) + 0

   1  ˜ Λ2 (s, t ) yn (s) − B 0 yn (s) + B 1 y˜ n (s) + G s, xn (s), yn (s) ds,

t yn+1 (t ) = yn (t ) +



0

where Λ1 (s, t ) = diag(λ11 (s, t ), λ12 (s, t ), . . . , λ1n1 (s, t )), Λ2 (s, t ) = diag(λ21 (s, t ), λ22 (s, t ), . . . , λ2n2 (s, t )), in which λ1i (s, t ), λ2 j (s, t ), i = 1, 2, . . . , n1 , j = 1, 2, . . . , n2 , are general Lagrange multipliers and x˜ n , y˜ n , F˜ , G˜ denote the restrictive variations, i.e., δ x˜ n = δ y˜ n = δ F˜ = δ G˜ = 0. Thus, we have

t

  Λ1 (s, t ) δ xn (s) − A 0 δ xn (s) ds,

δ xn+1 (t ) = δ xn (t ) + 0

t δ yn+1 (t ) = δ yn (t ) +

 B0 Λ2 (s, t ) δ yn (s) − δ yn (s) ds,



0

and the stationary conditions are obtained as



∂Λ1 (s, t ) + A 0 Λ 1 ( s , t ) = 0, ∂s ∂Λ1 (s, t ) B 0 + Λ 2 ( s , t ) = 0. ∂s 

1 + Λ1 (s, t ) s=t = 0,



1 + Λ2 (s, t ) s=t = 0,

Moreover, the general Lagrange multipliers can be readily identified by

 B0 Λ2 (s, t ) = − exp − (s − t ) .



 Λ1 (s, t ) = − exp − A 0 (s − t ) ,



Therefore, the variational iteration formula can be written as

t xn+1 (t ) = xn (t ) −







e (− A 0 (s−t )) xn (s) − Axn (s) − F s, xn (s), yn (s)

ds,

(2.11)

0

t yn+1 (t ) = yn (t ) −

e (−

B0

 (s−t ))



yn (s) −

1









B yn (s) + G s, xn (s), yn (s)

ds.

(2.12)

0

The following theorem shows that the sequences {xn (t )}n∞=1 , { yn (t )}n∞=1 defined by (2.11)–(2.12) with x0 (0) = x0 , y 0 (0) = y 0 converge to the solution of (2.10). Theorem 2. Let x(t ), xi (t ) ∈ (C 1 [0, T ])n1 , y (t ), y i (t ) ∈ (C 1 [0, T ])n2 , i = 0, 1, . . . . The sequences defined by (2.11)–(2.12) with x0 (0) = x0 , y 0 (0) = y 0 converge to the solutions of (2.10). Proof. By a similar process to the proof of Theorem 1, we can easily obtain



 E n x(t )  E n y (t )



 e ( ) √

Tn n!



 A 1  + l1 k1



T e /0 n 0 )( n )

e (

2π n(1 + O (1/n))

l2

n

 B 1 +k2





max0s T  E 0 x(s)



max0s T  E 0 y (s)

O(0 ) O(0 ) O(1) O(1)



max0s T  E 0 x(s) max0s T  E 0 y (s)

 ,

(2.13)

B0 

where e ( ) = max0st ,0t  T (e − A 0 (s−t ) , e −  (s−t ) ), li = max0s T li (s), ki = max0s T ki (s), i = 1, 2. Noting that 0 , T ,  A 1 ,  B 1 , e (0 ), li , ki , i = 1, 2, max0s T  E 0 x(s) and max0s T  E 0 y (s) are constants, we can derive from (2.13) that ( E n x(t ),  E n y (t )) T → 0 as n → ∞. 2 3. Illustrative examples In this section, some illustrative examples are given to show the efficiency of the VIM for solving SPIVPs. Example 1. Consider the linear SPIVP (cf. [19])



 y  (t ) = (t − 1) y (t ), 0 <   1, 0  t  0.9, y (0) = 0.2.

Using the VIM in the previous section, we construct the following correction functional as Case 2

(3.1)

Y. Zhao, A. Xiao / Computer Physics Communications 181 (2010) 947–956

t yn+1 (t ) = yn (t ) +

 s−1  λ(s, t ) yn (s) − yn (s) ds.



951

(3.2)

0

To find the optimal value λ, we have

t δ yn+1 (t ) = δ yn (t ) + δ

 s−1  λ(s, t ) yn (s) − yn (s) ds,



(3.3)

0

which leads to

δ yn+1 (t ) = δ yn (t ) + λ(s, t )δ yn (s) s=t −

t

 ∂λ(s, t ) s − 1 + λ(s, t ) δ yn (s) ds. ∂s 

(3.4)

0

Moreover, the stationary conditions are as follows



∂λ(s, t ) s − 1 + λ(s, t ) = 0. ∂s 

1 + λ(s, t ) s=t = 0,

(3.5)

Therefore, the general Lagrange multiplier can be readily identified by

λ(s, t ) = − exp

(2 − s − t )(s − t )



2

(3.6)

.

Substituting (3.6) into (3.2) results in the following iteration formula

t yn+1 (t ) = yn (t ) −

e

(2−s−t )(s−t )



2

s−1

yn (s) −



 yn (s) ds.

(3.7)

0

To get iterate sequence, we start with an initial approximation y 0 (t ) = y (0) = 0.2. By means of the formula (3.7), we have



t2

t



t

y 1 (t ) = 0.2 − 0.2 −e 2 + e  e −  = 0.2e −(1−t /2)t / , which is the exact solution of (3.1). Example 2. Consider the nonlinear SPIVP (cf. [19])



 y  (t ) = y (t ) − y 3 (t ), 0 <   1, t  0, y (0) = 12 .

(3.8)

Using the VIM in the previous section, we construct the following correction functional as Case 1

t yn+1 (t ) = yn (t ) +

  1 λ(s, t ) yn (s) + − y˜ n (s) + y˜ n3 (s) ds.



(3.9)

0

To find the optimal value λ(s, t ), we have

t δ yn+1 (t ) = δ yn (t ) + δ

  1 λ(s, t ) yn (s) + − y˜ n (s) + y˜ n3 (s) ds,



(3.10)

0

which leads to

t δ yn+1 (t ) = δ yn (t ) + λ(s, t )δ yn (s)|s=t −

∂λ(s, t ) δ yn (s) ds. ∂s

(3.11)

0

Moreover, the stationary conditions are as follow



1 + λ(s, t ) s=t = 0,

∂λ(s, t ) = 0. ∂s

(3.12)

Therefore, the general Lagrange multiplier can be readily identified by

λ(s, t ) = −1. Substituting (3.13) into (3.9) results in the following iteration formula

(3.13)

952

Y. Zhao, A. Xiao / Computer Physics Communications 181 (2010) 947–956

t

yn+1 (t ) = yn (t ) −



y n (s) +

1

− y n (s) +







yn3 (s)

ds.

(3.14)

0

To get iterate sequence, we start with an initial approximation y 0 (t ) = y (0) = 12 . By means of formula (3.14), we have

1

y 1 (t ) =

+

2 1

y 2 (t ) =

+

2 1

y 3 (t ) =

+

2

+

y 5 (t ) =

···

,

8 3t 8 3t

− −

27 t 4

−

2048  4 63 t 4 4

8 2048  t 13 19683

−

9 t3 128  3 17 t 3 3

+ +

3 t2 64  2 3 t2 2

, +

256  64  t 12 6561

+

27 t 5 5

+

567 t 6

2560  65536  6 t 11 45927

+

11669149696  13 2147483648  12 2952790016  11 t 10 729 1377 t 9 1701 t 8 1917 t 7

− y 4 (t ) =

3t

67108864  10 1 3t 3 t2

+

2 1

+

2

+

8 3t 8

+

−

+

−

64  2 3 t2

256  3 17 t 3

−

64  2

+

4194304  9 2097152  8 1835008  7

5 3 4 17 t 125 t t

256  3

−

4096  4 125 t 4

−

4096  4

+O +

,

,

5

721 t 5 81920  5

+O

t6

 ,

6

In fact, (3.8) is a Bernoulli equation and the exact solution is



y (t ) = 1



1 + 3 exp −



2t

.





By executing the command series(1/ 1 + 3 exp(− 2t  ), t = 0) in Maple, and comparing with the iterate sequence, we can obtain that



lim yn (t ) = 1



1 + 3 exp −

n→∞

2t



 .

(3.15)

Example 3. Consider the linear SPIVP (cf. [23])

⎧  ⎨  y (t ) = 2x(t ) − y (t ), 0 <   1, x (t ) = y (t ) − 2x(t ), t  0, ⎩ y (0) = 2.3, x(0) = 1.1.

(3.16)

The exact solution of the system is

y (t ) = x(t ) =



2 2 + 1 1 2 + 1





 y (0) + x(0) +

 2 + 1





y (0) − 2x(0) exp −(2 + 1)

t



, 

   t . y (0) − 2x(0) exp −(2 + 1)

2 + 1



 y (0) + x(0) −



1



Using the VIM in the previous section, we construct the following correction functionals as Case 2



t yn+1 (t ) = yn (t ) −

exp

s−t





y n (s) +



1



  yn (s) − 2xn (s) ds,

(3.17)

0

t xn+1 (t ) = xn (t ) −







exp 2(s − t ) xn (s) − yn (s) + 2xn (s) ds.

(3.18)

0

To get iterate sequences, we start with the initial approximations y 0 (t ) = y (0) = 2.3, x0 (t ) = x(0) = 1.1. By means of formulas (3.17)– (3.18), we have



t



y 1 (t ) = 2.3 − 0.1 1 − e −  , x1 (t ) = 1.1 −

0.1 (e −2t − e −t / ) 2 − 1



y 2 (t ) = 2.3 − 0.1 1 − e

− t



, t

t

t

−e −  +  e −2t − te −  + 2t  e −  + 0.2 , 4 2 − 4 + 1

Y. Zhao, A. Xiao / Computer Physics Communications 181 (2010) 947–956

953

Fig. 3.1. Results for Problem 3.

x2 (t ) = 1.1 −

0.1 (e −2t − e −t / ) 2 − 1 t

t

t

t

t

−  e −  − t  e −2t − te −  + 2t  e −  + 2t  e −2t , 8 3 − 12 2 + 6 − 1 t t t  t  −e −  +  e −2t − te −  + 2t  e −  y 3 (t ) = 2.3 − 0.1 1 − e −  + 0.2 4 2 − 4 + 1  2 −t t t + 0.2 6 e  + 4t 2  2 e −  − 6 2 e −2t − 8t  2 e −  + 2t  e −2t − 4t  2 e −2t  t t t   − 4t 2  e −  + t 2 e −  + 4t  e −  / 16 4 − 32 3 + 24 2 − 8 + 1 , + 0.2

x3 (t ) = 1.1 −

2 e −2t

0.1 (e −2t − e −t / )

+ 0.2

2 − 1 2 e −2t

t

−  e −  − t  e −2t − te −  + 2t  e −  + 2t  e −2t 8 3 − 12 2 + 6 − 1

 t + 0.2 −12 2 e −2t − t 2 e −2t − 12t  2 e −  − 4t 2  2 e −2t + 6t  e −2t t

t

t

t

− 4t 2  e −  + 6t  e −  − 12t  2 e −2t + t 2 e −  + 12 2 e −  + 4t 2  e −2t  t   + 4t 2  2 e −  / 32 5 − 80 4 + 80 3 + 40 2 + 10 − 1 . The exact and approximate solutions are plotted in Fig. 3.1. Here, we use y 3 , x3 as the approximate solution, and the parameter

 = 10−3 , 10−4 , 10−5 , respectively. Fig. 3.1 shows that the method gives very good approximation of the exact solution. Example 4. Consider the problem of thermal decomposition of ozone

⎧    ⎪ ⎨  x (t ) = − x(t ) + x(t ) y (t ) /1000 + 294 y (2), t  0,   y  (t ) = x(t ) 1 − y (t ) /98 − 3 y (t ), 0 <   1, ⎪ ⎩ x(0) = 1, y (0) = 0.

(3.19)

Using the VIM in the previous section, we construct the following correction functionals as Case 2



t xn+1 (t ) = xn (t ) −

exp

s−t





xn (s) +



1



xn (s)

 1000

+

xn (s) yn (s) 1000

 − 0.294 yn (s) ds,

(3.20)

0

t yn+1 (t ) = yn (t ) − 0





exp(s − t ) yn (s) − xn (1 − yn ) + 3 yn ds.

(3.21)

954

Y. Zhao, A. Xiao / Computer Physics Communications 181 (2010) 947–956

Fig. 3.2. Results for Problem 4.

To get iterate sequences, we start with the initial approximations x0 (t ) = x(0) = 1, y 0 (t ) = y (0) = 0 and formulas (3.20)–(3.21), we have

 = 0.001. By means of

x1 (t ) = e −t , 1 −t 1 −3t y 1 (t ) = e − e , 196 196 x2 (t ) ≈ 0.24659864e −t + 0.052721088e −2t + 0.75e −3t − 0.0017006803e −4t + 1.5te −t , y 2 (t ) ≈ −0.0025683743e −t + 0.00011731577e −2t + 0.0024987596e −3t

+ 0.0076530612te −t − 0.000078092461te −2t + 0.0076527956te −3t     + O e −4t + O te −4t , x3 (t ) ≈ 1.9119428e −t + 0.056447893e −2t − 0.92977056e −3t + 0.0076126967e −4t

− 0.7551020te −t + 0.043953040te −2t − 1.1249216te −3t + 0.013995562te −4t + 1.125t 2 e −t + 0.011479592t 2 e −2t − 5.8569346 ∗ 10−5t 2 e −3t + 0.0038263978t 2 e −4t       + O e −5t + O te −5t + O t 2 e −5t , y 3 (t ) ≈ 0.014550989e −t + 0.00079910065e −2t − 0.015471811e −3t + 0.00015224102e −4t

− 0.0095923574te −t − 0.00065890238te −2t − 0.0094895121te −3t + 0.00015679384te −4t + 0.0057397959t 2 e −t + 0.00046915242t 2 e −2t − 0.0057379512t 2 e −3t + 0.00010631823t 2 e −4t − 8.7854019 ∗ 10−5 t 3 e −2t − 1.8925469 ∗ 10−6t 3 e −3t + 8.7837248 ∗ 10−5 t 3 e −4t         + O e −5t + O te −5t + O t 2 e −5t + O t 3 e −5t . The approximate solutions and the numerical solutions computed by parallel multistep hybrid method (PMHM) (cf. [23]) are plotted in Fig. 3.2. Here, we use x3 , y 3 as the approximate solution, and the parameter  = 10−3 . Example 5. Consider the nonlinear problem from enzyme kinetics (cf. [19])

⎧      ⎪ ⎨ x (t ) = − x(t ) + x(t ) + κ − λ y (t ) , t  0,    y  (t ) = x(t ) − x(t ) + κ y (t ), 0 <   1, ⎪ ⎩ x(0) = 1, y (0) = 0, where κ , λ are bounded positive constants which satisfy κ < λ, the parameter  ∼ 10−6 . Using the VIM in the previous section, we construct the following correction functionals as Case 2

(3.22)

Y. Zhao, A. Xiao / Computer Physics Communications 181 (2010) 947–956

955

Fig. 3.3. Results for Problem 5.

t xn+1 (t ) = xn (t ) −

    exp(s − t ) xn (s) + xn (s) − xn (s) + κ − λ yn (s) ds,

(3.23)

0

t yn+1 (t ) = yn (t ) −



κ (s − t ) exp 





y n (s) −

1



  xn (s) − xn (s) + κ yn (s) ds. 



(3.24)

0

To get iterate sequences, we start with the initial approximations x0 (t ) = x(0) = 1, y 0 (t ) = y (0) = 0 and (3.23)–(3.24), we have for κ = 19, λ = 20

x1 (t ) = e −t , y 1 (t ) =

1000000 −t 1000000 −19000000t , e − e 18999999 18999999

 = 10−6 . By means of formulas

956

Y. Zhao, A. Xiao / Computer Physics Communications 181 (2010) 947–956

  x2 (t ) ≈ e −t − 0.05263158172e −2t + O e −19000000t , y 2 (t ) ≈ 0.05540166526e −t − 0.005685960965e −2t + 0.000145793885e −3t

  − 0.002770083395te −t + 0.0001457938782te −2t + O e −19000000t ,

and for

κ = 14, λ = 15

x1 (t ) = e −t , 1000000 −t 1000000 −14000000t y 1 (t ) = , e − e 13999999 13999999

  x2 (t ) ≈ e −t − 0.07142857653e −2t + O e −14000000t , y 2 (t ) ≈ 0.07653061844e −t − 0.01056851544e −2t + 0.0003644316170e −3t

  − 0.05102041545te −t + 0.000364431591te −2t + O e −14000000t .

The approximate solutions for κ = 19, λ = 20 and solutions, and the parameter  = 10−6 .

κ = 14, λ = 15 are plotted in Fig. 3.3. Here, we use x2 , y 2 as the approximate

4. Conclusion The VIM used in this paper is the variational iteration algorithm-I, there are also variational iteration algorithm-II and variational iteration algorithm-III [12]. In this paper, we apply the VIM to obtain the analytical or approximate analytical solutions of SPIVPs. The convergence results of VIM for solving SPIVPs are given. The illustrative examples show the efficiency of the method. When considering the system (2.10), the choice of correction functionals of Case 1 or Case 2 rely on the practical problems and this choice will result in the difference of the speed of convergence. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20]

A.M. Wazwaz, The variational iteration method for analytic treatment for linear and nonlinear ODEs, Appl. Math. Comput. 212 (2009) 120–134. A. Saadatmandi, M. Dehghan, Variational iteration method for solving a generalized pantograph equation, Comput. Math. Appl., in press, doi:10.1016/j.camwa.2009.03.017. D.K. Salkuyeh, Convergence of the variational iteration method for solving linear systems of ODEs with constant coefficients, Comput. Math. Appl. 56 (2008) 2027–2033. E.M. Dejager, F.R. Jiang, The Theory of Singular Perturbation, Elsevier Science B.V., The Netherlands, 1996. J.H. He, X.H. Wu, Variational iteration method: New development and applications, Comput. Math. Appl. 54 (2007) 881–894. J.H. He, Variational iteration method: some recent results and new interpretations, J. Comput. Appl. Math. 207 (2007) 3–17. J.H. He, X.H. Wu, Variational iteration method for autonomous ordinary differential systems, Appl. Math. Comput. 114 (2000) 115–123. J.H. He, Variational iteration method for delay differential equations, Commun. Nonlinear Sci. Numer. Simul. 2 (4) (1997) 235–236. J.H. He, Some applications of nonlinear fractional differential equation and their approximations, Bull. Sci. Technol. 15 (12) (1999) 86–90. J.H. He, A new approach to linear partial differential equations, Commun. Nonlinear Sci. Numer. Simul. 2 (4) (1997) 230–235. J.H. He, Approximate solution of nonlinear differential equation with convolution product nonlinearities, Comput. Methods Appl. Mech. Engrg. 167 (1998) 69–73. J.H. He, G.C. Wu, F. Austin, The variational iteration method which should be followed, Nonlinear Sci. Lett. A 1 (1) (2010) 1–30. J.F. Lu, Variational iteration method for solving two point boundary value problems, J. Comput. Appl. Math. 207 (2008) 92–95. L. Xu, Variational iteration method for solving integral equations, Comput. Math. Appl. 54 (2007) 1071–1078. M. Rafei, D.D. Ganji, H. Daniali, H. Pashaei, The variational iteration method for nonlinear oscillators with discontinuities, J. Sound Vibration 30 (6) (1992) 1338–1360. M.T. Darvishi, F. Khani, A.A. Soliman, The numerical simulation for stiff systems of ordinary differential equations, Comput. Math. Appl. 54 (2007) 1055–1063. M. Tatari, M. Dehghan, On the convergence of He’s variational iteration method, J. Comput. Appl. Math. 207 (2007) 121–128. M. Mamode, Variational iteration method and initial value problems, Appl. Math. Comput., in press, doi:10.1016/j.amc.2009.05.008. R.E. O’Malley Jr., Singular Perturbation Methods for Ordinary Differential Equations, Springer-Verlag, New York, 1990. R. Saadati, M. Dehghan, S.M. Vaezpour, M. Saravi, The convergence of He’s variational iteration method for solving integral equations, Comput. Math. Appl., in press, doi:10.1016/j.camwa.2009.03.008. [21] Z.H. Yu, Variational iteration method for solving the multi-pantograph delay equation, Phys. Lett. A 372 (2008) 6475–6479. [22] V. Marinca, N. Herisanu, C. Bota, Application of the variational iteration method to some nonlinear one dimensional oscillations, Meccanica 43 (2008) 75–79. [23] A.G. Xiao, Y.X. Zhao, Convergence of parallel multistep hybrid methods for singular perturbation problems, Appl. Math. Comput. 215 (2009) 2139–2148.