On asymptotics for difference equations

(1)

On asymptotics for difference equations

Proefschrift

ter verkrijging van de graad van doctor aan de Technische Universiteit Delft,

op gezag van de Rector Magnificus prof. ir. K.C.A.M. Luyben, voorzitter van het College voor Promoties,

in het openbaar te verdedigen op donderdag 2 februari 2012 om 15.00 uur

door

Morteza RAFEI,

Master of Science in Mechanical Engineering, Mazandaran University, Iran geboren te Tehran, Iran.

(2)

Dit proefschrift is goedgekeurd door de promotor: Prof. dr. ir. A. W. Heemink

Copromotor:

Dr. ir. W. T. van Horssen

Samenstelling promotiecommissie:

Rector Magnificus voorzitter

Prof. dr. ir. A. W. Heemink Technische Universiteit Delft, promotor Dr. ir. W. T. van Horssen Technische Universiteit Delft, copromotor Prof. A. K. Abramian Russian Academy of Sciences, Russia Prof. I. V. Andrianov RWTH Aachen University, Germany Prof. dr. ir. R. E. Kooij Technische Universiteit Delft Prof. dr. ir. H. X. Lin Technische Universiteit Delft Prof. dr. ir. C. W. Oosterlee Technische Universiteit Delft

This thesis has been completed in fulfillment of the requirements of the Delft University of Technology for the award of the Ph.D. degree. The research described in this thesis was carried out in the Delft Institute of Applied Mathematics, Delft University of Technology, The Netherlands.

ISBN 978-94-6186-023-1

Copyright c 2012 by M. Rafei

All rights reserved. No part of the material protected by this copyright notice may be re-produced or utilized in any form or by any means, electronic or mechanical, including pho-tocopying, recording, or by any information storage and retrieval system, without written permission from the author.

Cover design by Hanike Studio.

(3)

(4)

(5)

Chapter 1 Introduction

Be happy for this moment. This moment is your life. Omar Khayyam

E

uler, in the period 1732–1734, developed the fundamental concept of how to make asingle first order, ordinary differential equation (ODE) exact by means of integrating factors. It should be remarked that before 1734 integrating factors were already known and used by, for instance, Leibniz and Bernoulli. Euler, however, showed that all integrating fac-tors for a single first order ODE have to satisfy a first order partial differential equation. Also nth-order, linear ODEs with constant coefficients were solved by Euler by means of an inte-grating factor. Lagrange extended Euler’s method later to nth-order linear ODEs with non-constant coefficients. When the existence of solutions for a system of n first order ODEs has been established on a time-scale it is well-known that this system of n first order ODEs has n, and can not have more than n, functionally independent first integrals on this time-scale. A complete proof for this statement can be found in ((Forsyth 1912), pp. 637–643). Similar remarks on the existence of n functionally independent first integrals can also be found in for instance (Arnold 1978), (Forsyth 1956), (Ince and Sneddon 1987), and (Kaplan 1958). The fundamental concept to obtain these n first integrals by means of integrating vectors was recently presented in (van Horssen 1997) and (van Horssen 1999a). Based on this fundamental concept a perturbation method has been developed by Van Horssen in (van Horssen 1999b), and (van Horssen 2001) to study classes of weakly nonlinear, regularly or singularly perturbed ODEs. Approximations of first integrals for strongly nonlinear oscilla-tors based on the aforementioned concept were also constructed by Waluya in (Waluya and van Horssen 2003b) and (Waluya and van Horssen 2003a). When approximations of inte-grating vectors have been obtained an approximation of a first integral can be given. Also an error-estimate for this approximation of a first integral can be given on a time-scale. Like most methods for differential equations there is an analogous method for difference equations (see, for instance, (Agarwal 1992) and (Elaydi 2005)). Recently, first integrals, invariants and Lie group theory for ordinary difference equations (O∆Es) obtained a lot of attention in the literature (see for instance the list of references in (van Horssen 2002a)). Also recently, the fundamental concept of invariance factors for O∆Es to obtain invariants (or first integrals) for O∆Es has been presented in (van Horssen 2002a). It has been shown in (van Horssen 2002a) that in finding invariants for a system of first order difference

(8)

equa-2 1. Introduction

tions all invariance factors have to satisfy a functional equation (for more information on functional equations we refer the reader to ((Kuczma 1968), (Kuczma et al. 1990), (Aczel and Dhombres 1989), (Risteski and Covachev 2001), (Bowman 2002) and (Small 2007))). The main goal of the study as presented in this thesis is to investigate how the perturbation method based on invariance vectors can be extended to systems of first order O∆Es, and to what classes of nonlinear O∆Es the method is applicable. The problems that will be studied in this thesis are systems of weakly nonlinear, regularly perturbed O∆Es with a Van der Pol type, and with a Rayleigh type of nonlinearity. Also we study a general system of weakly nonlinear, regularly perturbed O∆Es which is applicable, for instance, in the stability analy-sis of a linear single degree of freedom oscillator with a time-varying mass. Throughout the study we are always interested in constructing accurate approximations of the solutions by applying an improved version of the multiple scales perturbation method to our systems. In chapter 2 of this thesis we construct the solutions for the general functional equation

Z(x, y, n) = Z(a11x+ a12y, a21x+ a22y, n + 1), (1.1)

where a11, a12, a21, and a22 are real constants. The functional equation (1.1) is related to a

system of two first order, linear difference equations. Linear transformations and an adapted version of the method of separation of variables will be used to construct the general solu-tion of this funcsolu-tional equasolu-tion. The solusolu-tion of this funcsolu-tional equasolu-tion does not seem to be available in the literature. In chapter 3 we use these solutions to construct asymptotic approximations of first integrals for second order difference equations such as the second order, weakly nonlinear, regularly perturbed O∆E with a Van der Pol type of nonlinearity, that is,

xn+2− 2 cos(θ0)xn+1+ xn= ε(1 − x2n+1)(xn+2− xn), (1.2)

where θ0is a constant, and ε is a small parameter. As mentioned before one of the main

goals of the study as presented in this thesis is to construct accurate approximations of first integrals for systems of O∆Es. Therefore, in chapter 4 a perturbation method based on invari-ance 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 iteration-scales, that is, on iteration-scales of order ε−1, where ε is a small parameter. To show how this perturbation method works, the method is applied to, for instance, a Rayleigh equation, that is,

xn+2− 2 cos(θ )xn+1+ xn+ εah2x2n+1= εh 1 −1 3 x_n+1− x_n h 2 (xn+1− xn), (1.3)

where θ and a are constants, and where h is the discretisation time step. Finally, in chap-ter 5 we apply the multiple scales perturbation method to a generalized system of weakly nonlinear, regularly perturbed ordinary difference equations, that is,

(9)

3 xn+1= f1,0(xn, yn) + ε f1,1(xn, yn) + O(ε2), yn+1= f2,0(xn, yn) + ε f2,1(xn, yn) + O(ε2), (1.4) where,                  f1,0(xn, yn) = a11xn+ a12yn, f2,0(xn, yn) = a21xn+ a22yn, f1,1(xn, yn) = b10+ b11xn+ b12yn+ b13x2n+ b14xnyn+ b15y2n +b16x3n+ b17x2nyn+ b18xny2n+ b19y3n, f2,1(xn, yn) = b20+ b21xn+ b22yn+ b23x2n+ b24xnyn+ b25y2n +b₂₆x3_n+ b27x2nyn+ b28xny2n+ b29y3n, (1.5)

and where ai j and bi jare constants, and 0 < ε 1. Such systems arise as a result of the

discretization of a system of nonlinear differential equations, or as a result in the stabil-ity analysis of nonlinear oscillations. In our procedure, asymptotic approximations of the solutions of the difference equations will be constructed which are valid on long iteration scales. The obtained results help us in a better understanding of the behavior of nonlinear oscillations. In particular, all kinds of bifurcations can be studied in detail.

(10)

(11)

Submitted as: M. Rafei and W.T. Van Horssen – ”On constructing solutions for the functional equation Z(x, y, n) =Z(a11x+ a12y, a21x+ a22y, n + 1)”, Journal of Mathematical Analysis and Applications, December

2009.

Chapter 2 On constructing solutions for the functional

equation Z(x, y, n) =Z(a

₁₁

x

+ a

₁₂

y, a

₂₁

x

+ a

₂₂

y, n + 1)

The question is not what you look at, but what you see. Henry David Thoreau

Abstract

The concept of invariance factors to obtain first integrals for difference equations was presented by (van Horssen 2002a). It was shown that all invariance factors have to satisfy a functional equation. One of the main difficulties in finding first integrals for a system of first order difference equations is solving the aforementioned functional equa-tion. In this chapter, we consider a functional equation which is related to a system of two first order, linear difference equations. Linear transformations and an adapted version of the method of separation of variables will be used to construct the general solution of this functional equation. The solution of this functional equation does not seem to be available in the literature.

2.1 Introduction

D

ifference equations play an important role in applications and in numerical analysis (see (Elaydi 2005), (Agarwal 1992), (Agarwal et al. 2000), (Gy¨ori and Ladas 1991) and (Mickens 1990)). The concept of invariance factors and invariance vectors to obtain invariants (or first integrals) for difference equations was presented by (van Horssen 2002a). It has been shown in (van Horssen 2002a) and (van Horssen 2007) that all invariance factors have to satisfy a functional equation. This concept turns out to be analogous to the concept of integrating factors and integrating vectors for ordinary differential equations. As it has been shown in (van Horssen 1999a), all integrating factors for a system of n first order, ordinary differential equations have to satisfy a system of1₂n(n + 1) first order, linear partial differential equations. Recently, first integrals, invariants and Lie group theory for ordinary difference equations (O∆Es) obtained a lot of attention in the literature (see the references in (van Horssen 2002a)). For instance, (Maeda 1987) showed that autonomous systems of first order O∆Es can be simplified or solved using an extension of Lie’s method. It has also been

