NUMERICAL ALGORITHMS FOR STATIONARY STATISTICAL PROPERTIES OF DISSIPATIVE DYNAMICAL SYSTEMS

. It is well-known that physical laws for large chaotic dynamical systems are revealed statistically. The main concern of this manuscript is nu- merical methods for dissipative chaotic inﬁnite dimensional dynamical systems that are able to capture the stationary statistical properties of the underlying dynamical systems. We ﬁrst survey results on temporal and spatial approximations that enjoy the desired properties. We then present a new result on fully discretized approximations of inﬁnite dimensional dissipative chaotic dynamical systems that are able to capture asymptotically the stationary statistical properties. The main ingredients in ensuring the convergence of the long time statistical properties of the numerical schemes are: (1) uniform dissipativity of the scheme in the sense that the union of the global attractors of the numerical approximations is pre-compact in the phase space; (2) convergence of the solutions of the numerical scheme to the solution of the continuous system on the unit time interval [0 , 1] modulo an initial layer, uniformly with respect to initial data from the union of the global attractors. The two conditions are reminiscent of the Lax equivalence theorem where stability and consistency are needed for the convergence of a numerical scheme. Applications to the complex Ginzburg-Landau equation and the two-dimensional Navier-Stokes equations in a periodic box are discussed.

1. Introduction. The long-time dynamics of many infinite-dimensional dynamical systems are very complex with abundant chaotic/turbulent behaviors. The observed complex behaviors are not necessarily related to the possible loss of regularity of the solution (say the three dimensional Navier-Stokes equations). Even the simple logistic map, the Lorenz 63 and the Lorenz 96 models possess intrinsic chaotic behavior which renders long-time approximation of a generic single trajectory extremely difficult. On the other hand, it is well-known that statistical properties of these kind of systems are much more important, physically relevant and stable than single trajectories [11,19,21,27,28,43]. Indeed, much of the classical turbulence theories, such as the famous Kolmogorov U 3 L scaling law of the energy dissipation 4600 XIAOMING WANG rate per unit mass as well as the Kolmogorov k − 5 3 energy spectrum in the inertial range in three dimensional homogeneous isotropic turbulence, are presented in the statistical forms [28,11]. Therefore, for complex physical processes, due to the intrinsic stochasticity, it is necessary to consider statistical properties (averaged quantities) of the system instead of properties of individual orbit (see for instance [28,11,43,27,21]).
For a given abstract autonomous continuous in time dynamical system determined by a semi-group {S(t), t ≥ 0} on a separable Banach space H, we recall that if the system reaches a statistical equilibrium in the sense that the statistics are time independent (stationary statistical properties), the probability measure µ on H that describes the stationary statistical properties can be characterized via either the strong (pull-back) or weak (push-forward) formulation [11,21,27,43,44].
where B(H) represents the σ-algebra of all Borel sets on H. Equivalently, the invariant measure µ can be characterized through the following push-forward weak invariance formulation for all bounded continuous test functionals Φ. Invariant measure (stationary statistical solution) for a discrete dynamical system generated by a map S discrete on a Banach space H is defined in a similar fashion with the continuous time t replaced by discrete time n = 0, 1, 2, . . . .
A closely related object associated with the long-time behavior of a dynamical system is the global attractor which we recall for convenience [11,14,38].
Here, dist H denotes the Hausdorff semi-distance in H between two subsets which is defined as where · H = · denotes the norm on H.
The global attractor for a discrete dynamical system induced by a map S discrete on a Banach space H is defined in a similar fashion with the continuous time t replaced by discrete time n = 0, 1, 2, . . . .
A dynamical system is called dissipative if it possesses a global attractor. It is easy to see, thanks to the invariance and the attracting property, that the global attractor, when it exists, is unique [14,38]. The reader is cautioned that our definition of dissipativity may be slightly different (weaker) from the traditional notation [14,38].
Notice that the global attractor is a set in the phase space. Knowledge of the global attractor only usually provides very little information on the dynamics. On the other hand, knowledge of an invariant measure would allow us to calculate various statistical quantities such as the moments.
We are usually interested in H Φ(u) dµ(u) (statistical average) for various test functionals Φ. These test functionals are also called observables in physics literature. One approach to estimate these observables is to estimate the invariant measure µ directly. This is the so-called directly approach [33]. In the finite dimensional case, one can try to approximate the probability density function (pdf) p associated with the invariant measure by solving the Liouville equation [21] ∂ ∂t p(u, t) + ∇ · (F (u)p(u, t)) = 0 (6) where the forcing term F (u) defines the dynamical system in the sense that D dt S(t)u = F (S(t)u). However, computing the invariant measure (or the associated pdf in the finite dimensional case) is usually very difficult and costly if the spatial dimension is high. One of the commonly used alternative methods in calculating the statistical quantity is to substitute spatial average by long time average under Boltzmann's assumption of ergodicity ( [11,21,27,44]) This is usually termed indirect method. Although the above relationship is true for each ergodic invariant measure µ and almost all initial data with respect to µ, the relationship is in general false for non-ergodic invariant measure since the long time average which exists for almost all initial data (with respect to the given invariant measure) may depend on the initial data and hence may not be a constant (the spatial average) ( [21,44]). One way to circumvent this difficulty is to replace the long time limit by Banach (generalized) limits ( [11], [22], section 4.2) which are bounded linear functionals on the space of bounded functions that agree with the usual long time limit on those functions whenever the long time limit exists. One may show via the so-called Bogliubov-Krylov argument that these generalized long time averages over trajectory lead to invariant measures (may depend on the chosen Banach limit and initial datum u) of the system for appropriate dissipative dynamical systems, and the spatial and temporal averages are equivalent (see for instance [11]  It is usually impossible to derive analytical formula for long-time statistical properties for most of the physically interesting systems. Therefore, we need to resort to numerical methods in generic situations. Even under the ergodicity assumption, it is not at all clear that classical numerical schemes which provide accurate approximation on finite time intervals will remain meaningful for stationary statistical properties (long time properties) since small error will be amplified and may accumulate over long time. (A noticeable exception is when the underlying dynamics is asymptotically stable where statistical approach is not necessary since there is no chaos. See for instance [12,15,20].) Indeed, let S k be the solution operator of a one-step scheme with time step k = ∆t, and assume that the scheme is of order m so that the following type of error estimate holds S(nk)u − S n k u H ≤ C exp(αnk)k m where C > 0, α are constants. We then have on a time interval [0, T ], an a priori error bound on the long time average of the order of k m exp(αT )−exp(αk) exp(αk)−1 which diverges as T approaches infinity for a positive α. The positivity of α follows from the existence of a positive Lyapunov exponent (the existence of chaotic behavior). Even if the long time averages of the scheme converge, the limit is not necessarily that of the original dynamical system under approximation since the two limits where T → ∞ and the time-step approaches zero are not commutable in general. Extra work is needed to verify the limit is the desired one.
Therefore, it is of great importance and a challenge to design and analyze numerical methods that are able to capture stationary statistical properties of infinite dimensional complex dynamical system. We will focus on dissipative systems for simplicity. Addressing issues like this is of great importance in the numerical study of climate change since the climate is customary estimated via long time integration of the system.
We will demonstrate below that the central idea in the design of numerical algorithms that are capable of capturing the long-time statistical properties is the faithfulness to the underlying (dissipative) dynamical system. More specifically, we will illustrate below that the key ingredient in algorithms that are able to capture the long-time statistical properties is the uniform dissipativity and the uniform convergence on the unit time interval (modulo an initial layer) for initial data coming out of a compact subset of the phase space. It is easy to see that the assumptions are natural. Since the underlying dynamical system is dissipative, it is natural to require that the numerical scheme inherit the dissipativity of the continuous in time system so that the scheme is uniformly dissipative (for small time steps). The uniform convergence of the numerical scheme for initial data from the global attractor on the unit time interval is also expected for most reasonable numerical schemes. Another ingredient in the theory, satisfied by almost all well behaved continuous in time dissipative systems, is the strong continuity of the underlying dynamical system on the unit time interval uniformly with respect to initial data from the union of the global attractors. Once the desired natural conditions are discovered, the proof of the main result is relatively straightforward.
In practice, the design of numerical schemes that are efficient and uniformly dissipative is the primary challenge. The verification of uniform convergence on the unit time interval, although tedious in many cases, is quite standard in the sense that it is usually a refinement of the standard convergence analysis of the numerical scheme under investigation. These points will be illustrated via two examples later.
These two requirements are reminiscent of the conditions in the Lax-Richtmyer equivalence theorem [23,22]. Here the uniform dissipativity plays the role of stability. The uniform convergence on the unit time interval for appropriate initial data can be interpreted as consistency. This is because a point in the phase space for statistical solutions is analogous to a point in the physical space for standard solutions. A statistical solution that emanates from a single point in the phase space is a solution to the underlying dynamical system (generated by PDEs for instance) in the usual sense. And the unit time interval (or any finite time interval) for long-time behavior plays the role of a single time step in classical numerical scheme.
There are abundant works on how to approximate the global attractors or some other invariant structures for various dissipative systems. See for instance [31,32,16,17,24,25,29,34, 51] among many others. There has been a lot of work on temporal approximation of dissipative dynamical systems such as the two dimensional incompressible Navier-Stokes system and the one-dimensional Kuramoto-Sivashinsky equation (see [12,16,18,31,32,41,9,10] among others) that preserves the dissipative property in some sense. Despite the importance of long-time statistical properties, there are relatively few results on numerical schemes designed to capture these features (see [30,33,42] for the case of Hamiltonian ODEs and related topics, and [3,4,13,39,40,48,49,50] for various case studies for dissipative PDEs and for general theory on temporal and spatial discretization that captures the long time statistical properties asymptotically). The primary new contribution of this manuscript is a new set of criteria on fully discretized one-step schemes that guarantees the convergence of the long time statistical properties of the algorithms under investigation. As a by-product, we also present an abstract result on criteria that ensure the convergence of the global attractors of the numerical schemes.
The rest of the manuscript is organized in the follows. In section 2, we present our abstract theory on numerical algorithms that are able to reproduce the long time statistical behavior of the underlying dynamical system asymptotically. Old results on temporal and spatial approximations, as well as a new result on fully discretized approximation will be presented. We illustrate the application of the abstract results on complex Ginzburg-Landau equation and the two-dimensional Navier-Stokes equation in a periodic box in section 3. Concluding remarks will be presented in the last section.
2. Abstract results. The purpose of this section is to present a few abstract results on different types of approximation (temporal, spatial or fully discretized) of a continuous in time dissipative dynamical system on an infinite dimensional phase space, which capture the long-time statistical properties asymptotically. The overriding theme is to be faithful to the original system in the sense that the approximations must be uniformly dissipative and uniformly convergent on the unit time interval (modulo an initial layer).
We first consider temporal discretization since long time approximation seems to be one of the key issues involved. The following abstract result on time discretizations that guarantees the convergence of long-time statistical property can be found in [49].
Theorem 1 (Temporal discretization). Let {S(t), t ≥ 0} be a continuous semigroup on a separable Hilbert space H which generates a continuous dissipative dynamical system (in the sense of possessing a compact global attractor A) on H. Let {S k , 0 < k ≤ k 0 } be a family of continuous maps on H which generates a family of discrete dissipative dynamical system (with global attractor A k ) on H. Suppose that the following two conditions are satisfied.
is pre-compact in H. TH2: [ Uniform convergence on the unit time interval] S k uniformly converges to S on the unit time interval (modulo an initial layer) and uniformly for initial data from the global attractor of S k in the sense that for any t 0 ∈ (0, 1) Moreover, we assume that S(t) is uniformly continuous on K over the unit time interval in the sense that for any T * ∈ [0, 1] Then, the invariant measures of the discrete dynamical system {S k , 0 < k ≤ k 0 } converge to invariant measures of the continuous dynamical system S. More precisely, let µ k ∈ IM k where IM k denotes the set of all invariant measures of S k . There must exist a subsequence, still denoted {µ k }, and µ ∈ IM (an invariant measure of S(t)), such that µ k weakly converges to µ, i.e., Moreover, extremal statistics converge in an upper-semi-continuous fashion in the sense that for any bounded continuous functional Φ on the phase space H, there exist ergodic invariant measures µ k ∈ IM k and an ergodic invariant measure µ ∈ IM, such that sup u0∈H lim sup lim sup Remark. The separable Hilbert space can be replaced by a ball in the same space. Such a set-up may be useful in certain applications [13].
Next, we consider spatial discretization. The following abstract result on spatial discretization that ensures convergence of long-time statistical properties can be found in [48].
Then for any sequence of invariant measures µ N ∈ IM N of the dynamical system S N , there must exists a subsequence, still denoted {µ N }, and an invariant measure where E * N is the lift operator induced by the continuous embedding E N in the sense that for any bounded continuous test functional ϕ on H Practical implementation of the schemes requires both temporal and spatial discretization. Our next goal here is to extend our criterion on semi-discrete in time or space schemes (discrete in time but continuous in space, or discrete in space but continuous in time dynamical systems) with convergent stationary statistical properties to the case of fully discrete schemes (finite dimensional in space and discrete in time dynamical systems, or maps on finite dimensional spaces). We anticipate that the guiding principle in selecting schemes that are able to capture stationary statistical properties, i.e., be faithful to the underlying dynamical system, remain essentially the same within the dynamical system approach. In particular, we can show that the main ingredient remains to preserve the stability (in the sense of uniform dissipativity) and consistency (in the sense of pathwise convergence on the unit time interval), reminiscent of the Lax-Richtmyer equivalence theorem. Just as in the semi-discrete in space case, we need a continuous linear map (embedding operator) E h that maps the finite dimensional space into the underlying infinite dimensional phase space.
Moreover, we assume that S(t) is uniformly continuous on K over the unit time interval in the sense that for any t ∈ [0, 1] Then for any sequence of invariant measures µ h,k ∈ IM h,k , the set of all invariant measures of S h,k , of the fully discrete dynamical system S h,k , there must exists a subsequence, still denoted {µ h,k }, and an invariant measure µ ∈ IM of S(t) such that where E * h is the lift operator induced by the continuous embedding E h in the sense that for any bounded continuous test functional ϕ on H Moreover, extremal statistics converge in an upper-semi-continuous fashion in the sense that for any bounded continuous functional Φ on the phase space H, there exist ergodic invariant measures µ h,k ∈ IM h,k and an ergodic invariant measure µ ∈ IM, such that sup u0∈H lim sup lim sup h,k→0 Proof. Thanks to the definition of the lift of invariant measures of S h,k (17), the continuity of the embedding operator E h , and the Kakutani-Riesz representation theorem [22], we see that {E * h µ h,k } must be a family of Borel probability measures on H. Moreover, since K = 0<k≤k1,0<h≤h1 E h A h,k is pre-compact in H by (H1), and since all invariant measures are supported on the global attractor [11,47] and µ h,k ∈ IM h,k , we see that {E * h µ h,k } is tight in the space of all Borel probability measures on H thanks to Prokhorov's theorem [1,22,11]. Hence, it must have a convergent subsequence, still denoted {E * h µ h,k }, and a Borel probability measure µ on H so that for all bounded and continuous functionals ϕ on H.
Our goal is to show that µ is invariant under S(t), i.e., µ ∈ IM. Now we fix a t ∈ (0, 1] and let n k = t k be the floor of t k (the largest integer dominated by t k ), and let ϕ be any smooth (C 1 ) test functional with compact support. We have where we have used the weak convergence of E * h µ h,k to µ, the boundedness and continuity of ϕ, ϕ • S(t), ϕ • E h , the invariance of µ h,k under S h,k , and the notation Thanks to the uniform continuity of S(t) (with τ = n k k in (16)), the assumption that ϕ ∈ C 1 with compact support, the mean value theorem, and the fact that the support of µ h,k is contained in the global attractor A h,k [48], and the precompactness of K, we have Likewise, thanks to the mean value theorem, the assumption that ϕ ∈ C 1 with compact support, and the uniform convergence of S h,k (with t 0 = t/2 in (15), and k ≤ t/3) over the unit time interval (modulo an initial layer), we have which is exactly the weak invariance (2) for the smooth (C 1 ) test functional with compact support and t ∈ (0, 1]. For a general bounded continuous test functional ϕ, we can first approximate it by a finite dimensional test functional of the form ϕ•P m where P m is the orthogonal projection onto the m-dimensional subspace spanned by the first m elements of a given (fixed) orthonormal basis of H [11]. (This is where the assumption that H is a separable Hilbert space is used.) We can then approximate ϕ • P m by smooth test functionals with compact support using mollifiers and truncation [22] since only the value of ϕ • P m on the compact global attractor is relevant for the statistics. This proves short time weak invariance (2) for any bounded continuous test functional ϕ and t ∈ (0, 1]. Now for a general t > 1, there exists a unique positive integer n and t * ∈ (0, 1] such that t = n + t * . Hence where we have used the semi-group property of S(t), the strong continuity of S(t), the short time weak invariance that we proved above with t = 1 n times, and t = t * one time.
This ends the proof of the convergence of the invariant measures. The first half of the results on the extremal statistics is a consequence of the fact that extremal statistics are saturated by ergodic invariant measures (see for instance [44], or [47] Theorem 5). We sketch the proof for convenience.
Recall that for any fixed bounded continuous test functional Φ and initial data u 0 , it is possible to choose a special Banach limit LIM which agrees with the lim sup on the orbit [22,47] and hence there exists invariant measures µ u0 ∈ IM such that Likewise, since Φ • E h is a bounded continuous test functional on H h , there exists a µ h,k,u0 ∈ IM h,k such that On the other hand, extremal points of the set of invariant measures are ergodic ( [44] or [47] Theorem 3). Hence, there exist ergodic invariant measures µ h,k ∈ IM h,k , µ ∈ IM such that Combining the above two sets of equations together with the tightness of IM k , IM, we arrive at (18,19).
As for the upper semi-convergence of the extremal statistics stated in (20), we have, thanks to the uniform dissipativity and Prokhorov's theorem, there exists a subsequence (still denoted {µ h,k }) and ν ∈ IM such that Since IM is tight in the space of Borel probability measures on H, there exists an ergodic invariant measure ν max ∈ IM such that supμ ∈IM H Φ(u) dμ = H Φ(u) dν max [47]. Therefore, This completes the proof of the theorem.

