Backward error analysis of a full discretization scheme for a class of semilinear parabolic partial differential equations

Nonlinear Analysis 52 (2003) 805 – 826

www.elsevier.com/locate/na

Backward error analysis of a full discretization scheme for a class of semilinear parabolic partial di%erential equations Karsten Matthies Institut fur Mathematik I, Freie Universitat Berlin, Arnimallee 2-6, 14195 Berlin, Germany Received 27 March 2001; accepted 2 October 2001

Abstract A numerical method for reaction–di%usion equations with analytic nonlinearity is presented, for which a rigorous backward error analysis is possible. We construct a modi1ed equation, which describes the behavior of the full discretization scheme up to exponentially small errors in the step size. In the construction the numerical scheme is 1rst exactly embedded into a nonautonomous equation. This equation is then averaged with only exponentially small remainder terms. The long-time behavior near hyperbolic equilibria, the persistence of homoclinic orbits and regularity properties are analyzed. ? 2002 Elsevier Science Ltd. All rights reserved. Keywords: Backward error analysis; Parabolic partial di%erential equations; Rapid periodic forcing; Exponential averaging

1. Introduction For an illustration of backward error analysis of di%erential equations, we consider 1rst an ordinary di%erential equation with analytic right-hand side x˙ = f(x);

x(0) = x0 ∈ Rm :

(1.1)

We want to compare the =ow (t; u) with numerical methods, e.g. some Runge–Kutta method with step-size h, which de1nes a local di%eomorphism xj+1 = (h; xj ): E-mail address: [email protected] (K. Matthies). 0362-546X/02/$ - see front matter ? 2002 Elsevier Science Ltd. All rights reserved. PII: S 0 3 6 2 - 5 4 6 X ( 0 2 ) 0 0 1 3 4 - 7

(1.2)

806

K. Matthies / Nonlinear Analysis 52 (2003) 805 – 826