(12)

6 2. On constructing solutions for a functional equation

shown in (Maeda 1987) that the linearized symmetry condition for such O∆Es leads to a set of functional equations, which are (in general) hard to solve. Extensions to nonautonomous systems and higher order O∆Es can be found in (Hydon 2000).

Finding an invariance vector for a given system of first order difference equations is a difficult and usually impossible task. If we consider a system of two first order, linear difference equations

xn+1= a11xn+ a12yn= f1(xn, yn, n) ,

yn+1= a21xn+ a22yn= f2(xn, yn, n) ,

(2.1)

where a11, a12, a21, and a22are real constants, then invariance factors µ1(xn, yn, n) and µ2(xn, yn, n)

have to satisfy the functional equation (see (van Horssen 2002a))

µ1(xn, yn, n) xn+ µ2(xn, yn, n) yn= µ1( f1, f2, n + 1) f1+ µ2( f1, f2, n + 1) f2, (2.2)

or equivalently,

Z(x, y, n) = Z(a11x+ a12y, a21x+ a22y, n + 1), (2.3)

and an invariant Z(x, y, n) = constant for the system (2.1) is given by

Z(x, y, n) = µ1(x, y, n) x + µ2(x, y, n) y. (2.4)

The purpose of this chapter is to find the general solution of (2.3) which can help us in obtaining the invariants for second order weakly perturbed linear problems. It is obvious that this solution plays an important role in perturbation theory. Moreover, we use the concept of separation of variables for partial difference equations, which was presented by (van Horssen 2002b). For the references concerning functional equations in a single variable, or in several variables, and complex vector functional equations, we refer the reader to Refs. ((Kuczma 1968), (Kuczma et al. 1990), (Aczel and Dhombres 1989), (Risteski and Covachev 2001), (Bowman 2002) and (Small 2007)).

This chapter is organised as follows. In section 2 and in section 3 the general solution of (2.3) will be constructed. Linear transformations and an adapted version of the method of separation of variables will be used to construct this general solution of the functional equation. Finally, in section 4 some conclusions will be drawn and some future directions for research will be indicated.

(13)

2.2. On solving the functional equation 7

In this section, we are going to construct the general solution of the functional equation (2.3). To obtain the solution of (2.3), we apply the method of separation of variables by substituting K(x, y)N(n) for Z(x, y, n) into (2.3), yielding

K(x, y)N(n) = K(k1, k2)N(n + 1), (2.5)

where k1= a11x+ a12y, and k2= a21x+ a22y. So,

K(x, y) K(k1, k2)

=N(n + 1)

N(n) = γ, (2.6)

where γ is a (complex-valued) separation constant. From (2.6) it follows that

K(x, y) = γK(k1, k2), and N(n + 1) = γN(n). (2.7)

Now, the transformation K(x, y) = G(z, w), where z = ax + by, and w = cx + dy in which a, b, c, and d are (complex-valued) constants, will be used to have

G(z, w) = γG(ak1+ bk2, ck1+ dk2), (2.8)

or equivalently,

G(z, w) = γG(g11x+ g12y, g21x+ g22y), (2.9)

where g11= aa11+ ba21, g12= aa12+ ba22, g21= ca11+ da21, g22= ca12+ da22. Now, let

g11x+ g12y= λ1z, and g21x+ g22y= λ2w. Therefore,

G(z, w) = γG(λ1z, λ2w), (2.10)

where λ1and λ2are (complex-valued) constants. From (2.9) and (2.10), we have

       aa11+ ba21= λ1a, aa12+ ba22= λ1b, ca11+ da21= λ2c, ca12+ da22= λ2d. (2.11)

To have a non-trivial solution for a, b, c, and d, it follows that

λ1,2=

a11+ a22±p(a11+ a22)2− 4(a11a22− a12a21)

2 , (2.12)

and we observe that λ1,2 are the eigenvalues of the matrix corresponding to the system of

linear, ordinary difference equations (2.1). The transformations leading to (2.10) are nonde-generate (for instance the transformation (x, y) ↔ (z, w) is bijective) when the eigenvalues λ1 and λ2 are distinct and nonzero, or when the dimension of the eigenspace is two for

(14)

coinciding eigenvalues. When the eigenvalues are coinciding and the dimension of the cor-responding eigenspace is one, or when at least one of the eigenvalues is zero, then a slightly different transformation should be introduced to avoid degeneracy. This case will be dis-cussed in the next section of this chapter. Now we introduce real constants r1, θ1, r2, and θ2

such that

λ1= r1exp(iθ1), and λ2= r2exp(iθ2), (2.13)

where ri> 0, and 0 6 θi< 2π for i = 1, 2. Then the method of separation of variables is

applied to (2.10) by putting G(z, w) = H(z)F(w), and therefore, H(z)

H(λ1z)

=γ F (λ2w)

F(w) = α, (2.14)

where α is a (complex-valued) separation constant. So, we now have two functional equa-tions with complex variables z and w. First, we solve the first one, that is,

H(z) = αH(λ1z), (2.15)

and then by a simple transformation we obtain the solution of the second one. Equation (2.15) is a q-difference equation. For more information on q-difference equations we refer the reader to Refs. ((Polyanin and Manzhirov 2007), (Carmichael 1912), (Heittokangas et al. 2000), (Barnett et al. 2007) and (Aczel 1966)). To obtain the general solution of (2.15) we use (2.13) and introduce

α = α1+ iα2, z= r exp(iθ ), H(z) = J(r, θ ), and J(r, θ ) = J1(r, θ ) + iJ2(r, θ ),

(2.16) where α1and α2are real-valued constants, and where J1and J2are real-valued functions.

Then, (2.15) becomes

J1(r, θ ) + iJ2(r, θ ) = (α1+ iα2)(J1(rr1, θ + θ1) + iJ2(rr1, θ + θ1)), (2.17)

and so, we obtain from (2.17)

J1(r, θ ) = α1J1(rr1, θ + θ1) − α2J2(rr1, θ + θ1), (2.18)

and,

J2(r, θ ) = α2J1(rr1, θ + θ1) + α1J2(rr1, θ + θ1). (2.19)

Now, we have to consider two cases: the eigenvalues λ1and λ2are real, or complex valued.

In subsection 2.1 we will first consider the case of real-valued eigenvalues, and in subsection 2.2 the case of complex-valued eigenvalues will be considered.

(15)

2.2.1 Real-valued eigenvalues

If λ1,2are real valued, then it follows from (2.13) that

λ1= r1, λ2= r2, and θ1= θ2= 0. (2.20)

We now consider two cases for α2: α2= 0, or α26= 0. For α2= 0 it follows from (2.18) and

(2.19) that

J1(r, θ ) = α1J1(rr1, θ ), and J2(r, θ ) = α1J2(rr1, θ ). (2.21)

Now we consider two cases for r1: r1= 1, or r16= 1.

Case 1. r1= 1 From (2.21) it follows that there are only nontrivial solutions for J1(r, θ )

when α1= 1. Now we consider two cases for r2: r2= 1, or r26= 1.

Case 1.1. r1= 1, and r2= 1 From (2.10) it then follows that there are only nontrivial

solutions for G(z, w) when γ = 1. Then, from (2.6) it follows that Z(x, y, n) = A(x, y), where Ais an arbitrary function.

Case 1.2. r1= 1, and r26= 1 From (2.10)-(2.21) it follows in this case that

G(z, w) = γG(z, r2w), (2.22)

where γ is a real-valued constant, and where z, and w are real-valued variables. If we intro-duce

G(z, w) = U (z, u), (2.23)

where u =ln(|w|)_ln(r

2), then (2.22) becomes

U(z, u) = γU(z, u + 1), (2.24)

by multiplying both sides in (2.24) by γu, it follows that

S(z, u) = S(z, u + 1), (2.25)

where S(z, u) = γuU(z, u). So, the general solution for G(z, w) is given by

G(z, w) = γ−

ln(|w|)