Remark 1.
The results remain the same if the separable Hilbert space H is replaced by a ball within the same space as long as the global attractors are contained in this ball, and the ball is invariant under the dynamics.
In application, the discrete dynamical systems {S h,k } are usually generated by numerical schemes, with finite dimensional phase space H h of the original infinite dimensional dynamical system (generated by a time-dependent PDE). In another word, S n h,k (u) is the solution to the numerical scheme. For the case of conformal spatial discretisation, the embedding operator can be taken as the natural inclusion operator. However, the case with non-conformal spatial discretization such as finite difference discretisation is more challenging. Appropriate interpolation operators can be used in the finite difference case (see the next section). The uniform dissipativity of the numerical scheme can be established via the existence of a uniform (in mesh size) absorbing ball in discrete form in another separable Hilbert space V which is compactly imbedded in H in the case of strongly dissipative system (see the next section for two examples). The finite time uniform convergence comes with classical numerical analysis for reasonable schemes and the smoothing property of the underlying dissipative system.
A by-product of the convergence analysis of the invariant measures presented here is the convergence of the global attractors of the scheme to that of the underlying system. This is also within expectation since the global attractors carry the support of the invariant measures [48]. The convergence of the global attractors under discretization has been discussed for the two dimensional Navier-Stokes system, reaction-diffusion equation, and for finite-dimensional dynamical systems (see [31,34,16,24,51,25] among others). Therefore our result on the convergence of global attractors may be viewed as a generalization and abstraction of these results. However, we would like to point out that the convergence of the global attractors is established under much weaker assumption: one only needs the uniform boundedness of the union of the global attractors (K), instead of the pre-compactness (plus finite time uniform convergence for data from K). Because of this important distinction, it is possible to have schemes that are able to capture the global attractor asymptotically but not necessarily the stationary statistical properties (invariant measures). There are also interesting works on persistence under approximation of various invariant sets (such as steady state, time periodic orbit, inertial manifold etc) both for PDEs and ODEs under appropriate assumptions (such as spectral gap condition that is usually associated with inertial manifold theory, see [34,17,37,38] and the references therein). We also notice that the convergence of invariant sets and the convergence of stationary statistical properties are two related but very different issues associated with the long-time behavior. It is easy to construct two dynamical systems with exactly the same global attractor or inertial manifolds as sets, but with totally different dynamics or stationary statistical properties.
Next, we show the convergence of the global attractors under weaker assumptions, namely the uniform boundedness of K (the union of the global attractors), and uniform convergence on finite time interval (modulo an arbitrary initial layer).
Theorem 4 (Convergence of Global Attractors). Let {S(t), t ≥ 0} be a continuous semi-group on a Banach space H which generates a dissipative dynamical system (in the sense of possessing a compact global attractor A) on H.
is bounded in H. H4: [ Finite time uniform convergence] S h,k uniformly converges to S on any finite time interval (modulo an initial layer) and uniformly for initial data from the global attractor of the scheme in the sense that there exists t 0 > 0 such that for any t > t 0 > 0 Then the global attractors converge in the sense of Hausdorff semi-distance, i.e.
Proof. Since K is bounded, for any given > 0, there exists a T ε > t 0 > 0 such that dist H (S(t)K, A) < ε 2 , ∀t ≥ T ε because the global attractor A attracts all bounded set, in particular K. Now for an arbitrary element u h,k ∈ A h,k in the global attractor of S h,k , we have that there exists a v h,k ∈ A h,k such that Thanks to the uniform convergence on [T ε , T ε + 1] and the fact that v h,k ∈ A h,k , we have, there exist k ε > 0, h ε > 0 such that This implies that This completes the proof for the convergence of the global attractors.
We would like to reiterate the point that the uniform boundedness assumption H3 is much weaker than the uniform dissipativity assumption H1 for infinite dimensional systems although they are equivalent for finite dimensional systems. This is an important difference and hence it is theoretically possible to have schemes that are able to capture the global attractor asymptotically but not the invariant measures necessarily for infinite dimensional systems. Conditions H4 and H2 are almost the same. H2 is slightly stronger than H4 in some sense since H2 requires the uniform convergence of the scheme on [t 0 , 1], ∀t 0 ∈ (0, 1) while H4 only requires the uniform convergence of the scheme on [t 0 , t] (t ≥ 1) with one t 0 ∈ (0, 1). On the other hand, H4 is slightly stronger than H2 in some other sense since only t = 1 is needed in H2. They are usually valid for the same type of reasonable numerical schemes and hence normally makes no difference. Therefore convergence of the global attractors is usually easier to establish than the convergence of the invariant measures (stationary statistical properties). Related results on convergence of global attractors can be found in [14,34,29] among others.