Our main intention is to compare this dynamical system with the =ow of some modi1ed ˜ equation x˙ = f(x), with f˜ inside the same class as f, rather than compare it to the =ow of the original Eq. (1.1). For ordinary di%erential equations there is a strong relationship between backward analysis of numerical schemes and the averaging of nonautonomous equations. One-step methods for Eq. (1.1) with step-size h can be interpreted as the exact time-h-map of a rapidly forced ordinary di%erential equation as
t  x˙ = f(x) + hp g x; ; h ; x(0) = x0 ∈ Rm : (1.3) h with g(:; ; :) periodic in = t=h with period 1, see [6]. For such rapidly forced equations like (1.3) Neishtadt [15] proved, that the equation can be transformed to
t  ˜ (1.4) y˙ = f(y) +  y; ; h h with  ¡ c2 exp(−c1 =h) and f˜ − f ¡ c3 hp on a bounded domain. Hence there is a modi1ed autonomous vector 1eld near f, which describes very accurately the behavior of the numerical method. Instead of going via the embedding into nonautonomous vector 1elds and using the Neishtadt theorem, there are more recent direct methods. For real analytic f bounded on some poly-disc around x0 , Hairer and ˜ x) Lubich [8, Theorem 1] construct directly a f(h; ˜ such that
c  2 (1.5) x1 − x(h) ˜ 6 c1 exp − h for h ¡ h0 with explicit estimates of c1 ; c2 and h0 . If x(t) remains inside some compact set K, where f is also real analytic and bounded, then such an extremely good approximation also holds for 1nite times T = O(h−1 )

c  for jh 6 T: (1.6) xj − x(jh) ˜ = O exp − h These exponential estimates can be used to explain, why several qualitative e%ects are well replicated by numerical schemes. Some of these properties are the local behavior near steady-states, stable periodic orbits, homoclinic orbits [6] and energy conservation for symplectic integrators [1,17]. There is a vast literature on di%erent aspects of backward error analysis for ordinary di%erential equations, see [8,16] and the review article [9] and references therein. Some early ideas for backward error analysis of partial differential equations can be found in [4,20,21]. However, no rigorous theorems on the existence of modi1ed equations are given. 1.1. Main results We will consider semilinear parabolic partial di%erential equations, see e.g. [10]. As a model take a system of reaction–di%usion equations with U = (u1 ; : : : ; un ),

K. Matthies / Nonlinear Analysis 52 (2003) 805 – 826

807

D = diag(d1 ; : : : ; dn ) with di ¿ 0, and F = (f1 ; : : : ; fn ) @ U (x; t) − DIU (x; t) = F(U (x; t)) @t

(1.7)

with periodic boundary conditions on  = [0; 2]d ; d = 1; 2; 3 and initial conditions s U (0) = U0 in the Sobolev space Hper (; Rn ) ⊂ C 0 (; Rn ), i.e. s ¿ d=2. For a de1nition n s of Hper (; R ) for noninteger s see e.g. [19, p. 48]. The nonlinearity F is assumed to s be an entire function Rn → Rn . Eq. (1.7) de1nes a semi=ow on Hper (; Rn ) = D(As=2 ) with the di%erential operator A = −D$. The solutions of this equation are extremely regular, they are in Gevrey classes G&s=2 (; Rn ) as shown in [5]. The Fourier coeKcients Vj ; j ∈ Zd of a function V in the Gevrey class G&s=2 (; Rn ) decay like j−s exp(−&j). A norm of G&s=2 (; Rn ) is given by  1=2 n   j |vk |2 (1 + dk j2 )s exp(2&j) ; (1.8) |V |G&s=2 =  k=1 j∈Zd

where vkj is the kth component of the Fourier coeKcient of the spatial Fourier expansion into exp(ijx) of the periodic functions V . Alternatively, we can de1ne Gevrey classes as G&s=2 (; Rn ) = D(As=2 exp(&(−$)1=2 ))

(1.9)

with A=−D$. The spatial Fourier modes decay exponentially fast, where the exponent depends on time. In [5] it is shown for the solution S(t; U0 ) of (1.7), that S(t; U0 ) ∈ Gts=2

for 0 6 t 6 T ∗ ;

S(t; U0 ) ∈ GTs=2∗

for t ¿ T ∗

(1.10)

holds for some maximal Gevrey exponent T ∗ ¿ 0. Furthermore the Gevrey norm remains bounded |S(t; U0 )|Gs=2

min(t; T ∗ )

6 C0

s 6 M0 , i.e. if only a suitable Sobolev norm of the solution as long as |S(t; U0 )|Hper remains bounded, see [5]. We consider a special discretization scheme of (1.7). Modi1ed equations will be constructed, which describe the behavior of the scheme up to exponentially small errors. We have to stress that the analyzed scheme is hardly used in practice due to very strict requirements on the size of the time steps. When for an explicit method standardly one would choose a step-size It ≈ C=(Ix)2 , we choose It ≈ C=(Ix)2(1+,) , which gives a bad, but still algebraic complexity of the method. Anyway, this scheme will give a 1rst example of a numerical scheme in parabolic partial di%erential equations for which a rigorous backward error analysis is possible, thus we can then make accurate predictions about its qualitative behavior. As outlined below it is impossible to give nearby modi1ed equations for standard time-discretization methods as they do

808

K. Matthies / Nonlinear Analysis 52 (2003) 805 – 826

not describe the initial transient from Sobolev to Gevrey class regularity correctly. More eKcient methods with this feature are still to be developed. A scheme of full discretization of (1.7) in space and time is analyzed. In space we use a Galerkin spectral method, which de1nes a system of analytic ordinary di%erential equations for  j UN (x; t) = UN (t) exp(ijx) j6N

with multi-index j ∈ Zd . As we consider periodic boundary conditions, the eigenfunctions for the Galerkin space are cosine and sine functions. Or when using the above complex notation they have the form exp(ijx). We de1ne the approximation space HN in this complex notation, but UN (x) ∈ Rn for all x ∈  and UN ∈ HN :      HN = (1.11) U j exp(ijx)|U j ∈ Cn ; U j = U −j :   d j∈Z ;j6N

The equation in the Galerkin space is given by U˙ N = PN IUN + PN F(UN ); UN (0) = PN U0 ;

(1.12)

where PN denotes the L2 projection on the real Galerkin approximation space HN , which is spanned by eigenfunctions of A. The ordinary di%erential equation (1.12) de1nes a local =ow on HN denoted by SN (t; UN ). Let 0N be the largest eigenvalue of −Id$ restricted to HN . Eq. (1.12) is then discretized in time by some analytic one-step method (e.g. Runge–Kutta) with step size h and order k. The full discretization de1nes a local discrete semigroup on HN which is analytic in U0 and its jth iterate is denoted by h;j N (U0 ); j = 0; 1; 2; : : : : We will couple the step size h and the number of Galerkin modes via h0N1+, ≈ 1

(1.13)

for , ¿ 0. This means, we choose N for given h, such that |h0N1+, − 1| is minimal. Our main theorem for this kind of discretization is the in1nite-dimensional counterpart of (1.6). Theorem 1. Consider Eq. (1.7); where F : Rn → Rn is an entire function. Suppose; s 6 M1 remains inside the that the local semi4ow S(t; U0 ) of initial data U0 with |U0 |Hper ˜ 6 ∞. s (M0 ) for times t with 0 ¡ t ¡ T ball BHper Let the full discretization h; N —with N (h) de8ned by (1.3)—be real analytic in U and be of order p for a complex extended domain. We assume there is C : R+ → R+ ,

K. Matthies / Nonlinear Analysis 52 (2003) 805 – 826

809

which is non-decreasing and is independent of &, such that |h; N (U0 ) − SN (h; U0 )|G&s=2 6 C(|PN U0 |G&s=2 )hp+1

(1.14)

for U0 in a complex extension of a G&s=2 -ball. s (M0 ) Then there exist h∗ ; t ∗ ; T and a modi8ed equation on the H s -ball BHper N UN ; h; N (h)); UN˙ + AUN = F(UN ) + F(

(1.15)

UN (0) = U0 :

(1.16)

This modi8ed equation de8nes a local semi4ow SNh; N (t; U0 ), that is exponentially close to the numerical method for 0 ¡ h 6 h∗ and N (h) de8ned by (1.13):  |h;j N (U0 ) − SNh; N (j · h; U0 )| 6 C(exp(−c1 min(t; t ∗ ) 0N ) + exp(−c2 h−,=(1+,) ))

(1.17)

for times 0 ¡ t = jh ¡ T with  T˜ if T˜ ¡ ∞; T=  O(h−,=(1+,) ; 0N ) if T˜ = ∞: The constants C; c1 ; c2 ; t ∗ are independent of h and N , when choosing N = N (h) as above. The possibly nonlocal perturbation s s × [0; h∗ ] × N → Hper FN : Hper

is uniformly bounded in h → 0; N = N (h) → ∞ for U ∈ BH s (M0 ): N )|H s 6 C1 hp : |F(U

(1.18)

Remark 1. (1) The exponential estimates are fully developed after a transient time t ∗ . But we have to stress; that exponential estimates hold for any 1xed positive time N (2) The correction FN is nonlocal, i.e. to evaluate F(U; h)(x) at a point x ∈  one has N to evaluate F(U; h) in the Sobolev space 1rst. While to evaluate the local nonlinearity F(U )(x) it suKces to compute F(U (x)). (3) The assumptions on h; N are ful1lled, when we are discretizing the equation in time by some implicit discretization schemes (e.g. higher order Radau methods). But the order k, that these methods have as integrators of ordinary di%erential equations cannot be achieved, see Section 2. (4) The constant t ∗ is given by t ∗ =min(T ∗ ; 1=2L), where T ∗ is the maximal Gevrey exponent as in (1.10). L is the Lipschitz constant of F + FN on the ball of radius M1 s in Hper . This is up to an error O(hp ) close to the Lipschitz constant of F. (5) The essential diKculty to consider other equations is to prove Gevrey regularity. Such results hold e.g. for Navier–Stokes equations with periodic boundary conditions [7] and for reaction di%usion equations on the sphere S 2 [2].

810

K. Matthies / Nonlinear Analysis 52 (2003) 805 – 826

(6) For 1xed N and h, one could try to apply backward error analysis techniques of [8]. But it is not directly obvious, that for N → ∞ the equation keeps the form of a semilinear parabolic equation as in (1.15). Furthermore our method also gives the existence of nonautonomous equations (see Section 3), whose solutions exactly coincide with the numerical method, which is important for the regularity results in Section 5. The rest of the paper is organized as follows. We will prove the theorem in Section 4. A main ingredient will be the embedding of the numerical scheme into a nonautonomous equation in Section 3. In the rest of this section we will explain, why standard time-discretization schemes cannot be described by modi1ed equations, which are only small perturbations of the original equation. In Section 2 we give examples of integrators which ful1ll (1.14). In the last section we give some examples, how the backward error results can be used to understand the qualitative behavior of the numerical scheme. 1.2. Analyticity and standard discretization schemes When applying standard discretization schemes, one will run into trouble capturing the transition from regularity of Sobolev type to the Gevrey regularity in (1.10). Functions in the Gevrey class G&s=2 for & ¿ 0 are real analytic in the space variable x, whereas the results of discretization schemes remain often in Sobolev spaces. For example, when we consider the linear heat equation @ U (x; t) − DIU (x; t) = 0 @t

(1.19)

and choose a backwards Euler scheme Uj+1 = (I − hD$)−1 Uj as a pure time discretization with step-size h has after time 1 only a regular solution in H 2=h , whereas solutions of (1.19) are analytic after any arbitrary small time. So we cannot expect, that the numerical solutions ful1ll a partial di%erential equation, which is close to (1.19). However, it is possible to write an exact modi1ed equation in this example: @ ˜ U = log[(I − hD$)−1 ]U˜ : @t But this equation is not a small or even bounded perturbation of (1.19), unless we restrict ourselves to a subspace such that hD$ is small. Extreme convergence results for pure spatial discretization by spectral Galerkin methods are known, but one has to assume that the initial condition U0 is already analytic, see [3]. Again, the initial transient cannot be described by a pure Galerkin method. For a Galerkin approximation of the linear heat Eq. (1.19) @ U = DIU; @t U (0) = PN U0 ; however, the decay of the truncated modes is particularly high. So the error is exp(−0N t)(I − PN )U0 and the error does not in=uence the lower modes, such that

K. Matthies / Nonlinear Analysis 52 (2003) 805 – 826

811

we already have the modi1ed equation. For semilinear equation as in the theorem Gevrey regularity G&s=2 gives also exponential bounds on the error of the Galerkin approximation of a single trajectory V (t), which has the regularity in (1.10)  |V (t) − PN V (t)|H s 6 C exp(− 0N min(t; c)): But then the error has an in=uence on the lower modes and thus the error in initial transient has an in=uence too, which again cannot be captured by the Galerkin method. 2. Time integrators Although assumption (1.14) is standard for ordinary di%erential equations, it does not hold in general for parabolic partial di%erential equations. We give examples of discretization schemes, which ful1ll the approximation assumption (1.14) of the theorem for some p ¿ 0. It assumes the convergence of the method independent of N and independent of the Banach space G&s=2 . In the proof of Theorem 2 we furthermore need estimates when extending HN ∩ BG&s=2 to a complex neighborhood, i.e. we have to complexify G&s=2 for & ¿ 0. Error estimates for strongly A(3) stable Runge–Kutta methods were given by Lubich and Ostermann in [11] and [12] in a general setting of semilinear parabolic equations on some Banach space X d U + AU = F(U ); dt

U (0) = U0 ∈ X  :

(2.1)

Here A is a sectorial operator, X  = D(A ) a fractional power space and f : X  → X a Lipschitz nonlinearity, see [10]. Eq. (1.7) is a semilinear parabolic equation on Gevrey classes, because A is sectorial on G&s=2 by the following lemma and because the nonlinearity is a Lipschitz map from G&s=2 to itself by, e.g. [14, Lemma 2]. Thus we can take initial values in X 0 = X = G&s=2 for  = 0 in the abstract equation. Lemma 2. The operator A = −diag(d1 ; : : : ; dn )$ is sectorial on G&s=2 (; Cn ) with D(A) = G&s=2+1 (; Cn ). Proof. We use De1nition 1.3.1 of Henry [10]. A is densely de1ned on G&s=2 ; as it is G s=2

G s=2

& & u and Auk → v. We de1ned on 1nite Fourier series. To check; that A is closed; let uk → s=2 then have to show u ∈ D(A) and Au = v. The convergence of Auk in G& gives; that uk converges in G&s=2+1 ; hence u ∈ D(A). As A is closed on L2 ; we have by |v − Au|L2 = 0; that the Fourier series of v and Au coincide and thus Au = v in G&s=2 . The operator A is sectorial on L2 (), its resolvent can be estimated M |(0 − A)−1 |L(L2 ; L2 ) 6 M (0) = |0 + 1|

812

K. Matthies / Nonlinear Analysis 52 (2003) 805 – 826

for 0 ∈ {0|5 6 |arg(0 + 1)| 6 } for some 5 ∈ (0; =2) with M depending on 5; see, e.g. [10]. As each Fourier monomial exp(ijx) is an eigenfunction of (0 − A)−1 , we have |(0 − A)−1 exp(ijx)|G&s=2 6 M (0)|exp(ijx)|G&s=2 : Thus for all v ∈ G&s=2 (; Cn ) we get the same resolvent estimate as in L2 |(0 − A)−1 v|G&s=2 6

M |v| s=2 |0 + 1| G&

for 0 in a sector {0|5 6 |arg(0 + 1)| 6 } for some 5 ∈ (0; 2 ). In the setting of abstract semilinear parabolic equations (2.1) there are two di%erent types of error estimates. First there are nonsmooth error estimates for general initial values as in [12]. Uniform estimates for 1nite times can only be given after some transient: The constant in front of the expected order hk will depend on a transient time and tend to in1nity for small times. Whereas in ordinary di%erential equations we get as in (1.14) an estimate of order hk+1 for the 1rst time step. Furthermore the order of convergence is even after the transient usually small. So they are not applicable to assure (1.14). Secondly there are smooth error estimates. For higher order estimates we need smoothness in time starting at t = 0. This cannot be assumed in general for arbitrary initial conditions. Even then the full order of convergence known from the theory of ordinary di%erential equations is not realized in general. We will use these smooth error estimates, but instead of assuming smoothness we use the coupling of 0N and h. We use results of Lubich and Ostermann [11, Theorem 4.1]. Their result holds for a large class of semidiscretization in time (strongly A(3) stable Runge–Kutta methods with internal stage order q and order k for k ¿ q + 1) and it is also true uniformly in N for the full discretization.    h (q+1) (q+1) |h; N (U0 ) − SN (h; U0 )|G&s=2 6 C |UN (0)|G&s=2 + |UN (t)|G&s=2 dt hq+1 ; 0

where UN(q+1) (t) denotes the (q+1)th time derivative, which can be bounded by C0Nq+1 . Using 0N1+, h 6 1 gives |h; N (U0 ) − SN (h; U0 )|G&s=2 6 C(0N h)q+1 = Ch,(q+1)=(1+,) : Hence assumption (1.14) holds for the full discretization with Galerkin discretization in space and some higher order strongly A(3) stable Runge–Kutta methods like higher order Radau methods in time, when they are coupled by 0N1+, h 6 C with appropriate , ¿ 0. For example, to apply Theorem 1 to a 3-stage RadauIIa method with q = 3, we choose , = 12 . Then this yields to the exponential estimates C exp(−ch−1=3 ) and N 6 Ch1=3 with p = 1 in the theorem. |F| 3 For linear parabolic equations with periodic boundary conditions with semigroup S(t; U ) and semidiscretization  of step size h, Lubich and Ostermann showed in [13]

K. Matthies / Nonlinear Analysis 52 (2003) 805 – 826

813

that the full order of convergence p can be achieved for 1nite times t = nh 6 T : |S(nh; U0 ) − n (U0 )| 6 Chp : So perhaps the above error estimates are not optimal, because we do not use the periodic boundary conditions. Convergence results for discretization in space by spectral Galerkin methods can be found e.g. in [18] and [3]. In [3] the Gevrey regularity is used to show exponential convergence of the Galerkin approximation, if one already starts with Gevrey regular initial values. In contrast we consider in this paper general initial values and take also the e%ect of time discretization into account. 3. Embedding into nonautonomous equations An essential tool in the proof of the theorem is the embedding into a nonautonomous equation. We rewrite Eq. (1.7) as @ U = −AU + F(U ); @t s U (0) = U0 ∈ Hper (; Rn )

(3.1)

with A = −diag(d1 $; : : : ; dn $) on  = [0; 2]d with d = 1; 2; 3, periodic boundary conditions and s ¿ d=2. To compare the numerical scheme h; N with equations of type (3.1) we are going to embed h; N into a rapidly forced ordinary di%erential equation
 @ t UN = −AUN + PN F(UN ) + hp G UN ; ; h; N ; @t h UN (0) = PN U0 :

(3.2)

The discretization scheme coincides exactly with the time-h map 7(h; 0; :) of a nonautonomous equation like in (3.2) h;j N (U0 ) = 7(jh; 0; PN U0 ; h; N ): In Proposition 3 we will construct a time-periodic perturbation G of (3.2), such that its evolution 7(t; 0; UN ; h; N ) has the desired property. The nonautonomous ordinary di%erential equation (3.2) can easily be extended to a partial di%erential equation with nonlinearities, which are nonlocal
 @ t U = −AU + PN F(PN U ) + hp G PN U; ; h; N ; @t h U (0) = PN U0 :

(3.3)

We apply a theorem like the Neishtadt theorem for ordinary di%erential equations for our class of equations as in [14]. This will complete the proof of the theorem in Section 4.

814

K. Matthies / Nonlinear Analysis 52 (2003) 805 – 826

We consider the ordinary di%erential equation (1.12) with associated =ow SN (t; U ). As an adaption of [6, Proposition 2.1] and its re1nement in [22] we construct a nonautonomous perturbation G as in (3.2) for our special discretization scheme. In particular, we get uniform estimates on the Gevrey norms for h → 0; N (h) → ∞. Proposition 3. Let the assumptions of Theorem 1 on the nonlinearity F and the full discretization scheme h; N hold: The nonlinearity F : Rn → Rn is entire in U and the discretization h; N is real analytic in U and h for 8xed N and ful8lls (1.14). Furthermore N and h are coupled as in (1.13) by setting 0N1+, h ≈ 1 with , ¿ 0. Then there exist h∗ ¿ 0, a nonautonomous vector 8eld G = G(UN ; 3; h; N ), a nonincreasing continuous function 8 : [0; h∗ ] → (0; ∞] with 8(0) = ∞ and a continuous non-decreasing function r : R+ → R+ such that: (i) G(UN ; 3; h; N ) ∈ HN is de8ned for all real 3 and 0 6 h 6 h∗ and UN ∈ HN with s 6 8(h). |UN |Hper (ii) For 8xed N the perturbation G is C ∞ -smooth in 3 and analytic in h; UN . (iii) The vector 8eld G is periodic in 3 with period 1. (iv) h; N (h; U0 ) = 7(h; 0; PN U0 ; h; N ), where 7 is the evolution of (3.2). (v) |G(UN ; 3; h; N )|G&s=2 6 r(|UN |G&s=2 ) independent of &; 3; h and N . (vi) For DN = HN ∩ BG&s=2 (R) ⊂ Rk the perturbation G can be extended to a complex ;-neighborhood taking distances in the G&s=2 norm: sup

UN ∈DN +;;3∈R

|G(UN ; 3; h; N )|G&s=2 6 B2 :

Proof. The basic idea of the construction is to interpolate for 1xed N between the identity and the numerical scheme  by a family t → 7(t; 0; U ; h; N ) of di%eomorphisms with 7(0; 0; U ; h; N ) = Id; 7(h; 0; U ; h; N ) = h; N (h; U ): An evolution is then de1ned, omitting the other arguments, by 7(t; s) = 7(t; 0) ◦ 7(s; 0)−1 : This can be extended to t ∈ [0; h] by iteration. We repeat the construction of [6]. This interpolation is done in step 1. In step 2 we will construct G, such that −A + PN F + G is the time derivative of 7. The properties (i) – (iv) will be checked in step 3. The estimates in (v) and (vi) will be derived in step 4. Step 1: Interpolation. Choose < : R → [0; 1] a C ∞ cut-o% function analytic in

= 0; 1 with  1 for 6 0; <( ) = 0 for ¿ 1:

K. Matthies / Nonlinear Analysis 52 (2003) 805 – 826

815

Then we de1ne 7 by interpolation for t ∈ [0; h]
t

t  SN (t; UN ) + 1 − < SN (t − h; h; N (h; UN )): (3.4) 7(t; 0; UN ; h; N ):=< h h As SN is only a local =ow in time, we need the restriction |UN |H s 6 8(h), then 7(0; 0; UN ; h; N ) = SN (0; UN ) = UN ;

(3.5)

7(h; 0; UN ; h; N ) = SN (0; h; N (UN )) = h; N (UN ):

(3.6)

We can then extend the de1nition to all t ¿ 0 by
t  7(t; 0; UN ; h; N ):=7 t − h; 0; h;[t=h] (U ); h; N ; N N h

(3.7)

where [ ] denotes the largest integer less than . Thus 7 is an interpolation of (jh; UN ); j ∈ Z. Then 7(t; 0; 7(jh; 0; UN ; h; N ); h; N ) = 7(t + jh; 0; UN ; h; N ) and the curve t → 7(t; 0; UN ; h; N ) is clearly analytic for t=h ∈ N. By the recursion formula above, it is enough to check C ∞ -smoothness at t = h to prove that t → m 7(t; 0; UN ; h; N ) is C ∞ . We denote by D± the mth t-derivative from above (+) and from below (−). Then we get for arbitrary m, when dropping the dependence on h; N : m D+ |t=h 7(t; 0; UN ) m = D+ |t=0 7(t; 0; h; N (UN ))

t

t   m SN (t; h; N (UN )) + 1 − < SN (t − h; h;2 N (UN )) = D+ |t=0 < h h m m |t=0 SN (t; (h; UN )) = D− |t=0 SN (t; (h; UN )) = D+

t

t   m SN (t; UN ) + 1 − < SN (t − h; h; N (UN )) |t=h < = D− h h m |t=h 7(t; 0; UN ): = D−

This proves C ∞ -regularity of 7 in t. Step 2: Construction of G. To construct the perturbation, we rescale time by letting 3 = t=h and writing ˜ 0; UN ; h; N ) = 7(3h; 0; UN ; h; N ): 7(3;

(3.8)

Thus we have for 3 ∈ [0; 1]: ˜ 0; UN ; h; N ) = <(3)SN (h3; UN ) + (1 − <(3))SN (h3 − h; h; N (h; UN )): 7(3; We de1ne ˜ 0; VN ; h; N ) − (−AUN + PN F(UN ))); G(UN ; 3; h; N ):=h−p (h−1 D3 7(3;

(3.9)

816

K. Matthies / Nonlinear Analysis 52 (2003) 805 – 826

where VN = VN (UN ; 3; h; N ) is de1ned implicitly by ˜ 0; VN ; h; N ): UN = 7(3;

(3.10)

This can be solved for VN by the implicit function theorem. 7˜ is analytic in h; VN and smooth in =. The trivial solution is VN = VN (UN ; 3; 0; N ) = UN for h = 0 by (3.5). ˜ 0; VN ; 0; N ) = Id is invertible. Hence we can solve (3.10) for The linearization DVN 7(3; VN = VN (UN ; 3; h; N ) uniformly for 0 6 h 6 h∗ in a ball |UN | 6 8(h) and 7˜ is de1ned. Thus VN will be de1ned on a reduced domain 0 6 3 6 2; 0 6 h 6 h∗ ; |UN | 6 8(h). Step 3: Basic properties of G. First we extend the domain of G to all real 3 by 1-periodic extension to show (iii). Thus we have to check periodicity for 0 6 3 6 2. Let 0 6 3 6 1; h ¿ 0 then we get from (3.7) ˜ 0; h; N (VN (UN ; 3 + 1; h; N )); h; N ) 7(3; ˜ + 1; 0; VN (UN ; 3 + 1; h; N ); h; N ) = UN ; =7(3 ˜ 0)−1 on both sides which implies by applying 7(3; h; N (VN (UN ; 3 + 1; h; N )) = VN (UN ; 3; h; N ): Hence, we have when dropping arguments h and N : hp+1 (G(UN ; 3 + 1) − G(UN ; 3)) ˜ + 1; 0; VN (3 + 1)) − D3 7(3; ˜ 0; VN (3)) = D3 7(3 ˜ 0; h; N VN (3 + 1)) − D3 7(3; ˜ 0; VN (3)) = D3 7(3; ˜ 0; VN (3)) − D3 7(3; ˜ 0; VN (3)) = 0 = D3 7(3; proving periodicity for h ¿ 0. Next we will extend G down to h = 0. Let dj denote various remainder terms real analytic in 0 6 h 6 h∗ ; |UN | 6 8(h) and C ∞ in 3. By the pth order approximation property we have h; N (UN ) = SN (h; UN ) + hp+1 d0 : This implies ˜ UN ) = <(3)SN (h3; UN ) + (1 − <(3))SN (h3 − h; h; N (h; UN )) 7(3; = <(3)SN (h3; UN ) + (1 − <(3))SN (h3 − h; SN (h; UN ))   + (1 − <(3)) SN (h3 − h; h; N (UN )) − SN (h3 − h; SN (h; UN )) = SN (h3; UN ) + hp+1 d1 :

(3.11)

˜ 0; VN (3; UN )) = UN , this yields to When replacing UN by VN and using 7(3; VN (3; UN ) = SN (−h3; UN ) + hp+1 d2 :

(3.12)

K. Matthies / Nonlinear Analysis 52 (2003) 805 – 826

817

From (3.11), (3.12) we obtain ˜ VN (3; UN )) = −ASN (h3; VN ) + PN F(SN (h3; VN )) + h−1 D3 hp+1 d1 h−1 D3 7(3; = −AUN + PN F(UN ) + hp+1 d3 + hp D3 d1 ; where hp+1 d3 = (−A + PN F)(SN (h3; SN (−h3; UN ) + hp+1 d2 )) − (−AUN + PN F(UN )):

(3.13)

Setting G = hd3 + D3 d1 gives the di%erentiability requirements in (ii) and by continuous extension to h = 0 also (iii) can be proved. Then G is de1ned on [0; h∗ ] for 3 ∈ R and |UN | 6 8(h), thus proving (i). The interpolation property (iv) is ful1lled by construction of 7. Step 4: Estimates uniform in h; N (h). Properties (v) and (vi) are the essential new ingredients to [6, Proposition 2.1]. We show 1rst the estimates for G = hd3 + D3 d1 in (v). Starting with D3 d1 , we obtain by (3.11): hp+1 D3 d1 = D3 [(1 − <(3))(SN (h3 − h; h; N (UN )) − SN (h3 − h; SN (h; UN )))] = −< (3)(SN (h3 − h; h; N (UN )) − SN (h3 − h; SN (h; UN ))) + (1 − <(3))h(−ASN (h3 − h; h; N (UN )) + PN F(SN (h3 − h; h; N (UN ))) + ASN (h3 − h; SN (h; UN )) − PN F(SN (h3 − h; SN (h; UN )))): Using the mean value theorem we get, where K = {V ∈ > h; N (UN ) + (1 − >)SN (h; UN ); > ∈ [0; 1]}; hp+1 |D3 d1 |G&s=2

 6 |h; N (UN ) − SN (h; UN )|G&s=2 |< (3)|sup|DSN (h3 − h; V )|L(G&s=2 ;G&s=2 ) K

 +h|1 − <(3)|sup|D(ASN (h3 − h; V ) + PN F(SN (h3 − h; V ))|L(G&s=2 ;G&s=2 ) : K

So, we need to estimate |DV SN (h3 − h; V )|L(G&s=2 ;G&s=2 ) . To do this, we use the variational equation, which is ful1lled by W (t) = DV SN (t; V (t)) Wt + APN W = DV PN F(PN V (t))W: Thus, we get |DV SN (h3 − h; V )|G&s=2 6 |exp(−A(h3 − h))PN |G&s=2

818

K. Matthies / Nonlinear Analysis 52 (2003) 805 – 826

 +|

h3−h

0

exp(−A(h3 − h − #))

×DV PN F(PN V (#)) d#|G&s=2 :

(3.14)

To bound this expression, we 1rst have to analyze  t V (t) = exp(−At)PN V0 + exp(−A(t − #))PN F(V (#)) d#: 0

The linear part can be bounded using h0N1+, ≈ 1. PN F is locally Lipschitz on HN with constants L(VN ) uniformly bounded for N → ∞ in the G&s=2 -norm. Thus V (t) is bounded by a constant C for |t| 6 h and N 6 N (h) by the Gronwall inequality  t |V (t)|G&s=2 6 C + C|PN F(V (#))| d# 0

 6C +

0

t

CL|V (#)|G&s=2 d# 6 C:

Then (3.14) is uniformly bounded when using h0N1+, ≈ 1. Using the bound (1.14) for |h; N (UN ) − SN (h; UN )|G&s=2 we have the estimate hp+1 |D3 d1 |G&s=2 6 (c + hC0N )C(|UN |H&s=2 )hp+1 6 CC(|UN |G&s=2 )hp+1

(3.15)

as 0N h ¡ 0N1+, h ≈ 1 is bounded. Next, we deal with hp+1 d3 = (−A + PN F)(SN (h3; SN (−h3; UN ) + hp+1 d2 )) − (−AUN + PN F(UN )); hence we get an estimate depending on d2 , where K = {V ∈ SN (−h3; UN ) + >hp+1 d2 ; > ∈ [0; 1]} hp+1 |d3 |G&s=2 6 sup|D(ASN (h3; V ) + PN F(SN (h3; V ))|L(G&s=2 ;G&s=2 ) hp+1 |d2 |G&s=2 K

6 |0N + L|Chp+1 |d2 |G&s=2 6 Chp |d2 |G&s=2 : It remains to estimate d2 given by (3.12) hp+1 d2 (UN ; 3) = VN (3; UN ) − SN (−h3; UN ): From (3.11) we get ˜ UN ) − hp+1 d1 (UN )) = UN : SN (−h3; 7(3; Then replacing UN by VN (UN ) yields to SN (−h3; UN − hp+1 d1 (VN )) = VN :

(3.16)

K. Matthies / Nonlinear Analysis 52 (2003) 805 – 826

819

Hence, hp+1 d2 (UN ; 3) = SN (−h3; UN − d1 (VN )) − SN (−h3; UN ) with hp+1 |d2 (UN ; 3)|G&s=2 6 Chp+1 |d1 (VN ; 3)|G&s=2 : ˜ VN (UN ; 3); h) = UN , where 7˜ is a near-identity di%eoVN (UN ; 3) is de1ned by 7(3; morphism. So VN remains uniformly bounded per construction. Hence we can estimate d1 in the same way as D3 d1 above, and we get hp+1 |d1 (VN ; 3)|G&s=2 6 Chp+1 . Thus

hp+1 |d3 (UN ; 3)|G&s=2 6 Chp for UN ∈ G&s=2 bounded. Hence for G = hd3 + D3 d1 we obtain |G(UN ; 3; h; N )|G&s=2 6 r(|UN |G&s=2 )

and (v) is shown. It remains to check (vi). We need uniform bounds when extending G to a complex domain DN + ;. First we check that SN (h3; UN ) does not leave DN + ;0 for 3 ∈ [0; 1] and UN ∈ DN + ;. For the imaginary part of UN we get the following equation: d Im(UN ) = −A Im(UN ) + Im(F(UN )): dt As Im(F(UN )) = 0 for Im(UN ) = 0, there is a linear factor of Im(UN ) in Im(F(UN )). Thus we can estimate using the Gronwall inequality  h3 −h3A Im(UN (0)) + e−(h3− )A Im(F(UN ( ))) d |G&s=2 |Im(SN (h3; UN ))|G&s=2 6 |e  6;+

0

0

h3

C Im(UN ( )) d 6 ;eh3C 6 ;eh

∗

C

= ;0

and SN (h3; UN ) does not leave a slightly enlarged complex extended domain. Now all estimates for (v) still hold in the extended domain, essentially because the error estimate |h; N (U0 ) − SN (h; U0 )|G&s=2 6 C(|UN (0)|G&s=2 )hp+1 given in (1.14) does extend to the complex extended domain by assumption. This completes the proof of the proposition.

4. Proof of Theorem 1 We are now in the position to prove Theorem 1. By Proposition 3 the numerical scheme h; N coincides exactly with the time-h map of the ordinary di%erential equation
 t U˙ N + AUN = PN F(UN ) + hp G UN ; ; h; N : h We extend this equation to a partial di%erential equation 
t U˜˙ + AU˜ = PN F(PN U˜ ) + hp G PN U˜ ; ; h; N (h) h

(4.1)

820

K. Matthies / Nonlinear Analysis 52 (2003) 805 – 826

˜ 0; :), such that the numerical scheme and S(t; ˜ 0; PN U0 ) with a local time evolution S(t; coincide. We want to apply [14] to (4.1). We assume two groups of hypotheses. The 1rst group implies spatial regularity of the solutions and the second ensures analyticity in the nonlinearities. We consider abstract equations
t  @ U − DIU = F(U ) + hp G U; ; h : @t h The 1rst hypothesis (H.1) is to assure, that the solutions are Gevrey regular s with s ¿ d=2 U0 ∈ Hper

F(:); G(:; t) : BG&s=2 (R; 0) → G&s=2 di%erentiable in U s F  (U ); G  (U ) ∈ L(L2 (; Rn ); L2 (; Rn )) for U bounded in Hper

G continuous in t; h F(:); G(:) bounded with constants not depending on &; t :

(H.1)

|F(U )|G&s=2 6 Cs a(Cs |U |G&s=2 ) |G(U; t)|G&s=2 6 Cs b(Cs |U |G&s=2 ) with a; b monotone increasing on [0; R=Cs ) The second hypothesis (H.2) describes the needed analytic properties of F and G, when they are restricted to Galerkin approximation spaces. PN F; PN G : DN → HN is real analytic for all N ∈ N sup

UN ∈DN +;

sup

|PN F(UN )|G&s=2 6 B1

UN ∈DN +;; t∈R

(H.2)

|PN G(UN ; t; h)|G&s=2 6 B2

Then the following version of [14, Theorem 2 with (3.17), (3.22), (3.23), Remark 9] holds: Theorem 2. Assume that hypotheses (H.1) and (H.2) hold. Let G be periodic in =t=h with period 1 and 8x the perturbation order p ¿ 0. Suppose for initial values U0 with |U0 |H s 6 M1 ; that the forward orbit remains bounded; i.e. |S(U0 ; t)|H s 6 M0 ¡ R for t ∈ [0; T ). Then there exists a real analytic and time-periodic change of coordinates for 0 ¡ h ¡ h1
t  U = V + hW V; ; h (4.2) h s s . Its image W (Hper ; t=h; h) is 8nite with the following properties: W is bounded on Hper dimensional for 8xed h and W (:; 0; :) = 0. The transformed nonautonomous terms are exponentially small after a transient; but the equation may contain additional

K. Matthies / Nonlinear Analysis 52 (2003) 805 – 826

821

nonlocal terms FN = (fN 1 ; : : : ; fN n )


 
@ N h) +  V; t ; h + b1 V; t ; N; h V − DIV = F(V ) + F(V; h h @t
t  + b2 V; ; N; h : h

(4.3)

The index of the approximation space N (h) is given by 0N1+, h ≈ 1. Then the correction terms are given by
 t  b1 = (I − PN ) F(V + hW (PN V )) − F(V ) + hp G V + hPN W (V ); ; h ; h (4.4)  b2 = I + h

−1  @ PN F(V + hW (PN V )) − F(PN V + hW (PN V )) W (PN V ) @PN V

t  t  (4.5) + hp G V + hW (PN V ); ; h − hp G PN V + hW (PN V ); ; h h h

or can be estimated by sup

s ¡M0 |V |Hper

sup

N )|H s 6 c3 hp ; |F(V per

s ¡M1 |V |Hper

p −,=(1+,) s 6 c2 h exp(−c1 h |(V )|Hper ):

Here the constants c1 ; c2 ; c3 ; h1 only depend on B1 and B2 . Furthermore FN and  are smooth in h. We apply this theorem to (4.1). So we only have to replace F(:) by PN F(PN :), where N is given by the discretization scheme. First we check the assumptions of the theorem. By [14, Remark 6] F and therefore PN F ful1ll (H.1) and (H.2). The perturbation term G satis1es (H.1) and (H.2) by Proposition 3. As the nonlinearities in (4.1) have their range in HN only, the form (4.3) in Theorem 2 can be improved. We obtain b1 ≡ b2 ≡ 0 when using the same N for the discretization and Theorem 2, see (4.4) and (4.5). The correction b1 is zero, because PN F(PN :) and G have their range in HN = ker(I − PN ). The other correction b2 is zero, because components outside HN do not change the values of PN F(PN :) and G. So we get by the above theorem an averaged equation
 N h) +  V; t ; h; N (h) ; V˙ + AV = PN F(PN V ) + F(V; h V (0) = PN U0 with local evolution S˜V (t; 0; PN U0 ) and strong estimates on the nonautonomous part: sup

p −,=(1+,) s 6 c2 h exp(−c1 h ||Hper ):

s (M1 ;0)×R+ ×[0;h1 ) V (0)∈BHper

822

K. Matthies / Nonlinear Analysis 52 (2003) 805 – 826

We also have S˜V (nh; 0; PN U0 ) = h;n N (U0 );

(4.6)

since the coordinate change W is zero on multiples of the period h. The modi1ed Eq. (1.15) given in Theorem 1 is then N UN ; h; N (h)); UN˙ + AUN = F(UN ) + F( UN (0) = U0 ;

(4.7) (4.8)

which de1nes a local semi=ow SNh; N (t; U0 ) depending on N and h. The estimates on FN in (1.18) then also follow from Theorem 2. Estimating the di%erence between UN (t) = SNh; N (t; U0 ) and V (t) = S˜V (t; 0; PN U0 ) will prove the desired exponential estimate (1.17) because of (4.6). We use PN V ( ) = V ( ) and that exp(−A(t − ))U ( ) ∈ Gts=2 , if U ( ) ∈ G s=2 . The time T ∗ denotes the maximal Gevrey exponent as in (1.10). |UN (t) − V (t)|H s =|SNh; N (t; U0 ) − S˜V (t; 0; PN U0 )|H s  t   N UN ( )) =|exp(−At)(Id − PN )U0 |H s +  exp(−A(t − )) F(UN ( )) + F( 0


   N ( )) −  V ( ); ; h; N (h) − PN F(PN V ( )) − F(V d  h Hs adding inside the integral 0 = −PN F(PN UN ( )) + PN F(PN UN ( ))   t     N N s exp(−0N +1 t) + 6 |U0 |Hper  exp(−A(t − )){F(U ( )) − PN F(PN U ( ))} d   +L

0

0

t

Hs

|V ( ) − UN ( )|H s d + tc2 hp exp(−c1 h−,=(1+,) )

 s exp(−0N +1 t) + t2B1 exp(− 6 |U0 |Hper 0N +1 min(t; T ∗ )) + tc2 hp exp(−c1 h−,=(1+,) ) + L



t

0

|V ( ) − UN ( )|H s d :

First we consider small times t 6 min(1=2L; T ∗ ) only, where L is the Lipschitz constant N + PN F(PN :). This Lipschitz constant is close to the Lipschitz constant of F of F(:) N is analytic on a complex s (M0 ). As F as a Nemitski operator on a suitable ball BHper p N extension of a bounded domain in HN , we get DF 6 Ch using the Cauchy estimate, as in [14]. The last inequality implies for t 6 t ∗ = min(1=2L; T ∗ ) y(t) = max |UN ( ) − V ( )|H s 6 C(t) 06 6t

 s exp(−0N +1 t) + t2B1 exp(− := 2(|U0 |Hper 0N +1 min(t; T ∗ ))

K. Matthies / Nonlinear Analysis 52 (2003) 805 – 826

823

+ tc2 hp exp(−c1 h−,=(1+,) ))  6 C(exp(− 0N +1 t) + exp(−c1 h−,=(1+,) )): For large times t ¿ t ∗ = min(1=2L; T ∗ ), we use a simple Gronwall inequality and get exponential estimates: |UN (t) − V (t)|H s 6 C(t ∗ )(exp(L(t − t ∗ )))  6 C(exp(− 0N +1 t ∗ ) + exp(−c1 h−,=(1+,) )) exp(L(t − t ∗ )): This (1.17) on a time-scale O(h−,=(1+,) ;  implies the desired exponential estimate ∗ 0N (h) ), where the constants C; c1 ; c2 ; t do not depend on h and N = N (h). The theorem is proved. 5. Global errors In this section, we use the modi1ed equation to describe some aspects of the long-time behavior of the analyzed numerical scheme. Several dynamical features of semilinear parabolic equations are replicated by most numerical methods usually with an error of the order h, for a review see [18]. Using the modi1ed equation we can show that these features are replicated with an error of exponential small order—when compared with the modi1ed equation. Furthermore the exact embedding (3.3) gives an analyticity result. Proposition 4. The numerical solutions h;j N (U0 ) are analytic in space uniformly in

s=2 h; N (h). They have the Gevrey regularity as in (1.10): h;j N (U0 ) ∈ Gmin( jh; T ∗ ) as long j as |h; N (U0 )|H s remains bounded.

Proof. The numerical method h;n N (U0 ) coincides with the solution of (3.3) with initial conditions PN (U0 ). This di%erential equation ful1lls (H.1) and thus has Gevrey regular solutions with (1.10) for small times; see [5;14]. Furthermore as long as the Sobolev norm |h;j N (U0 )|H s remains bounded; the argument can be iterated and |h;j N (U0 )|G&s=2 is bounded too. Thus the solutions are analytic in x. For the geometric analysis, consider for example an equilibrium U0 of the semigroup de1ned by (1.7), i.e. S(t; U0 ) = U0 ; ∀t ¿ 0. It is supposed to be a hyperbolic, thus for any 1xed time t ¿ 0 the spectrum & of the linearization DS(t; U0 ) at U0 splits into two parts by the unit circle: there exists C ∈ (0; 1) such that & = &s ∪ &u &s = {0 ∈ &0| ¡ C}

and

&u = {0 ∈ &0| ¿ C−1 }:

Then in [18, Chapter 4] it is shown, that several numerical schemes with step size h have a hyperbolic equilibrium Uh with |Uh − U0 | 6 Ch. The respective phase portraits are close of order h, too. We compare the numerical scheme with the modi1ed equation. We can show, that our discretization scheme h; N replicates the dynamics of a slightly modi1ed partial di%erential equation near an equilibrium, with an only exponentially small error.

824

K. Matthies / Nonlinear Analysis 52 (2003) 805 – 826

Proposition 5. Suppose the original Eq. (1.7) possesses a hyperbolic equilibrium at U0 . Then the numerical scheme h; N with a coupling like in (1.13) with , = 12 has a hyperbolic 8xed point: h; N (Uh ) = Uh ; which is exponentially close to a hyperbolic equilibrium UN h of the modi8ed Eq. (1.15) with: |Uh − UN h |H s 6 C exp(−ch−1=3 );

(5.1)

|UN h − U0 |H s 6 Ch:

(5.2)

Moreover the local unstable manifolds of Uh and UN h are exponentially close on BH s (8; U0 ) and their local stable manifolds are exponentially close on BH s (8; 0) in∗ tersected with the image of W s; loc (Uh ) under the regularizing map h;[TN =h] (:). Proof. This proposition is given for the comparison of rapidly forced equations and averaged equations in [14; Proposition 11]: The solution for projected initial conditions of the rapidly forced equation coincides exactly with the numerical scheme at multiples of the step size h; this also gives the needed C 1 error estimates for the persistence of invariant manifolds. The averaged equation is replaced by the modi1ed Eq. (1.15). After having established the persistence of the equilibrium and the 1xed point; which is suspended to a periodic orbit of period h; it is then not diKcult to show exponential closeness using (1.17). Hence the special full discretization scheme, discussed in this paper replicates the dynamical behavior near hyperbolic equilibria with an exponentially small error comparing to a modi1ed equation. Using the modi1ed equation more geometric theory of the considered discretization scheme can be developed, e.g. the analysis of the in=uence of rapid forcing on homoclinic orbits [14, Proposition 14] holds in the same way for the discretization scheme. Homoclinic orbits are codimension 1 phenomenon in dissipative dynamical systems with continuous time. For systems with discrete time, a homoclinic orbit can become transversal—the homoclinic orbit ‘splits’. Thus it can persist in open sets of parameter space. Here we show, that for our discretization this set is in fact at most exponentially small in h. Proposition 6. Suppose the original Eq. (1.7) depends on a real parameter 0 and possesses a hyperbolic equilibrium at U0 and generic homoclinic orbit to U0 for parameter 00 . Then the numerical scheme h; N with a coupling like in (1.13) with , = 12 has a homoclinic orbit to h; N (Uh ) = Uh only in an exponentially small wedge of in parameter space (0; h). Furthermore the size of the splitting near the local stable manifold of hyperbolic 8xed point can be estimated: d(h) 6 C exp(−ch−1=3 )

(5.3)

Proof. The existence of the homoclinic orbit for the discretized system; which is a small perturbation of the original equation; is due to the nondegeneracy and transvers × R. The result on exponential small splittings sality in the extended phase space Hper follows with Proposition 5. The local stable and unstable manifolds are exponentially

K. Matthies / Nonlinear Analysis 52 (2003) 805 – 826

825

close to those of the averaged system; where the manifolds do not split at all. When the local unstable manifold evolves forward in time; it only needs some 1nite time to reach the local stable manifold near the homoclinic orbit. In this 1nite time the forward orbit of the local unstable manifolds of UN h and Uh remain exponentially close. So the stable and unstable manifolds can only intersect; when the distance between the local stable and the forward orbit of the local unstable manifold of the averaged system is exponentially small. This holds only in an exponentially small wedge in parameter space due to transverse unfolding. The size of the splitting can be at most exponentially small due to the shadowing of exponential order of stable and unstable manifold of the discretized system by the continuous time modi1ed system; for which there is no splitting at all. A completely rigorous proof is given in similar setting in [14; Proposition 14]. Using the modi1ed equation one should be also able to analyze the in=uence of discretization on other invariant objects like periodic orbits, inertial manifolds, attractors, etc. more accurately. Acknowledgements I would like to thank my supervisor Prof. Dr. Bernold Fiedler for discussions and encouragement. Furthermore I am grateful to Prof. Dr. Christian Lubich and Dr. Claudia Wul% for inspiring discussions. This work was supported by the DFG-Schwerpunktprogramm “Analysis, Ergodentheorie und eKziente Simulation dynamischer Systeme.” References [1] G. Benettin, A. Giogilli, On the Hamiltonian interpolation of near to the identity symplectic mappings with application to symplectic integration algorithms, J. Stat. Phys. 74 (1994) 1117–1143. [2] C. Cao, M. Rammaha, E. Titi, Gevrey regularity for nonlinear analytic parabolic equations on the sphere, J. Dyn. Di%erential Equations 12 (2000) 411–433. [3] A. Doelman, E. Titi, Regularity of solutions and the convergence of the Galerkin method in the Ginzburg-Landau equation, Numer. Funct. Anal. Opt. 14 (1993) 299–321. [4] K. Eriksson, D. Estep, P. Hansbo, C. Johnson, Introduction to adaptive methods for di%erential equations, in: A. Iserles (Ed.), Acta Numerica 1995, Cambridge University Press, Cambridge, 1995, pp. 105–158. [5] A. Ferrari, E. Titi, Gevrey regularity for nonlinear analytic parabolic equations, Commun. Partial Di%erential Equations 23 (1998) 1–16. [6] B. Fiedler, J. Scheurle, Discretization of homoclinic orbits, rapid forcing and “invisible” Chaos, Mem. AMS 570 (1996). [7] C. Foias, R. Temam, Gevrey class regularity for the solution of the Navier Stokes Equations, J. Funct. Anal. 87 (1989) 359–369. [8] E. Hairer, C. Lubich, The life-span of backward error analysis for numerical integrators, Numer. Math. 76 (1997) 441–462. [9] E. Hairer, C. Lubich, Asymptotic expansions and backward error analysis for numerical integrators, in: R. de la Llave, L.R. Petzold, J. Lorenz (Eds.), Dynamics of Algorithms, IMA, Vol. 118, Springer, Berlin, 2000, pp. 91–106. [10] D. Henry, Geometric Theory of Semilinear Parabolic Equations, Lecture Notes in Mathematics, Vol. 840, Springer, New York, 1983.

826

K. Matthies / Nonlinear Analysis 52 (2003) 805 – 826

[11] C. Lubich, A. Ostermann, Runge–Kutta methods for parabolic equations and convolution quadrature, Math. Comput. 60 (1993) 105–131. [12] C. Lubich, A. Ostermann, Interior estimates for time discretization of parabolic equations, Appl. Numer. Math. 18 (1995) 241–251. [13] C. Lubich, A. Ostermann, Runge–Kutta time discretization of reaction–di%usion and Navier–Stokes equations: nonsmooth-data error estimates and long time behaviour, Appl. Numer. Math. 22 (1996) 279–292. [14] K. Matthies, Time-averaging under fast periodic forcing of parabolic partial di%erential equations: exponential estimates, J. Di%erential Equations 174 (2001) 133–180. [15] A. Neishtadt, On the separation of motions in systems with rapidly rotating phase, J. Appl. Math. Mech. 48 (1984) 134–139. [16] S. Reich, Backward error analysis for numerical integrators, SIAM Numer. Anal. 36 (1999) 1549–1570. [17] J.M. Sanz-Serna, Symplectic integrators for Hamiltonian problems: an overview, Acta Numer. 1 (1992) 243–286. [18] A. Stuart, Perturbation theory for in1nite dimensional dynamical systems, in: M. Ainsworth, J. Levesley, W.A. Light, M. Merletta (Eds.), Theory and Numerics of Ordinary and Partial Di%erential Equations, Advances in Numerical Analysis, Vol. IV, Clarendon Press, Oxford, 1995, pp. 181–290. [19] R. Temam, In1nite-Dimensional Dynamical Dystems in Mechanics and Physics, Applied Mathematical Sciences, Vol. 68, Springer, New York, 1988. [20] L.N. Trefethen, Group velocity in 1nite di%erence schemes, SIAM Rev. 24 (1982) 113–136. [21] R. Warming, B. Hyett, The modi1ed equation approach to the stability and accuracy analysis of 1nite-di%erence methods, J. Comput. Phys. 14 (1974) 159–179. [22] Y. Zou, W. Beyn, Invariant manifolds for nonautonomous systems with applications to one-step methods, J. Dyn. Di%erential Equations 10 (1998) 379–407.