ln(r2)g(z,ln(|w|) ln(r2)

(16)

where g is an arbitrary function which is 1-periodic in its second argument. Since K(x, y) = G(z, w), Z(x, y, n) = K(x, y)N(n), and N(n) = c0γn, where c0 is a constant, it then easily

follows that the general solution of Z(x, y, n) is given by

Z(x, y, n) = A(a21x+ (a22− r2)y, ((1 − a11)x − a12y)r₂−n, (1 − a11)x − a12y), (2.27)

or equivalently,

Z(x, y, n) = B((a11− r2)x + a12y, (a21x+ (a22− 1)y)r2−n, a21x+ (a22− 1)y), (2.28)

Since, following (2.25), the function S(z, u) is a 1-periodic function in u and we are only interested in two independent invariants and not in equilibrium solutions we exclude the last argument in the general solution for Z(x, y, n). Other possible solutions are a combination of the first two arguments. As a result we consider the ”general” solution for Z(x, y, n) as

Z(x, y, n) = A(a21x+ (a22− r2)y, ((1 − a11)x − a12y)r−n₂ ), (2.29)

or,

Z(x, y, n) = B((a11− r2)x + a12y, (a21x+ (a22− 1)y)r−n2 ), (2.30)

where A and B are arbitrary functions. In the subsequent cases, we also cancel out for the same reason the 1-periodic arguments in the general solutions.

Case 2. r16= 1 In this case, the general solutions for (2.21) can readily be obtained (see

(Polyanin and Manzhirov 2007)), and are given by J1(r, θ ) = α −ln(r) ln(r1) 1 j1(θ ), and J2(r, θ ) = α −ln(r) ln(r1) 1 j2(θ ), (2.31)

where j1(θ ) and j2(θ ) are arbitrary functions. Now, from (2.16) we obtain

H(z) = α− ln(|z|) ln(r1)_h z |z| , (2.32)

where h is a complex-valued function. Now, from (2.14) we have F(w) = γ

αF(r2w), (2.33)

which can be solved by some simple transformations α → γ α, z→ w, r1→ r2, and H→ F. (2.34) Therefore, we have F(w) =γ α −ln(|w|)_ln(r2) f w |w| , (2.35)

where f is an arbitrary function. Now we consider two cases for r2: the first case is r2= r1,

(17)

Case 2.1. r16= 1, and r1= r2 From (2.35) we obtain

F(w) =γ α −ln(|w|) ln(r1) _f w |w| . (2.36)

In this case, z and w are real. Therefore, after some calculations, since we have G(z, w) = H(z)F(w), we obtain G(z, w) = γ− ln(|w|) ln(r1)_gz w . (2.37) From (2.7), we have N(n) = c0γn, (2.38)

where c0 is a constant. Since K(x, y) = G(z, w) and Z(x, y, n) = K(x, y)N(n), the general

solution for Z(x, y, n) is given by

Z(x, y, n) = A(xr−n₁ , yr−n₁ ). (2.39)

Case 2.2. r16= 1, r26= r1, and r26= 1 Since G(z, w) = H(z)F(w), from (2.32) and (2.35)

we obtain G(z, w) = γ− ln(|w|) ln(r2)g ln(|z|) ln(r1) −ln(|w|) ln(r2) . (2.40)

since K(x, y) = G(z, w) and Z(x, y, n) = K(x, y)N(n), the general solution for Z(x, y, n) is given by

Z(x, y, n) = A((a21x+ (a22− r2)y)r−n₁ , ((r1− a11)x − a12y)r−n₂ ), (2.41)

or,

Z(x, y, n) = B(((a11− r2)x + a12y)r−n1 , (a21x+ (a22− r1)y)r−n2 ), (2.42)

where A and B are arbitrary functions. To verify whether (2.41) leads to the invariants of (2.1) we have to prove that

∆((a21x+ (a22− r2)y)r−n1 ) = 0, (2.43)

and,

∆(((r1− a11)x − a12y)r−n2 ) = 0, (2.44)

where ∆( f (xn, yn, n)) = f (xn+1, yn+1, n + 1) − f (xn, yn, n). By taking apart in (2.43) and

(2.44) the coefficients of x, and y, and setting these coefficients equal to zero, we will have r2_1,2− (a11+ a22)r1,2+ a11a22− a12a21= 0, (2.45)

(18)

which is in fact (2.12), and therefore (2.41) indeed leads to the invariants of (2.1). For α26= 0,

from (2.18) we obtain J2(r, θ ) = α1 α2 J1(r, θ ) − 1 α2 J1 r r1 , θ , (2.46)

then, by substituting (2.46) into (2.19) we obtain

α12+ α22 J1(rr1, θ ) − 2α1J1(r, θ ) + J1 r r1 , θ = 0. (2.47)

Now we consider two cases for r1: r1= 1, or r16= 1.

Case 1. r1= 1 From (2.47), it follows that

α12+ α22 J1(r, θ ) − 2α1J1(r, θ ) + J1(r, θ ) = 0, (2.48)

which has no nontrivial solution. Also, from (2.15), it follows that if r1= 1, then α2= 0,

which is a contradiction. Therefore, there is no solution for Z(x, y, n) in this case.

Case 2. r16= 1 We introduce

J1(r, θ ) = L1(p, θ ), (2.49)

where p =_ln(rln(r)

1). Therefore, (2.47) becomes

α₁2+ α₂2 L1(p + 1, θ ) − 2α1L1(p, θ ) + L1(p − 1, θ ) = 0. (2.50)

Substituting L1(p, θ ) = V (p)W (θ ) into (2.50) yields

α₁2+ α₂2 V (p + 1) − 2α1V(p) +V (p − 1) = 0. (2.51)

By introducing α = r3exp(iθ3), two independent real-valued solutions for V (p) can be

writ-ten as V(p) = c1r₃−pcos(pθ3), (2.52) and, V(p) = c2r −p 3 sin(pθ3), (2.53)

where c1and c2are arbitary constants. When we use the first one, L1(p, θ ) is given by

(19)

and consequently, from (2.16) and (2.49), we have

J1(r, θ ) = r −ln(r) ln(r1) 3 h1(θ ) cos ln(r) ln(r1) θ3 , (2.55)

where h1is an arbitary function. From (2.46) we obtain J2. Then, by using (2.16) we obtain

H(z) = α− ln(|z|) ln(r1)_h z |z| , (2.56)

where h is an arbitary function. By using (2.53) instead of (2.52) we obtain the same solution for H(z). Now, from (2.14) we have

F(w) = γ

αF(r1w), (2.57)

which can be solved by some simple transformations α → γ α, z→ w, r1→ r2, and H→ F. (2.58) Therefore, we have F(w) =γ α −ln(|w|)_ln(r2) f w |w| . (2.59)

Now we consider two cases for r2: the first case is r2= r1, and the second one is r26= r1,

r26= 1. Since the general solution for F(w) is the same as (2.35), the general solution for

Z(x, y, n) in the first case is as presented in (2.39), and the general solution in the second case is as presented in (2.41), or (2.42).

2.2.2 Complex-valued eigenvalues

If λ1,2are complex valued, then it follows from (2.12) and (2.13) that

r1= r2=

√

a11a22− a12a21, a11+ a22− 2r1cos(θ1) = 0, and θ2= −θ1, (2.60)

where 0 < θ1< 2π. We now consider two cases for α2: α2= 0, or α26= 0.

For α2= 0, from (2.18) and (2.19) we obtain

J1(r, θ ) = α1J1(rr1, θ + θ1), and J2(r, θ ) = α1J2(rr1, θ + θ1), (2.61)

(20)

Case 1. r1= 1 From (2.61) we have

J1(r, θ ) = α1J1(r, θ + θ1), (2.62)

and

J2(r, θ ) = α1J2(r, θ + θ1). (2.63)

The general solutions of (2.62) and (2.63) can readily be obtained, and are given by

J1(r, θ ) = α −θ θ1 1 j1(r), and J2(r, θ ) = α −θ θ1 1 j2(r), (2.64)

where j1(r) and j2(r) are arbitrary functions. Now, from (2.16) we obtain

H(z) = α i θ1ln z |z| h(|z|), (2.65)

αF(λ2w), (2.66)

which can be solved by some simple transformations

α → γ α, λ1→ λ2, z→ w, and H→ F. (2.67) Therefore, we have F(w) =γ α −_θ1i ln w |w| f(|w|). (2.68)

After some calculations, since we have G(z, w) = H(z)F(w), we obtain

G(z, w) = γ− i θ1ln w |w| g |z| , |w| , zw |zw| . (2.69) From (2.7) we have N(n) = c₀γn, (2.70)

where c0 is a constant. Since K(x, y) = G(z, w) and Z(x, y, n) = K(x, y)N(n), the general

Z(x, y, n) = A((a21x+ (a22− exp(−iθ1))y) exp(−inθ1),

((exp(iθ1) − a11)x − a12y) exp(inθ1)).

(21)

J1(r, θ ) = α1J1(rr1, θ + θ1). (2.72) Now we introduce J1(r, θ ) = L1(p, q), (2.73) where p =_ln(rln(r) 1) and q = θ θ1. Therefore, (2.72) becomes L1(p, q) = α1L1(p + 1, q + 1). (2.74)

Now we apply the method of separation of variables. Substituting L1(p, q) = V (p)W (q), in

which V (p) and W (q) are real-valued functions, into (2.74), yields α1V(p + 1)

V(p) =

W(q)

W(q + 1) = β , (2.75)

where β is an arbitary separation constant. The equation can easily be solved, yielding

V(p) = c1 β α1 p , and W(q) = c2 1 β q , (2.76)

where c1and c2are arbitary constants. Therefore,

L1(p, q) = α₁−pR(p − q), (2.77) and so we have J1(r, θ ) = α −ln(r) ln(r1) 1 j1 ln(r) ln(r1) − θ θ1 , (2.78)

and similarly we can find from (2.61) that

J2(r, θ ) = α −ln(r) ln(r1) 1 j2 ln(r) ln(r1) − θ θ1 , (2.79)

where j1(r) and j2(r) are arbitrary functions. Now, from (2.16) we obtain

H(z) = α− ln(|z|) ln(r1)_h ln(|z|) ln(r1) + i θ1 ln z |z| , (2.80)

αF(λ2w), (2.81)

which can be solved by some simple transformations

α → γ

(22)

16 2. On constructing solutions for a functional equation Therefore, we have F(w) =γ α −ln(|w|)_ln(r1) f ln(|w|) ln(r1) − i θ1 ln w |w| . (2.83)

After some calculations, since we have G(z, w) = H(z)F(w), we obtain

G(z, w) = γ− ln(|w|) ln(r1)_g ln(|z|) ln(r1) + i θ1 ln z |z| ,ln(|w|) ln(r1) − i θ1 ln w |w| , z w . (2.84) From (2.7) we have N(n) = c0γn. (2.85)

Since K(x, y) = G(z, w) and Z(x, y, n) = K(x, y)N(n), the general solution for Z(x, y, n) is given by

Z(x, y, n) = A((a₂₁x+ (a₂₂− r1exp(−iθ1))y)r−n1 exp(−inθ1),

((r1exp(iθ1) − a11)x − a12y)r−n₁ exp(inθ1)).

(2.86) For α26= 0, from (2.18) we obtain

J2(rr1, θ + θ1) = α1 α2 J1(rr1, θ + θ1) − 1 α2 J1(r, θ ), (2.87) and therefore, J2(r, θ ) = α1 α2 J1(r, θ ) − 1 α2 J1 r r1 , θ − θ1 . (2.88)

Then, by substituting (2.88) into (2.19) we obtain

α₁2+ α₂2 J1(rr1, θ + θ1) − 2α1J1(r, θ ) + J1 r r1 , θ − θ1 = 0. (2.89)

Now, we consider two cases for r1: r1= 1, or r16= 1.

α₁2+ α₂2 J1(r, θ + θ1) − 2α1J1(r, θ ) + J1(r, θ − θ1) = 0. (2.90) Now we introduce J1(r, θ ) = L1(r, q), (2.91) where q = θ θ1. Therefore, (2.90) becomes α12+ α22 L1(r, q + 1) − 2α1L1(r, q) + L1(r, q − 1) = 0. (2.92)

(23)

Now, we apply the method of separation of variables. Substituting V (r)W (q) for L1(r, q), in

which V (r) and W (q) are real-valued functions, into (2.92) yields

α₁2+ α22 W (q + 1) − 2α1W(q) +W (q − 1) = 0, (2.93)

By introducing α = r3exp(iθ3), two independent real-valued solutions for W (q) can be

writ-ten as

W(q) = c1r₃−qcos(qθ3), (2.94)

and,

W(q) = c2r₃−qsin(qθ3), (2.95)

where c1and c2are arbitary constants. When we use the first one, L1(r, q) is given by

L1(r, q) = r −q

3 cos(qθ3)R(r), (2.96)

and consequently, a solution for J1(r, θ ) is given by

J1(r, θ ) = r −θ θ1 3 cos θ θ3 θ1 j1(r), (2.97)

where j1is an arbitary function. From (2.88) we can obtain similarly J2. Then, by using

(2.16) we have H(z) = α i θ1ln z |z| h(|z|), (2.98)

where h is a complex-valued function. By using (2.95) instead of (2.94) we obtain the same solution for H(z). Since H(z) is the same as H(z) in case 1 (when α2= 0), and also r1= 1,

we have the same result for Z(x, y, n), that is,

((exp(iθ1) − a11)x − a12y) exp(inθ1)). (2.99) Case 2. r16= 1 We introduce in (2.89) J1(r, θ ) = L1(p, q), (2.100) where p = _ln(rln(r) 1), and q = θ