Remark 2.
We also point out that the semi-discretization in time results derived in [49] are special cases of the fully discretized results presented here. In order to recover the semi-discretization in time case, we simply set H h = H and E h = I (the identity operator) in our fully discretized results above. Likewise, the semidiscretization in space case presented in [48] is a special case of the fully discretized results proved above. In order to recover the semi-discretization in space case, we simply set studied as a model equation for chaotic infinite dimensional dynamical system (see [38] and references therein). Numerical analysis and simulations on the model are also available (see for instance [24,25] for some of the numerical analysis works among many others). The model takes the form together with periodic boundary condition on the interval Ω = [0, 1]. Here R > 0, ν, µ are real parameters. The equation generates a continuous dissipative infinite dimensional dynamical system on H = L 2 per (Ω) [38]. There is a natural orthonormal basis {Ψ l = e 2πilx , l = 0, ±1, ±2, . . . } which are the eigenfunctions of the principal linear operator ∂ 2 ∂x 2 with the associated periodic boundary condition.
For simplicity, we consider backward Euler in time and centered difference in space for discretization. More specifically, let k be the time-step size, h = 1/J be the grid size. Let u n j be our approximation to u(jh, nk) with periodicity in j of period J. Introducing the standard forward (right) difference operator D + , the backward (left) difference operator D − , and the centered difference operator D 2 = D + D − so that (D 2 u n ) j = 1 h 2 (u n j+1 − 2u n j + u n j−1 ), the scheme that we consider then takes the form u n+1 Here the unknown at each time step is u n = (u n 0 , · · · , u n J−1 ) together with periodic extension when the index goes beyond 0 or J − 1. Therefore, the phase space for the fully discrete numerical algorithm is where the subscript per denotes periodicity in j with period J. The space is equipped with the discrete L 2 norm H h has a natural orthonormal basis ψ l = (1, e 2πilh , . . . , e 2πil(J−1)h ), l = −J 2 , . . . , 0, . . . , J 2 that consists of the eigenvectors of D 2 . A discrete H 1 space, denoted H 1 h , can then be defined via the following discrete H 1 norm It is known that the numerical scheme (25) is uniquely solvable and dissipative under appropriate time-step restriction [24].
A natural embedding operator from H h to H is given by the following interpolation operator (Notice that our extension operator E h is denoted Θ −1 in [24].) It is easy to see that E h is an isometry. Furthermore, thanks to Lemma 2.3 [24], the discrete and continuous H 1 norms are equivalent in the sense that 2 This H 1 norm equivalence together with Theorem 3.4 of [24] implies that the scheme generates a dissipative dynamical system with uniform dissipativity since per and hence pre-compact in H. The uniform convergence of the scheme over the unit time interval for initial data from K, i.e., the verification of (15), is a standard albeit tedious (see Lemma 2.5 [24] for a related result in the semi-discrete case). The uniform continuity of the solution semigroup associated with the complex Ginzburg-Landau equation (24)) is straightforward. We leave the details to the interested reader.
Therefore, the conditions in the main theorem can be verified, and hence the results are valid. In particular, long time statistical properties of the numerical scheme (25) converge to those of the equation (24).

