- Email: [email protected]
Invariant tori of dissipatively perturbed Hamiltonian systems under symplectic discretization
~
APPLIED
"~q~ NUMERICAL MATHEMATICS ELSEVIER
Applied Numerical Mathematics 29 (1999) 57-71
Invariant toil of dissipatively perturbed Hamiltonian systems under symplectic discretization Ernst Hairer a, 1, Christian Lubich b,, a Section de Mathdmatiques, Universitd de GenOve, CH-1211 GenOve 24, Switzerland b Mathematisches Institut, Universitiit Tiibingen, Aufder Morgenstelle 10, D-72076 Tiibingen, Germany
Abstract In a recent paper, Stoffer showed that, under a very weak restriction on the step size, weakly attractive invariant tori of dissipative perturbations of integrable Hamiltonian systems persist under symplectic numerical discretizations. Stoffer's proof works directly with the discrete scheme. Here, we show how such a result, together with approximation estimates, can be obtained by combining Hamiltonian perturbation theory and backward error analysis of numerical integrators. In addition, we extend Stoffer's result to dissipative perturbations of nonintegrable Hamiltonian systems in the neighborhood of a KAM toms. © 1999 Elsevier Science B.V. and IMACS. All rights reserved.
Keywords: Weakly damped nonlinear oscillations; Symplectic numerical integrator; Hamiltonian perturbation theory; Backward error analysis; Invariant manifold; Long-time approximation
I. Introduction
We study the long-time behavior of numerical discretizations of perturbed Hamiltonian systems OH --(p,q) +sP(p,q), Oq OH gl = + ~ p (P, q) + sQ(p, q),
-
0 < s < < 1.
(1)
A classical example is Van der Pol's equation, where the unperturbed system is the harmonic oscillator, and where the perturbation creates a limit cycle; see, for example, [3], where this equation and many others are studied by the method of averaging. In higher dimensions, systems (1) often possess weakly attractive invariant toil, as shown in [8, Section 12] for perturbations of systems of non-resonant harmonic oscillators, and in [4] in a far-reaching framework. More recently, Stoffer [12] studied the * Corresponding author. E-mail: [email protected]. I E-mail: [email protected]. 0168-9274/99/$19.00 © 1999 Elsevier Science B.V. and IMACS. All rights reserved. PII: S0168-9274(98)00029-4
E. Hairer, C. Lubich /Applied Numerical Mathematics 29 (1999) 57-71
58
existence and properties of attractive invariant tori for dissipative perturbations of integrable Hamiltonian systems. Here, we assume, more generally, that the unperturbed Hamiltonian system has an invariant Lagrangian toms with linear flow, such as a KAM toms, and that the perturbation is dissipative, ensuring that (1) has a weakly attractive invariant toms. The precise assumptions and proof of existence are given in Section 3. The proof relies on averaging transformations of [11,12], and on an invariant manifold theorem of [7,10]. Does such a weakly attractive invariant toms persist under numerical discretization? For dissipative perturbations of integrable systems, Stoffer [12] shows that symplectic discretization methods preserve the invariant toms under the very weak restriction that the step size is logarithmically small in the perturbation parameter e. Stoffer's proof works directly with the discrete scheme, rebuilding the familiar constructions for continuous problems in the more complicated discrete situation, requiring great technical skill and effort. In a different approach, for the unperturbed case e = 0 of a Hamiltonian system in the neighborhood of a KAM toms, in [6] the authors show that a symplectic discretization has a nearby toms which inherits from the KAM toms the property of being 'very sticky' (a notion of Perry and Wiggins [11]): every numerical solution starting near the toms remains within twice the initial distance from the toms over a time span exponentially long in the distance and the step size. This result uses the backward error analysis of symplectic integrators, which interprets a step of the discrete numerical method as the almost-exact flow of a modified Hamiltonian system, up to an error exponentially small in the step size. This estimate is proved independently in [2,6]. Here, we show how the existence and closeness of an attractive invariant toms of a symplectic discretization of (1) can be derived using the backward error analysis result. We extend Stoffer's discrete result to situations where the unperturbed Hamiltonian is not integrable. In Section 2, for illustration, we begin with a numerical experiment with Van der Pol's equation. In Section 3, we formulate the precise assumptions about the Hamiltonian system and its perturbation, and we show the existence of a weakly attractive invariant toms. We review the required results from backward error analysis in Section 4. The existence and approximation properties of a weakly attractive invariant toms of the symplectic discretization are studied in Section 5. Finally, in Appendix A, we give two useful coordinate transformations of Perry and Wiggins [11] and of Stoffer [12], needed in the proofs.
2. Formal analysis and numerical experiment We assume that, after a suitable coordinate transformation, the system (1) is
OF h = --~-~o(a, qg) + eA(a, ~p), OF (o = w + - 7 - ( a , ~o) + e ~ ( a , qg),
(2)
06l
where a 6 •d, q9 6 ]~d, O) is a constant vector and the functions F, A, • are 2zr-periodic in the variables ~0y. For motivation, we give a non-rigorous explanation without stating the technical assumptions and without going into detail. The rigorous analysis will be presented in subsequent sections.
E. Hairer, C. Lubich /Applied Numerical Mathematics 29 (1999) 57-71
59
I. Averaging principle. Replace the first equation in (2) by the so-called averaged system _
h = eA(a),
--
1
/"
A(a) - (2zr)~ J A(a, ~o) d~0,
(3)
-g'd
with ~d = ~d/27rzd the standard toms. (Due to the periodicity of F(a, ~o) the average of OF/O~o vanishes.) The claim is that, under suitable assumptions, the averaged system gives a good approximation to the original one. Here, we consider the situation where (3) has an asymptotically stable fixed point a*, i.e., A(a*) = 0 and the eigenvalues of the Jacobian A'(a*) each have a negative real part. Hence, (3) together with the second equation of (2) has the invariant toms a = a*. Blind use of the "averaging principle" then hints at the existence of an attractive invariant toms of (2). II. B a c k w a r d error analysis. The idea is to interpret formally the numerical solution (a,, ~0n) of a pth order one-step method applied to (2) as the exact solution (a(tn), ~(tn)) at t~ = nh of the so-called modified equations a=
(4) =
+
where
OF G(a, ~o) -= - - - ( a ,
~o) + hP G p+l (a, ~o) + h p+I G p+z(a, tp) + . . . ,
A(a, qg) = A(a, qg) d- hP Ap+l (a, qg) -k- h p+I Ap+z(a, ~o) + . . . , and similarly f o r / : and qS. This allows us to study properties of the numerical solution by examining the modified equations. Assuming that (2) has an attractive invariant toms, what can be said about its numerical approximation? Applying the averaging principle to (4) yields = h~G~÷,(~) + . . . + e A ( ~ )
+eh~A,,÷,(~) + - . - .
(5)
In general, the average Gp+l (a*) does not vanish, and (5) can have an asymptotically stable fixed point (close to a*) only if h p is small compared to e. The distance of the fixed points for Eqs. (3) and (5) is O(hP/e). In contrast, applying a symplectic integrator, i.e., a method for which (an, ~0n) ~ (an+l, ~0n+l) is a symplectic transformation in the Hamiltonian case e = 0, the e-independent terms of the modified equation are Hamiltonian, i.e.,
G j(a, ~o) --
aFj
--
(a, qg).
The corresponding averages vanish identically, so that Gj (a) = 0 for all j. Consequently, Eq. (5) has an asymptotically stable fixed point if h p is small compared to 1, and its distance to a* is O(hP). Symplectic integrators not only show the correct qualitative behavior for much larger step sizes, but the distance to the exact invariant torus is smaller by a factor e.
E. Hairer, C. Lubich /Applied Numerical Mathematics 29 (1999) 57-71
60
This explanation is non-rigorous, because the averaging principle requires further assumptions, and the series in the modified equations do not converge.
Numerical experiment. Consider Van der Pol's equation in the form [~=-q + e(1-q2)P'
(6)
4=p, with small e, which is a perturbation of the harmonic oscillator. A symplectic change to polar coordinates (q = v~asin99, p = x/2acos99) puts the system in the form (2) with co = 1, F(a, 99) = 0, A(a, 99) = 2a cos 2 99(1 - 2a sin 2 99) and qb(a, 99) -- cos 99sin 99(1 - 2a sin 2 qg). The angular average of A (a, 99) is A(a) = a(1 - a/2), so that with a* = 2 we have A(a*) = 0, A'(a*) < O. We solve (6) numerically with two initial values, (P0, q0) = (0, 1.4) and (P0, q0) = (0, 2.6), and with three methods: the non-symplectic explicit and implicit Euler methods, and the symplectic Euler method; all have order 1. The results are displayed in Fig. 1. For large step sizes (compared to the perturbation parameter e), the non-symplectic methods give a completely wrong numerical solution, whereas the symplectic method is qualitatively correct. For smaller stepsizes, the non-symplectic numerical methods also show a limit cycle. We explain these observations by the averaging principle and a backward error analysis argument. For a differential equation :~ = f (y) + eg (y), the explicit Euler discretization is Yn+l = Yn + h ( f (yn) + eg (y~)). Here, the modified differential equation is
y= f(y) +eg(y)-~
h
f ' ( y ) f ( y ) +O(h 2 +eh).
(7)
In the above coordinates, we get a system (4) with G2(a, 99) = a. Therefore, the right-hand side of (5) is approximately ha + ea(1 - a/2), which has a = 2 + 2h/e as positive root. Hence, the limit cycle of the numerical solution of the explicit Euler method has approximate radius 2 + 2h/e (Fig. l) which is far from the correct value. The implicit Euler discretization is adjoint to the explicit Euler method. Therefore, its modified differential equation is (7) with h replaced by - h . In this case, the radius of the limit cycle is approximately 2 - 2h/e (for h < e), which again agrees very well with Fig. 1. For the symplectic Euler method
h / OH ) Pn+l = Pn + ~--~-a (Pn+l, qn) + eP(Pn+l, qn) , OH (Pn+l,qn) + eQ(Pn+l, qn)), qn+l =qn + h /~-~p the modified differential equations for Van der Pol's equation are p=--q'We(1 _~2)~+
h 2~+O(h
2_+_eh),
h
q = p - ~ q +O(h2-F eh). Since the method is symplectic, all e-independent terms in (5) vanish. Therefore, the radius of the limit cycle is of size 2 + O(h); see Fig. 1.
E. Hairer, C. Lubich /Applied Numerical Mathematics 29 (1999) 57-71 I
61
i
implicit Euler, h = 0.18
symplectic Euler, h = 0.18
explicit Euler, h = 0.18
%F i
i
implicit Euler, h = 0.06
symplectic Euler, h = 0.06
implicit Eule'r, h = 0.02
V
i
llll
explicit Euler, h --- 0.06
M
Fig. 1. Numerical experiments with Van der Pol's equation, e = 0.05.
3. Dissipatively perturbed Hamiltonian systems with a weakly attractive invariant torus 3.1. Assumptions about the unperturbed Hamiltonian system We assume there is a real-analytic symplectic coordinate transform tp : (a, ~o) ~ (p, q) which takes the Hamiltonian H0 = H o tp into the form H0(a, ~o) = w , a + F(a, ~o).
(8)
Here, o9 6 •d is a constant vector of frequencies and o9. a :- OgTa. For positive constants 7, v, we assume w satisfies the strong non-resonance condition Iog " kl >1 71kl-",
k ~ Zd\{0}.
(9)
E. Hairer, C. Lubich /Applied Numerical Mathematics 29 (1999) 57-71
62
The function F is assumed 2zr-periodic in the variables ~0j and real-analytic in the complex neighbourhood
D(r, s) = { (a, ~o) e C d x (ca/2zrza): laj[ < r, IIm~ojl < s for all j} of the toms {a = 0} c ~d x "IFd. Further, we assume that with a constant M and for all (a, ~o) e D(r, s),
IF(a, ~0)1 <~M [a[2.
(10)
The equations of motion of the unperturbed system with Hamiltonian H0 are & = O(la12),
~b= o~ + O(lal).
Hence, the toms a = 0 is invariant under the flow of this Hamiltonian system, and the flow is linear on this toms. Note, (8) with (10) is precisely the form delivered by the KAM theorem, in the neighborhood of a KAM toms (see, for example, [1]). Condition (9) is also required by the KAM theorem. Hamiltonian systems satisfying the very assumptions (8)-(10) are studied in [11].
3.2. Assumptions about the perturbation In the coordinates (a, ~o), the equations for the perturbed system become aF
= --g-z(a, ~o) + eA(a, ~0), oW
~F
(11)
(o = w + ff-aa (a, ~o) + cO(a, ~o). We assume that A and • (which may also depend on e) are real-analytic on D(r, s) such that, for all (a, ¢p) e D(r, s),
[A(a, ~0)[ ~< MA,
]O(a, ~o)] <<,M~,
(12)
where the bounds are independent of e. Further, we assume that A has a small angular average A(0), and we require dissipativity of the perturbation by assuming that A'(0) has eigenvalues with only strictly negative real parts. More precisely, we assume ]A-(0)[ ~< M0[ log~[ -(v+3),
(13)
and that there is a norm on R d (independent of e) such that, with respect to the induced matrix norm,
[le'T(°ll ~ e t#
forallt>0,
with/~ < 0 .
(14)
The constants M0 and # are assumed independent of s.
3.3. Coordinate transformations to a normal form By [11, Theorem 5] (see the variant Theorem A.1 in Appendix A), the ~0-dependence of F can be averaged out up to terms exponentially small in the distance to the toms. Using Theorem A.1 with Y = 0 and 8 = p2, we obtain: for every positive p <~r/4, there is a real-analytic symplectic coordinate transformation such that, in the new variables, the Hamiltonian H0 is
Hi(b, ~) = ~o. b + Z(b) q- R(b, ~p),
(15)
E. Hairer, C. Lubich / Applied Numerical Mathematics 29 (1999) 57-71
63
where, for all (b, lp) E D(2p, 2tr) (with a = s/4),
IZ(b)[ <~Cp 2,
JR(b, ~)1 <~CP 2e-KIps,
(16)
with ot = 1/(v + 2) and positive constants C and K independent of p. Now, we return to the perturbed system (11), rewriting it in the coordinates (b, ~), where we choose p so small that e -K/p~ = e 3, or equivalently,
p = (K/3)"+2lloge[ -w+2).
(17)
With the gradient f2(b) = Z'(b), the perturbed system is [9 = g B ( b , lp ) -'[- 0(E3p2),
(18)
~t = o) + ~2(b) + ekV(b, ~ ) + 0(83/9), where the O(.) terms arise from the gradient of R(b, ~p), using (16) and Cauchy's estimates for the derivatives. The angular average of B(b, gr) is
1 B (b) -- (2zr)~ [ B(b, ¢) de Td 1 (Ob(a ' Ob O~ (27r) d f \ o a ~o)A(a, ~o)+ ~-~ (a, ~o)~(a, ~0)) d e t - ~ - ( a , ~o) dqg, 7d t /
where a = a(b, ~) is considered as a function of b and q9 via the relation 9 = ~0(b, ap). Since, by Theorem A.1, a - b = O ( p 2 ) , Ob/Oa = I + O(p), ab/8~o = O(/92) and 0ap/0~0 = I + O(p), we obtain
[-B(b) -A(b)] ~ Cp
for Ibl ~< p.
(19)
From the remark after Theorem A.1 applied with /t = 0, for a = O(p 2) Cauchy's estimate yields 8b/Oa = I + O(p 2) and 8~/09 = I + O(p2); hence, for b = 0, we have the sharper estimates
[B(0) - A(0)[ <~ Cp 2,
[B-7(0) - AV(0) I ~< Cp.
(20)
Stoffer's transformation [12], Theorem A.2 of the Appendix, averages out the ~-dependence of the perturbation functions B and ¢P in the leading powers of e. We apply Theorem A.2 with m = 2 and 3 = 1 to (18) without the O(.) terms, then, by a real-analytic change of coordinates which is O(e)-close to the identity,
~;= 8 ( F' (c) -~-0(82/p2)), 0=
(21)
+ O(c) +
where F and O are bounded by Cp in D(p/2, a/2), and the constant C is independent of e. Bounds for the Lipschitz constants of the O(.)-terms can be obtained by Cauchy's estimates.
3.4. The attractive invariant torus The following result extends Theorem 1 of Stoffer [12], who considered dissipative perturbations of
integrable Hamiltonian systems.
64
E. Hairer, C. Lubich /Applied Numerical Mathematics 29 (1999) 57-71
Theorem 1. Under the assumptions of Sections 3.1 and 3.2, there exists to > 0 such that, f o r 0 < e <~ to, the system (11) (or, equivalently, (1)) has an attractive invariant torus, which attracts an O(p)-neighborhood (with p of(17)) with an exponential rate proportional to e.
Proof. Like Stoffer [12], we verify that the time-T flow map of (21) satisfies the conditions of the invariant manifold theorem of [7,10]. Stoffer shows this for the time-e 3 flow map of an equation like (21). Since below we use that the conditions are satisfied by the time-T flow map for arbitrary fixed T > 0, and because our assumptions differ from [12], we include the proof. The estimates (20) together with the closeness to the identity of the coordinate transformation in Theorem A.2 yield
IF(o) - ~(o)1 < Ir(o) - ~(o)1 + I~(o) - ~(o)1 = o ( t / p ) + O(p~),
Ir'(o) - ~(o) 1 ~< [r'(o) - ~(o) l + [~(o)
- A--7(0) [ =
O(e/p 2) + O(p).
Therefore, Q(c) := F ( c ) - -A(O) - A--~(O) • c is bounded by O(p 2) and Lipschitz bounded by O(p) for [cl < p / 2 . With J := A'(0), (21) reads (22)
= t J c + t-A(O) + t Q ( c ) + O ( t 3 / p 2 ) ,
0 = ~o+ O(c) + o(t2/p). For an arbitrary T, 0 < T ~< 1, consider the flow of the differential equation over the interval 0 ~< t ~< T, where the flow is O(e)-close to linear. Applying the variation of constants formula to (22) yields c(t) = eteJc(O) + j - I (e,eS _ I)A(0) + O ( t t p 2 ) .
By conditions (13) and (14), this implies
Ic(t)[
~< e'~[c(O)l + c t e l l o g e l - l p + o(tep2),
showing that, for sufficiently small e, the time-T flow of (21) maps the strip Bp/4 × ~fd, with Bp/4 = {c E ~J: Icl ~< p/4}, into itself. Similarly, we obtain bounds for the derivatives of the solution with respect to the initial values: Oc(T) Oc(T) 0-~ <~ Lcc = e ~ + O(ep), 00(0) ~ Lco = O(e3/p2), OO(T) ~ Loc = 0(1),
00(7) 00(0)
I
<<.Loo = O ( e 2 / p ) .
Thus, for sufficiently small e and hence small p, Lcc + Loo + 2
L~.oL%c<~e T~'/2 < l.
(23)
Hence, the invariant manifold theorem [7, Theorem 1; 10, Theorem 3] shows that the time-T map has an attractive invariant torus {(c*(0), 0): 0 ~ Td}, where c* :qU - - + np/4 is Lipschitz bounded by <<.2Leo~(1 - L , . - Loo) = o ( e e / p 2 ) . This invariant toms attracts orbits of the time-T flow map in Bp/4 × 7f d with attractivity factor Lcc + )~Loc <~ e T~u/2. Standard arguments [10] show that the toms is invariant under the flow of (21). [] This proof shows that the invariant torus of (11) is O(p) close (with p of (17)) to the invariant toms a = 0 of the unperturbed Hamiltonian system. This estimate could be improved under more stringent
E. Hairer, C. Lubich / Applied Numerical Mathematics 29 (1999) 57-71
65
conditions on A(0) than (13). Note, the attractive neighborhood and the attractivity rate depend only on the constants in Sections 3.1 and 3.2.
4. Symplectic integrators and backward error analysis For the numerical solution of a perturbed differential equation
= f ( y ) + sg(y),
y(O) = Yo,
(24)
with real-analytic f and g and small parameter s, consider applying a Runge-Kutta method of order p ~> 1 with step size h > 0, S
yl=Yo+hEbiY/
with
i=1 S
Yi'=f(Yi)+sg(Yi),
Yi=Yo+hZaijY
~,
i = 1 . . . . . s.
(25)
j=l
Backward error analysis interprets the numerical solution y~ as the almost exact solution y(h) of a modified differential equation y = f ( y ) + s~,(y),
y(0) = Y0,
(26)
with suitably truncated series
f ( Y ) = f ( Y ) + ht'fp+l (Y) + " " + hU-I fN(Y),
(27)
g(y) = g(y) + hPgp+l(Y) + . . . + hN-l gN(y), where the fj are independent of s, h, N, but the gj may depend on s. Expanding Yl and ~(h) in Taylor series and comparing equal powers of h yields recursions for the coefficient functions such that the numerical solution y~ is equal to ~(h) up to an error bounded by CN hN, for arbitrary N. However, in general, the series diverge as N --+ ~ . Benettin and Giorgilli [2] and Hairer and Lubich [6] show that there is a truncation index N -----N(h) such that the difference between Yl and the solution ~(h) of the modified equation (26) with (27) is exponentially small in the step size h. This extends to the following result with optimal bounds between the original and the modified functions.
Theorem 2. Let f and g be real-analytic and bounded by M in a polydisc BR(YO) = {y ~ cd: [yj -Yojl < R for all j} of radius R around Yo. Then there exist positive constants ho, C, K (depending only on M, R and the size of the coefficients of the Runge-Kutta method) such that for 0 < h <<,ho and 0 < s <<.1, with N = N(h) equal to the largest integer satisfying hN <<,ho, the difference between Yl and the solution of(26) with (27) is bounded by ]YI - y(h)] ~< Ce -K/h.
The functions f and ~ of(27) are real-analytic in BR(YO) with If(y)-f(y)l<~ChP,
] g ( y ) - g ( y ) ] < ~ C h p foryEBR/2(Yo).
E. Hairer, C. Lubich /Applied Numerical Mathematics 29 (1999) 57-71
66
Remark. This result is true for arbitrary one-step methods defining an analytic transformation. For partitioned Runge-Kutta methods, such as the symplectic Euler method used in Section 2, this follows from [6]; for other methods, one must apply analogous results of [2]. Proof. The estimate for Yl - y(h) is from [6, Theorem 1], applied to the differential equation (24), with a truncation index N proportional to R/(2Mh) independent of e, because f + eg is uniformly bounded by 2M on BR(yo) for lel ~< 1. The O(h p) bound for f ( y ) - f ( y ) follows from bounds for the coefficient functions fj given by [6, Lemma 11] applied to ~ = f ( y ) . By applying that lemma to (24), a bound of the same type is obtained for
( f ( y ) + e~(y, e)) - ( f (y) + eg(y)), uniformly for all complex e with lel ~< 1, where we have indicated the dependence on e of ~(y) = ~,(y, e) for clarity. For any fixed y ~ Blc/2(yo), 1 ([()7(y) + e'g,(y,e)) - (f(y)-t- eg(y))] -- If(y) - f ( y ) ] ) "~(y, e) - g(y) = -~ is an analytic function of e in the complex unit disk lel ~< 1, bounded by 2Ch p for lel = 1. By the maximum principle, it is bounded by 2Ch p for lel ~< 1. [] The Runge-Kutta method is symplectic if its one-step map is symplectic whenever the differential equation is Hamiltonian. From [5] (cf. also [9, Section 3(e)] and [2] for related, but weaker results), the corresponding modified differential equations are Hamiltonian, with modified Hamiltonian
:- H + hP Hp+l + . . . + h N - I H N ,
(28)
where the Hamiltonian corrections Hj are composed of derivatives of H and are thus analytic (and single-valued) on the domain of H. Therefore, a symplectic Runge-Kutta method applied to the perturbed Hamiltonian system (1) has the modified differential equations pq=
- - ( / ~ , ~ ) +e/5(fi, q), 0q OH
(29)
+
where f f =- P + hP Pp+l -q- " " q- h N - 1 P N ,
(30)
0.= Q + hp Qp+I q- " " -~- h N - I Q N , with e-dependent perturbation functions Pj, Qj. With a suitable truncation index N = N(h), which can be chosen the same for all (p, q) in a bounded complex neighborhood U of the invariant torus of (1), Theorem 2 gives
IP, - ,~(h)l + Iql - q(h)l ~< Ce-K/h.
(31)
E. Hairer, C. Lubich/Applied Numerical Mathematics 29 (1999) 57-71
67
Here, (Pl, ql) is the result of a Runge-Kutta step for (1) and (,E(h), ~(h)) is the solution of the modified equations (29) with initial value (P0, q0) e U. Further, uniformly for (p, q) ~ U, [/t(P' q) - H ( p ' q ) I <~ChP' IP(P, q ) - P(P,q)I + I Q ( P , q ) -
(32) Q(P,q)l <~Cht"
5. The attractive invariant torus of the discretization The following result shows that a symplectic (possibly partitioned) Runge-Kutta method preserves the qualitative behavior of the perturbed Hamiltonian system under a very mild restriction on the step size. For perturbed integrable systems, the existence of an attractive invariant torus of the discretization was shown previously by Stoffer [12] in a very technical proof working directly with the discrete numerical scheme and making no use of backward error analysis. Theorem 3. Let (1)fulfill the assumptions of Sections 3.1 and 3.2, and consider its numerical solution using a symplectic Runge-Kutta method of order p. Let t¢ ~ (v + 2)/p and 3x > 1. There exists eo > 0 such that, for 0 < e <~eo and 0 < h ~< Iloge[ -3~,
(33)
the Runge-Kutta map has an attractive invariant torus ~,h close to the invariant torus T~ of (1): for the Hausdorff distance of the tori,
diStH(~,/,, ~ ) ~< ChPl loge[ ~. The constant C is independent of e, h. The torus attracts an O(p )-neighborhood (with p of(17)) with an exponential rate proportional to e, uniformly in h.
Remark. We do not know if the distance estimate of the tori can be improved to Ch p. For the case where the unperturbed Hamiltonian system is integrable, such an improved estimate is stated without proof by Stoffer [12, Theorem 5]. However, it appears that his techniques also yield only an estimate with a positive power of I log el. Proof. We show that the modified differential equations (29) have an invariant toms (part (a)) with the stated properties (part (b)), and then infer the result for the Runge-Kutta map by a simple perturbation argument (part (c)). In this proof, C denotes a generic constant independent of e and h which takes on different values on different occurrences. (a) Consider the modified Hamiltonian/4 of (28) rewritten in the variables (b, lp) of (15). This is an O (h P)-perturbation of Hi on D (2p, 2a), with p of (17). By Theorem A. 1 (with Z of (15) in the role of Y, with the O(h p) perturbation in the role of Q, and without considering the remainder R(b, ~p)), there is a real-analytic symplectic coordinate transformation (b, lp) = q~0(b, ap) in terms of which the modified Hamiltonian becomes Hi(b, g r ) = w - b + Z(b') + R(b, ~),
E. Hairer, C. Lubich/Applied Numerical Mathematics 29 (1999) 57-71
68
where Z and/~ satisfy the same bounds as Z and R of (16). Moreover, Theorem A.1 (with ~ = Ch p) shows that
[~-bl <<.Ch",
[~/-~[ <~ChP/p.
(34)
This implies further I Z ( b ) - Z(b)l ~< Ch p,
(35)
IR(b, ~p) - R(b('b, ~), ~p('b, ~))1 ~< Chpe3"
(36)
In the coordinates (b, ~p), the modified differential equations (29) become *,.%
b
_
V,)
(37)
= w + ~ (b') + etP (b, ~p) + O(e3), with ~ = Z'. By (35) and Cauchy's estimate, I~(b) - 12(b)1 ~< ChP/p, and by (32) and (34), IB(b, ~) - B(b, ~)1 <<-ChP/P ,
I~(b, 1/i) - kO(b, ap) I ~< ChP/p.
Further, by (36) and Cauchy's estimate, the O(e 3) remainder terms, which are formed by the gradient of/~, are O(hpe3/p) close to those of (18). To (37) we apply Stoffer's transformation which took (18) to (21), yielding a system of the form (21) with additional terms of size eChP/p in the first components, and of size ChP/p in the second components. By a further Stoffer transformation (applied with m = 2 and 6 = hP/p by disregarding the O(.) terms in (21)), the 0-dependence of these perturbations can be eliminated in the leading powers of e, yielding c = e (F ( ~ + O ( d / p2) ), _
(38)
o=+gJ(c-3
where/~ and O differ from F and O in (21) by additional O(hP/p) terms, and the remainders given by the O(.) terms differ from those of (21) by additional o(eZ6/p) terms. For an arbitrary T, 0 < T ~< 1, we consider the time-T flow map of (38). Since h p <~Cp 3 and hence 6 ~< Cp 2 by our restriction on the stepsize (33), the influence of the perturbation terms is negligible, and this flow map satisfies Lipschitz conditions with constants as in the proof of Theorem 1. This implies the existence of an attractive invariant torus ~.h for the modified differential equations. (b) From the estimates of part (a), the right-hand sides of (21) and (38) differ at most by eChP/p in the first components, and by C h P / p in the second components. Thus, the difference of the time-T flow maps satisfies
I (T) - ¢(T)I
eChP/p,
IO(T) - O(T)[
ChP/p.
By [10, Corollary 4], the functions c*, c-* : ']i"d ---> N d whose graphs define the invariant tori of (21) and (38), respectively, differ by no more than 1
[5*(0) -- c*(O)[ < 1 -- e rEu/2 (e ChP/p + )~ChP/p) < ChP/p.
E. Hairer, C. Lubich /Applied Numerical Mathematics 29 (1999) 57-71
69
Here, X = 2Lco/(1 - Lcc - Loo) = O ( e 2 / p 2) is a Lipschitz constant of 5*. Transforming back to the original coordinates and using (34) yields (39)
diStH(~.h, T~) <<.C h P / p .
(c) Finally, by Theorem 2 and Cauchy's estimate for the derivative, the Runge-Kutta map is an exponentially small (in h) Lipschitz perturbation of the time-h flow map of the modified differential equation. Therefore, by [ 10, Theorem 3 and Corollary 4], it has an attractive invariant torus ~.h satisfying Ce-K/h
distn(~,h, ~',h) ~ 1 -- e h~/2 ~ ChP' because 3x > 1 in (33). Together with (39), this completes the proof.
(40) []
Remark. As in Corollaries 7 and 8 of [6], the following estimate holds near the toms, under the step size restriction of Theorem 3: the difference between the numerical solution of the Runge-Kutta method and the solution of the modified differential equation remains exponentially small in h for exponentially long time, uniformly in e.
Appendix A. Averaging transformations We put together two coordinate transformations based on the method of averaging which are needed in the proofs of the previous sections. Here, IIf [Ip,~ denotes the supremum norm of a function f on the complex domain D ( p , a ) defined in Section 3.1. For a function f that does not depend on the angular variables, we omit the subscript a and write IIf [I;. A. 1. Perry and Wiggins' transformation
The following result arises in the course of the proof of [11, Theorem 5].
Theorem A.1. Consider a Hamiltonian Ho(a, ~o) = w . a + F(a, ~o)
with F(a, ~o) = Y ( a ) + Q(a, ~o),
where o) E R d satisfies the strong non-resonance condition (9), and where Y and Q are real-analytic functions on D(4p, 4 a ) , f o r some positive p and a, with
IIYIhp ~< MP 2,
IIQll4p,4~ <~ M r ,
f o r a positive 6 <<.p2. There exists a real-analytic symplectic coordinate transformation 4)o : D ( 2 p , 2a) --+ D ( 4 p , 4a),
~b0: (b, 7z) v-~ (a, ~0)
taking the Hamiltonian Hi = 14o o 4)0 into the form
H~ (b, ~p) = co. b + Z ( b ) + R(b, ~p). There are positive constants C, K (depending only on d, y, v, a, M ) such that, with ot = 1/(v + 2),
IIZII2p <~ CP 2,
IIRII2R,2~~ C6e -K/p".
E. Hairer, C. Lubich /Applied Numerical Mathematics 29 (1999) 57-71
70
Moreover, if(a, ~o) = ~Po(b, aP) with (b, aP) ~ D(p, ~r), then I b - a l <~C&
l a P - f l <.. C&lp.
R e m a r k . If, instead of IIQII4p,4~ ~< M~, we assume that [Q( a, ¢P)l ~< M(lal 2 +¢S)
for (a, q)) ~ D(4p,4cr),
then we have the estimates
laP -~1 <<.C(Ibl2 + & + P3)/P,
Ib - al ~< C(Ibl 2 + ~ + p3),
(A.1)
for (b, aP) ~ D(p, or). R e m a r k . The bounds for Ib - al, laP - 91 are stronger than in [11]. They are obtained by tracing the construction in [11, Section 4.2]; in particular, see formula (33) with the present 8, (r, 2(r, 2p, 2p in the roles of s, &, ~r, p, r, respectively. For the estimates (A.1), apply the iterative lemma once with ~t = or/4 (in the notation of [11], Eq. (61)) and s times with 3 = (r/4s. A.2. Stoffer's transformation The following variant of L e m m a 2 of Stoffer [12] is obtained from its proof. T h e o r e m A.2. Consider a perturbed integrable Hamiltonian system b = sB(b, aP), ~t = w + S2(b) + sqJ (b, aP), where w ~ •d satisfies the strong non-resonance condition (9), and ~ , B, !P are real-analytic on D (p, ~r), for some positive p and or, bounded by
liBiip,~, ~< M, I1 ' -
1112]]p ~< i p ,
liB - BIIp
< M6,
Jl~lip,~,~< M, U6,
for a positive 6 <<.1. Here, B and qJ are the angular averages o r B and qJ, respectively. For each integer m >~ 1, there exist positive constants Po, Co, C, K (depending only on d, y, v, a, M, m) such that, for 0 < p <~Po and 0 < s <~ Cop, there exists a real-analytic coordinate transformation ¢~.,n :D(p/2, or/2) -+ D(p, or),
dP~,m:(c, 0) w-> (b, aP),
which puts the differential system into the form
= s(ro(c) + smr,,(c, O) + rR(c, o)), = 09 "31-(~o(C) "~- emO)m(C, O) "~- E (~R(C, 0),
where, with fl = 1/(v + 1),
lie0-
llp/2
Moreover, if (b, aP) = Ic - bl <~Ce&
~)e.m
IIFm IIp/2,~/2 <~C(~/P) m,
IIFR Ilpl2,,~/2 <~C~e-KIJ,
IIO,, Ilpn ~/2 ~< C ( ~ I P ) m-I ,
IIOR IIp/2,~/2 ~ CSe-KIJ.
(C, 0), then IO - aPl <<.Ce&.
E. Hairer, C. Lubich / Applied Numerical Mathematics 29 (1999) 57-71
71
References [1] G. Benettin, L. Galgani, A. Giorgilli and J.M. Strelcyn, A proof of Kolmogorov's theorem on invariant tori using canonical transformations defined by the Lie method, Nuovo Cimento 79B (1984) 201-223. [2] G. Benettin and A. Giorgilli, On the Hamiltonian interpolation of near to the identity symplectic mappings with application to symplectic integration algorithms, J. Statist. Phys. 74 (1994) 1117-1143. [3] N.N. Bogoliubov and Yu.A. Mitropolsky, Asymptotic Methods in the Theory of Non-Linear Oscillations (Gordon and Breach, New York, 1961). [4] N.N. Bogoliubov, Yu.A. Mitropoliskii and A.M. Samoilenko, Methods of Accelerated Convergence in Nonlinear Mechanics (Springer, Berlin, 1976). [5] E. Hairer, Backward analysis of numerical integrators and symplectic methods, Ann. Numer. Math. 1 (1994) 107-132. [6] E. Hairer and C. Lubich, The life-span of backward error analysis for numerical integrators, Numer. Math. 74 (1997) 441-462. [7] U. Kirchgraber, E Lasagni, K. Nipp and D. Stoffer, On the application of invariant manifold theory, in particular to numerical analysis, in: Intemational Series of Numerical Mathematics 97 (Birkh~iuser, Basel, 1991) 189-197. [8] U. Kirchgraber and E. Stiefel, Methoden derAnalytischen St6rungsrechnung und lhre Anwendungen (Teubner, Stuttgart, 1978). [9] J. Moser, Lectures on Hamiltonian systems, Mem. Amer. Math. Soc. 81 (1968) 1-60. [10] K. Nipp and D. Stoffer, Attractive invariant manifolds for maps: existence, smoothness and continuous dependence on the map, Report 92-11, SAM, ETH ZUrich (1992). [11] A.D. Perry and S. Wiggins, KAM tori are very sticky: rigorous lower bounds on the time to move away from an invariant Lagrangian torus with linear flow, Physica D 71 (1994) 102-121. [12] D. Stoffer, On the qualitative behaviour of symplectic integrators. Part III: Perturbed integrable systems, Preprint, Department of Mathematics, ETH ZUrich (1996).
APPLIED
"~q~ NUMERICAL MATHEMATICS ELSEVIER
Applied Numerical Mathematics 29 (1999) 57-71
Invariant toil of dissipatively perturbed Hamiltonian systems under symplectic discretization Ernst Hairer a, 1, Christian Lubich b,, a Section de Mathdmatiques, Universitd de GenOve, CH-1211 GenOve 24, Switzerland b Mathematisches Institut, Universitiit Tiibingen, Aufder Morgenstelle 10, D-72076 Tiibingen, Germany
Abstract In a recent paper, Stoffer showed that, under a very weak restriction on the step size, weakly attractive invariant tori of dissipative perturbations of integrable Hamiltonian systems persist under symplectic numerical discretizations. Stoffer's proof works directly with the discrete scheme. Here, we show how such a result, together with approximation estimates, can be obtained by combining Hamiltonian perturbation theory and backward error analysis of numerical integrators. In addition, we extend Stoffer's result to dissipative perturbations of nonintegrable Hamiltonian systems in the neighborhood of a KAM toms. © 1999 Elsevier Science B.V. and IMACS. All rights reserved.
Keywords: Weakly damped nonlinear oscillations; Symplectic numerical integrator; Hamiltonian perturbation theory; Backward error analysis; Invariant manifold; Long-time approximation
I. Introduction
We study the long-time behavior of numerical discretizations of perturbed Hamiltonian systems OH --(p,q) +sP(p,q), Oq OH gl = + ~ p (P, q) + sQ(p, q),
-
0 < s < < 1.
(1)
A classical example is Van der Pol's equation, where the unperturbed system is the harmonic oscillator, and where the perturbation creates a limit cycle; see, for example, [3], where this equation and many others are studied by the method of averaging. In higher dimensions, systems (1) often possess weakly attractive invariant toil, as shown in [8, Section 12] for perturbations of systems of non-resonant harmonic oscillators, and in [4] in a far-reaching framework. More recently, Stoffer [12] studied the * Corresponding author. E-mail: [email protected]. I E-mail: [email protected]. 0168-9274/99/$19.00 © 1999 Elsevier Science B.V. and IMACS. All rights reserved. PII: S0168-9274(98)00029-4
E. Hairer, C. Lubich /Applied Numerical Mathematics 29 (1999) 57-71
58
existence and properties of attractive invariant tori for dissipative perturbations of integrable Hamiltonian systems. Here, we assume, more generally, that the unperturbed Hamiltonian system has an invariant Lagrangian toms with linear flow, such as a KAM toms, and that the perturbation is dissipative, ensuring that (1) has a weakly attractive invariant toms. The precise assumptions and proof of existence are given in Section 3. The proof relies on averaging transformations of [11,12], and on an invariant manifold theorem of [7,10]. Does such a weakly attractive invariant toms persist under numerical discretization? For dissipative perturbations of integrable systems, Stoffer [12] shows that symplectic discretization methods preserve the invariant toms under the very weak restriction that the step size is logarithmically small in the perturbation parameter e. Stoffer's proof works directly with the discrete scheme, rebuilding the familiar constructions for continuous problems in the more complicated discrete situation, requiring great technical skill and effort. In a different approach, for the unperturbed case e = 0 of a Hamiltonian system in the neighborhood of a KAM toms, in [6] the authors show that a symplectic discretization has a nearby toms which inherits from the KAM toms the property of being 'very sticky' (a notion of Perry and Wiggins [11]): every numerical solution starting near the toms remains within twice the initial distance from the toms over a time span exponentially long in the distance and the step size. This result uses the backward error analysis of symplectic integrators, which interprets a step of the discrete numerical method as the almost-exact flow of a modified Hamiltonian system, up to an error exponentially small in the step size. This estimate is proved independently in [2,6]. Here, we show how the existence and closeness of an attractive invariant toms of a symplectic discretization of (1) can be derived using the backward error analysis result. We extend Stoffer's discrete result to situations where the unperturbed Hamiltonian is not integrable. In Section 2, for illustration, we begin with a numerical experiment with Van der Pol's equation. In Section 3, we formulate the precise assumptions about the Hamiltonian system and its perturbation, and we show the existence of a weakly attractive invariant toms. We review the required results from backward error analysis in Section 4. The existence and approximation properties of a weakly attractive invariant toms of the symplectic discretization are studied in Section 5. Finally, in Appendix A, we give two useful coordinate transformations of Perry and Wiggins [11] and of Stoffer [12], needed in the proofs.
2. Formal analysis and numerical experiment We assume that, after a suitable coordinate transformation, the system (1) is
OF h = --~-~o(a, qg) + eA(a, ~p), OF (o = w + - 7 - ( a , ~o) + e ~ ( a , qg),
(2)
06l
where a 6 •d, q9 6 ]~d, O) is a constant vector and the functions F, A, • are 2zr-periodic in the variables ~0y. For motivation, we give a non-rigorous explanation without stating the technical assumptions and without going into detail. The rigorous analysis will be presented in subsequent sections.
E. Hairer, C. Lubich /Applied Numerical Mathematics 29 (1999) 57-71
59
I. Averaging principle. Replace the first equation in (2) by the so-called averaged system _
h = eA(a),
--
1
/"
A(a) - (2zr)~ J A(a, ~o) d~0,
(3)
-g'd
with ~d = ~d/27rzd the standard toms. (Due to the periodicity of F(a, ~o) the average of OF/O~o vanishes.) The claim is that, under suitable assumptions, the averaged system gives a good approximation to the original one. Here, we consider the situation where (3) has an asymptotically stable fixed point a*, i.e., A(a*) = 0 and the eigenvalues of the Jacobian A'(a*) each have a negative real part. Hence, (3) together with the second equation of (2) has the invariant toms a = a*. Blind use of the "averaging principle" then hints at the existence of an attractive invariant toms of (2). II. B a c k w a r d error analysis. The idea is to interpret formally the numerical solution (a,, ~0n) of a pth order one-step method applied to (2) as the exact solution (a(tn), ~(tn)) at t~ = nh of the so-called modified equations a=
(4) =
+
where
OF G(a, ~o) -= - - - ( a ,
~o) + hP G p+l (a, ~o) + h p+I G p+z(a, tp) + . . . ,
A(a, qg) = A(a, qg) d- hP Ap+l (a, qg) -k- h p+I Ap+z(a, ~o) + . . . , and similarly f o r / : and qS. This allows us to study properties of the numerical solution by examining the modified equations. Assuming that (2) has an attractive invariant toms, what can be said about its numerical approximation? Applying the averaging principle to (4) yields = h~G~÷,(~) + . . . + e A ( ~ )
+eh~A,,÷,(~) + - . - .
(5)
In general, the average Gp+l (a*) does not vanish, and (5) can have an asymptotically stable fixed point (close to a*) only if h p is small compared to e. The distance of the fixed points for Eqs. (3) and (5) is O(hP/e). In contrast, applying a symplectic integrator, i.e., a method for which (an, ~0n) ~ (an+l, ~0n+l) is a symplectic transformation in the Hamiltonian case e = 0, the e-independent terms of the modified equation are Hamiltonian, i.e.,
G j(a, ~o) --
aFj
--
(a, qg).
The corresponding averages vanish identically, so that Gj (a) = 0 for all j. Consequently, Eq. (5) has an asymptotically stable fixed point if h p is small compared to 1, and its distance to a* is O(hP). Symplectic integrators not only show the correct qualitative behavior for much larger step sizes, but the distance to the exact invariant torus is smaller by a factor e.
E. Hairer, C. Lubich /Applied Numerical Mathematics 29 (1999) 57-71
60
This explanation is non-rigorous, because the averaging principle requires further assumptions, and the series in the modified equations do not converge.
Numerical experiment. Consider Van der Pol's equation in the form [~=-q + e(1-q2)P'
(6)
4=p, with small e, which is a perturbation of the harmonic oscillator. A symplectic change to polar coordinates (q = v~asin99, p = x/2acos99) puts the system in the form (2) with co = 1, F(a, 99) = 0, A(a, 99) = 2a cos 2 99(1 - 2a sin 2 99) and qb(a, 99) -- cos 99sin 99(1 - 2a sin 2 qg). The angular average of A (a, 99) is A(a) = a(1 - a/2), so that with a* = 2 we have A(a*) = 0, A'(a*) < O. We solve (6) numerically with two initial values, (P0, q0) = (0, 1.4) and (P0, q0) = (0, 2.6), and with three methods: the non-symplectic explicit and implicit Euler methods, and the symplectic Euler method; all have order 1. The results are displayed in Fig. 1. For large step sizes (compared to the perturbation parameter e), the non-symplectic methods give a completely wrong numerical solution, whereas the symplectic method is qualitatively correct. For smaller stepsizes, the non-symplectic numerical methods also show a limit cycle. We explain these observations by the averaging principle and a backward error analysis argument. For a differential equation :~ = f (y) + eg (y), the explicit Euler discretization is Yn+l = Yn + h ( f (yn) + eg (y~)). Here, the modified differential equation is
y= f(y) +eg(y)-~
h
f ' ( y ) f ( y ) +O(h 2 +eh).
(7)
In the above coordinates, we get a system (4) with G2(a, 99) = a. Therefore, the right-hand side of (5) is approximately ha + ea(1 - a/2), which has a = 2 + 2h/e as positive root. Hence, the limit cycle of the numerical solution of the explicit Euler method has approximate radius 2 + 2h/e (Fig. l) which is far from the correct value. The implicit Euler discretization is adjoint to the explicit Euler method. Therefore, its modified differential equation is (7) with h replaced by - h . In this case, the radius of the limit cycle is approximately 2 - 2h/e (for h < e), which again agrees very well with Fig. 1. For the symplectic Euler method
h / OH ) Pn+l = Pn + ~--~-a (Pn+l, qn) + eP(Pn+l, qn) , OH (Pn+l,qn) + eQ(Pn+l, qn)), qn+l =qn + h /~-~p the modified differential equations for Van der Pol's equation are p=--q'We(1 _~2)~+
h 2~+O(h
2_+_eh),
h
q = p - ~ q +O(h2-F eh). Since the method is symplectic, all e-independent terms in (5) vanish. Therefore, the radius of the limit cycle is of size 2 + O(h); see Fig. 1.
E. Hairer, C. Lubich /Applied Numerical Mathematics 29 (1999) 57-71 I
61
i
implicit Euler, h = 0.18
symplectic Euler, h = 0.18
explicit Euler, h = 0.18
%F i
i
implicit Euler, h = 0.06
symplectic Euler, h = 0.06
implicit Eule'r, h = 0.02
V
i
llll
explicit Euler, h --- 0.06
M
Fig. 1. Numerical experiments with Van der Pol's equation, e = 0.05.
3. Dissipatively perturbed Hamiltonian systems with a weakly attractive invariant torus 3.1. Assumptions about the unperturbed Hamiltonian system We assume there is a real-analytic symplectic coordinate transform tp : (a, ~o) ~ (p, q) which takes the Hamiltonian H0 = H o tp into the form H0(a, ~o) = w , a + F(a, ~o).
(8)
Here, o9 6 •d is a constant vector of frequencies and o9. a :- OgTa. For positive constants 7, v, we assume w satisfies the strong non-resonance condition Iog " kl >1 71kl-",
k ~ Zd\{0}.
(9)
E. Hairer, C. Lubich /Applied Numerical Mathematics 29 (1999) 57-71
62
The function F is assumed 2zr-periodic in the variables ~0j and real-analytic in the complex neighbourhood
D(r, s) = { (a, ~o) e C d x (ca/2zrza): laj[ < r, IIm~ojl < s for all j} of the toms {a = 0} c ~d x "IFd. Further, we assume that with a constant M and for all (a, ~o) e D(r, s),
IF(a, ~0)1 <~M [a[2.
(10)
The equations of motion of the unperturbed system with Hamiltonian H0 are & = O(la12),
~b= o~ + O(lal).
Hence, the toms a = 0 is invariant under the flow of this Hamiltonian system, and the flow is linear on this toms. Note, (8) with (10) is precisely the form delivered by the KAM theorem, in the neighborhood of a KAM toms (see, for example, [1]). Condition (9) is also required by the KAM theorem. Hamiltonian systems satisfying the very assumptions (8)-(10) are studied in [11].
3.2. Assumptions about the perturbation In the coordinates (a, ~o), the equations for the perturbed system become aF
= --g-z(a, ~o) + eA(a, ~0), oW
~F
(11)
(o = w + ff-aa (a, ~o) + cO(a, ~o). We assume that A and • (which may also depend on e) are real-analytic on D(r, s) such that, for all (a, ¢p) e D(r, s),
[A(a, ~0)[ ~< MA,
]O(a, ~o)] <<,M~,
(12)
where the bounds are independent of e. Further, we assume that A has a small angular average A(0), and we require dissipativity of the perturbation by assuming that A'(0) has eigenvalues with only strictly negative real parts. More precisely, we assume ]A-(0)[ ~< M0[ log~[ -(v+3),
(13)
and that there is a norm on R d (independent of e) such that, with respect to the induced matrix norm,
[le'T(°ll ~ e t#
forallt>0,
with/~ < 0 .
(14)
The constants M0 and # are assumed independent of s.
3.3. Coordinate transformations to a normal form By [11, Theorem 5] (see the variant Theorem A.1 in Appendix A), the ~0-dependence of F can be averaged out up to terms exponentially small in the distance to the toms. Using Theorem A.1 with Y = 0 and 8 = p2, we obtain: for every positive p <~r/4, there is a real-analytic symplectic coordinate transformation such that, in the new variables, the Hamiltonian H0 is
Hi(b, ~) = ~o. b + Z(b) q- R(b, ~p),
(15)
E. Hairer, C. Lubich / Applied Numerical Mathematics 29 (1999) 57-71
63
where, for all (b, lp) E D(2p, 2tr) (with a = s/4),
IZ(b)[ <~Cp 2,
JR(b, ~)1 <~CP 2e-KIps,
(16)
with ot = 1/(v + 2) and positive constants C and K independent of p. Now, we return to the perturbed system (11), rewriting it in the coordinates (b, ~), where we choose p so small that e -K/p~ = e 3, or equivalently,
p = (K/3)"+2lloge[ -w+2).
(17)
With the gradient f2(b) = Z'(b), the perturbed system is [9 = g B ( b , lp ) -'[- 0(E3p2),
(18)
~t = o) + ~2(b) + ekV(b, ~ ) + 0(83/9), where the O(.) terms arise from the gradient of R(b, ~p), using (16) and Cauchy's estimates for the derivatives. The angular average of B(b, gr) is
1 B (b) -- (2zr)~ [ B(b, ¢) de Td 1 (Ob(a ' Ob O~ (27r) d f \ o a ~o)A(a, ~o)+ ~-~ (a, ~o)~(a, ~0)) d e t - ~ - ( a , ~o) dqg, 7d t /
where a = a(b, ~) is considered as a function of b and q9 via the relation 9 = ~0(b, ap). Since, by Theorem A.1, a - b = O ( p 2 ) , Ob/Oa = I + O(p), ab/8~o = O(/92) and 0ap/0~0 = I + O(p), we obtain
[-B(b) -A(b)] ~ Cp
for Ibl ~< p.
(19)
From the remark after Theorem A.1 applied with /t = 0, for a = O(p 2) Cauchy's estimate yields 8b/Oa = I + O(p 2) and 8~/09 = I + O(p2); hence, for b = 0, we have the sharper estimates
[B(0) - A(0)[ <~ Cp 2,
[B-7(0) - AV(0) I ~< Cp.
(20)
Stoffer's transformation [12], Theorem A.2 of the Appendix, averages out the ~-dependence of the perturbation functions B and ¢P in the leading powers of e. We apply Theorem A.2 with m = 2 and 3 = 1 to (18) without the O(.) terms, then, by a real-analytic change of coordinates which is O(e)-close to the identity,
~;= 8 ( F' (c) -~-0(82/p2)), 0=
(21)
+ O(c) +
where F and O are bounded by Cp in D(p/2, a/2), and the constant C is independent of e. Bounds for the Lipschitz constants of the O(.)-terms can be obtained by Cauchy's estimates.
3.4. The attractive invariant torus The following result extends Theorem 1 of Stoffer [12], who considered dissipative perturbations of
integrable Hamiltonian systems.
64
E. Hairer, C. Lubich /Applied Numerical Mathematics 29 (1999) 57-71
Theorem 1. Under the assumptions of Sections 3.1 and 3.2, there exists to > 0 such that, f o r 0 < e <~ to, the system (11) (or, equivalently, (1)) has an attractive invariant torus, which attracts an O(p)-neighborhood (with p of(17)) with an exponential rate proportional to e.
Proof. Like Stoffer [12], we verify that the time-T flow map of (21) satisfies the conditions of the invariant manifold theorem of [7,10]. Stoffer shows this for the time-e 3 flow map of an equation like (21). Since below we use that the conditions are satisfied by the time-T flow map for arbitrary fixed T > 0, and because our assumptions differ from [12], we include the proof. The estimates (20) together with the closeness to the identity of the coordinate transformation in Theorem A.2 yield
IF(o) - ~(o)1 < Ir(o) - ~(o)1 + I~(o) - ~(o)1 = o ( t / p ) + O(p~),
Ir'(o) - ~(o) 1 ~< [r'(o) - ~(o) l + [~(o)
- A--7(0) [ =
O(e/p 2) + O(p).
Therefore, Q(c) := F ( c ) - -A(O) - A--~(O) • c is bounded by O(p 2) and Lipschitz bounded by O(p) for [cl < p / 2 . With J := A'(0), (21) reads (22)
= t J c + t-A(O) + t Q ( c ) + O ( t 3 / p 2 ) ,
0 = ~o+ O(c) + o(t2/p). For an arbitrary T, 0 < T ~< 1, consider the flow of the differential equation over the interval 0 ~< t ~< T, where the flow is O(e)-close to linear. Applying the variation of constants formula to (22) yields c(t) = eteJc(O) + j - I (e,eS _ I)A(0) + O ( t t p 2 ) .
By conditions (13) and (14), this implies
Ic(t)[
~< e'~[c(O)l + c t e l l o g e l - l p + o(tep2),
showing that, for sufficiently small e, the time-T flow of (21) maps the strip Bp/4 × ~fd, with Bp/4 = {c E ~J: Icl ~< p/4}, into itself. Similarly, we obtain bounds for the derivatives of the solution with respect to the initial values: Oc(T) Oc(T) 0-~ <~ Lcc = e ~ + O(ep), 00(0) ~ Lco = O(e3/p2), OO(T) ~ Loc = 0(1),
00(7) 00(0)
I
<<.Loo = O ( e 2 / p ) .
Thus, for sufficiently small e and hence small p, Lcc + Loo + 2
L~.oL%c<~e T~'/2 < l.
(23)
Hence, the invariant manifold theorem [7, Theorem 1; 10, Theorem 3] shows that the time-T map has an attractive invariant torus {(c*(0), 0): 0 ~ Td}, where c* :qU - - + np/4 is Lipschitz bounded by <<.2Leo~(1 - L , . - Loo) = o ( e e / p 2 ) . This invariant toms attracts orbits of the time-T flow map in Bp/4 × 7f d with attractivity factor Lcc + )~Loc <~ e T~u/2. Standard arguments [10] show that the toms is invariant under the flow of (21). [] This proof shows that the invariant torus of (11) is O(p) close (with p of (17)) to the invariant toms a = 0 of the unperturbed Hamiltonian system. This estimate could be improved under more stringent
E. Hairer, C. Lubich / Applied Numerical Mathematics 29 (1999) 57-71
65
conditions on A(0) than (13). Note, the attractive neighborhood and the attractivity rate depend only on the constants in Sections 3.1 and 3.2.
4. Symplectic integrators and backward error analysis For the numerical solution of a perturbed differential equation
= f ( y ) + sg(y),
y(O) = Yo,
(24)
with real-analytic f and g and small parameter s, consider applying a Runge-Kutta method of order p ~> 1 with step size h > 0, S
yl=Yo+hEbiY/
with
i=1 S
Yi'=f(Yi)+sg(Yi),
Yi=Yo+hZaijY
~,
i = 1 . . . . . s.
(25)
j=l
Backward error analysis interprets the numerical solution y~ as the almost exact solution y(h) of a modified differential equation y = f ( y ) + s~,(y),
y(0) = Y0,
(26)
with suitably truncated series
f ( Y ) = f ( Y ) + ht'fp+l (Y) + " " + hU-I fN(Y),
(27)
g(y) = g(y) + hPgp+l(Y) + . . . + hN-l gN(y), where the fj are independent of s, h, N, but the gj may depend on s. Expanding Yl and ~(h) in Taylor series and comparing equal powers of h yields recursions for the coefficient functions such that the numerical solution y~ is equal to ~(h) up to an error bounded by CN hN, for arbitrary N. However, in general, the series diverge as N --+ ~ . Benettin and Giorgilli [2] and Hairer and Lubich [6] show that there is a truncation index N -----N(h) such that the difference between Yl and the solution ~(h) of the modified equation (26) with (27) is exponentially small in the step size h. This extends to the following result with optimal bounds between the original and the modified functions.
Theorem 2. Let f and g be real-analytic and bounded by M in a polydisc BR(YO) = {y ~ cd: [yj -Yojl < R for all j} of radius R around Yo. Then there exist positive constants ho, C, K (depending only on M, R and the size of the coefficients of the Runge-Kutta method) such that for 0 < h <<,ho and 0 < s <<.1, with N = N(h) equal to the largest integer satisfying hN <<,ho, the difference between Yl and the solution of(26) with (27) is bounded by ]YI - y(h)] ~< Ce -K/h.
The functions f and ~ of(27) are real-analytic in BR(YO) with If(y)-f(y)l<~ChP,
] g ( y ) - g ( y ) ] < ~ C h p foryEBR/2(Yo).
E. Hairer, C. Lubich /Applied Numerical Mathematics 29 (1999) 57-71
66
Remark. This result is true for arbitrary one-step methods defining an analytic transformation. For partitioned Runge-Kutta methods, such as the symplectic Euler method used in Section 2, this follows from [6]; for other methods, one must apply analogous results of [2]. Proof. The estimate for Yl - y(h) is from [6, Theorem 1], applied to the differential equation (24), with a truncation index N proportional to R/(2Mh) independent of e, because f + eg is uniformly bounded by 2M on BR(yo) for lel ~< 1. The O(h p) bound for f ( y ) - f ( y ) follows from bounds for the coefficient functions fj given by [6, Lemma 11] applied to ~ = f ( y ) . By applying that lemma to (24), a bound of the same type is obtained for
( f ( y ) + e~(y, e)) - ( f (y) + eg(y)), uniformly for all complex e with lel ~< 1, where we have indicated the dependence on e of ~(y) = ~,(y, e) for clarity. For any fixed y ~ Blc/2(yo), 1 ([()7(y) + e'g,(y,e)) - (f(y)-t- eg(y))] -- If(y) - f ( y ) ] ) "~(y, e) - g(y) = -~ is an analytic function of e in the complex unit disk lel ~< 1, bounded by 2Ch p for lel = 1. By the maximum principle, it is bounded by 2Ch p for lel ~< 1. [] The Runge-Kutta method is symplectic if its one-step map is symplectic whenever the differential equation is Hamiltonian. From [5] (cf. also [9, Section 3(e)] and [2] for related, but weaker results), the corresponding modified differential equations are Hamiltonian, with modified Hamiltonian
:- H + hP Hp+l + . . . + h N - I H N ,
(28)
where the Hamiltonian corrections Hj are composed of derivatives of H and are thus analytic (and single-valued) on the domain of H. Therefore, a symplectic Runge-Kutta method applied to the perturbed Hamiltonian system (1) has the modified differential equations pq=
- - ( / ~ , ~ ) +e/5(fi, q), 0q OH
(29)
+
where f f =- P + hP Pp+l -q- " " q- h N - 1 P N ,
(30)
0.= Q + hp Qp+I q- " " -~- h N - I Q N , with e-dependent perturbation functions Pj, Qj. With a suitable truncation index N = N(h), which can be chosen the same for all (p, q) in a bounded complex neighborhood U of the invariant torus of (1), Theorem 2 gives
IP, - ,~(h)l + Iql - q(h)l ~< Ce-K/h.
(31)
E. Hairer, C. Lubich/Applied Numerical Mathematics 29 (1999) 57-71
67
Here, (Pl, ql) is the result of a Runge-Kutta step for (1) and (,E(h), ~(h)) is the solution of the modified equations (29) with initial value (P0, q0) e U. Further, uniformly for (p, q) ~ U, [/t(P' q) - H ( p ' q ) I <~ChP' IP(P, q ) - P(P,q)I + I Q ( P , q ) -
(32) Q(P,q)l <~Cht"
5. The attractive invariant torus of the discretization The following result shows that a symplectic (possibly partitioned) Runge-Kutta method preserves the qualitative behavior of the perturbed Hamiltonian system under a very mild restriction on the step size. For perturbed integrable systems, the existence of an attractive invariant torus of the discretization was shown previously by Stoffer [12] in a very technical proof working directly with the discrete numerical scheme and making no use of backward error analysis. Theorem 3. Let (1)fulfill the assumptions of Sections 3.1 and 3.2, and consider its numerical solution using a symplectic Runge-Kutta method of order p. Let t¢ ~ (v + 2)/p and 3x > 1. There exists eo > 0 such that, for 0 < e <~eo and 0 < h ~< Iloge[ -3~,
(33)
the Runge-Kutta map has an attractive invariant torus ~,h close to the invariant torus T~ of (1): for the Hausdorff distance of the tori,
diStH(~,/,, ~ ) ~< ChPl loge[ ~. The constant C is independent of e, h. The torus attracts an O(p )-neighborhood (with p of(17)) with an exponential rate proportional to e, uniformly in h.
Remark. We do not know if the distance estimate of the tori can be improved to Ch p. For the case where the unperturbed Hamiltonian system is integrable, such an improved estimate is stated without proof by Stoffer [12, Theorem 5]. However, it appears that his techniques also yield only an estimate with a positive power of I log el. Proof. We show that the modified differential equations (29) have an invariant toms (part (a)) with the stated properties (part (b)), and then infer the result for the Runge-Kutta map by a simple perturbation argument (part (c)). In this proof, C denotes a generic constant independent of e and h which takes on different values on different occurrences. (a) Consider the modified Hamiltonian/4 of (28) rewritten in the variables (b, lp) of (15). This is an O (h P)-perturbation of Hi on D (2p, 2a), with p of (17). By Theorem A. 1 (with Z of (15) in the role of Y, with the O(h p) perturbation in the role of Q, and without considering the remainder R(b, ~p)), there is a real-analytic symplectic coordinate transformation (b, lp) = q~0(b, ap) in terms of which the modified Hamiltonian becomes Hi(b, g r ) = w - b + Z(b') + R(b, ~),
E. Hairer, C. Lubich/Applied Numerical Mathematics 29 (1999) 57-71
68
where Z and/~ satisfy the same bounds as Z and R of (16). Moreover, Theorem A.1 (with ~ = Ch p) shows that
[~-bl <<.Ch",
[~/-~[ <~ChP/p.
(34)
This implies further I Z ( b ) - Z(b)l ~< Ch p,
(35)
IR(b, ~p) - R(b('b, ~), ~p('b, ~))1 ~< Chpe3"
(36)
In the coordinates (b, ~p), the modified differential equations (29) become *,.%
b
_
V,)
(37)
= w + ~ (b') + etP (b, ~p) + O(e3), with ~ = Z'. By (35) and Cauchy's estimate, I~(b) - 12(b)1 ~< ChP/p, and by (32) and (34), IB(b, ~) - B(b, ~)1 <<-ChP/P ,
I~(b, 1/i) - kO(b, ap) I ~< ChP/p.
Further, by (36) and Cauchy's estimate, the O(e 3) remainder terms, which are formed by the gradient of/~, are O(hpe3/p) close to those of (18). To (37) we apply Stoffer's transformation which took (18) to (21), yielding a system of the form (21) with additional terms of size eChP/p in the first components, and of size ChP/p in the second components. By a further Stoffer transformation (applied with m = 2 and 6 = hP/p by disregarding the O(.) terms in (21)), the 0-dependence of these perturbations can be eliminated in the leading powers of e, yielding c = e (F ( ~ + O ( d / p2) ), _
(38)
o=+gJ(c-3
where/~ and O differ from F and O in (21) by additional O(hP/p) terms, and the remainders given by the O(.) terms differ from those of (21) by additional o(eZ6/p) terms. For an arbitrary T, 0 < T ~< 1, we consider the time-T flow map of (38). Since h p <~Cp 3 and hence 6 ~< Cp 2 by our restriction on the stepsize (33), the influence of the perturbation terms is negligible, and this flow map satisfies Lipschitz conditions with constants as in the proof of Theorem 1. This implies the existence of an attractive invariant torus ~.h for the modified differential equations. (b) From the estimates of part (a), the right-hand sides of (21) and (38) differ at most by eChP/p in the first components, and by C h P / p in the second components. Thus, the difference of the time-T flow maps satisfies
I (T) - ¢(T)I
eChP/p,
IO(T) - O(T)[
ChP/p.
By [10, Corollary 4], the functions c*, c-* : ']i"d ---> N d whose graphs define the invariant tori of (21) and (38), respectively, differ by no more than 1
[5*(0) -- c*(O)[ < 1 -- e rEu/2 (e ChP/p + )~ChP/p) < ChP/p.
E. Hairer, C. Lubich /Applied Numerical Mathematics 29 (1999) 57-71
69
Here, X = 2Lco/(1 - Lcc - Loo) = O ( e 2 / p 2) is a Lipschitz constant of 5*. Transforming back to the original coordinates and using (34) yields (39)
diStH(~.h, T~) <<.C h P / p .
(c) Finally, by Theorem 2 and Cauchy's estimate for the derivative, the Runge-Kutta map is an exponentially small (in h) Lipschitz perturbation of the time-h flow map of the modified differential equation. Therefore, by [ 10, Theorem 3 and Corollary 4], it has an attractive invariant torus ~.h satisfying Ce-K/h
distn(~,h, ~',h) ~ 1 -- e h~/2 ~ ChP' because 3x > 1 in (33). Together with (39), this completes the proof.
(40) []
Remark. As in Corollaries 7 and 8 of [6], the following estimate holds near the toms, under the step size restriction of Theorem 3: the difference between the numerical solution of the Runge-Kutta method and the solution of the modified differential equation remains exponentially small in h for exponentially long time, uniformly in e.
Appendix A. Averaging transformations We put together two coordinate transformations based on the method of averaging which are needed in the proofs of the previous sections. Here, IIf [Ip,~ denotes the supremum norm of a function f on the complex domain D ( p , a ) defined in Section 3.1. For a function f that does not depend on the angular variables, we omit the subscript a and write IIf [I;. A. 1. Perry and Wiggins' transformation
The following result arises in the course of the proof of [11, Theorem 5].
Theorem A.1. Consider a Hamiltonian Ho(a, ~o) = w . a + F(a, ~o)
with F(a, ~o) = Y ( a ) + Q(a, ~o),
where o) E R d satisfies the strong non-resonance condition (9), and where Y and Q are real-analytic functions on D(4p, 4 a ) , f o r some positive p and a, with
IIYIhp ~< MP 2,
IIQll4p,4~ <~ M r ,
f o r a positive 6 <<.p2. There exists a real-analytic symplectic coordinate transformation 4)o : D ( 2 p , 2a) --+ D ( 4 p , 4a),
~b0: (b, 7z) v-~ (a, ~0)
taking the Hamiltonian Hi = 14o o 4)0 into the form
H~ (b, ~p) = co. b + Z ( b ) + R(b, ~p). There are positive constants C, K (depending only on d, y, v, a, M ) such that, with ot = 1/(v + 2),
IIZII2p <~ CP 2,
IIRII2R,2~~ C6e -K/p".
E. Hairer, C. Lubich /Applied Numerical Mathematics 29 (1999) 57-71
70
Moreover, if(a, ~o) = ~Po(b, aP) with (b, aP) ~ D(p, ~r), then I b - a l <~C&
l a P - f l <.. C&lp.
R e m a r k . If, instead of IIQII4p,4~ ~< M~, we assume that [Q( a, ¢P)l ~< M(lal 2 +¢S)
for (a, q)) ~ D(4p,4cr),
then we have the estimates
laP -~1 <<.C(Ibl2 + & + P3)/P,
Ib - al ~< C(Ibl 2 + ~ + p3),
(A.1)
for (b, aP) ~ D(p, or). R e m a r k . The bounds for Ib - al, laP - 91 are stronger than in [11]. They are obtained by tracing the construction in [11, Section 4.2]; in particular, see formula (33) with the present 8, (r, 2(r, 2p, 2p in the roles of s, &, ~r, p, r, respectively. For the estimates (A.1), apply the iterative lemma once with ~t = or/4 (in the notation of [11], Eq. (61)) and s times with 3 = (r/4s. A.2. Stoffer's transformation The following variant of L e m m a 2 of Stoffer [12] is obtained from its proof. T h e o r e m A.2. Consider a perturbed integrable Hamiltonian system b = sB(b, aP), ~t = w + S2(b) + sqJ (b, aP), where w ~ •d satisfies the strong non-resonance condition (9), and ~ , B, !P are real-analytic on D (p, ~r), for some positive p and or, bounded by
liBiip,~, ~< M, I1 ' -
1112]]p ~< i p ,
liB - BIIp
< M6,
Jl~lip,~,~< M, U6,
for a positive 6 <<.1. Here, B and qJ are the angular averages o r B and qJ, respectively. For each integer m >~ 1, there exist positive constants Po, Co, C, K (depending only on d, y, v, a, M, m) such that, for 0 < p <~Po and 0 < s <~ Cop, there exists a real-analytic coordinate transformation ¢~.,n :D(p/2, or/2) -+ D(p, or),
dP~,m:(c, 0) w-> (b, aP),
which puts the differential system into the form
= s(ro(c) + smr,,(c, O) + rR(c, o)), = 09 "31-(~o(C) "~- emO)m(C, O) "~- E (~R(C, 0),
where, with fl = 1/(v + 1),
lie0-
llp/2
Moreover, if (b, aP) = Ic - bl <~Ce&
~)e.m
IIFm IIp/2,~/2 <~C(~/P) m,
IIFR Ilpl2,,~/2 <~C~e-KIJ,
IIO,, Ilpn ~/2 ~< C ( ~ I P ) m-I ,
IIOR IIp/2,~/2 ~ CSe-KIJ.
(C, 0), then IO - aPl <<.Ce&.
E. Hairer, C. Lubich / Applied Numerical Mathematics 29 (1999) 57-71
71
References [1] G. Benettin, L. Galgani, A. Giorgilli and J.M. Strelcyn, A proof of Kolmogorov's theorem on invariant tori using canonical transformations defined by the Lie method, Nuovo Cimento 79B (1984) 201-223. [2] G. Benettin and A. Giorgilli, On the Hamiltonian interpolation of near to the identity symplectic mappings with application to symplectic integration algorithms, J. Statist. Phys. 74 (1994) 1117-1143. [3] N.N. Bogoliubov and Yu.A. Mitropolsky, Asymptotic Methods in the Theory of Non-Linear Oscillations (Gordon and Breach, New York, 1961). [4] N.N. Bogoliubov, Yu.A. Mitropoliskii and A.M. Samoilenko, Methods of Accelerated Convergence in Nonlinear Mechanics (Springer, Berlin, 1976). [5] E. Hairer, Backward analysis of numerical integrators and symplectic methods, Ann. Numer. Math. 1 (1994) 107-132. [6] E. Hairer and C. Lubich, The life-span of backward error analysis for numerical integrators, Numer. Math. 74 (1997) 441-462. [7] U. Kirchgraber, E Lasagni, K. Nipp and D. Stoffer, On the application of invariant manifold theory, in particular to numerical analysis, in: Intemational Series of Numerical Mathematics 97 (Birkh~iuser, Basel, 1991) 189-197. [8] U. Kirchgraber and E. Stiefel, Methoden derAnalytischen St6rungsrechnung und lhre Anwendungen (Teubner, Stuttgart, 1978). [9] J. Moser, Lectures on Hamiltonian systems, Mem. Amer. Math. Soc. 81 (1968) 1-60. [10] K. Nipp and D. Stoffer, Attractive invariant manifolds for maps: existence, smoothness and continuous dependence on the map, Report 92-11, SAM, ETH ZUrich (1992). [11] A.D. Perry and S. Wiggins, KAM tori are very sticky: rigorous lower bounds on the time to move away from an invariant Lagrangian torus with linear flow, Physica D 71 (1994) 102-121. [12] D. Stoffer, On the qualitative behaviour of symplectic integrators. Part III: Perturbed integrable systems, Preprint, Department of Mathematics, ETH ZUrich (1996).