θ1. To solve the partial difference equation (2.89) we apply the

adapted version of the method of separation of variables as presented in (van Horssen 2002b). Substituting L1(p, q) = V (p)W (q) into (2.89), and dividing by W (q + 1) yields

α₁2+ α₂2 V (p + 1) − 2α1V(p)

W(q)

W(q + 1)+V (p − 1)

W(q − 1)

(24)

Obviously, (2.101) is not directly separable. When we, however, apply the difference opera-tor ∆q= Eq− 1 to (2.101) we obtain −2α1V(p)∆q _W_(q) W(q + 1) +V (p − 1)∆q W (q − 1) W(q + 1) = 0, (2.102)

which is a separable equation if we divide (2.102) by ∆q(W (q − 1)/W (q + 1)). From (2.102)

it follows that V(p) = βV (p − 1), and − 2α1β ∆q W(q) W(q + 1) + ∆q W (q − 1) W(q + 1) = 0, (2.103)

where β is an arbitrary separation constant. Substituting (2.103) into (2.101) and noticing that V (p + 1) = β2V(p − 1) we obtain

α₁2+ α22

β2W(q + 1) − 2α1βW (q) + W (q − 1) = 0. (2.104)

By introducing α = r3exp(iθ3), two independent real-valued solutions for W (q) can be

writ-ten as W(q) = c1 1 β r3 q cos(qθ3), (2.105) and, W(q) = c2 1 β r3 q sin(qθ3), (2.106)

where c1and c2are arbitary constants. From (2.103) it follows that

V(p) = c₃βp, (2.107)

and, therefore a solution for L1(p, q) is given by

L1(p, q) = r −q

3 cos(qθ3)R(p − q), (2.108)

where R(p − q) is an arbitrary function, and consequently, a solution for J1(r, θ ) is given by

J1(r, θ ) = r −θ θ1 3 cos θ θ3 θ1 j1 ln(r) ln(r1) − θ θ1 , (2.109)

where j1is an arbitary function. From (2.88) we obtain J2. Then, by using (2.16) we obtain

H(z) = α i θ1ln z |z| h ln(|z|) ln(r1) + i θ1 ln z |z| , (2.110)

where h is an arbitary function. Now, from (2.14) we have F(w) = γ

(25)

2.3. Degenerate Cases 19

which can be solved by some simple transformations α → γ

α, λ1→ λ2, z→ w, and H→ F. (2.112)

Now, by considering (2.60) and (2.111), we have

F(w) =γ α −_θ1i ln w |w| f ln(|w|) ln(r1) − i θ1 ln w |w| . (2.113)

After some calculations, since we have G(z, w) = H(z)F(w), we obtain G(z, w) = γ− i θ1ln w |w| gln(|z|)_ln(r 1)+ i θ1ln z |z| ,ln(|w|)_ln(r 1) − i θ1ln w |w| ,_|zw|zw. (2.114) By using (2.106) instead of (2.105) we obtain the same solution for G(z, w). From (2.7) we have

N(n) = c0γn. (2.115)

Since K(x, y) = G(z, w) and Z(x, y, n) = K(x, y)N(n), the general solution for Z(x, y, n) is given by

Z(x, y, n) = A((a₂₁x+ (a₂₂− r1exp(−iθ1))y)r1−nexp(−inθ1),

(2.116) where the constants r1and θ1can be obtained from (2.12) and (2.13).

2.3 Degenerate Cases

In this section, we discuss two degenerate cases. When the eigenvalues are coinciding, and the dimension of the corresponding eigenspace is one, or when at least one of the eigenvalues is zero. For these cases a slightly different transformation for (z, w) should be introduced to avoid degeneracy in the transformation. We introduce the transformation Z(x, y, n) = P(z, w, n), where

z= a(n)x + b(n)y,

w= c(n)x + d(n)y. (2.117)

By applying the same procedure as presented before a separable form P(z, w, n) = P(µ1z, µ2w, n+

1) can be obtained when a(n), b(n), c(n), and d(n) satisfy       

a(n + 1)a11+ b(n + 1)a21= µ1a(n),

a(n + 1)a12+ b(n + 1)a22= µ1b(n),

c(n + 1)a11+ d(n + 1)a21= µ2c(n),

c(n + 1)a12+ d(n + 1)a22= µ2d(n),

(26)

where µ1, and µ2are constants.

Case 1. r1= r26= 0, and the dimension of the corresponding eigenspace is one In this

case we have r1= r2=a11+a₂ 22, and we can choose µ1= µ2= r1. The solution of the system

(2.118), where a(0) = 1, b(0) = 0, c(0) = 0, and d(0) = 1, is given by        a(n) = r₁+ n(r₁− a11), b(n) = −na12, c(n) = −na21, d(n) = r1+ n(r1− a22). (2.119)

Since the solution for P(z, w, n) is A0(zr−n1 , wr −n

1 ), where A0 is an arbitrary function, the

Z(x, y, n) = A(((r1+ n(r1− a11))x − na12y)r1−n, (na21x− (r1+ n(r1− a22))y)r1−n), (2.120)

where A is an arbitrary function.

Case 2.1. r1= 0, and r26= 0 Since one of the eigenvalues is zero, a11a22− a12a21= 0, and

r2= a11+ a22. From (2.11) we obtain

a= a₂₁, b = −a₁₁, c = a₁₁, and d= a₁₂, (2.121) and by applying the same procedure as presented before, the solution for Z(x, y, n) is given by Z(x, y, n) = A 0 n a21x− a11y , (a₁₁x+ a₁₂y)r−n₂ , (2.122) or, Z(x, y, n) = B 0 n a22x− a12y , (a21x+ a22y)r−n2 , (2.123)

where 00_{is defined to be 1, and where A and B are arbitrary functions.}

Case 2.2. r1= r2= 0 Since both of the eigenvalues are zero, a11a22− a12a21= 0, and

a11+ a22= 0. If the dimension of the corresponding eigenspace is two, then the solution for

Z(x, y, n) is given by Z(x, y, n) = A0 n x, 0n y , (2.124)

(27)

2.4. Conclusions 21

and if the dimension of the eigenspace is one, then, from (2.117) and (2.118), the solution for Z(x, y, n) is given by Z(x, y, n) = A 0 n+1 (0 − n2_a 11)x − n2a12y , 0 n+1 n2_a 21x− (0 − n2a22)y , (2.125)

where A is an arbitrary function.

2.4 Conclusions

In this chapter the solution of a functional equation, which has its origin in the process of finding first integrals for a system of two first order ordinary difference equations, has been presented. The results obtained in this chapter can be summerised as follows.

THEOREM. Consider a system of two first order linear difference equations,

xn+1= a11xn+ a12yn= f1(xn, yn, n) ,

yn+1= a21xn+ a22yn= f2(xn, yn, n) ,

(2.126)

where a11, a12, a21, and a22are real constants, and the eigenvalues of the system are

λ1= r1exp(iθ1), and λ2= r2exp(iθ2). (2.127)

Then, invariance factors µ1(xn, yn, n) and µ2(xn, yn, n) have to satisfy the functional equation

µ1(xn, yn, n) xn+ µ2(xn, yn, n) yn= µ1( f1, f2, n + 1) f1+ µ2( f1, f2, n + 1) f2, (2.128)

or equivalently,

Z(x, y, n) = Z(a₁₁x+ a₁₂y, a₂₁x+ a₂₂y, n + 1), (2.129) and an invariant Z(x, y, n) = constant for the system of difference equations (2.126) is given by

Z(x, y, n) = µ1(x, y, n) x + µ2(x, y, n) y. (2.130)

If λ1,2 are real valued, and r1= r2= 1, then Z(x, y, n) = A(x, y), where A is an arbitrary

function. If λ1,2are real valued, and also r1= 1 and r26= 1, then a solution of (2.129) is

Z(x, y, n) = A(a21x+ (a22− r2)y, ((1 − a11)x − a12y)r−n2 ), (2.131)

or,

(28)

where A and B are arbitrary functions. If λ1,2are real valued, and r1= r26= 1, then a solution

of (2.129) is

Z(x, y, n) = A(xr−n₁ , yr−n₁ ). (2.133) If λ1,2are real valued, and also r16= 1, r26= 1 and r16= r2, then a solution of (2.129) is

Z(x, y, n) = A((a21x+ (a22− r2)y)r1−n, ((r1− a11)x − a12y)r−n2 ), (2.134)

or,

Z(x, y, n) = B(((a11− r2)x + a12y)r₁−n, (a21x+ (a22− r1)y)r−n₂ ). (2.135)

If λ1,2are complex valued, and r1= r2= 1, then a solution of (2.129) is

((exp(iθ1) − a11)x − a12y) exp(inθ1)).

