On constructing accurate approximations of first integrals for difference equations

Commun Nonlinear Sci Numer Simulat 18 (2013) 835–850

Contents lists available at SciVerse ScienceDirect

Commun Nonlinear Sci Numer Simulat journal homepage: www.elsevier.com/locate/cnsns

On constructing accurate approximations of first integrals for difference equations M. Rafei ⇑, W.T. Van Horssen Delft Institute of Applied Mathematics (DIAM), Delft University of Technology, Mekelweg 4, 2628 CD Delft, The Netherlands

a r t i c l e

i n f o

Article history: Received 28 March 2011 Received in revised form 11 September 2012 Accepted 11 September 2012 Available online 24 September 2012 Keywords: Invariance vector First integrals Functional equation Nonlinear difference equation Multiple scales perturbation method

a b s t r a c t In this paper, a perturbation method based on invariance factors and multiple scales will be presented for weakly nonlinear, regularly perturbed systems of ordinary difference equations. Asymptotic approximations of first integrals will be constructed on long iterationscales, that is, on iteration-scales of order e1 , where e is a small parameter. It will be shown that all invariance factors have to satisfy a functional equation. To show how this perturbation method works, the method is applied to a Van der Pol equation, and a Rayleigh equation. It will be explicitly shown for the first time in the literature how these multiple scales should be introduced for systems of difference equations to obtain very accurate approximations of first integrals on long iteration-scales. Ó 2012 Elsevier B.V. All rights reserved.

1. Introduction Nowadays, the multiple time-scales perturbation method for differential equations is a well developed, well accepted, and a very popular method to approximate solutions of weakly nonlinear differential equations. This method was developed in the period 1935–1962 by Krylov and Bogoliubov, Kuzmak, Kevorkian, and Cole, Cochran, and Mahony. In the early 1970s, Nayfeh popularized this method by writing many papers and books on this subject (see for instance [1]). More recent books on this method and its historical development are, for instance, the books by Andrianov and Manevitch [2], Holmes [3], Kevorkian and Cole [4], Murdock [5], and Verhulst [6]. For scientists and engineers the analysis of nonlinear dynamical systems is an important field of research since the solutions of these systems can exhibit counterintuitive and sometimes unexpected behavior. To obtain useful information from these systems the construction of (approximations of) first integrals by means of computing (approximate) integrating factors can play an important role. In [7] a perturbation method based on integrating vectors and multiple scales was presented for regularly perturbed differential equations. Recently, first integrals, invariants and Lie group theory for ordinary difference equations (ODEs) obtained a lot of attention in the literature (see for instance the list of references in [8]). Van Horssen [9] developed a perturbation method for a single first order difference equation based on invariance factors. Also recently, the fundamental concept of invariance factors for ODEs to obtain invariants (or first integrals) for ODEs has been presented in [8]. In [10] a perturbation method based on invariance factors has been presented for regularly perturbed systems of ordinary difference equations (ODEs). Using straightforward (but naive) approximations for the invariance vectors it has been shown how approximations of first integrals (including error-estimates) can be obtained. It turned out that some of the approximations for the first integrals are valid ⇑ Corresponding author. E-mail addresses: [email protected] (M. Rafei), [email protected] (W.T. Van Horssen). 1007-5704/$ - see front matter Ó 2012 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.cnsns.2012.09.010

836

M. Rafei, W.T. Van Horssen / Commun Nonlinear Sci Numer Simulat 18 (2013) 835–850

on the iteration scales for which the original problem for the system of ODEs has been proved to be well-posed. However, some of the approximations are only valid on a much smaller iteration scale. This is caused by the use of straightforward, naive approximations for the invariance factors, and as a well-known consequence in perturbation theory, the so-called secular terms occur in the approximations. To avoid secular terms in the approximations for the solution of a difference equation, the well-known perturbation method based on multiple (iteration-) scales can sometimes be used successfully. The development of the multiple scales perturbation method for ordinary difference equations (ODEs) started in 1960 with the work of Torng [11]. From the results obtained in [11], it is clear that the solution of a weakly perturbed (non) linear ODEs behaves differently on different iteration scales. In 1977, Hoppensteadt and Miranker introduced in [12] the multiple scales perturbation method for ODEs. It is interesting to notice that for difference equations a formulation of the multiple-scales methods completely in terms of difference operators has also recently been developed in [13]. To apply this method, an additional, new (iteration-) scale is usually introduced in the approximation for the solution. For instance sn ¼ en, where n is the original fast iteration-scale, and where sn is the new slow iteration-scale. In the multiple (iteration-) scales perturbation method these variables n and sn will be treated as independent variables. Now it should be observed that snþ1 ¼ sn þ e, and so by introducing a new (iteration-) scale, one adds in fact a new difference equation to the original system of ODEs (note that sn will now be treated as a variable dependent on n). This simple observation led in [13] to a perturbation method based on invariance vectors and multiple scales. It turned out that this perturbation method is consistent and straightforward. Moreover, it turned out that secular terms can be avoided and that error-estimates can easily be obtained on long iteration-scales. It has been shown in [8] that in finding invariants for a system of first order difference equations all invariance factors have to satisfy a functional equation (for more information on functional equations and how to solve some of them we refer the reader to [14–19]). The aim of this paper is to construct asymptotic approximations of first integrals for a system of first order ODEs. After presenting the concepts, we will explicitly show how highly accurate approximations of invariants for a second order, weakly nonlinear difference equation (with a Van der Pol type of nonlinearity), and also for a weakly nonlinear Rayleigh equation can be constructed by using the multiple scales perturbation method. Finally, it should be remarked that the advantage of the method of invariance factors is that the method can be applied to problems which cannot be integrated or be solved completely when the small parameter is zero. For such perturbed, essentially nonlinear problems one can find approximations of some of the first integrals by applying the perturbation method based on invariance factors. Classical perturbation methods can usually only be applied to problems which are linear when the small parameter is zero. So, the perturbation method as presented in this paper can be applied to a much larger class of problems as compared to the classical perturbation methods. The outline of this paper is as follows. In Section 2 the concept of invariance factors for a system of first order ODEs will be given, and in Section 3 some of the first results in the development of a multiple scales perturbation method for a system of ODEs based on invariance factors will be presented. In Sections 4 and 5 approximations of first integrals for systems of weakly nonlinear ODEs of Van der Pol type, and of Rayleigh type, respectively, will be constructed. Finally, in Section 6 of this paper some conclusions will be drawn and some future directions for research will be indicated. 2. On invariance factors for ODEs In this section we are going to present an overview of the concept of invariance factors (or vectors) for a system of k first order ODEs (where k is fixed and k 2 N), for a reference concerning this concept, we refer the reader to [9]. Consider

xnþ1 ¼ f ðxn ; nÞ;

ð1Þ

for n ¼ 0; 1; 2; . . ., and where xn ¼ ðx1;n ; x2;n ; . . . ; xk;n ÞT ; f ¼ ðf1 ; f2 ; . . . ; fk ÞT in which the superscript indicates the transpose, and where fi ¼ fi ðxn ; nÞ ¼ fi ðx1;n ; x2;n ; . . . ; xk;n ; nÞ are sufficiently smooth functions (for i ¼ 1; 2; . . . ; k). We also assume that an invariant for (1) can be represented by

Iðxnþ1 ; n þ 1Þ ¼ Iðxn ; nÞ ¼ constant () DIðxn ; nÞ ¼ 0: We now try to find an invariant for (1). By multiplying each ith equation, in (1), with a factor li ðf ðxn ; nÞ; n þ 1Þ for i ¼ 1; 2; . . . ; k and by adding the resulting equations we obtain

lðxnþ1 ; n þ 1Þ  xnþ1 ¼ lðf ðxn ; nÞ; n þ 1Þ  f ðxn ; nÞ;

ð2Þ

li ðxnþ1 ; n þ 1Þ ¼ ð3Þ

T

where l ¼ ðl1 ; l2 ; . . . ; lk Þ , and li ¼ li ðxn ; nÞ ¼ li ðx1;n ; x2;n ; . . . ; xk;n ; nÞ for i ¼ 1; 2; . . . ; k. In fact, when l is an invariance vector we now obtain an exact difference equation (2). The relationship between I and l follows from the equivalence of (2) and (3), yielding