3.2.
Two-dimensional Navier-Stokes equations. Here we consider the application of the main abstract result to the well-known two-dimensional Navier-Stokes system in a periodic box Ω in stream-function -vorticity formulation [26,6,37].
where ω is the vorticity, ψ is the stream fucntion, ν is the viscosity, and f is the (time-independent) external body force. We first recall a few notations and facts from a recent work [13] (section 4) for convenience.
Consider a 2D periodic box Ω = (0, 1) × (0, 1). For a given positive integer n, we define N = N x = N y = 2n + 1, we set h x = 1/N x = 1/N y = h y = h. All variables are evaluated at the regular numerical grid (x i , y j ), with x i = ih, y j = jh, Let H h = {all periodic function over the given 2D numerical grid with grid size h}. For f ∈ H h , assume its discrete Fourier expansion is given by The collocation interpolation operator E h is then defined as

4614
XIAOMING WANG E h maps H h into H = L 2 per (Ω). As a result, its collocation Fourier spectral approximations to first and second order partial derivatives (in x direction) are given by The corresponding collocation spectral differentiations in y directions can be defined in the same way. In turn, the discrete Laplacian, gradient and divergence operators can be defined as at the point-wise level. Moreover, given any periodic grid functions f and g (over the 2-D numerical grid), the discrete L 2 norm and inner product are introduced as i,j=0 Meanwhile, such a discrete L 2 inner product can also be viewed in the Fourier space rather than in physical space, with the help of Parseval identity: in which (f N c ) k1,l1 , (ĝ N c ) k1,l1 are the Fourier interpolation coefficients of the grid functions f and g in the expansion as in (28). Furthermore, a detailed calculation shows that the following formulas of summation by parts are also valid at the discrete level: This allows us to define the discrete H 1 norm as It is known that the embedding operator defined in (29) has the following desired continuity property (Lemma 4.1 [13] for instance).
where C is a constant independent of ϕ. A simple interpolation argument implies that the above inequality remains valid for fractional order Sobolev spaces H δ , δ ∈ (0, 1) as well. Next, we introduce our fully discrete algorithm where we treat the nonlinear convection term explicitly for the sake of numerical convenience, and the diffusion term implicitly to preserve the L 2 stability.
where we have utilized a now standard technique in ensuring the skew-symmetry of the nonlinear term [36]. It is observed that the numerical velocity u n+1 = ∇ ⊥ N ψ n+1 is automatically divergence-free: at any time step. Meanwhile, note that the nonlinear term is a spectral approximation to 1 2 u n ·∇ω and 1 2 ∇ · (uω) at time step t n . Furthermore, a careful application of summation by parts formula (35) gives In other words, the nonlinear convection term appearing in the numerical scheme (37), so-called skew symmetric form, makes the nonlinear term orthogonal to the vorticity field in the L 2 space, without considering the temporal discretization. This property is crucial in the stability analysis for the Fourier collocation spectral scheme (37)-(39) as presented in section 4 in [13]. We now recall that sections 4.3, 4.4 and 4.5 in [13] imply that a sufficiently large ball in H δ h , denoted B δ (C 1 ), is invariant under the discrete dynamics determined by the scheme (37) with a mild time-step restriction that is independent of N . Moreover, the scheme generates a dissipative dynamical system on B δ (C 1 ) for sufficiently small time-step size k. Furthermore, the attractors A h,k of the scheme are uniformly bounded in H 1 per independent of h and k as long as these parameters are sufficiently small. This, together with the continuity of the embedding operator E h , implies the uniform dissipativity of the scheme.
The uniform continuity of the solution semigroup to the 2D Navier-Stokes equation in a periodic box with data from a ball in H 1 per is well-known [7]. The uniform convergence over the unit time interval for initial data from a ball in H 1 per is straightforward although tedious. We leave the details to the interested reader.
Combining the above we see that the abstract result with H h (= H δ h ) replaced by the ball B δ (C 1 ) ⊂ H h is applicable to the efficient scheme (37) for the 2D Navier-Stokes equation in a periodic box. Therefore, the long-time statistical properties of the scheme (37) converge to those of the 2D NSE.
The result above can be easily generalized to the barotropic quasi-geostrophic equation with dissipation and external forcing [27]. ∂q ∂t + ∇ ⊥ ψ · ∇q = Dψ + f, q = −∆ψ + F ψ + βy, F ≥ 0, β > 0, Here ψ(x, y; t) represents the stream-function of the barotropic flow, q denotes the potential vorticity, F > 0 is the F-plane constant related to the stratification of the fluid, β > 0 is the beta-plane constant, d 1 is the Ekman damping coefficient, d 2 is the eddy/Newtonian viscosity coefficient, d j , j ≥ 3 are the coefficients of various hyper-viscosity, and f represents external forcing. The convergence of long-time statistical properties for the Fourier spectral discretization of this equation was discussed in [48]. 4. Conclusions and remarks. We have presented several general/abstract results on the convergence of stationary statistical properties of various approximations of infinite dimensional dissipative dynamical systems. Temporal, spatial as well as combined temporal-spatial discretizations are all considered. The two natural conditions that are crucial to the convergence of the stationary statistical properties are uniform dissipativity of the scheme; and some kind of uniform convergence of the scheme on [0, 1] modulo an initial layer. The conditions are natural and reminiscent of those in the Lax equivalence theory. It is easy to understand the necessity of the uniform dissipativity, at least heuristically, since the invariant measures on supported on the global attractors. It is also hard to image any scheme to be able to capture the long-time statistics without being able to approximate the short time behavior well, i.e., violating the short-time convergence result. Hence, these conditions are also heuristically necessary.
Applications to the complex Ginzburg-Landau equation and the two-dimensional Navier-Stokes equation in a periodic box are sketched. One of the immediate goals is to implement the numerical scheme for the 2D NSE for the study of 2D turbulence. We will report our numerical results elsewhere.
There are many questions that merit further investigation. For instance, the convergence presented here is weak convergence without any rate. It would be very desirable to have algorithms with a rate of convergence, at least for a few physically important observables, similar to the case of stochastic differential equations [35,8]. It would also be desirable to consider higher order in time schemes so that we could take relatively large time-step in order to reach statistical equilibrium quickly. Some special cases are considered in [40,50] but a general theory is still missing.
We hope that this work will stimulate further work on numerical schemes that are able to capture stationary statistical properties of infinite dimensional dissipative systems.