(2.136)

If λ1,2are complex valued, and r1= r26= 1, then a solution of (2.129) is

Z(x, y, n) = A((a21x+ (a22− r1exp(−iθ1))y)r₁−nexp(−inθ1),

(2.137)

If r1= r26= 0, and the dimension of the corresponding eigenspace is one, then a solution of

(2.129) is

Z(x, y, n) = A(((r1+ n(r1− a11))x − na12y)r1−n, (na21x− (r1+ n(r1− a22))y)r1−n). (2.138)

If r1= 0, and r26= 0, then a solution of (2.129) is

Z(x, y, n) = A 0 n a21x− a11y , (a11x+ a12y)r−n₂ , (2.139) or, Z(x, y, n) = B 0 n a22x− a12y , (a21x+ a22y)r−n2 , (2.140) where 00_{is defined to be 1. If r}

1= r2= 0, and the dimension of the eigenspace is two, then

a solution of (2.129) is Z(x, y, n) = A0 n x, 0n y , (2.141)

and if the dimension of the eigenspace is one, then a solution of (2.129) is given by

Z(x, y, n) = A 0 n+1 (0 − n2_a 11)x − n2a12y , 0 n+1 n2_a 21x− (0 − n2a22)y . (2.142)

When for a system of k first order linear or nonlinear O∆Es k functionally independent first integrals have been constructed, the solution for this system of O∆Es can be computed from

(29)

2.4. Conclusions 23

the k first integrals. As it has been shown in this chapter, the solution Z(x, y, n) of the linear functional equation corresponding to a linear system of O∆Es can be linear in x and y, and also n-dependent.

It is to be expected that the presented procedure can also be applied to the corresponding functional equation of a system of more than two linear, first order ordinary difference equa-tions, or a system of weakly perturbed nonlinear, first order ordinary difference equations. Of course, these extensions will be interesting subjects for future research in the field of perturbation theory and applied mathematics.

(30)

(31)

Published as: M. Rafei and W.T. Van Horssen – ”On asymptotic approximations of first integrals for second order difference equations”, Nonlinear Dynamics, vol. 61, no. 3, pp. 535–551, February 2010.

Chapter 3 On asymptotic approximations of first integrals

for second order difference equations

Shall I refuse my dinner because I do not fully understand the process of digestion?

O. Heaviside

Abstract

In this chapter, the concept of invariance factors for second order difference equations to obtain first integrals or invariants will be presented. It will be shown that all invariance factors have to satisfy a functional equation. (van Horssen 2007) developed a perturba-tion method for a single first order difference equaperturba-tion based on invariance factors. This perturbation method will be reviewed shortly, and will be extended to second order dif-ference equations. Also in this chapter we will construct approximations of first integrals for second order linear, and weakly nonlinear difference equations.

3.1 Introduction

F

or scientists and engineers the analysis of nonlinear dynamical systems is an impor-tant 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) in-tegrating factors can play an important role.

The fundamental concept of how to make a single first order ordinary differential equation (ODE) exact by means of integrating factors was discovered by Euler in the period 1732– 1734. Euler showed that for a first order ODE all of integrating factors have to satisfy a single, first order, linear partial differential equation. Finding an integrating factor for a given first order ODE was and still is a difficult and usually impossible task. Euler, how-ever, used special types of integrating factors obtaining (and so solving) classes of first order ODEs. In the latter part of the 19th century, Sophus Lie introduced the notion of continuous groups (currently known as Lie groups) to unify and to extend various solution methods for ODEs. Lie also showed that the existence of an integrating factor is equivalent to having a

(32)

26 3. On asymptotic approximations of first integrals

Lie group invariance or symmetry of the given ODE. To find a symmetry for a given (system of) ODE(s) is in general, however, extremely difficult and usually an impossible task. In Van Horssen’s papers ((van Horssen 1997), (van Horssen 1999a) and (van Horssen 1999b)) the fundamental concept of how to make second and higher order ODEs as well as systems of first order ODEs exact by means of integrating factors and integrating vectors has been presented. Like most methods for differential equations there is an analogous method for difference equations (see for instance (Agarwal 1992) and (Elaydi 2005)). Recently, first in-tegrals, invariants and Lie group theory for ordinary difference equations (O∆Es) obtained a lot of attention in the literature (see for instance the list of references in (van Horssen 2002a)). Also recently, the fundamental concept of invariance factors for O∆Es to obtain invariants (or first integrals) for O∆Es has been presented in (van Horssen 2002a).

It has been shown in (van Horssen 2002a) 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 ((Kuczma 1968), (Kuczma et al. 1990), (Aczel and Dhombres 1989), (Risteski and Covachev 2001), (Bowman 2002) and (Small 2007))). The aim of this chapter is to con-struct asymptotic approximations of first integrals for a system of first order O∆Es. After presenting the concepts, we will explicitly show how invariants for a linear, second order difference equation, and also for a second order weakly nonlinear difference equation (with a Van der Pol type of nonlinearity) can be constructed. To obtain highly accurate approxima-tions of invariants it also will be proposed to include the method of multiple-scales into the perturbation method based on invariance factors. 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 (van Horssen and ter Brake 2009).

The outline of this chapter is as follows. In section 2 the concept of invariance factors for a system of first order O∆Es will be given, and in section 3 some of the first results in the development of a perturbation method for a system of O∆Es based on invariance factors will be presented. In section 4 and 5 approximations of first integrals for systems of linear and weakly nonlinear O∆Es, respectively, will be constructed. Finally, in section 6 of this chapter some conclusions will be drawn and some future directions for research will be indicated.

3.2 On invariance factors for O∆Es

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 O∆Es (where k is fixed and k ∈ N), for the reference

(33)

3.2. On invariance factors for O∆Es 27

concerning this concept, we refer the reader to (van Horssen 2007). Consider

x_n+1= f (x_n, n), (3.1)