(

Iðxnþ1 ; n þ 1Þ ¼ lðxnþ1 ; n þ 1Þ  xnþ1 ; Iðxn ; nÞ ¼ lðf ðxn ; nÞ; n þ 1Þ  f ðxn ; nÞ:

ð4Þ

By reducing the index n þ 1 by 1 in the first part of (4) I can be eliminated from (4), and then it follows that all invariance vectors for the system of difference equations (1) have to satisfy the functional equation

lðxn ; nÞ  xn ¼ lðf ðxn ; nÞ; n þ 1Þ  f ðxn ; nÞ:

ð5Þ

M. Rafei, W.T. Van Horssen / Commun Nonlinear Sci Numer Simulat 18 (2013) 835–850

837

When an invariance vector has been determined from the functional equation (5) an invariant for (1) easily follows from (4), yielding

Iðxn ; nÞ ¼ lðxn ; nÞ  xn :

ð6Þ

Finding an invariance vector for a given system of first order difference equations is a difficult and usually impossible task. On the other hand, we can use invariance vectors of some special form, and so we can obtain invariants for special classes of systems of k first order difference equations. Examples of this approach have been given in [8]. 3. A perturbation method based on invariance factors and multiple scales for systems of first order ODEs In [10] a perturbation method based on invariance factors has been presented for regularly perturbed problems. Using straightforward, but naive approximations for the invariance factors, it has been shown how approximations of first integrals can be obtained. It turns out that some of the approximations for the invariance factors contain so-called secular terms, and as a consequence some of the approximations for the first integrals are valid only on an iteration-scale that is much smaller than the iteration-scale for which the original problem has been proved to be well-posed. For a regularly perturbed problem it is well known that the solution of the problem can behave differently on different iteration-scales. For instance the solution can behave differently on iteration-scales of order 1 and of order e1 . In general no secular terms may occur in the approximations to describe this behavior correctly on long iteration-scales. To avoid those secular terms one can sometimes introduce successfully an additional new (slow) iteration-scale sn ¼ en. This is a short description of the well-known method of two iteration-scales. In some sense this method is not consistent since n and sn will be treated as independent variables. Now it should be observed that snþ1 ¼ sn þ e. When the new iteration-scale sn is treated as being given by a difference equation, we use the variable sn in a consistent way, that is, sn depends on n. In fact by introducing one (or more) new (iteration-) scale(s) we add one (or more) new difference equation(s) to the original system of ODEs. This simple observation will lead to a new, consistent, and straightforward perturbation method based on invariance factors and multiple scales. In this section some of the first results in the development of a perturbation method for difference equation based on invariance factors and multiple scales will be presented. The analysis in the first part of this section will be restricted to the following single, first order ODE (see also [9]):



xnþ1 ¼ f ðxn ; n; eÞ;

ð7Þ

snþ1 ¼ sn þ e; where

e is a small parameter, and sn ¼ en. In most applications the function f has the form f ðxn ; n; eÞ ¼ f0 ðxn ; nÞ þ ef1 ðxn ; nÞ þ e2 f2 ðxn ; nÞ þ   

ð8Þ

It is obvious from (5) that the functional equation for system (7) becomes

l1 ðxn ; n; sn ; eÞxn þ l2 ðxn ; n; sn ; eÞsn ¼ l1 ðxnþ1 ; n þ 1; sn þ e; eÞxnþ1 þ l2 ðxnþ1 ; n þ 1; sn þ e; eÞsnþ1 :

ð9Þ

Since we are interested in simple invarinace factors, we assume that

l2 ðxn ; n; sn ; eÞ ¼

1

sn

ð10Þ

:

Therefore, functional Eq. (9) becomes

l1 ðxn ; n; sn ; eÞxn ¼ l1 ðxnþ1 ; n þ 1; sn þ e; eÞxnþ1 : Now it will also be assumed that

ð11Þ

l1 can be expanded in power series in e, that is,

l1 ðxn ; n; sn ; eÞ ¼ l1;0 ðxn ; n; sn Þ þ el1;1 ðxn ; n; sn Þ þ e2 l1;2 ðxn ; n; sn Þ þ   

ð12Þ

The expansions (8) and (12) are then substituted into (11), yielding

fl1;0 ðxn ; n; sn Þ þ el1;1 ðxn ; n; sn Þ þ   gxn ¼ fl1;0 ðf0 þ ef1 þ    ; n þ 1; sn þ eÞ þ el1;1 ðf0 þ ef1 þ    ; n þ 1; sn þ eÞ þ   gðf0 þ ef1 þ   Þ:

ð13Þ

Now it should be observed that (for i ¼ 0; 1; 2; . . .)

l1;i ðf0 þ ef1 þ e2 f2 þ    ; n þ 1; sn þ eÞ ¼ l1;i ðf0 ; n þ 1; sn Þ þ e½D1 ðl1;i Þðf0 ; n þ 1; sn Þf1 þ D3 ðl1;i Þðf0 ; n þ 1; sn Þ

" 1 1 þ e2 D1 ðl1;i Þðf0 ; n þ 1; sn Þf2 þ D1;1 ðl1;i Þðf0 ; n þ 1; sn Þf12 þ D3;3 ðl1;i Þðf0 ; n þ 1; sn Þ 2 2 # þ D1;3 ðl1;i Þðf0 ; n þ 1; sn Þf1 þ e3 ð  Þ þ   

ð14Þ

838

M. Rafei, W.T. Van Horssen / Commun Nonlinear Sci Numer Simulat 18 (2013) 835–850

where Dk ¼ @X@ is differentiation with respect to the kth argument of k

2

l1;i , and Dk;m ¼ @X k@@Xm is differentiation with respect to

the kth, and mth argument of l1;i , for k ¼ 1; 3; m ¼ 1; 3, and i ¼ 0; 1; 2; . . .; Then, by using (14) and by taking apart in (13) the O(1)-terms, the OðeÞ-terms, the Oðe2 Þ-terms, and so on, it follows that the O(1)-problem becomes

l1;0 ðxn ; n; sn Þxn ¼ l1;0 ðf0 ; n þ 1; sn Þf0 ;

ð15Þ

that the OðeÞ-problem becomes

l1;1 ðxn ; n; sn Þxn ¼ l1;1 ðf0 ; n þ 1; sn Þf0 þ l1;0 ðf0 ; n þ 1; sn Þf1 þ D1 ðl1;0 Þðf0 ; n þ 1; sn Þf0 f1 þ D3 ðl1;0 Þðf0 ; n þ 1; sn Þf0 ;

ð16Þ

and that the Oðe2 Þ-problem becomes

1 2 þ 1; sn Þf0 f12 þ l1;1 ðf0 ; n þ 1; sn Þf1 þ D1 ðl1;1 Þðf0 ; n þ 1; sn Þf0 f1 þ D3 ðl1;0 Þðf0 ; n þ 1; sn Þf1

l1;2 ðxn ; n; sn Þxn ¼ l1;2 ðf0 ; n þ 1; sn Þf0 þ l1;0 ðf0 ; n þ 1; sn Þf2 þ D1 ðl1;0 Þðf0 ; n þ 1; sn Þ½f12 þ f0 f2  þ D1;1 ðl1;0 Þðf0 ; n 1 þ D3;3 ðl1;0 Þðf0 ; n þ 1; sn Þf0 þ D3 ðl1;1 Þðf0 ; n þ 1; sn Þf0 þ D1;3 ðl1;0 Þðf0 ; n þ 1; sn Þf0 f1 : 2

ð17Þ

For a given function f ðxn ; n; eÞ the Oð1Þ-problem, the OðeÞ-problem, and so on now have to be solved. For some examples the reader is referred to [9]. For a system of two first order ODEs,

8 2 > < xnþ1 ¼ f1 ðxn ; yn ; n; eÞ ¼ f1;0 ðxn ; yn ; nÞ þ ef1;1 ðxn ; yn ; nÞ þ e f1;2 ðxn ; yn ; nÞ þ    ynþ1 ¼ f2 ðxn ; yn ; n; eÞ ¼ f2;0 ðxn ; yn ; nÞ þ ef2;1 ðxn ; yn ; nÞ þ e2 f2;2 ðxn ; yn ; nÞ þ    > : snþ1 ¼ sn þ e;

ð18Þ

the same procedure can be followed. An invariance vector lðxn ; yn ; n; sn ; eÞ for (18) has to satisfy (5). Now it will also be assumed that l1 and l2 can be expanded in power series in e, that is,

(

l1 ðxn ; yn ; n; sn ; eÞ ¼ l1;0 ðxn ; yn ; n; sn Þ þ el1;1 ðxn ; yn ; n; sn Þ þ e2 l1;2 ðxn ; yn ; n; sn Þ þ    l2 ðxn ; yn ; n; sn ; eÞ ¼ l2;0 ðxn ; yn ; n; sn Þ þ el2;1 ðxn ; yn ; n; sn Þ þ e2 l2;2 ðxn ; yn ; n; sn Þ þ   

ð19Þ

It is obvious from (5) that the functional equation for a system of three first order ODEs becomes

l1 ðxn ; yn ; n; sn ; eÞxn þ l2 ðxn ; yn ; n; sn ; eÞyn þ l3 ðxn ; yn ; n; sn ; eÞsn ¼ l1 ðxnþ1 ; ynþ1 ; n þ 1; sn þ e; eÞxnþ1 þ l2 ðxnþ1 ; ynþ1 ; n þ 1; sn þ e; eÞynþ1 þ l3 ðxnþ1 ; ynþ1 ; n þ 1; sn þ e; eÞsnþ1 :

ð20Þ

We assume that

l3 ðxn ; yn ; n; sn ; eÞ ¼

1

sn

ð21Þ

:

Therefore, functional Eq. (20) becomes

l1 ðxn ; yn ; n; sn ; eÞxn þ l2 ðxn ; yn ; n; sn ; eÞyn ¼ l1 ðxnþ1 ; ynþ1 ; n þ 1; sn þ e; eÞxnþ1 þ l2 ðxnþ1 ; ynþ1 ; n þ 1; sn þ e; eÞynþ1 : ð22Þ After substitution of (18) and of (19) into (22) we will have

li;j ðf1;0 þ ef1;1 þ e2 f1;2 þ    ; f2;0 þ ef2;1 þ e2 f2;2 þ    ; n þ 1; sn þ eÞ ¼ li;j ðf1;0 ; f2;0 ; n þ 1; sn Þ þ e½D1 ðli;j Þðf1;0 ; f2;0 ; n þ 1; sn Þf1;1 þ D2 ðli;j Þðf1;0 ; f2;0 ; n þ 1; sn Þf2;1 þ D4 ðli;j Þðf1;0 ; f2;0 ; n þ 1; sn Þ "

2

þ e D2 ðli;j Þðf1;0 ; f2;0 ; n þ 1; sn Þf2;2 þ D1 ðli;j Þðf1;0 ; f2;0 ; n þ 1; sn Þf1;2 þ D1;2 ðli;j Þðf1;0 ; f2;0 ; n þ 1; sn Þf1;1 f2;1 1 1 2 2 þ D1;1 ðli;j Þðf1;0 ; f2;0 ; n þ 1; sn Þf1;1 þ D2;2 ðli;j Þðf1;0 ; f2;0 ; n þ 1; sn Þf2;1 þ D1;4 ðli;j Þðf1;0 ; f2;0 ; n þ 1; sn Þf1;1 2 2 # 1 þ D2;4 ðli;j Þðf1;0 ; f2;0 ; n þ 1; sn Þf2;1 þ D4;4 ðli;j Þðf1;0 ; f2;0 ; n þ 1; sn Þ þ    2 where Dk ¼ @X@ is differentiation with respect to the kth argument of k

the kth, and mth argument of lem becomes

ð23Þ

2

li;j , and Dk;m ¼ @Xk@@Xm is differentiation with respect to

li;j , for k ¼ 1; 2; 4; m ¼ 1; 2; 4; i ¼ 1; 2, and j ¼ 0; 1; 2; . . .; Then, it follows that the Oð1Þ-prob-

l1;0 ðxn ; yn ; n; sn Þxn þ l2;0 ðxn ; yn ; n; sn Þyn ¼ l1;0 ðf1;0 ; f2;0 ; n þ 1; sn Þf1;0 þ l2;0 ðf1;0 ; f2;0 ; n þ 1; sn Þf2;0

ð24Þ

M. Rafei, W.T. Van Horssen / Commun Nonlinear Sci Numer Simulat 18 (2013) 835–850

839

and the OðeÞ-problem becomes

l1;1 ðxn ; yn ; n; sn Þxn þ l2;1 ðxn ; yn ; n; sn Þyn ¼ l1;0 ðf1;0 ; f2;0 ; n þ 1; sn Þf1;1 þ D1 ðl1;0 Þðf1;0 ; f2;0 ; n þ 1; sn Þf1;0 f1;1 þ D2 ðl1;0 Þ  ðf1;0 ; f2;0 ; n þ 1; sn Þf1;0 f2;1 þ l1;1 ðf1;0 ; f2;0 ; n þ 1; sn Þf1;0 þ l2;0 ðf1;0 ; f2;0 ; n þ 1; sn Þf2;1 þ D1 ðl2;0 Þðf1;0 ; f2;0 ; n þ 1; sn Þf1;1 f2;0 þ D2 ðl2;0 Þðf1;0 ; f2;0 ; n þ 1; sn Þf2;0 f2;1 þ l2;1 ðf1;0 ; f2;0 ; n þ 1; sn Þf2;0 þ D4 ðl1;0 Þðf1;0 ; f2;0 ; n þ 1; sn Þf1;0 þ D4 ðl2;0 Þðf1;0 ; f2;0 ; n þ 1; sn Þf2;0 : ð25Þ Now we define new functions

(

Z 0 ðxn ; yn ; n; sn Þ ¼ l1;0 ðxn ; yn ; n; sn Þxn þ l2;0 ðxn ; yn ; n; sn Þyn ; Z 1 ðxn ; yn ; n; sn Þ ¼ l1;1 ðxn ; yn ; n; sn Þxn þ l2;1 ðxn ; yn ; n; sn Þyn

ð26Þ

and by using (26), the Oð1Þ-problem becomes

Z 0 ðxn ; yn ; n; sn Þ ¼ Z 0 ðf1;0 ; f2;0 ; n þ 1; sn Þ;

ð27Þ

and the OðeÞ-problem becomes

Z 1 ðxn ; yn ; n; sn Þ ¼ Z 1 ðf1;0 ; f2;0 ; n þ 1; sn Þ þ D1 ðZ 0 Þðf1;0 ; f2;0 ; n þ 1; sn Þf1;1 þ D2 ðZ 0 Þðf1;0 ; f2;0 ; n þ 1; sn Þf2;1 þ D4 ðZ 0 Þ  ðf1;0 ; f2;0 ; n þ 1; sn Þ:

ð28Þ

In the next two sections it will be shown how the perturbation method based on invariance factors and multiple scales can be applied to a system of two first order ODEs. In Section 4 a second order, weakly nonlinear, regularly perturbed ODEs with a Van der Pol type of nonlinearity will be considered, and in Section 5 a special case of nonlinear Rayleigh oscillator will be studied. 4. A weakly nonlinear, regularly perturbed system of two ODEs In this section approximations of first integrals for a second order, weakly nonlinear, regularly perturbed ODE with a Van der Pol type of nonlinearity will be considered. The Van der Pol equation [1] corresponds to a nonlinear oscillatory system that has both input and output sources of energy. This equation is given by

_ €x þ x20 x ¼ lð1  x2 Þx;

ð29Þ

where x ¼ xðtÞ; l is a non-negative small parameter, and where x0 is a bounded constant. This leads to the following difference equation, when a central difference scheme [20] is used to discretize Eq. (29):

xnþ1  2xn þ xn1 h

2

xnþ1  xn1 ; 2h

ð30Þ

ð1  x2n Þðxnþ1  xn1 Þ:

ð31Þ

þ x20 xn ¼ lð1  x2n Þ

or equivalently, 2

xnþ1 þ ðx20 h  2Þxn þ xn1 ¼

lh 2

In fact (31) can be considered as a central finite difference approximation of (29). In this case xn is an approximation of xðtn Þ at tn ¼ nh, where h is the discretization time step. Now let

(

2

h x20  2 ¼ 2 cosðhÞ; lh 2

¼ e:

ð32Þ

In this paper ordinary differential equations are discretized to obtain difference equations, and since the main goal of this paper is to develop perturbation methods for difference equations, it will be assumed that e tends to zero (that is, l tends to zero), and that h is small but fixed. Substitution of these new constants h and e from (32) into (31) and shifting the index by 1, yields

xnþ2  2 cosðhÞxnþ1 þ xn ¼ eð1  x2nþ1 Þðxnþ2  xn Þ;

ð33Þ

where e is a small parameter, that is, 0 < e  1, and where h is constant (which is related to the stepsize in making the continuous Van der Pol equation discrete). Now (33) is transformed into a system of first order difference equations,

(

xnþ1 ¼ yn ; 1 2 ynþ1 ¼ 1eð1y 2 Þ ð2yn cosðh0 Þ  xn  exn ð1  yn ÞÞ:

ð34Þ

n

Two (functionally independent) approximations of first integrals for system (34) have been constructed in [10]. For one of these approximations, the straightforward perturbation approach yields secular terms. For that reason, we now introduce

840

M. Rafei, W.T. Van Horssen / Commun Nonlinear Sci Numer Simulat 18 (2013) 835–850

multiple scales approach to avoid these secular terms. By introducing a new (slow) iteration-scale sn ¼ en, we add (as mentioned earlier) a new difference equation snþ1 ¼ sn þ e to the original system of ODEs. It is also more convenient to work in polar coordinates



xn ¼ r n sinðun þ hÞ;

ð35Þ

yn ¼ r n sinðun Þ: In polar coordinates the Van der Pol system (34) based on multiple scales becomes

8   1 2 > < r nþ1 ¼ r n 1 þ 4 ½r nðcosð4uÞ  1Þ þ 4ðcosð2uÞ þ 1Þe þ    ;  unþ1 ¼ un  h þ  14 r2n sinð4un Þ þ 12 r2n sinð2un Þ  sinð2un Þ e þ    ; > : snþ1 ¼ sn þ e:

ð36Þ

Then we substitute (36) into (18) and (27). So, the Oð1Þ-problem will be

Z 0 ðr n ; un ; n; sn Þ ¼ Z 0 ðrn ; un  h; n þ 1; sn Þ:

ð37Þ

The general solution of (37) is given by

Z 0 ðr n ; un ; n; sn Þ ¼ A0 ðr n ; un þ nh; sn Þ;

ð38Þ

where A0 is an arbitrary function. Since we are dealing with a second order difference Eq. (33), we need two functionally independent approximations of the invariants. The function A0 is still arbitrary, and will now be chosen to be as simple as possible to obtain relatively simple approximations of the invariants. 4.1. Case 1 First we take

Z 0 ðr n ; un ; n; sn Þ ¼ B0 ðrn ; sn Þ;

ð39Þ

where B0 is an arbitrary function. By substituting (39) into (28), the OðeÞ-problem then becomes

1 Z 1 ðr n ; un ; n; sn Þ ¼ Z 1 ðrn ; un  h; n þ 1; sn Þ þ D1 ðB0 Þðr n ; sn Þrn cosð2un Þ þ D1 ðB0 Þðr n ; sn Þr3n cosð4un Þ þ D1 ðB0 Þ 4 1 3  ðr n ; sn Þr n  D1 ðB0 Þðr n ; sn Þrn þ D2 ðB0 Þðr n ; sn Þ: 4

ð40Þ

It is quite obvious that Z 1 should be bounded in un in order to avoid secular terms. From (40) it then follows that B0 should satisfy

1 D1 ðB0 Þðr n ; sn Þr n  D1 ðB0 Þðrn ; sn Þr 3n þ D2 ðB0 Þðr n ; sn Þ ¼ 0: 4

ð41Þ

The general solution of the first order partial differential equation for B0 can readily be obtained, yielding

B0 ðr n ; sn Þ ¼ C 0

   4 1  2 expð2sn Þ ; rn

ð42Þ

where C 0 is an arbitrary function. The OðeÞ-problem (40) can now be solved completely. Since we are interested in secularfree, simple solutions we now take

B0 ðr n ; sn Þ ¼

 1

 4 expð2sn Þ; 2 rn

ð43Þ

we then have

Z 1 ðr n ; un ; n; sn Þ ¼ Z 1 ðrn ; un  h; n þ 1; sn Þ þ

2 expð2sn Þ ð4 cosð2un Þ þ r 2n cosð4un ÞÞ: r2n

ð44Þ

Now, we are going to construct a particular solution of (44). To do this, we look for a particular solution in the form

Z 1;p ðrn ; un ; n; sn Þ ¼

2 expð2sn Þ ða1 ðnÞ sinð2un Þ þ a2 ðnÞ cosð2un Þ þ a3 ðnÞr2n sinð4un Þ þ a4 ðnÞr 2n cosð4un ÞÞ: r 2n

ð45Þ

By substituting (45) into (44) and by taking apart in the resulting equation the coefficients of sinð2un Þ, cosð2un Þ; r 2n sinð4un Þ, and r 2n cosð4un Þ, and setting these coefficients equal to zero, we will have the following system

8 a1 ðn þ 1Þ ¼ cosð2hÞaðnÞ  sinð2hÞbðnÞ þ 4 sinð2hÞ; > > > < a2 ðn þ 1Þ ¼ sinð2hÞaðnÞ þ cosð2hÞbðnÞ  4 cosð2hÞ; > > a3 ðn þ 1Þ ¼ cosð4hÞcðnÞ  sinð4hÞdðnÞ þ sinð4hÞ; > : a4 ðn þ 1Þ ¼ sinð4hÞcðnÞ þ cosð4hÞdðnÞ  cosð4hÞ:

ð46Þ

M. Rafei, W.T. Van Horssen / Commun Nonlinear Sci Numer Simulat 18 (2013) 835–850

841

The eigenvalues of the homogeneous system related to (46) are cosð2hÞ  i sinð2hÞ, and cosð4hÞ  i sinð4hÞ, and a particular solution of (46) is given by

8 2 a1 ðnÞ ¼ 1cosð2hÞ ðsinð2hÞ þ sinð2nhÞ  sinð2nh þ 2hÞÞ; > > > > > < a2 ðnÞ ¼ 2 ðsinð2hÞ  sinð2nhÞ  sinð2nh þ 2hÞÞ; sinð2hÞ 1 > a3 ðnÞ ¼ 2ðcosð3hÞcosðhÞÞ ðsinðhÞ  sinð3hÞ  sinð4nh þ hÞ þ sinð4nh þ 3hÞÞ; > > > > : a ðnÞ ¼ 1 ðsinð4hÞ  sinð4nhÞ  sinð4nh þ 4hÞÞ: 4

ð47Þ

2 sinð4hÞ

The general solution Z 1;h of the equation corresponding to the homogeneous equation (44) is

Z 1;h ðrn ; un ; n; sn Þ ¼ A1 ðr n ; un þ nh; sn Þ;

ð48Þ

where A1 is an arbitrary function, which can be used to avoid secular terms at the O (e2 ) level. At this moment, however, we are not interested in the higher order approximations. For that reason, we take A1  0. An approximation IA ðxn ; yn ; nÞ of an invariant Iðxn ; yn ; nÞ ¼ constant for (36) is now constructed, and is given by

  4 2e expð2sn Þ IA ðr n ; un ; n; sn Þ ¼ 1  2 expð2sn Þ þ ða1 ðnÞ sinð2un Þ þ a2 ðnÞ cosð2un Þ þ a3 ðnÞr 2n sinð4un Þ rn r 2n þ a4 ðnÞr 2n cosð4un ÞÞ;

ð49Þ

or equivalently,

IA ðr nþ1 ; unþ1 ; n þ 1; snþ1 Þ ¼ IA ðr n ; un ; n; sn Þ þ e2 Rðr n ; un ; n; sn ; eÞ;

ð50Þ

where 2

e

! !!   8f 1;1 4 4 4 expð2sn Þ Rðrn ; un ; n; sn ; eÞ ¼ 1  2 expð2ðsn þ eÞÞ  1  2 expð2sn Þ  e þ2 1 2 3 f1;0 f1 f1;0 f1;0 " expð2ðsn þ eÞÞ ða1 ðn þ 1Þ sinð2f 2 Þ þ a2 ðn þ 1Þ cosð2f 2 Þ þ a3 ðn þ 1Þf12 sinð4f 2 Þ þ 2e f12 þ a4 ðn þ 1Þf12 cosð4f 2 ÞÞ 

expð2ðsn ÞÞ ða1 ðn þ 1Þ sinð2f 2;0 Þ þ a2 ðn þ 1Þ cosð2f 2;0 Þ 2 f1;0 #

2 2 sinð4f 2;0 Þ þ a4 ðn þ 1Þf1;0 cosð4f 2;0 ÞÞ ; þa3 ðn þ 1Þf1;0

ð51Þ

where f1 ; f2 ; f1;0 ; f2;0 and f1;1 are given by (18) and (36). Since the Van der Pol system (36) is autonomous, and a1 ðnÞ; a2 ðnÞ; a3 ðnÞ, and a4 ðnÞ are bounded, then R is bounded. From (50) it follows that

IA ðr n ; un ; n; sn Þ ¼ IA ðr 0 ; u0 ; 0; 0Þ þ e2

n1 X Rðri ; ui ; i; si ; eÞ:

ð52Þ

i¼0

From (50) and (52), it can be shown that

IA ðr n ; un ; n; sn Þ ¼ constant þ Oðe2 nÞ;

ð53Þ

it then follows that

(

IA ðr n ; un ; n; sn Þ ¼ constant þ Oðe2 Þ for n ¼ Oð1Þ; IA ðr n ; un ; n; sn Þ ¼ constant þ OðeÞ for n ¼ Oð1e Þ:

ð54Þ

From (49) and (54) it then follows that

  4 1  2 expð2sn Þ ¼ constant þ OðeÞ; rn

ð55Þ

on an iteration-scale of order e1 . So far only one approximation of a first integral has been determined. Another (functionally independent) approximation of a first integral can also be obtained in a similar and straightforward way as follows from the next subsection. 4.2. Case 2 In this case we take in (38)

Z 0 ðr n ; un ; n; sn Þ ¼ un þ nh:

ð56Þ

842

M. Rafei, W.T. Van Horssen / Commun Nonlinear Sci Numer Simulat 18 (2013) 835–850

It turns out in [10] that no secular terms occur in the approximation of a first integral that coincides with un þ nh ¼ constant when e ¼ 0. From (28), the OðeÞ-problem then becomes

1 1 Z 1 ðr n ; un ; n; sn Þ ¼ Z 1 ðrn ; un  h; n þ 1; sn Þ  sinð2un Þ þ r2n sinð2un Þ  r 2n sinð4un Þ: 2 4

ð57Þ

Now, we are going to construct a particular solution of (57) in the form

Z 1;p ðrn ; un ; n; sn Þ ¼ b1 ðnÞ sinð2un Þ þ b2 ðnÞ cosð2un Þ þ b3 ðnÞr 2n sinð2un Þ þ b4 ðnÞr2n cosð2un Þ þ b5 ðnÞr 2n sinð4un Þ þ b6 ðnÞr2n cosð4un Þ:

ð58Þ

By substituting (58) into (57), and by taking apart in the resulting equation the coefficients of sinð2un Þ; cosð2un Þ; r2n sin ð2un Þ; r2n cosð2un Þ; r 2n sinð4un Þ, and r2n cosð4un Þ, and setting these coefficients equal to zero, we will have the following system

8 b1 ðn þ 1Þ ¼ cosð2hÞaðnÞ  sinð2hÞbðnÞ þ cosð2hÞ; > > > > > b2 ðn þ 1Þ ¼ sinð2hÞaðnÞ þ cosð2hÞbðnÞ þ sinð2hÞ; > > > < b ðn þ 1Þ ¼ cosð2hÞcðnÞ  sinð2hÞdðnÞ  1 cosð2hÞ; 3 2 > b4 ðn þ 1Þ ¼ sinð2hÞcðnÞ þ cosð2hÞdðnÞ  12 sinð2hÞ; > > > > 1 > > b5 ðn þ 1Þ ¼ cosð4hÞeðnÞ  sinð4hÞf ðnÞ þ 4 cosð4hÞ; > : b6 ðn þ 1Þ ¼ sinð4hÞeðnÞ þ cosð4hÞf ðnÞ þ 14 sinð4hÞ:

ð59Þ

The eigenvalues of the homogeneous system related to (59) are cosð2hÞ  i sinð2hÞ; cosð2hÞ  i sinð2hÞ, and cosð4hÞ  i sinð4hÞ, and a particular solution of (59) is given by

8 1 b1 ðnÞ ¼ 2 sinð2hÞ ðsinð2nh þ 2hÞ þ sinð2nhÞ  sinð2hÞÞ; > > > > > > 1 > ðsinð2nh þ 2hÞ  sinð2nhÞ  sinð2hÞÞ; > b2 ðnÞ ¼ 2ðcosð2hÞ1Þ > > > > > 1 > < b3 ðnÞ ¼ 4 sinð2hÞ ðsinð2hÞ  sinð2nh þ 2hÞ  sinð2nhÞÞ; 1 > > ðsinð2hÞ  sinð2nh þ 2hÞ þ sinð2nhÞÞ; b4 ðnÞ ¼ 4ðcosð2hÞ1Þ > > > > > > 1 > > > b5 ðnÞ ¼ 8 sinð4hÞ ðsinð4nh þ 4hÞ þ sinð4nhÞ  sinð4hÞÞ; > > > : b ðnÞ ¼ 1 ðsinð4nh þ 3hÞ  sinð4nh þ hÞ  sinð3hÞ þ sinðhÞÞ: 6

ð60Þ

8ðcosð3hÞcosðhÞÞ

The general solution Z 1;h of the equation corresponding to the homogeneous Eq. (57) is

Z 1;h ðrn ; un ; n; sn Þ ¼ A1 ðr n ; un þ nh; sn Þ;

ð61Þ

where A1 is an arbitrary function, which can be used to avoid secular terms at the O (e2 ) level. At this moment, however, we are not interested in the higher order approximations. For that reason, we take A1  0. An approximation IA ðr n ; un ; n; sn Þ of an invariant Iðr n ; un ; n; sn Þ ¼ constant for (36) is now constructed, and is given by

IA ðr n ; un ; n; sn Þ ¼ un þ nh þ eðb1 ðnÞ sinð2un Þ þ b2 ðnÞ cosð2un Þ þ b3 ðnÞr 2n sinð2un Þ þ b4 ðnÞr 2n cosð2un Þ þ b5 ðnÞr 2n  sinð4un Þ þ b6 ðnÞr 2n cosð4un ÞÞ;

ð62Þ

where b1 ðnÞ; b2 ðnÞ; . . ., and b6 ðnÞ are given by (60), and

IA ðr nþ1 ; unþ1 ; n þ 1; snþ1 Þ ¼ IA ðrn ; un ; n; sn Þ þ e2 Rðr n ; un ; n; sn ; eÞ;

ð63Þ

where

e2 Rðrn ; un ; n; sn ; eÞ ¼ f2  f2;0  ef2;1 þ eðb1 ðn þ 1Þ sinð2f 2 Þ þ b2 ðn þ 1Þ cosð2f 2 Þ þ b3 ðn þ 1Þf12 sinð2f 2 Þ þ b4 ðn þ 1Þf12 cosð2f 2 Þ þ b5 ðn þ 1Þf12 sinð4f 2 Þ þ b6 ðn þ 1Þf12 cosð4f 2 Þ  ðb1 ðn þ 1Þ sinð2f 2;0 Þ 2 2 þ b2 ðn þ 1Þ cosð2f 2;0 Þ þ b3 ðn þ 1Þf1;0 sinð2f 2;0 Þ þ b4 ðn þ 1Þf1;0 cosð2f 2;0 Þ 2 2 sinð4f 2;0 Þ þ b6 ðn þ 1Þf1;0 cosð4f 2;0 ÞÞÞ; þ b5 ðn þ 1Þf1;0

ð64Þ

where f1 ; f2 ; f1;0 ; f2;0 and f2;1 are given by (18) and (36). Since the Van der Pol system (36) is autonomous, and b1 ðnÞ; b2 ðnÞ; b3 ðnÞ; b4 ðnÞ; b5 ðnÞ, and b6 ðnÞ are bounded, then R is bounded. From (63) it follows that

IA ðr n ; un ; n; sn Þ ¼ IA ðr 0 ; u0 ; 0; 0Þ þ e2

n1 X Rðr i ; ui ; i; si ; eÞ:

ð65Þ

i¼0

From (63) and (65), it can be shown that

IA ðr n ; un ; n; sn Þ ¼ constant þ Oðe2 nÞ;

ð66Þ

M. Rafei, W.T. Van Horssen / Commun Nonlinear Sci Numer Simulat 18 (2013) 835–850

843

it then follows that

(

IA ðr n ; un ; n; sn Þ ¼ constant þ Oðe2 Þ for n ¼ Oð1Þ; IA ðr n ; un ; n; sn Þ ¼ constant þ OðeÞ for n ¼ Oð1e Þ:

ð67Þ

From (62) and (67) it then follows that

un þ nh ¼ constant þ OðeÞ; on an iteration-scale of order

ð68Þ

e . 1

5. A nonlinear Rayleigh oscillator In this section approximations of first integrals for a second order, weakly nonlinear, ODE with a Rayleigh type of nonlinearity [21] will be considered. The mathematical model that describes the flow-induced vibrations of the oscillator with one degree of freedom in a uniform wind-field can be given by

  x_ 2 _ €x þ x þ bx2 ¼ c 1  x; 3

ð69Þ

where x ¼ xðtÞ, and b and c are constants. In [21] Waluya approximated first integral and periodic solutions of a Rayleigh equation when b ¼ OðeÞ and c ¼ Oðe2 Þ, and when b ¼ Oð1Þ and c ¼ OðeÞ. In this section we consider

  x_ 2 _ €x þ x þ eax2 ¼ e 1  x; 3

ð70Þ

where a is a non-negative constant. This leads to the following difference equation, when a central difference scheme [20] is used to discretize Eq. (70):

xnþ1  2xn þ xn1 h

2

  1 xn  xn1 2 xn  xn1 þ xn þ eax2n ¼ e 1  ; 3 h h

ð71Þ

or equivalently,

  1 xn  xn1 2 2 2 xnþ1 þ ðh  2Þxn þ xn1 þ eah x2n ¼ eh 1  ðxn  xn1 Þ: 3 h

ð72Þ

In fact (72) can be considered as a central finite difference approximation of (70). In this case xn is an approximation of xðtn Þ at tn ¼ nh, where h is the discretization time step. Now let 2

h  2 ¼ 2 cosðhÞ:

ð73Þ

Substitution of this new constant h from (73) into (72) and shifting the index by 1, yields

  1 xnþ1  xn 2 2 ðxnþ1  xn Þ; xnþ2  2 cosðhÞxnþ1 þ xn þ eah x2nþ1 ¼ eh 1  3 h

ð74Þ

where e is a small parameter, that is, 0 < e  1, and where h is constant (which is related to the stepsize in making the continuous Rayleigh equation discrete). Now (74) is transformed into a system of first order difference equations,

(

xnþ1 ¼ yn ; 2

ynþ1 ¼ 2yn cosðhÞ  xn  eah y2n þ eh

2



1  13

yn xn 2  yn xn h

h

:

ð75Þ

In the further analysis it will be assumed that e is the perturbation parameter which tends to zero, and the discretization parameter h is small, but fixed (and independent of e). For the Rayleigh system (75), the straightforward perturbation approach yields secular terms. For that reason, we now introduce multiple scales approach to avoid these secular terms. By introducing a new (slow) iteration-scale sn ¼ en, we add a new difference equation snþ1 ¼ sn þ e to the original system of ODEs. It is also more convenient to work in polar coordinates



xn ¼ rn sinðun þ hÞ; yn ¼ r n sinðun Þ:

ð76Þ

In polar coordinates the Rayleigh system (70) based on multiple scales becomes

8 2 3 er n 2 > > < r nþ1 ¼ rn  24hðcosð2hÞ1Þ ð12h G1 þ 6ar n h G2 þ rn G3 Þ þ    2 3 e unþ1 ¼ un  h  24hð34 cosð2hÞþcosð4hÞÞ ð12h H1 þ 6ar n h H2 þ r 2n H3 Þ þ    > > : snþ1 ¼ sn þ e;

ð77Þ

844

M. Rafei, W.T. Van Horssen / Commun Nonlinear Sci Numer Simulat 18 (2013) 835–850

where

G1 ¼  cosð2un  hÞ þ cosð2un Þ þ cosð2un þ hÞ  cosð2un þ 2hÞ  cosð2hÞ þ 1; G2 ¼  sinðun  hÞ þ sinðun þ hÞ þ sinð3un  hÞ  sinð3un þ hÞ; G3 ¼ 3 cosð2un  2hÞ þ 8 cosð2un  hÞ  3 cosð2un Þ  8 cosð2un þ hÞ þ 7 cosð2un þ 2hÞ  cosð2un þ 4hÞ  cosð4un  hÞ þ 3 cosð4un Þ  2 cosð4un þ hÞ  2 cosð4un þ 2hÞ þ 3 cosð4un þ 3hÞ  cosð4un þ 4hÞ þ 3 cosðhÞ þ 6 cosð2hÞ  3 cosð3hÞ  6; H1 ¼ sinð2un  3hÞ  sinð2un  2hÞ  3 sinð2un  hÞ þ 3 sinð2un Þ þ 3 sinð2un þ hÞ  3 sinð2un þ 2hÞ  sinð2un þ 3hÞ þ sinð2un þ 4hÞ  6 sinðhÞ þ 2 sinð2hÞ þ 2 sinð3hÞ  sinð4hÞ;

ð78Þ

H2 ¼ 3 cosðun  3hÞ þ 9 cosðun  hÞ  9 cosðun þ hÞ þ 3 cosðun þ 3hÞ þ cosð3un  3hÞ  3 cosð3un  hÞ þ 3 cosð3un þ hÞ  cosð3un þ 3hÞ; H3 ¼ 3 sinð2un  4hÞ  10 sinð2un  3hÞ þ 3 sinð2un  2hÞ þ 24 sinð2un  hÞ  26 sinð2un Þ  12 sinð2un þ hÞ þ 30 sinð2un þ 2hÞ  8 sinð2un þ 3hÞ  9 sinð2un þ 4hÞ þ 6 sinð2un þ 5hÞ  sinð2un þ 6hÞ þ sinð4un  3hÞ  3 sinð4un  2hÞ þ 8 sinð4un Þ  6 sinð4un þ hÞ  6 sinð4un þ 2hÞ þ 8 sinð4un þ 3hÞ  3 sinð4un þ 5hÞ þ sinð4un þ 6hÞ þ 42 sinðhÞ  24 sinð2hÞ  9 sinð3hÞ þ 12 sinð4hÞ  3 sinð5hÞ: Then we substitute (77) into (18) and (27). So, the Oð1Þ-problem will be

Z 0 ðr n ; un ; n; sn Þ ¼ Z 0 ðrn ; un  h; n þ 1; sn Þ:

ð79Þ

The general solution of (79) is given by

Z 0 ðr n ; un ; n; sn Þ ¼ K 0 ðrn ; kn ; sn Þ;

ð80Þ

where kn ¼ un þ nh is just defined for simplicity, and K 0 is an arbitrary function. By substituting (80) into (28), the OðeÞ-problem then becomes

Z 1 ðr n ; un ; n; sn Þ ¼ Z 1 ðrn ; un  h; n þ 1; sn Þ þ

1 2 f3r n D1 ðK 0 Þðrn ; kn ; sn Þðr 2n L1 þ 4h L2 Þ 24hð3  4 cosð2hÞ þ cosð4hÞÞ 2

þ 3D2 ðK 0 Þðrn ; kn ; sn Þðr 2n L3 þ 4h L4 Þ þ 24hD3 ðK 0 Þðr n ; kn ; sn ÞL5 þ r n D1 ðK 0 Þðr n ; kn ; sn Þðr2n L6 3

2

3

2

þ 6ar n h L7 þ 12h L8 Þ þ D2 ðK 0 Þðrn ; kn ; sn Þðr 2n L9 þ 6ar n h L10 þ 12h L11 Þg;

ð81Þ

where

L1 ¼ 2 cosðhÞ þ 8 cosð2hÞ  3 cosð3hÞ  2 cosð4hÞ þ cosð5hÞ  6; L2 ¼ 4 cosð2hÞ þ cosð4hÞ þ 3; L3 ¼ 14 sinðhÞ þ 8 sinð2hÞ þ 3 sinð3hÞ  4 sinð4hÞ þ sinð5hÞ; L4 ¼ 6 sinðhÞ  2 sinð2hÞ  2 sinð3hÞ þ sinð4hÞ; L5 ¼ 4 cosð2hÞ þ cosð4hÞ þ 3; L6 ¼ 3 cosð2un  4hÞ  8 cosð2un  3hÞ  3 cosð2un  2hÞ þ 24 cosð2un  hÞ  10 cosð2un Þ  24 cosð2un þ hÞ þ 18 cosð2un þ 2hÞ þ 8 cosð2un þ 3hÞ þ cosð2un þ 6hÞ þ cosð4un  3hÞ  3 cosð4un  2hÞ þ 8 cosð4un Þ  6 cosð4un þ hÞ  6 cosð4un þ 2hÞ þ 8 cosð4un þ 3hÞ  9 cosð2un þ 4hÞ  3 cosð4un þ 5hÞ þ cosð4un þ 6hÞ; L7 ¼ sinðun  3hÞ þ 3 sinðun  hÞ þ 3 sinðun þ hÞ  sinðun þ 3hÞ  sinð3un  3hÞ þ 3 sinð3un  hÞ  3 sinð3un þ hÞ þ sinð3un þ 3hÞ; L8 ¼ cosð2un  3hÞ  cosð2un  2hÞ  3 cosð2un  hÞ þ 3 cosð2un Þ þ 3 cosð2un þ hÞ  3 cosð2un þ 2hÞ  cosð2un þ 3hÞ þ cosð2un þ 4hÞ; L9 ¼ 3 sinð2un  4hÞ þ 10 sinð2un  3hÞ  3 sinð2un  2hÞ  24 sinð2un  hÞ þ 26 sinð2un Þ þ 12 sinð2un þ hÞ  30 sinð2un þ 2hÞ þ 8 sinð2un þ 3hÞ þ 9 sinð2un þ 4hÞ  6 sinð2un þ 5hÞ þ sinð2un þ 6hÞ  sinð4un  3hÞ þ 3 sinð4un  2hÞ  8 sinð4un Þ þ 6 sinð4un þ hÞ þ 6 sinð4un þ 2hÞ  8 sinð4un þ 3hÞ þ 3 sinð4un þ 5hÞ  sinð4un þ 6hÞ;

ð82Þ

M. Rafei, W.T. Van Horssen / Commun Nonlinear Sci Numer Simulat 18 (2013) 835–850

845

L10 ¼ 3 cosðun  3hÞ  9 cosðun  hÞ þ 9 cosðun þ hÞ  3 cosðun þ 3hÞ  cosð3un  3hÞ þ 3 cosð3un  hÞ  3 cosð3un þ hÞ þ cosð3un þ 3hÞ; L11 ¼  sinð2un  3hÞ þ sinð2un  2hÞ þ 3 sinð2un  hÞ  3 sinð2un Þ  3 sinð2un þ hÞ þ 3 sinð2un þ 2hÞ þ sinð2un þ 3hÞ  sinð2un þ 4hÞ: It is quite obvious that Z 1 should be bounded in un in order to avoid secular terms. From (81) it then follows that K 0 should satisfy 2

2

3r n D1 ðK 0 Þðr n ; kn ; sn Þðr2n L1 þ 4h L2 Þ þ 3D2 ðK 0 Þðr n ; kn ; sn Þðr2n L3 þ 4h L4 Þ þ 24hD3 ðK 0 Þðr n ; kn ; sn ÞL5 ¼ 0:

ð83Þ

The general solution of the first order partial differential equation for K 0 can readily be obtained, yielding

K 0 ðr n ; kn ; sn Þ ¼ M0



 1 2 2 2 cosðhÞ1 ðr cosðhÞ  r þ 2h Þ expðh s Þ; r expðk sinðhÞÞ ; n n n n n r 2n

ð84Þ

where M 0 is an arbitrary function. Since we are dealing with a second order difference Eq. (74), we need two functionally independent approximations of the invariants. The function M 0 is still arbitrary, and will now be chosen to be as simple as possible to obtain relatively simple approximations of the invariants. 5.1. Case 1 Since we are interested in simple approximations of first integrals we now take (see also (80) and (84))

Z 0 ðr n ; un ; n; sn Þ ¼

1 2 2 ðr cosðhÞ  r 2n þ 2h Þ expðhsn Þ: r2n n

ð85Þ

The OðeÞ-problem (81) will then be

Z 1 ðr n ; un ; n; sn Þ ¼ Z 1 ðr n ; un  h; n þ 1; sn Þ þ

h expðhsn Þ 2 3 f12h N 1 þ 6ar n h N2 þ r 2n N3 g; 6r 2n ðcosð2hÞ  1Þ

ð86Þ

where

N1 ¼  cosð2un  hÞ þ cosð2un Þ þ cosð2un þ hÞ  cosð2un þ 2hÞ; N2 ¼  sinðun  hÞ þ sinðun þ hÞ þ sinð3un  hÞ  sinð3un þ hÞ; N3 ¼ 3 cosð2un  2hÞ þ 8 cosð2un  hÞ  3 cosð2un Þ  8 cosð2un þ hÞ þ 7 cosð2un þ 2hÞ  cosð2un þ 4hÞ  cosð4un  hÞ þ 3 cosð4un Þ  2 cosð4un þ hÞ  2 cosð4un þ 2hÞ þ 3 cosð4un þ 3hÞ  cosð4un þ 4hÞ:

ð87Þ

Now, we are going to construct a particular solution of (86) in the form

Z 1;p ðrn ; un ; n; sn Þ ¼

h expðhsn Þ 2 3 f12h ðc1 ðnÞ sinð2un Þ þ c2 ðnÞ cosð2un ÞÞ þ 6ar n h ðc3 ðnÞ sinðun Þ þ c4 ðnÞ cosðun Þ 6r 2n ðcosð2hÞ  1Þ þ c5 ðnÞ sinð3un Þ þ c6 ðnÞ cosð3un ÞÞ þ r 2n ðc7 ðnÞ sinð2un Þ þ c8 ðnÞ cosð2un Þ þ c9 ðnÞ sinð4un Þ þ c10 ðnÞ cosð4un ÞÞg:

ð88Þ

By substituting (88) into (86), and by taking apart in the resulting equation the coefficients of sinð2un Þ; cosð2un Þ; r n sinðun Þ; rn cosðun Þ; r n sinð3un Þ; r n cosð3un Þ; r2n sinð2un Þ, r 2n cosð2un Þ; r 2n sinð4un Þ; r 2n cosð4un Þ, and setting these coefficients equal to zero, we will have the following system

8 c1 ðn þ 1Þ ¼ cosð2hÞa1 ðnÞ  sinð2hÞa2 ðnÞ  sinðhÞ þ sinð2hÞ þ sinð3hÞ  sinð4hÞ; > > > > > c2 ðn þ 1Þ ¼ sinð2hÞa1 ðnÞ þ cosð2hÞa2 ðnÞ þ cosðhÞ  cosð2hÞ  cosð3hÞ þ cosð4hÞ; > > > > > c3 ðn þ 1Þ ¼ cosðhÞa3 ðnÞ  sinðhÞa4 ðnÞ  cosð2hÞ þ 1; > > > > > c > 4 ðn þ 1Þ ¼ sinðhÞa3 ðnÞ þ cosðhÞa4 ðnÞ  sinð2hÞ; > > > > c > 5 ðn þ 1Þ ¼ cosð3hÞa5 ðnÞ  sinð3hÞa6 ðnÞ  cosð2hÞ þ cosð4hÞ; > > > > c > 6 ðn þ 1Þ ¼ sinð3hÞa5 ðnÞ þ cosð3hÞa6 ðnÞ  sinð2hÞ þ sinð4hÞ; > > < c7 ðn þ 1Þ ¼ cosð2hÞa7 ðnÞ  sinð2hÞa8 ðnÞ þ 8 sinðhÞ  3 sinð2hÞ  8 sinð3hÞ > þ7 sinð4hÞ  sinð6hÞ; > > > > > ðn þ 1Þ ¼ sinð2hÞa c 7 ðnÞ þ cosð2hÞa8 ðnÞ  8 cosðhÞ þ 3 cosð2hÞ þ 8 cosð3hÞ > 8 > > > > 7 cosð4hÞ þ cosð6hÞ þ 3; > > > > > ðn þ 1Þ ¼ cosð4hÞa ðnÞ  sinð4hÞa10 ðnÞ  sinð3hÞ þ 3 sinð4hÞ  2 sinð5hÞ c 9 > 9 > > > > 2 sinð6hÞ þ 3 sinð7hÞ  sinð8hÞ; > > > > > ðn þ 1Þ ¼ sinð4hÞa ðnÞ þ cosð4hÞa10 ðnÞ þ cosð3hÞ  3 cosð4hÞ þ 2 cosð5hÞ c 10 9 > > : þ2 cosð6hÞ  3 cosð7hÞ þ cosð8hÞ:

ð89Þ

846

M. Rafei, W.T. Van Horssen / Commun Nonlinear Sci Numer Simulat 18 (2013) 835–850

The eigenvalues of the homogeneous system related to (89) are cosðhÞ  i sinðhÞ; cosð2hÞ  i sinð2hÞ, cosð2hÞ i sinð2hÞ; cosð3hÞ  i sinð3hÞ, and cosð4hÞ  i sinð4hÞ, and a particular solution of (89) is given by

8 1 c1 ðnÞ ¼ 2 sinð2hÞ ðcosð2nh  hÞ  cosð2nhÞ  cosð2nh þ 3hÞ þ cosð2nh þ 4hÞ  cosðhÞ þ cosð3hÞ  cosð4hÞ þ 1Þ; > > > > c ðnÞ ¼ 1 ðsinð2nh  hÞ  sinð2nhÞ  sinð2nh þ 3hÞ þ sinð2nh þ 4hÞ þ sinðhÞ þ sinð3hÞ  sinð4hÞÞ; > 2 > 2 sinð2hÞ > > > 1 > c3 ðnÞ ¼ 2 sinðhÞ ðsinðnh  hÞ þ sinðnhÞ  sinðnh þ hÞ  sinðnh þ 2hÞ þ 2 sinðhÞ þ sinð2hÞÞ; > > > > 1 > > ð cosðnh  hÞ  cosðnhÞ þ cosðnh þ hÞ þ cosðnh þ 2hÞ  cosð2hÞ þ 1Þ; c4 ðnÞ ¼ 2 sinðhÞ > > > > 1 > c ðnÞ ¼ ðcosð3nhÞ  cosð3nh þ hÞ  cosð3nh þ 2hÞ þ cosð3nh þ 3hÞ þ cosðhÞ þ cosð2hÞ  cosð3hÞ  1Þ; > 5 2ðcosð2hÞcosðhÞÞ > > > 1 > > c6 ðnÞ ¼ 2 sinð3hÞ ðcosð3nh  hÞ  cosð3nh þ hÞ þ cosð3nh þ 2hÞ  cosð3nh þ 4hÞ  cosð2hÞ þ cosð4hÞÞ; > > > > 1 > c7 ðnÞ ¼ 2 sinð2hÞ ð3 cosð2nh  2hÞ  8 cosð2nh  hÞ þ 6 cosð2nhÞ  4 cosð2nh þ 2hÞ þ 8 cosð2nh þ 3hÞ > > > < 6 cosð2nh þ 4hÞ þ cosð2nh þ 6hÞ þ 8 cosðhÞ þ cosð2hÞ  8 cosð3hÞ þ 6 cosð4hÞ  cosð6hÞ  6Þ; > > > > 1 > ðnÞ ¼ ð3 sinð2nh  2hÞ  8 sinð2nh  hÞ þ 6 sinð2nhÞ  4 sinð2nh þ 2hÞ þ 8 sinð2nh þ 3hÞ  6 sinð2nh þ 4hÞ c 8 > 2 sinð2hÞ > > > > þ sinð2nh þ 6hÞ  8 sinðhÞ þ 7 sinð2hÞ  8 sinð3hÞ þ 6 sinð4hÞ  sinð6hÞÞ; > > > 1 > > c9 ðnÞ ¼ 2 sinð4hÞ ðcosð4nh  hÞ  3 cosð4nhÞ þ 2 cosð4nh þ hÞ þ 2 cosð4nh þ 2hÞ  2 cosð4nh þ 3hÞ  2 cosð4nh þ 4hÞ > > > > > þ2 cosð4nh þ 5hÞ þ 2 cosð4nh þ 6hÞ  3 cosð4nh þ 7hÞ þ cosð4nh þ 8hÞ  3 cosðhÞ  2 cosð2hÞ þ 2 cosð3hÞ > > > > þ2 cosð4hÞ  2 cosð5hÞ  2 cosð6hÞ þ 3 cosð7hÞ  cosð8hÞ þ 3Þ; > > > > 1 > > c10 ðnÞ ¼ 2 sinð4hÞ ðsinð4nh  hÞ  3 sinð4nhÞ þ 2 sinð4nh þ hÞ þ 2 sinð4nh þ 2hÞ  2 sinð4nh þ 3hÞ  2 sinð4nh þ 4hÞ > > > > þ2 sinð4nh þ 5hÞ þ 2 sinð4nh þ 6hÞ  3 sinð4nh þ 7hÞ þ sinð4nh þ 8hÞ  sinðhÞ  2 sinð2hÞ þ 2 sinð3hÞ þ 2 sinð4hÞ > : 2 sinð5hÞ  2 sinð6hÞ þ 3 sinð7hÞ  sinð8hÞÞ: ð90Þ The general solution Z 1;h of the equation corresponding to the homogeneous Eq. (86) is

Z 1;h ðrn ; un ; n; sn Þ ¼ A1 ðr n ; un þ nh; sn Þ;

ð91Þ 2

where A1 is an arbitrary function, which can be used to avoid secular terms at the O (e ) level. At this moment, however, we are not interested in the higher order approximations. For that reason, we take A1  0. An approximation IA ðr n ; un ; n; sn Þ of an invariant Iðr n ; un ; n; sn Þ ¼ constant for (77) is now constructed, and is given by

IA ðr n ; un ; n; sn Þ ¼

1 2 2 ðr cosðhÞ  r 2n þ 2h Þ expðhsn Þ r2n n " he expðhsn Þ 2 3 12h ðc1 ðnÞ sinð2un Þ þ c2 ðnÞ cosð2un ÞÞ þ 6ar n h ðc3 ðnÞ sinðun Þ þ c4 ðnÞ cosðun Þ þ 2 6r n ðcosð2hÞ  1Þ þ c5 ðnÞ sinð3un Þ þ c6 ðnÞ cosð3un ÞÞ þ r 2n ðc7 ðnÞ sinð2un Þ þ c8 ðnÞ cosð2un Þ þ c9 ðnÞ sinð4un Þ # þ c10 ðnÞ cosð4un ÞÞ ;

ð92Þ

where c1 ðnÞ; c2 ðnÞ; . . ., and c10 ðnÞ are given by (90), and

IA ðr nþ1 ; unþ1 ; n þ 1; snþ1 Þ ¼ IA ðrn ; un ; n; sn Þ þ e2 Rðr n ; un ; n; sn ; eÞ;

ð93Þ

where 2

2

e2 Rðrn ; un ; n; sn ; eÞ ¼

2 2 f1;0 cosðhÞ  f1;0 þ 2h f12 cosðhÞ  f12 þ 2h expðhð s þ e ÞÞ  expðhsn Þ n 2 f12 f1;0 2

 

3 3 f1;0 cosðhÞ þ 2h f1;0  f1;0  4hf1;1 3 f1;0

he expðhsn Þ þ

he 6ðcosð2hÞ  1Þ

( expðhðsn þ eÞÞ 2 3 ½12h ðc1 ðn þ 1Þ sinð2f 2 Þ þ c2 ðn þ 1Þ cosð2f 2 ÞÞ þ 6ah f1 ðc3 ðn þ 1Þ sinðf2 Þ f12

þ c4 ðn þ 1Þ cosðf2 Þ þ c5 ðn þ 1Þ sinð3f 2 Þ þ c6 ðn þ 1Þ cosð3f 2 ÞÞ þ f12 ðc7 ðn þ 1Þ sinð2f 2 Þ þ c8 ðn þ 1Þ cosð2f 2 Þ þ c9 ðn þ 1Þ sinð4f 2 Þ þ c10 ðn þ 1Þ cosð4f 2 ÞÞ expðhsn Þ 2 3  ½12h ðc1 ðn þ 1Þ sinð2f 2;0 Þ þ c2 ðn þ 1Þ cosð2f 2;0 ÞÞ þ 6ah f1;0 ðc3 ðn þ 1Þ sinðf2;0 Þ 2 f1;0 2 ðc7 ðn þ 1Þ sinð2f 2;0 Þ þ c4 ðn þ 1Þ cosðf2;0 Þ þ c5 ðn þ 1Þ sinð3f 2;0 Þ þ c6 ðn þ 1Þ cosð3f 2;0 ÞÞ þ f1;0 )

þ c8 ðn þ 1Þ cosð2f 2;0 Þ þ c9 ðn þ 1Þ sinð4f 2;0 Þ þ c10 ðn þ 1Þ cosð4f 2;0 ÞÞ ;

ð94Þ

847

M. Rafei, W.T. Van Horssen / Commun Nonlinear Sci Numer Simulat 18 (2013) 835–850

where f1 ; f2 ; f1;0 ; f2;0 and f1;1 are given by (18) and (77). Since the Rayleigh system (77) is autonomous, and c1 ðnÞ; c2 ðnÞ; . . ., and c10 ðnÞ are bounded, then R is bounded. From (93) it follows that

IA ðr n ; un ; n; sn Þ ¼ IA ðr 0 ; u0 ; 0; 0Þ þ e2

n1 X Rðri ; ui ; i; si ; eÞ:

ð95Þ

i¼0

From (93) and (95), it can be shown that

IA ðr n ; un ; n; sn Þ ¼ constant þ Oðe2 nÞ;

ð96Þ

it then follows that

(

IA ðr n ; un ; n; sn Þ ¼ constant þ Oðe2 Þ for n ¼ Oð1Þ;

ð97Þ

IA ðr n ; un ; n; sn Þ ¼ constant þ OðeÞ for n ¼ Oð1e Þ: From (92) and (97) it then follows that

1 2 2 ðr cosðhÞ  r 2n þ 2h Þ expðhsn Þ ¼ constant þ OðeÞ; r2n n

ð98Þ

on an iteration-scale of order e1 . So far only one approximation of a first integral has been determined. Another (functionally independent) approximation of a first integral can also be obtained in a similar and straightforward way as follows from the next subsection. 5.2. Case 2 Since we are interested in simple approximations of first integrals we now take (see also (80) and (84))

Z 0 ðr n ; un ; n; sn Þ ¼ rncosðhÞ1 expðkn sinðhÞÞ;

ð99Þ

where kn ¼ un þ nh. The OðeÞ-problem (81) will then be cosðhÞ1

Z 1 ðr n ; un ; n; sn Þ ¼ Z 1 ðr n ; un  h; n þ 1; sn Þ þ

rn expððun þ nhÞ sinðhÞÞ 2 3 f12h P1 þ 6ar n h P2 þ r 2n P3 g; 24hð3  4 cosð2hÞ þ cosð4hÞÞ

ð100Þ

where

P1 ¼  cosð2un  3hÞ þ 2 cosð2un  2hÞ þ 2 cosð2un  hÞ  6 cosð2un Þ þ 6 cosð2un þ 2hÞ  2 cosð2un þ 3hÞ  2 cosð2un þ 4hÞ þ cosð2un þ 5hÞ; P2 ¼  sinðun  4hÞ  sinðun  3hÞ þ 5 sinðun  2hÞ þ 3 sinðun  hÞ  9 sinðun Þ  3 sinðun þ hÞ þ 7 sinðun þ 2hÞ þ sinðun þ 3hÞ  2 sinðun þ 4hÞ þ sinð3un  3hÞ  sinð3un  2hÞ  3 sinð3un  hÞ þ 3 sinð3un Þ þ 3 sinð3un þ hÞ  3 sinð3un þ 2hÞ  sinð3un þ 3hÞ þ sinð3un þ 4hÞ;

ð101Þ

P3 ¼ 2 cosð2un  4hÞ þ 8 cosð2un  3hÞ  6 cosð2un  2hÞ  16 cosð2un  hÞ þ 28 cosð2un Þ  28 cosð2un þ 2hÞ þ 16 cosð2un þ 3hÞ þ 6 cosð2un þ 4hÞ  8 cosð2un þ 5hÞ þ 2 cosð2un þ 6hÞ  cosð4un  3hÞ þ 4 cosð4un  2hÞ  3 cosð4un  hÞ  8 cosð4un Þ þ 14 cosð4un þ hÞ  14 cosð4un þ 3hÞ þ 8 cosð4un þ 4hÞ þ 3 cosð4un þ 5hÞ  4 cosð4un þ 6hÞ þ cosð4un þ 7hÞ: Now, we are going to construct a particular solution of (100) in the form cosðhÞ1

Z 1 ðr n ; un ; n; sn Þ ¼ Z 1 ðr n ; un  h; n þ 1; sn Þ þ

rn expððun þ nhÞ sinðhÞÞ 2 ½12h ðd1 ðnÞ sinð2un Þ þ d2 ðnÞ cosð2un ÞÞ 24hð3  4 cosð2hÞ þ cosð4hÞÞ

3

þ 6ar n h ðd3 ðnÞ sinðun Þ þ d4 ðnÞ cosðun Þ þ d5 ðnÞ sinð3un Þ þ d6 ðnÞ cosð3un ÞÞ þ r 2n ðd7 ðnÞ sinð2un Þ þ d8 ðnÞ cosð2un Þ þ d9 ðnÞ sinð4un Þ þ d10 ðnÞ cosð4un ÞÞ:

ð102Þ

By substituting (102) into (100), and by taking apart in the resulting equation the coefficients of sinð2un Þ, cosð2un Þ, r n sinðun Þ, rn cosðun Þ, rn sinð3un Þ, r n cosð3un Þ, r2n sinð2un Þ, r 2n cosð2un Þ, r 2n sinð4un Þ, r2n cosð4un Þ, and setting these coefficients equal to zero, we will have the following system

848

M. Rafei, W.T. Van Horssen / Commun Nonlinear Sci Numer Simulat 18 (2013) 835–850

8 d1 ðn þ 1Þ ¼ cosð2hÞa1 ðnÞ  sinð2hÞa2 ðnÞ þ 3sinðhÞ  6sinð2hÞ þ 6 sinð4hÞ  2sinð5hÞ  2 sinð6hÞ þ sinð7hÞ; > > > > d2 ðn þ 1Þ ¼ sinð2hÞa1 ðnÞ þ cosð2hÞa2 ðnÞ  cosðhÞ þ 6 cosð2hÞ  6cosð4hÞ þ 2cosð5hÞ þ 2 cosð6hÞ  cosð7hÞ  2; > > > > > d3 ðn þ 1Þ ¼ cosðhÞa3 ðnÞ  sinðhÞa4 ðnÞ þ 4 cosðhÞ þ 4 cosð2hÞ  6cosð3hÞ  cosð4hÞ þ 2cosð5hÞ  3; > > > > d4 ðn þ 1Þ ¼ sinðhÞa3 ðnÞ þ cosðhÞa4 ðnÞ þ 14 sinðhÞ þ 2sinð2hÞ  8sinð3hÞ  sinð4hÞ þ 2 sinð5hÞ; > > > > > d5 ðn þ 1Þ ¼ cosð3hÞa5 ðnÞ  sinð3hÞa6 ðnÞ þ cosðhÞ þ 3 cosð2hÞ  3cosð3hÞ  3cosð4hÞ > > > > þ3cosð5hÞ þ cosð6hÞ  cosð7hÞ  1; > > > > > d6 ðn þ 1Þ ¼ sinð3hÞa5 ðnÞ þ cosð3hÞa6 ðnÞ þ sinðhÞ þ 3 sinð2hÞ  3sinð3hÞ  3 sinð4hÞ > > < þ3sinð5hÞ þ sinð6hÞ  sinð7hÞ; > d7 ðn þ 1Þ ¼ cosð2hÞa7 ðnÞ  sinð2hÞa8 ðnÞ  24sinðhÞ þ 30 sinð2hÞ  28sinð4hÞ > > > > þ16sinð5hÞ þ 6sinð6hÞ  8sinð7hÞ þ 2sinð8hÞ; > > > > > ðn þ 1Þ ¼ sinð2hÞa7 ðnÞ þ cosð2hÞa8 ðnÞ þ 8cosðhÞ  26cosð2hÞ þ 28cosð4hÞ d 8 > > > > 16cosð5hÞ  6cosð6hÞ þ 8 cosð7hÞ  2cosð8hÞ þ 6; > > > > > d9 ðn þ 1Þ ¼ cosð4hÞa9 ðnÞ  sinð4hÞa10 ðnÞ  sinðhÞ þ 4sinð2hÞ  3sinð3hÞ  8sinð4hÞ > > > > þ14sinð5hÞ  14sinð7hÞ þ 8 sinð8hÞ þ 3sinð9hÞ  4 sinð10hÞ þ sinð11hÞ; > > > > > ðn þ 1Þ ¼ sinð4hÞa9 ðnÞ þ cosð4hÞa10 ðnÞ þ cosðhÞ  4cosð2hÞ þ 3 cosð3hÞ þ 8 cosð4hÞ d 10 > : 14cosð5hÞ þ 14cosð7hÞ  8 cosð8hÞ  3cosð9hÞ þ 4 cosð10hÞ  cosð11hÞ:

ð103Þ

The eigenvalues of the homogeneous system related to (103) are cosðhÞ  i sinðhÞ; cosð2hÞ  i sinð2hÞ, cosð2hÞ  i sinð2hÞ; cosð3hÞ  i sinð3hÞ, and cosð4hÞ  i sinð4hÞ, and a particular solution of (103) is given by

8 1 > d1 ðnÞ ¼ 2 sinð2hÞ ðcosð2nh  3hÞ  2 cosð2nh  2hÞ  cosð2nh  hÞ þ 4 cosð2nhÞ > > > > 2 cosð2nh þ hÞ þ 2 cosð2nh þ 3hÞ  4 cosð2nh þ 4hÞ þ cosð2nh þ 5hÞ > > > > > þ2 cosð2nh þ 6hÞ  cosð2nh þ 7hÞ þ 3 cosðhÞ þ 2 cosð2hÞ  3 cosð3hÞ > > > > þ4 cosð4hÞ  cosð5hÞ  2 cosð6hÞ þ cosð7hÞ  4Þ; > > > 1 > > d2 ðnÞ ¼ 2 sinð2hÞ ðsinð2nh  3hÞ  2 sinð2nh  2hÞ  sinð2nh  hÞ þ 4 sinð2nhÞ > > > > > 2 sinð2nh þ hÞ þ 2 sinð2nh þ 3hÞ  4 sinð2nh þ 4hÞ þ sinð2nh þ 5hÞ > > > > þ2 sinð2nh þ 6hÞ  sinð2nh þ 7hÞ þ sinðhÞ  2 sinð2hÞ  sinð3hÞ > > > > þ4 sinð4hÞ  sinð5hÞ  2 sinð6hÞ þ sinð7hÞÞ; > > > 1 > > > d3 ðnÞ ¼ 2 sinðhÞ ðsinðnh  4hÞ þ 2 sinðnh  3hÞ  4 sinðnh  2hÞ  8 sinðnh  hÞ > > > > þ6 sinðnhÞ þ 12 sinðnh þ hÞ  4 sinðnh þ 2hÞ  8 sinðnh þ 3hÞ > > > > þ sinðnh þ 4hÞ þ 2 sinðnh þ 5hÞ  20 sinðhÞ þ 10 sinð3hÞ  2 sinð5hÞÞ; > > > 1 > ð cosðnh  4hÞ  2 cosðnh  3hÞ þ 4 cosðnh  2hÞ þ 8 cosðnh  hÞ > d4 ðnÞ ¼ 2 sinðhÞ > > > > 6 cosðnhÞ  12 cosðnh þ hÞ þ 4 cosðnh þ 2hÞ þ 8 cosðnh þ 3hÞ > > > > >  cosðnh þ 4hÞ  2 cosðnh þ 5hÞ þ 4 cosðhÞ  8 cosð2hÞ  6 cosð3hÞ > > > > þ2 cosð4hÞ þ 2 cosð5hÞ þ 6Þ; > > > 1 > d5 ðnÞ ¼ 2 cosðhÞþ1 ðcosð3nh  hÞ  3 cosð3nh þ hÞ þ 3 cosð3nh þ 3hÞ  cosð3nh þ 5hÞ > > > > > þ2 cosðhÞ  3 cosð3hÞ þ cosð5hÞÞ; > > > > d ðnÞ ¼ 1 > ðsinð3nh  hÞ  3 sinð3nh þ hÞ þ 3 sinð3nh þ 3hÞ  sinð3nh þ 5hÞ 6 > 2 cosðhÞþ1 > > > < þ4 sinðhÞ  3 sinð3hÞ þ sinð5hÞÞ; 1 ðcosð2nh  4hÞ  4 cosð2nh  3hÞ þ 4 cosð2nh  2hÞ þ 4 cosð2nh  hÞ d7 ðnÞ ¼ sinð2hÞ > > > > 11 cosð2nhÞ þ 8 cosð2nh þ hÞ  8 cosð2nh þ 3hÞ þ 11 cosð2nh þ 4hÞ > > > > 4 cosð2nh þ 5hÞ  4 cosð2nh þ 6hÞ þ 4 cosð2nh þ 7hÞ  cosð2nh þ 8hÞ > > > > > 12 cosðhÞ  4 cosð2hÞ þ 12 cosð3hÞ  12 cosð4hÞ þ 4 cosð5hÞ þ 4 cosð6hÞ > > > > 4 cosð7hÞ þ cosð8hÞ þ 11Þ; > > > 1 > d ðnÞ ¼ ðsinð2nh  4hÞ  4 sinð2nh  3hÞ þ 4 sinð2nh  2hÞ þ 4 sinð2nh  hÞ > 8 > sinð2hÞ > > > 11 sinð2nhÞ þ 8 sinð2nh þ hÞ  8 sinð2nh þ 3hÞ þ 11 sinð2nh þ 4hÞ > > > > > 4 sinð2nh þ 5hÞ  4 sinð2nh þ 6hÞ þ 4 sinð2nh þ 7hÞ  sinð2nh þ 8hÞ > > > > 4 sinðhÞ þ 4 sinð2hÞ þ 4 sinð3hÞ  10 sinð4hÞ þ 4 sinð5hÞ þ 4 sinð6hÞ > > > > > 4 sinð7hÞ þ sinð8hÞÞ; > > > 1 > d ðnÞ ¼ ðcosð4nh  3hÞ  4 cosð4nh  2hÞ þ 3 cosð4nh  hÞ þ 8 cosð4nhÞ 9 > 2 sinð4hÞ > > > > 13 cosð4nh þ hÞ  4 cosð4nh þ 2hÞ þ 17 cosð4nh þ 3hÞ  17 cosð4nh þ 5hÞ > > > > þ4 cosð4nh þ 6hÞ þ 13 cosð4nh þ 7hÞ  8 cosð4nh þ 8hÞ  3 cosð4nh þ 9hÞ > > > > > þ4 cosð4nh þ 10hÞ  cosð4nh þ 11hÞ þ 10 cosðhÞ þ 8 cosð2hÞ  18 cosð3hÞ > > > > þ17 cosð5hÞ  4 cosð6hÞ  13 cosð7hÞ þ 8 cosð8hÞ þ 3 cosð9hÞ  4 cosð10hÞ cosð11hÞ  8Þ; > > > 1 > ðsinð4nh  3hÞ  4 sinð4nh  2hÞ þ 3 sinð4nh  hÞ þ 8 sinð4nhÞ > d10 ðnÞ ¼ 2 sinð4hÞ > > > > 13 sinð4nh þ hÞ  4 sinð4nh þ 2hÞ þ 17 sinð4nh þ 3hÞ  17 sinð4nh þ 5hÞ > > > > > þ4 sinð4nh þ 6hÞ þ 13 sinð4nh þ 7hÞ  8 sinð4nh þ 8hÞ  3 sinð4nh þ 9hÞ > > > > þ4 sinð4nh þ 10hÞ  sinð4nh þ 11hÞ þ 16 sinðhÞ  16 sinð3hÞ þ 17 sinð5hÞ > > : 4 sinð6hÞ  13 sinð7hÞ þ 8 sinð8hÞ þ 3 sinð9hÞ  4 sinð10hÞ þ sinð11hÞÞ:

ð104Þ

849

M. Rafei, W.T. Van Horssen / Commun Nonlinear Sci Numer Simulat 18 (2013) 835–850

The general solution Z 1;h of the equation corresponding to the homogeneous Eq. (100) is

Z 1;h ðrn ; un ; n; sn Þ ¼ A1 ðr n ; un þ nh; sn Þ;

ð105Þ

where A1 is an arbitrary function, which can be used to avoid secular terms at the O (e2 ) level. At this moment, however, we are not interested in the higher order approximations. For that reason, we take A1  0. An approximation IA ðr n ; un ; n; sn Þ of an invariant Iðrn ; un ; n; sn Þ ¼ constant for (77) is now constructed, and is given by cosðhÞ1

ern expððun þ nhÞ sinðhÞÞ IA ðr n ; un ; n; sn Þ ¼ r cosðhÞ1 expððun þ nhÞ sinðhÞÞ þ n 24hð3  4 cosð2hÞ þ cosð4hÞÞ 2 3  12h ðd1 ðnÞ sinð2un Þ þ d2 ðnÞ cosð2un ÞÞ þ 6ar n h ðd3 ðnÞ sinðun Þ þ d4 ðnÞ cosðun Þ þ d5 ðnÞ sinð3un Þ

þ d6 ðnÞ cosð3un ÞÞ þ r2n ðd7 ðnÞ sinð2un Þ þ d8 ðnÞ cosð2un Þ þ d9 ðnÞ sinð4un Þ þ d10 ðnÞ cosð4un ÞÞ ; ð106Þ where d1 ðnÞ; d2 ðnÞ; . . ., and d10 ðnÞ are given by (104), and so

IA ðr nþ1 ; unþ1 ; n þ 1; snþ1 Þ ¼ IA ðr n ; un ; n; sn Þ þ e2 Rðr n ; un ; n; sn ; eÞ;

ð107Þ

where

n

cosðhÞ1 cosðhÞ2 e2 Rðrn ; un ; n; sn ; eÞ ¼ f1cosðhÞ1 expðf2 sinðhÞÞ  f1;0 expðf2;0 sinðhÞÞ  ef1;0 expðf2;0 sinðhÞÞðf2;1 f1;0 sinðhÞ e

þ f1;1 cosðhÞ  f1;1 Þ þ 24hð3  4 cosð2hÞ þ cosð4hÞÞ h 2 3 cosðhÞ1  f1 expðf2 sinðhÞÞ12h ðd1 ðn þ 1Þ sinð2f 2 Þ þ d2 ðn þ 1Þ cosð2f 2 ÞÞ þ 6ah f1 ðd3 ðn þ 1Þ sinðf2 Þ

þ d4 ðn þ 1Þ cosðf2 Þ þ d5 ðn þ 1Þ sinð3f 2 Þ þ d6 ðn þ 1Þ cosð3f 2 ÞÞ þ f12 ðd7 ðn þ 1Þ sinð2f 2 Þ þ d8 ðn þ 1Þ cosð2f 2 Þ þ d9 ðn þ 1Þ sinð4f 2 Þ þ d10 ðn þ 1Þ cosð4f 2 ÞÞ cosðhÞ1

 f1;0

2

expðf2;0 sinðhÞÞ12h ðd1 ðn þ 1Þ sinð2f 2;0 Þ þ d2 ðn þ 1Þ cosð2f 2;0 ÞÞ

3

þ 6ah f1;0 ðd3 ðn þ 1Þ sinðf2;0 Þ þ d4 ðn þ 1Þ cosðf2;0 Þ þ d5 ðn þ 1Þ sinð3f 2;0 Þ þ d6 ðn þ 1Þ cosð3f 2;0 ÞÞ io 2 ðd7 ðn þ 1Þ sinð2f 2;0 Þ þ d8 ðn þ 1Þ cosð2f 2;0 Þ þ d9 ðn þ 1Þ sinð4f 2;0 Þ þ d10 ðn þ 1Þ cosð4f 2;0 ÞÞ þ f1;0  expððn þ 1Þh sinðhÞÞ;

ð108Þ

where f1 ; f2 ; f1;0 ; f2;0 ; f1;1 and f2;1 are given by (18) and (77). Since the Rayleigh system (77) is autonomous, and d1 ðnÞ; d2 ðnÞ; . . ., and d10 ðnÞ are bounded, then it follows from (108) that R is relatively bounded with respect to Z 0 (the principal part of IA ). From (107) it follows that

IA ðr n ; un ; n; sn Þ ¼ IA ðr 0 ; u0 ; 0; 0Þ þ e2

n1 X Rðri ; ui ; i; si ; eÞ:

ð109Þ

i¼0

From (107) and (109), it can be shown that

IA ðr n ; un ; n; sn Þ ¼ constant þ Oðe2 nÞ;

ð110Þ

it then follows that

(

IA ðr n ; un ; n; sn Þ ¼ constant þ Oðe2 Þ for n ¼ Oð1Þ; IA ðr n ; un ; n; sn Þ ¼ constant þ OðeÞ for n ¼ Oð1e Þ:

ð111Þ

From (106) and (111) it then follows that

rncosðhÞ1 ¼ constant expððun þ nhÞ sinðhÞÞ þ expððun þ nhÞ sinðhÞÞOðeÞ; on an iteration-scale of order

ð112Þ

e1 .

6. Conclusions and remarks In this paper a perturbation method based on invariance vectors and multiple scales has been presented for regularly perturbed weakly nonlinear systems of ODEs. It has been shown that all invariance factors have to satisfy a functional equation, and also it has been shown how secular-free asymptotic approximations for first integrals for such systems can be constructed on long iteration-scales. To show how this perturbation method actually works, the method has been applied to a Van der Pol difference equation, and to a Rayleigh equation. From the computations it follows that the method is elementary, consistent, and straightforward. In Sections 4 and 5 of this paper an asymptotic justification of the presented pertur-

850

M. Rafei, W.T. Van Horssen / Commun Nonlinear Sci Numer Simulat 18 (2013) 835–850

bation method has been given; that is, error estimates have been given on long iteration-scales. The examples in Sections 4 and 5 of this paper show that the presented perturbation method can be applied to a large class of problems to obtain similar asymptotic results on long iteration-scales. The goal of this paper is to present the multiple scales perturbation method, and the method is applied to (relatively) simple problems. In a forthcoming paper a strongly perturbed, nonlinear problem will be studied (that is, when e ¼ 0 the problem is nonlinear). Showing that the method can be applied to a large class of problems. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19]

Nayfeh AH. Perturbation methods. New York: Wiley-Interscience; 1973. Andrianov IV, Manevitch LI. Asymptotology: ideas, methods, and applications. Dordrecht: Kluwer Academic; 2002. Holmes MH. Introduction to perturbation methods. New York: Springer; 1995 [Dordrecht: Kluwer Academic; 2002]. Kevorkian J, Cole JD. Multiple scale and singular perturbation methods. New York: Springer; 1996. Murdock JA. Perturbations, theory and methods. New York: Wiley; 1991. Verhulst F. Methods and applications of singular perturbations. New York: Springer; 2005. Van Horssen WT. A perturbation method based on integrating vectors and multiple scales. SIAM, J Appl Math 1999(59):1444–67. Van Horssen WT. On invariance factors and invariance vectors for difference equations. J Differ Equ Appl 2002(8):1133–46. Van Horssen WT. On asymptotic approximations of first integrals for differential and difference equations. J Indonesian Math Soc 2007(13):1–15. Rafei M, Van Horssen WT. On asymptotic approximations of first integrals for second order difference equations. Nonlinear Dyn 2010(61):535–51. Torng HC. Second order non-linear difference equations containing small parameters. J Franklin Inst 1960(269):97–104. Hoppensteadt FC, Miranker WL. Multitime methods for systems of difference equations. Stud Appl Math 1977(56):273–89. Van Horssen WT, Ter Brake MC. On the multiple scales perturbation methods for difference equations. Nonlinear Dyn 2009(55):401–18. Aczel J. Lectures on functional equations and their applications. New York: Academic Press Inc.; 1966. Kuczma M. Functional equations in a single variable. Poland: Polish Scientific Publishers; 1968. Aczel J, Dhombres J. Functional equations in several variables. New York: Cambridge University Press; 1989. Kuczma M, Choczewski B, Ger R. Iterative functional equations. New York: Cambridge University Press; 1990. Small CG. Functional equations and how to solve them. Springer; 2007. Rafei M, Van Horssen WT. On constructing solutions for the functional equation Zðx; y; nÞ ¼ Zða11 x þ a12 y; a21 x þ a22 y; n þ 1Þ. J Math Anal Appl 2009; submitted for publication. [20] Agarwal RP. Difference equations and inequalities: theory, methods, and applications. New York: Marcel Dekker Inc.; 1992. [21] Waluya SB, Van Horssen WT. Asymptotic approximations of first integrals for a nonlinear oscillator. Nonlinear Anal Theory Methods Appl 2002:1327–46.