for n = 0, 1, 2, . . . , and where x_n= (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 (3.1) can be represented by

I(x_n+1, n + 1) = I(xn, n) = constant ⇔ ∆I(xn, n) = 0. (3.2)

We now try to find an invariant for (3.1). By multiplying each i−th equation, in (3.1), with a factor µi(xn+1, n + 1) = µi( f (xn, n), n + 1) for i = 1, 2, . . . , k and by adding the resulting

equations we obtain

µ (x_n+1, n + 1) · x_n+1= µ( f (x_n, n), n + 1) · f (x_n, n), (3.3) where µ = (µ1, µ2, . . . , µk)T, and µi= µi(xn, n) = µi(x1,n, x2,n, . . . , xk,n, n) for i = 1, 2, . . . , k.

In fact, when µ is an invariance vector we now obtain an exact difference equation (3.2). The relationship between I and µ follows from the equivalence of (3.2) and (3.3), yielding

(

I(x_n+1, n + 1) = µ(xn+1, n + 1) · xn+1,

I(x_n, n) = µ( f (x_n, n), n + 1) · f (x_n, n). (3.4) By reducing the index n + 1 by 1 in the first part of (3.4) I can be eliminated from (3.4), and then it follows that all invariance vectors for the system of difference equations (3.1) have to satisfy the functional equation

µ (x_n, n) · x_n= µ( f (x_n, n), n + 1) · f (x_n, n). (3.5) When an invariance vector has been determined from the functional equation (3.5) an invari-ant for (3.1) easily follows from (3.4), yielding

I(x_n, n) = µ(x_n, n) · x_n. (3.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 (van Horssen 2002a).

(34)

3.3 A perturbation method based on invariance factors for

systems of first order O∆Es

In this section some of the first results in the development of a perturbation method for difference equation based on invariance factors will be presented. The analysis in the first part of this section will be restricted to the following single, first order O∆E (see also (van Horssen 2007)):

xn+1= f (xn, n; ε), (3.7)

where ε is a small parameter. In most applications the function f has the form

f(x_n, n; ε) = f0(xn, n) + ε f1(xn, n) + ε2f2(xn, n) + · · · . (3.8)

An invariance factor µ(xn, n; ε) for (3.7) has to satisfy (3.5). Now it will also be assumed

that µ can be expanded in a power series in ε, that is,

µ (xn, n; ε) = µ0(xn, n) + ε µ1(xn, n) + ε2µ2(xn, n) + · · · . (3.9)

The expansions (3.8) and (3.9) are then substituted into (3.5), yielding {µ0(xn, n) + ε µ1(xn, n) + · · · }xn= {µ0( f0(xn, n) + ε f1(xn, n) + · · · , n + 1)

+ε µ1( f0(xn, n) + ε f1(xn, n) + · · · , n + 1) + · · · }( f0(xn, n) + ε f1(xn, n) + · · · ).

(3.10)

Now it should be observed that (for i = 0, 1, 2, . . . )

µi( f0+ ε f1+ ε2f2+ · · · , n + 1) = µi( f0, n + 1) + ε µi0( f0, n + 1) f1

+ε2[(µ_i0( f0, n + 1) f2+1₂µ_i00( f0, n + 1) f12]

+ε3_{(· · · ) + · · · ,}

(3.11)

where the prime0denotes differentiation with respect to the first argument of µi. Then, by

using (3.11) and by taking apart in (3.10) the O(1)-terms, the O(ε)-terms, the O(ε2_)-terms,

and so on, it follows that the O(1)-problem becomes

µ0(xn, n)xn= µ0( f0(xn, n), n + 1) f0(xn, n), (3.12)

that the O(ε)-problem becomes

µ1(xn, n)xn= µ1( f0(xn, n), n + 1) f0(xn, n) + µ0( f0(xn, n), n + 1) f1(xn, n)

+µ₀0( f₀(x_n, n), n + 1) f1(xn, n) f0(xn, n),

(35)

3.3. A perturbation method based on invariance factors 29

and that the O(ε2)-problem becomes

µ2(xn, n)xn= µ2( f0(xn, n), n + 1) f0(xn, n) + µ0( f0(xn, n), n + 1) f2(xn, n) +µ₀0( f0(xn, n), n + 1)[ f12(xn, n) + f2(xn, n) f0(xn, n)] +1 2µ 00 0( f0(xn, n), n + 1) f12(xn, n) f0(xn, n) +µ1( f0(xn, n), n + 1) f1(xn, n) +µ₁0( f0(xn, n), n + 1) f0(xn, n) f1(xn, n). (3.14)

For a given function f (xn, n; ε) the O(1)-problem, the O(ε)-problem, and so on now have to

be solved. For some examples the reader is refered to (van Horssen 2007). For a system of two first order O∆Es,

xn+1= f1(xn, yn, n; ε) = f1,0(xn, yn, n) + ε f1,1(xn, yn, n) + ε2f1,2(xn, yn, n) + · · · ,

yn+1= f2(xn, yn, n; ε) = f2,0(xn, yn, n) + ε f2,1(xn, yn, n) + ε2f2,2(xn, yn, n) + · · · ,

(3.15) the same procedure can be followed. An invariance vector µ(xn, yn, n; ε) for (3.15) has to

satisfy (3.5). Now it will also be assumed that µ = (µ1, µ2)T can be expanded in a power

series in ε, that is,

µ1(xn, yn, n; ε) = µ1,0(xn, yn, n) + ε µ1,1(xn, yn, n) + ε2µ1,2(xn, yn, n) + · · · ,

µ2(xn, yn, n; ε) = µ2,0(xn, yn, n) + ε µ2,1(xn, yn, n) + ε2µ2,2(xn, yn, n) + · · · .

(3.16)

It is obvious from (3.5) that the functional equation for a system of two first order O∆Es becomes

µ1(xn, yn, n; ε)xn+ µ2(xn, yn, n; ε)yn= µ1(xn+1, yn+1, n + 1; ε)xn+1

+µ2(xn+1, yn+1, n + 1; ε)yn+1.

(3.17)

After substitution of (3.15) and of (3.16) into (3.17) we will have µi, j( f1,0+ ε f1,1+ ε2f1,2+ · · · , f2,0+ ε f2,1+ ε2f2,2+ · · · , n + 1) = µi, j( f1,0, f2,0, n + 1) + ε[D1(µi, j)( f1,0, f2,0, n + 1) f1,1 +D2(µi, j)( f1,0, f2,0, n + 1) f2,1] + ε2[D2(µi, j)( f1,0, f2,0, n + 1) f2,2 +D1(µi, j)( f1,0, f2,0, n + 1) f1,2+ D1,2(µi, j)( f1,0, f2,0, n + 1) f1,1f2,1 +1₂D1,1(µi, j)( f1,0, f2,0, n + 1) f1,12 +12D2,2(µi, j)( f1,0, f2,0, n + 1) f 2 2,1] + · · · , (3.18) where the subscripts 1 and 2 of D are differentiation with respect to the first and second arguments of µi, j, respectively, for i = 1, 2 and j = 0, 1, 2, . . . ; Then, it follows that the

O(1)-problem becomes

µ1,0(xn, yn, n)xn+ µ2,0(xn, yn, n)yn=

µ1,0( f1,0(xn, yn, n), f2,0(xn, yn, n), n + 1) f1,0(xn, yn, n)

+µ2,0( f1,0(xn, yn, n), f2,0(xn, yn, n), n + 1) f2,0(xn, yn, n),

(36)

and the O(ε)-problem becomes

µ1,1(xn, yn, n)xn+ µ2,1(xn, yn, n)yn= f1,1(xn, yn, n)µ1,0( f1,0(xn, yn, n), f2,0(xn, yn, n), n + 1) + f_1,0(x_n, y_n, n) f_1,1(x_n, y_n, n)D₁(µ1,0)( f1,0(xn, yn, n), f2,0(xn, yn, n), n + 1) + f1,0(xn, yn, n) f2,1(xn, yn, n)D2(µ1,0)( f1,0(xn, yn, n), f2,0(xn, yn, n), n + 1) + f1,0(xn, yn, n)µ1,1( f1,0(xn, yn, n), f2,0(xn, yn, n), n + 1) + f2,1(xn, yn, n)µ2,0( f1,0(xn, yn, n), f2,0(xn, yn, n), n + 1) + f1,1(xn, yn, n) f2,0(xn, yn, n)D1(µ2,0)( f1,0(xn, yn, n), f2,0(xn, yn, n), n + 1) + f2,0(xn, yn, n) f2,1(xn, yn, n)D2(µ2,0)( f1,0(xn, yn, n), f2,0(xn, yn, n), n + 1) + f2,0(xn, yn, n)µ2,1( f1,0(xn, yn, n), f2,0(xn, yn, n), n + 1). (3.20) In the next two sections it will be shown how the perturbation method can be applied to a system of two first order O∆Es. In section 4 a system of two linear equations will be considered, and in section 5 a system of two weakly nonlinear O∆Es will be studied.

3.4 A system of two linear, first order O∆Es

To show how the method can be applied in practice we will consider the following differential equation:

¨

x+ µ ˙x + x = 0, (3.21)

where x = x(t), and where µ is a small, positive damping parameter. This leads to the following difference equation, when a central difference scheme is used to discretize the equation (3.21): xn+1− 2xn+ xn−1 h2 + µ xn+1− xn−1 2h + xn= 0, (3.22) or equivalently, xn+1+ (h2− 2)xn+ xn−1= µ h 2 (xn−1− xn+1). (3.23)

In fact (3.23) can be considered as a central finite difference approximation of (3.21). In this case xnis an approximation of x(tn) at tn= nh, where h is the discretisation time step. Now

let ( h2− 2 = −2 cos(θ0), µ h 2 = ε. (3.24) Substitution of these new constants θ0and ε from (3.24) into (3.23) and shifting the index

by 1, yields

(37)

3.4. A system of two linear, first order O∆Es 31

Now (3.25) is transformed into a system of first order difference equations, xn+1= yn, yn+1= 2yncos(θ0) − xn+ ε(xn− yn+1), (3.26) or equivalently, xn+1= yn, yn+1=_1+ε1 (2yncos(θ0) − xn+ εxn), (3.27) and after expanding (1 + ε)−1into 1 − ε + O(ε2) for ε small, it follows that (3.27) becomes

xn+1= yn,

yn+1= 2yncos(θ0) − xn+ 2ε(xn− yncos(θ0)) + O(ε2).

(3.28)

Then we substitute (3.28) into (3.15) and (3.19). So, the O(1)-problem will be µ1,0(xn, yn, n)xn+ µ2,0(xn, yn, n)yn= µ1,0(yn, 2yncos(θ0) − xn, n + 1)yn

+µ2,0(yn, 2yncos(θ0) − xn, n + 1)(2yncos(θ0) − xn).

(3.29) Now we define a new function

Z0(xn, yn, n) = µ1,0(xn, yn, n)xn+ µ2,0(xn, yn, n)yn, (3.30)

and by using (3.30), equation (3.29) becomes

Z0(xn, yn, n) = Z0(yn, 2yncos(θ0) − xn, n + 1). (3.31)

According to the results as obtained in Appendix A, the general solution of (3.31) is given by Z0(xn, yn, n) = A0 x2_n− 2x_nyncos(θ0) + y2n, n + 1 θ0 arctan ynsin(θ0) xn− yncos(θ0) , (3.32)

where A0is an arbitary function. For more information regarding how to solve the functional

equation (3.31), we refer the reader to Appendix A. According to (3.30) and (3.32), we now have

µ1,0(xn, yn, n)xn+ µ2,0(xn, yn, n)yn= A0 x2_n− 2xnyncos(θ0) + y2n, n +θ10arctan ynsin(θ0) xn−yncos(θ0) . (3.33) Since we are dealing with a second order difference equation (3.25), we need two function-ally independent approximations of the invariants. The function A0is still arbitrary, and will

now be chosen to be as simple as possible to obtain relatively simple approximations of the invariants.

(38)

3.4.1 Case 1

First we take A0≡ x2n− 2xnyncos(θ0) + y2n, µ1,0(xn, yn, n) = 0. (3.34) So, from (3.30) and (3.33), we can find

µ2,0=

1 yn

(x2_n− 2xnyncos(θ0) + y2n). (3.35)

By substituting the values of µ1,0and µ2,0into (3.20), the O(ε)-problem then becomes

µ1,1(xn, yn, n)xn+ µ2,1(xn, yn, n)yn= µ1,1(yn, 2yncos(θ0) − xn, n + 1)yn

+µ2,1(yn, 2yncos(θ0) − xn, n + 1)(2yncos(θ0) − xn) − 4(xn− yncos(θ0))2,

(3.36) or equivalently,

Z1(xn, yn, n) = Z1(yn, 2yncos(θ0) − xn, n + 1) − 4(xn− yncos(θ0))2, (3.37)

where Z1(xn, yn, n) = µ1,1(xn, yn, n)xn+ µ2,1(xn, yn, n)yn. When we put µ2,0= 0 instead of

µ1,0= 0, then µ1,0=_x1_n(xn2− 2xnyncos(θ0) + y2n). The interesting point is that the

O(ε)-problem then also leads to (3.37). According to (3.31) and (3.32), the general solution Z1,h

of the equation corresponding to the homogeneous equation (3.37) is

Z1,h(xn, yn, n) = A1 x2_n− 2xnyncos(θ0) + y2n, n + 1 θ0 arctan _y nsin(θ0) xn− yncos(θ0) . (3.38)

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

Z1,p(xn, yn, n) = a(n)x2n+ b(n)xnyn+ c(n)y2n. (3.39)

By substituting (3.39) into (3.37), we obtain

a(n)x2_n+ b(n)xnyn+ c(n)y2n= a(n + 1)y2n+ b(n + 1)yn(2yncos(θ0) − xn)

+c(n + 1)(2yncos(θ0) − xn)2− 4(xn− yncos(θ0))2,

{a(n + 1) − 4 cos2_(θ

0)a(n) − 2 cos(θ0)b(n) − c(n) − 4 cos2(θ0)}x2n

+{b(n + 1) + 4 cos(θ0)a(n) + b(n) + 8 cos(θ0)}xnyn+ {c(n + 1) − a(n) − 4}y2n= 0.

(39)

3.4. A system of two linear, first order O∆Es 33

Obviously the coefficients of x2_n, xnyn, and y2nhave to be zero, yielding

  

a(n + 1) = 4 cos2_(θ

0)a(n) + 2 cos(θ0)b(n) + c(n) + 4 cos2(θ0),

b(n + 1) = −4 cos(θ0)a(n) − b(n) − 8 cos(θ0),

c(n + 1) = a(n) + 4.

(3.42)

The eigenvalues of the homogeneous system related to (3.42) are 1, and cos(2θ0)±i sin(2θ0),

and a particular solution of (3.42) is given by 

 

a(n) = S1{(2n − 1) sin(2θ0) + sin(2nθ0) + sin(2(n + 1)θ0)},

b(n) = −4S1cos(θ0){n sin(2θ0) + sin(2nθ0)},

c(n) = S₁{(2n + 1) sin(2θ0) + sin(2nθ0) + sin(2(n − 1)θ0)},

(3.43)

where S1=_sin(2θ1 ₀₎. Therefore, according to (3.39) and (3.43), a particular solution of (3.37)

is given by

Z1,p(xn, yn, n) = S1{[(2n − 1) sin(2θ0) + sin(2nθ0) + sin(2(n + 1)θ0)]x2n

−4 cos(θ0)[n sin(2θ0) + sin(2nθ0)]xnyn+ [(2n + 1) sin(2θ0)

+ sin(2nθ0) + sin(2(n − 1)θ0)]y2n}.

(3.44) So, from (3.36)-(3.38), we have

µ1,1(xn, yn, n)xn+ µ2,1(xn, yn, n)yn= Z1,h(xn, yn, n) + Z1,p(xn, yn, n), (3.45)

where Z1,h(xn, yn, n) and Z1,p(xn, yn, n) are defined in (3.38) and (3.44), respectively. We will

now choose the invariant as simple as possible, that is, A1≡ 0 in (3.45), and so for this case

we have

µ1,1(xn, yn, n)xn+ µ2,1(xn, yn, n)yn= S1{[(2n − 1) sin(2θ0) + sin(2nθ0)

+ sin(2(n + 1)θ0)]x2n− 4 cos(θ0)[n sin(2θ0) + sin(2nθ0)]xnyn+ [(2n + 1) sin(2θ0)

+ sin(2nθ0) + sin(2(n − 1)θ0)]y2n}.

(3.46) An approximation IA(xn, yn, n) of an invariant I(xn, yn, n) = constant for (3.28) is now given

by (see also (3.6)):

IA(xn, yn, n) = µ1,0(xn, yn, n)xn+ µ2,0(xn, yn, n)yn

+ε{µ1,1(xn, yn, n)xn+ µ2,1(xn, yn, n)yn}.

(3.47)

Then by using (3.33), (3.34) and (3.46), we conclude that

IA(xn, yn, n) = x2n− 2xnyncos(θ0) + y2n+ εS1{[(2n − 1) sin(2θ0) + sin(2nθ0)

+ sin(2(n + 1)θ0)]x2n− 4 cos(θ0)[n sin(2θ0) + sin(2nθ0)]xnyn+ [(2n + 1) sin(2θ0)

+ sin(2nθ0) + sin(2(n − 1)θ0)]y2n},

(40)

34 3. On asymptotic approximations of first integrals and, IA(xn+1, yn+1, n + 1) = IA(xn, yn, n) + ε2R(xn, yn, n; ε), (3.49) where ε2R(xn, yn, n; ε) = −2 f1( f2− f2,0− ε f2,1) cos(θ0) + ( f22− f2,02 − 2ε f2,0f2,1) +ε{b(n + 1) f1( f2− f2,0) + c(n + 1)( f22− f2,02 )}, (3.50)

where f1, f2, f2,0 and f2,1 are given by (3.15) and (3.28). Since the system (3.28) is

n-independent, that is, f1 and f2are n-independent, and b(n) and c(n) are of O(n), then R

is unbounded in n and of O(n). From (3.49) it follows that

IA(xn, yn, n) = IA(x0, y0, 0) + ε2 n−1

∑

i=0

R(xi, yi, i; ε). (3.51)

From (3.48)-(3.51), it can be shown that

IA(xn, yn, n) = x02+ y20− 2x0y0cos(θ0) + O(ε2n2), (3.52)

it then follows that (

IA(xn, yn, n) = constant + O(ε2) for n= O(1),

IA(xn, yn, n) = constant + O(ε) for n= O(√1_ε). (3.53)

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 straight-forward way as follows from the next subsection.

3.4.2 Case 2

Now we take ( A0≡ n +_θ1₀arctan _y nsin(θ0) xn−yncos(θ0) , µ1,0(xn, yn, n) = 0. (3.54) So, from (3.30) and (3.33), we can find

µ2,0= 1 yn n+ 1 θ0 arctan _y nsin(θ0) xn− yncos(θ0) . (3.55)

By substituting the values of µ1,0and µ2,0into (3.20), the O(ε)-problem then becomes

µ1,1(xn, yn, n)xn+ µ2,1(xn, yn, n)yn= µ1,1(yn, 2yncos(θ0) − xn, n + 1)yn

+µ2,1(yn, 2yncos(θ0) − xn, n + 1)(2yncos(θ0) − xn) +2y_θn(xn−yncos(θ0)) sin(θ0)

0(x2n−2xnyncos(θ0)+y2n),

(41)

3.4. A system of two linear, first order O∆Es 35 or equivalently, Z1(xn, yn, n) = Z1(yn, 2yncos(θ0) − xn, n + 1) + 2yn(xn− yncos(θ0)) sin(θ0) θ0(x2n− 2xnyncos(θ0) + y2n) , (3.57)

where Z1(xn, yn, n) = µ1,1(xn, yn, n)xn+ µ2,1(xn, yn, n)yn. When we put µ2,0= 0 instead of

µ1,0= 0, then µ1,0= _x1_n n+ 1 θ0arctan _y nsin(θ0) xn−yncos(θ0)

. The interesting point is that again the O(ε)-problem leads to (3.57). According to (3.31) and (3.32), the general solution Z1,h

of the equation corresponding to the homogeneous equation (3.57) is

Z1,h(xn, yn, n) = A1 x2_n− 2xnyncos(θ0) + y2n, n + 1 θ0 arctan ynsin(θ0) xn− yncos(θ0) . (3.58)

Now, we are going to construct a particular solution of (3.57). Since 2 sin(θ0)

θ0(x2n−2xnyncos(θ0)+y2n)

satisfies the homogeneous equation related to (3.57), we consider a particular solution of (3.57) in the form Z1,p(xn, yn, n) = 2 sin(θ0) θ0(x2n− 2xnyncos(θ0) + y2n) {a(n)x2 n+ b(n)xnyn+ c(n)y2n}. (3.59)

By substituting (3.76) into (3.57), we obtain a(n)x2

n+ b(n)xnyn+ c(n)y2n= a(n + 1)yn2+ b(n + 1)yn(2yncos(θ0) − xn)

+c(n + 1)(2yncos(θ0) − xn)2+ yn(xn− yncos(θ0)),

{a(n + 1) − 4 cos2_(θ

0)a(n) − 2 cos(θ0)b(n) − c(n) + cos(θ0)}x2n

+{b(n + 1) + 4 cos(θ0)a(n) + b(n) − 1}xnyn+ {c(n + 1) − a(n)}y2n= 0.

(3.61)

Obviously the coefficients of x2

n, xnyn, and y2nhave to be zero, yielding

  

a(n + 1) = 4 cos2_(θ

0)a(n) + 2 cos(θ0)b(n) + c(n) − cos(θ0),

b(n + 1) = −4 cos(θ0)a(n) − b(n) + 1,

c(n + 1) = a(n).

(3.62)

The eigenvalues of the homogeneous system related to (3.62) are 1, and cos(2θ0)±i sin(2θ0),

and a particular solution of (3.62) is given by 

 

a(n) = S2{cos((2n + 1)θ0) − cos(θ0)},

b(n) = 4S2sin2(nθ0),

c(n) = S₂{cos((2n − 1)θ0) − cos(θ0)},

(42)

where S2=_{4 sin}12_(θ

0). Therefore, according to (3.76) and (3.63), a particular solution of (3.57)

is given by

Z1,p(xn, yn, n) =_θ 2S2sin(θ0)

0(x2n−2xnyncos(θ0)+y2n)

{x2

n[cos((2n + 1)θ0) − cos(θ0)]

+4xnynsin2(nθ0) + y2n[cos((2n − 1)θ0) − cos(θ0)]}.

(3.64)

So, from (3.56)-(3.58), we have

µ1,1(xn, yn, n)xn+ µ2,1(xn, yn, n)yn= A1 x2_n− 2xnyncos(θ0) + y2n, n +θ10arctan _y nsin(θ0) xn−yncos(θ0) + 2S2sin(θ0) θ0(x2n−2xnyncos(θ0)+y2n){x 2 n[cos((2n + 1)θ0) − cos(θ0)]

(3.65) The function A1 is still arbitrary but will be chosen to be as simple as possible: A1≡ 0,

yielding µ1,1(xn, yn, n)xn+ µ2,1(xn, yn, n)yn= 2S2sin(θ0) θ0(x2n−2xnyncos(θ0)+y2n){x 2 n[cos((2n + 1)θ0) − cos(θ0)]

(3.66) As for the case 1, we can construct an exact difference equation of (3.28) up to O(ε2), in an almost similar way, yielding

IA(xn+1, yn+1, n + 1) = IA(xn, yn, n) + ε2R(xn, yn, n; ε), (3.67) where IA(xn, yn, n; ε) = n +_θ1₀arctan _y nsin(θ0) xn−yncos(θ0) + 2εS2sin(θ0) θ0(x2n−2xnyncos(θ0)+y2n) {x2 n[cos((2n + 1)θ0) − cos(θ0)]

+4xnynsin2(nθ0) + y2n[cos((2n − 1)θ0) − cos(θ0)]},

(3.68) and, ε2R(xn, yn, n; ε) =_θ1 n arctan f2sin(θ0) f1− f2cos(θ0) − arctan f2,0sin(θ0) f1− f2,0cos(θ0) − f1f2,1sin(θ0) f2 1−2 f1f2,0cos(θ0)+ f2,02 ε o +2ε sin(θ0) θ0 n_{a(n+1) f}2 1+b(n+1) f1f2+c(n+1) f22 f2 1−2 f1f2cos(θ0)+ f22 −a(n+1) f 2 1+b(n+1) f1f2,0+c(n+1) f2,02 f2 1−2 f1f2,0cos(θ0)+ f2,02 o , (3.69) where f1, f2, f2,0 and f2,1 are given by (3.15) and (3.28). Since the system (3.28) is

n-independent, and a(n), b(n) and c(n) are bounded, then R is bounded. From (3.67), it follows that IA(xn, yn, n) = IA(x0, y0, 0) + ε2 n−1

∑

i=0 R(xi, yi, i; ε), (3.70)

(43)

3.5. A weakly nonlinear, regularly perturbed system of two O∆Es 37 where IA(x0, y0, 0) = _θ1 0arctan _y 0sin(θ0) x0−y0cos(θ0)

= constant. Finally, it follows from (3.67)-(3.70) that

IA(xn, yn, n) = constant + O(ε2n), (3.71)

it then follows that

IA(xn, yn, n) = constant + O(ε2) for n= O(1),

IA(xn, yn, n) = constant + O(ε) for n= O(1_ε).

(3.72)

The exact solution of the system of difference equations (3.27) is

xn= c1cos(nθ0) + c2sin(nθ0) − nε(c1cos(nθ0) + c2sin(nθ0)) + O(ε2n2),

yn= xn+1,

(3.73)

and therefore, it can readily be verified that IA(xn, yn, n) as given by (3.48) (or by (3.68))

satisfies IA(xn, yn, n) = constant + O(ε2n2) (or IA(xn, yn, n) = constant + O(ε2n)).

3.5 A weakly nonlinear, regularly perturbed system of two

O∆Es

In this section approximations of first integrals for a second order, weakly nonlinear, regu-larly perturbed O∆E with a Van der Pol type of nonlinearity will be considered. The Van der Pol equation ((Mickens 1996) and (Nayfeh 1973)) corresponds to a nonlinear oscillatory system that has both input and output sources of energy. This equation is given by

¨

x+ ω02x= µ(1 − x2) ˙x, (3.74)

where x = x(t), µ is a non-negative small parameter, and where ω0is a bounded constant.

This leads to the following difference equation, when a central difference scheme (Potts 1983) is used to discretize the equation (3.74):

xn+1− 2xn+ xn−1 h2 + ω 2 0xn= µ(1 − x2n) xn+1− xn−1 2h , (3.75) or equivalently, xn+1+ (ω02h2− 2)xn+ xn−1=µ h 2 (1 − x 2 n)(xn+1− xn−1). (3.76)

(44)

In fact (3.76) can be considered as a central finite difference approximation of (3.74). In this case xnis an approximation of x(tn) at tn= nh, where h is the discretisation time step. Now

let ( h2ω₀2− 2 = −2 cos(θ0), µ h 2 = ε. (3.77) Substitution of these new constants θ0and ε from (3.77) into (3.76) and shifting the index

by 1, yields

xn+2− 2 cos(θ0)xn+1+ xn= ε(1 − x2n+1)(xn+2− xn), (3.78)

where ε is a small parameter, that is, 0 < ε 1, and where θ0is constant (which is related

to the stepsize in making the continuous Van der Pol equation discrete). Now (3.78) is transformed into a system of first order difference equations,

xn+1= yn, yn+1= 2yncos(θ0) − xn+ ε(yn+1− xn)(1 − y2n), (3.79) or equivalently, ( xn+1= yn, yn+1=_1−ε(1−y1 2 n)(2yncos(θ0) − xn− εxn(1 − y 2 n)), (3.80)

and after expanding (1 − ε(1 − y2n))−1into 1 + ε(1 − y2n) + O(ε2), it follows that (3.80)

be-comes

xn+1= yn,

yn+1= 2yncos(θ0) − xn+ 2ε(xn− yncos(θ0))(y2n− 1) + O(ε2).

(3.81)

Then we substitute (3.81) into (3.15) and (3.19). So, the O(1)-problem will be µ1,0(xn, yn, n)xn+ µ2,0(xn, yn, n)yn= µ1,0(yn, 2yncos(θ0) − xn, n + 1)yn

+µ2,0(yn, 2yncos(θ0) − xn, n + 1)(2yncos(θ0) − xn).

(3.82) Now we define a new function

Z0(xn, yn, n) = µ1,0(xn, yn, n)xn+ µ2,0(xn, yn, n)yn, (3.83)

and by using (3.83), equation (3.82) becomes

Z0(xn, yn, n) = Z0(yn, 2yncos(θ0) − xn, n + 1). (3.84)

The obtained functional equation is the same as the functional equation (3.31). Therefore, the general solution of (3.84) is given by

Z0(xn, yn, n) = A0 x2_n− 2xnyncos(θ0) + y2n, n + 1 θ0 arctan _y nsin(θ0) xn− yncos(θ0) , (3.85)

On asymptotics for difference equations