Space-time FLAVORS: finite difference, multisymlectic, and pseudospectral integrators for multiscale PDEs
SSpace-time FLAVORS: finite difference, multisymlectic, andpseudospectral integrators for multiscale PDEs
Molei Tao, Houman Owhadi, Jerrold E. MarsdenOctober 23, 2018
Abstract
We present a new class of integrators for stiff PDEs. These integrators aregeneralizations of FLow AVeraging integratORS (FLAVORS) for stiff ODEs andSDEs introduced in [32] with the following properties: (i)
Multiscale : they are basedon flow averaging and have a computational cost determined by mesoscopic stepsin space and time instead of microscopic steps in space and time; (ii)
Versatile :the method is based on averaging the flows of the given PDEs (which may havehidden slow and fast processes). This bypasses the need for identifying explicitly(or numerically) the slow variables or reduced effective PDEs; (iii)
Nonintrusive : Apre-existing numerical scheme resolving the microscopic time scale can be used asa black box and easily turned into one of the integrators in this paper by turningthe large coefficients on over a microscopic timescale and off during a mesoscopictimescale; (iv)
Convergent over two scales : strongly over slow processes and in thesense of measures over fast ones; (v)
Structure-preserving : for stiff HamiltonianPDEs (possibly on manifolds), they can be made to be multi-symplectic, symmetry-preserving (symmetries are group actions that leave the system invariant) in allvariables and variational.
Multi-scale PDEs can be divided into two (possibly over-lapping) categories: PDEswith highly oscillating or rough coefficients and PDEs with large (or stiff) coefficients.Classical numerical methods are usually: (i) stable but arbitrarily inaccurate for theformer category (consider, for instance, a finite element method for the elliptic oper-ator − div( a ∇ ) with a rapidly changing coefficient a ∈ L ∞ ), or (ii) unstable for thelatter category. Accurate numerical methods for the former category, called numericalhomogenization methods, are, in absence of local ergodicity or scale separation, basedon the compactness of the solution space (we refer, for instance, to [26, 2, 27]). Numer-ical methods for the latter category are, in essence, based on the existence of slow andfast variables (or components) [14]. When fast variables converge toward Dirac (singlepoint) distributions, asymptotic-preserving schemes [15] allow for simulations with largetime steps. We also refer to [18, 25] for multi-scale transport equations and hyperbolicsystems of conservation laws with stiff diffusive relaxation. Well-identified slow variables1 a r X i v : . [ m a t h . NA ] A p r an be simulated with large time-steps using the two-scale structure of the original stiffPDEs (we refer to [1] and [12] for existing examples; slow variables satisfy a non-stiffPDE that can be identified in analogy to equations (A.9) and (A.13) of [32]; we alsorefer to [14] for a definition of slow variables).In this paper, we consider the second category of PDEs and propose a generaliza-tion of FLow AVeraging integratORS (FLAVORS) (introduced in [32] for stiff ODEsand SDEs) to stiff PDEs. Multi-scale integrators for stiff PDEs are obtained withoutthe identification of slow variables by turning on and off stiff coefficients in single-step(legacy) integrators (used as black boxes) and alternating microscopic and mesoscopictime steps (Subsection 2.2). We illustrate the generality of the proposed strategy by ap-plying it to finite difference methods in Section 2, multi-symplectic integrators in Section3, and pseudospectral methods in Section 4 (although we have not done so in this paper,the proposed strategy can also be applied to finite element methods or finite volumemethods). The convergence of the proposed strategy, after semi-discretization in space,is analyzed in Subsection 5.1, where a non-asymptotic error bound indicates the two-scale convergence ([32], i.e., strong with respect to hidden slow variables and weak withrespect to hidden fast variables) of PDE-FLAVORS. As illustrated by numerical (Figure2) and theoretical results (Section 5), an explicit tuning (( h/(cid:15) ) (cid:28) H (cid:28) h/(cid:15) ) betweenmicroscopic h and mesoscopic ( H ) time-steps and the stiff parameter 1 /(cid:15) is necessaryand sufficient for convergence. We also show in Section 6 that applying the FLAVORstrategy to characteristics leads to accurate approximations of solutions of stiff PDEs.These results, along with those of [32], diverge from the concept that, in situationswhere the slow variables are not linear functions of the original variables, multiscalealgorithms “do not work” “if the slow variables are not explicitly identified and madeuse of ” (page 2 of [13]). Consider a multiscale PDE: F (1 , (cid:15) − , x, t, u ( x, t ) , u x ( x, t ) , u t ( x, t ) , u xx ( x, t ) , u xt ( x, t ) , u tt ( x, t ) , . . . ) = 0 (1)where F is a given function (possibly nonlinear), (cid:15) is a small positive real parameter and x and t are spatial and temporal coordinates.To obtain a numerical solution of (1), the simplest single-scale finite difference ap-proach employs a uniform rectangular mesh with time step length h and space step length k , and approximates the solution u at its values at discrete grid points. Differential oper-ators will be approximated by finite differences; for instance, according to forward spaceforward time rules: u x ( ik, jh ) ≈ ( u i +1 ,j − u i,j ) /k and u t ( ik, jh ) ≈ ( u i,j +1 − u i,j ) /h , where u ij is the numerical solution at discrete grid point with space index i and time index j . After this discretization, the original PDE is approximated by a finite dimensionalalgebraic system, which can be solved to yield the numerical solution.2f course, a necessary condition for obtaining stability and accuracy in the numericalsolution is that h and k have to be small enough. A quantitative statement on howsmall they need to be will depend on the specific PDE and discretization. For 1Dlinear advection equations u x − au t = 0 and forward time forward space discretizations,the h < k/a CFL condition [11] has to be met to ensure stability, which is also aneccessary condition for accuracy [19]. Intuitively, the CFL condition guarantees thatinformation does not propagate faster than what the numerical integrator can handle.The Von Neumann stability analysis [9] helps determine analogous CFL conditions forlinear equations with arbitrary discretizations. The stability of numerical schemes forgeneral nonlinear equations remains a topic of study. We refer to [31] for additionaldiscussions on single-scale finite difference schemes. In general, the presence of a stiffcoefficient (cid:15) − in equation (1) requires h and k to scale with (cid:15) in order to guarantee thestability of numerical integration schemes. This makes the numerical approximation ofthe solution of (1) computationally untractable when (cid:15) is close to 0. FLAVORs are multiscale in the sense that they accelerate computation by adoptingboth larger time and space steps. A finite difference scheme can be FLAVORized byemploying two rules:Figure 1:
Mesh used by FLAVORS. A uniform mesoscopic space step is used and two alternating smalland mesoscopic time steps are used. Stiffness is turned on in red regions and turned off otherwise. K that does not scale with (cid:15) , and an alternating temporal grid corresponds to two time steps, microscopic h (scalingwith (cid:15) ) and mesoscopic H − h ( H independent from (cid:15) ). It is worth mentioning that whenusing this non-uniform mesh, grid sizes have to be taken into consideration when deriva-tives are approximated by finite differences. 1st-order derivatives are straightforward toobtain, and we refer to Section 3 for approximations of higher order derivatives.Second, the stiff parameter (cid:15) − should be temporarily set to be 0 (i.e., turned off)when the current time step is the mesoscopic H − h ; if the small time step h is usedinstead, the large value of (cid:15) − needs to be restored, or in other words, stiffness shouldbe turned on again.The rule of thumb is that k and h should be chosen such that the integration of (1)with these step sizes and stiffness turned on is stable and accurate. On the other hand,there is another pair of step size values such that the same integration with stiffnessturned off is stable and accurate, and K and H should be chosen to be an order ofmagnitude smaller than these values. FLAVORS does not require a microscopic k , butonly a mesoscopic space-step K , a microscopic time-step h , and a mesoscopic time-step H . The intuition is as follows: adopt the point of view of semi-discrete approach forPDE integration, in which space is discretized first and the PDE is approximated by asystem of ODEs. The integration (in the time) of the resulting finite dimensional ODEsystem can be accelerated by applying the FLAVOR strategy to any legacy scheme(used as a black box). Turning on and off stiff coefficients in the legacy scheme andalternating microscopic time steps (stiffness on) with mesoscopic time steps (stiffnesson) preserves the symmetries of that scheme and at the same time induces an averagingof the dynamic of (possibly hidden) slow variables with respect to the fast ones. Withthis strategy, the FLAVORized scheme advances in mesoscopic time steps without losingstability. The (possibly hidden) slow dynamic is captured in a strong sense, while the fastone is captured only in the (weak) sense of measures. A rigorous proof of convergenceof the proposed method relies on the assumption of existence of (possibly hidden) slowvariables and of local ergodicity of (possibly hidden) fast variables (we refer to Section 5).It is important to observe that the proposed method does not require the identificationof slow variables. Consider a specific stiff PDE: u t + f ( u ) x = (cid:15) − u (1 − u ) (2)in which f ( u ) = sin u and 0 < (cid:15) (cid:28)
1. Use the boundary condition of u ( x = 0 , t ) = u ( x = L, t ) and the initial condition of u ( x, t = 0) = sin( πx ). This system contains twoscales: the fast process corresponds to u quickly converging towards 1 or −
1, and theslow process corresponds to the front (with steep gradients) separating u > u < O (1) velocity. 4e will FLAVORize the following Lax-Friedrichs finite difference scheme: (cid:40) u i +1 ,j +1 = ¯ u i +1 ,j − h (cid:16) f u (¯ u i +1 ,j ) u i +2 ,j − u i,j k + (cid:15) − ¯ u i +1 ,j (1 − ¯ u i +1 ,j ) (cid:17) ¯ u i +1 ,j (cid:44) u i +2 ,j + u i,j (3)where u i,j = u i + L/k,j and u i, = sin ( π ( i − k ). If the domain of integration is restrictedto [0 , L ] × [0 , T ], then i = 1 , , . . . , (cid:98) L/k (cid:99) +1, and j = 1 , , . . . , (cid:98) T /h (cid:99) +1. We use h = 0 . (cid:15) and k = 0 . (cid:15) for our purposes, both of which we found numerically at the order of thestability limit. In our experiment, we chose (cid:15) = 2 · − , and therefore h = 0 . k = 0 . ˜ u i +1 ,j = ¯ u i +1 ,j − h (cid:16) f u (¯ u i +1 ,j ) u i +2 ,j − u i,j K + (cid:15) − ¯ u i +1 ,j (1 − ¯ u i +1 ,j ) (cid:17) ¯ u i +1 ,j (cid:44) ( u i +2 ,j + u i,j ) / u i +1 ,j +1 = ˜ u i +2 ,j +˜ u i,j − ( H − h ) (cid:16) f u ( ˜ u i +2 ,j +˜ u i,j ) ˜ u i +2 ,j − ˜ u i,j K (cid:17) (4)where u i,j = u i + L/K,j and u i, = sin ( π ( i − K ). If the domain of integration is re-stricted to [0 , L ] × [0 , T ], then i = 1 , , . . . , (cid:98) L/K (cid:99) + 1, and j = 1 , , . . . , (cid:98) T /H (cid:99) + 1. Weuse the same h as before, and choose H = 0 .
005 and K = 0 .
01, which ensures that thestability of the integration remains independent of (cid:15) .Errors of FLAVOR based on Lax-Friedrichs with different H and h values are com-puted by comparing the results to a benchmark Lax-Friedrichs integration with fine steps h = 0 . (cid:15) and k = 0 . (cid:15) . More precisely, we calculated the distance between two vectorsrespectively corresponding to FLAVOR and Lax-Friedrichs integrations, which containordered u ( x, t ) values on the intersection of FLAVOR and Lax-Friedrichs meshes (whichis in fact the FLAVOR mesh as long as H is a multiple of 0 . (cid:15) ). 1-norm is used andnormalized by the number of discrete points to mimic the L norm for the continuoussolution. Experimental settings are (cid:15) = 2 · − , L = 2 and T = 2. As we can seein Figure 2, FLAVOR is indeed uniformly convergent in the sense that the error scaleswith H , as long as h takes an appropriate value. This is not surprising, because we havealready proven in the ODE case that the error is bounded by a function of H (uniformlyin (cid:15) ) as long as (cid:0) h(cid:15) (cid:1) (cid:28) H (cid:28) h/(cid:15) , and this error can be made arbitrarily small as H ↓ H can still be much larger than (cid:15) as (cid:15) ↓ H = 0 .
005 and K = 0 .
01) in comparison to thebenchmark ( h = 0 . k = 0 . HK hk = 312 . (cid:15) →
0, and this statement will be true forall FLAVOR examples shown in this paper.5 ε L distance between u x,t by FLAVOR and benchmarkH e rr o r Figure 2:
Errors of FLAVOR based on Lax-Friedrichs as a function of H and h . H samples multiplesof 0 . (cid:15) , starting from 2x to 50x with 1x increment, and h ranges from 0 . (cid:15) to 3 (cid:15) with 0 . (cid:15) increment.Errors with magnitude bigger than 1 are not plotted, for they indicate unstable integrations. We refer to [7, 21, 22] for a discussion on the geometry of Hamiltonian PDEs (e.g.,multi-symplectic structure). We will now recall the Euclidean coordinate form of aHamiltonian PDE: M z t + K z x = ∇ z H ( z ) (5)where z ( x, t ) is a n-dimensional vector, M and K are arbitrary skew-symmetric matriceson R n , and H : R n → R is an arbitrary smooth function. The solution preserves themulti-symplectic structure in the following sense: ∂ t ι ( U, V ) + ∂ x κ ( U, V ) = 0 (6)where ι and κ are differential 2-forms defined by ι ( x, y ) = (cid:104)M x, y (cid:105) and κ ( x, y ) = (cid:104)K x, y (cid:105) (7)6 u ( x ,t ) u ( x ,t ) Figure 3:
Numerical solutions to (2) by Lax-Friedrichs (left, Eq. 3) and its FLAVORization (right,Eq. 4). and U and V are two arbitrary solutions to the variational equation (the solution isidentified with dz : R (cid:55)→ R n ): M dz t + K dz x = D zz H ( z ) dz, dz ( x, t ) ∈ R n (8)Preservation of multi-symplecticity can be partially and intuitively interpreted as a con-servation of infinitesimal volume in the jet bundle, which generalizes the conservation ofphase space volume in Hamiltonian ODE settings to field theories.A broad spectrum of PDEs fall in the class of Hamiltonian PDEs, including general-ized KdV, nonlinear Schr¨odinger models, nonlinear wave equations, atmospheric flows,fluid-structure interactions, etc. [4, 3, 6, 7]. We also refer to [8] and references thereinfor surveys on numerical recipes, and to [20] for an application to numerical nonlinearelastodynamics.Hamiltonian PDEs (5) can be viewed as Euler-Lagrange equations for field theories,which are obtained by applying Hamilton’s principle (i.e., a variational principle of δ S /δz = 0) to the following action: S ( z ( · , · )) = (cid:90) (cid:90) L ( z, z t , z x ) dt dx (9)where the Lagrangian density is given by L ( z, z t , z x ) = 12 (cid:104)M z t , z (cid:105) + 12 (cid:104)K z x , z (cid:105) − H ( z ) (10)This variational view of Hamiltonian PDEs will intrinsically guarantee the preserva-tion of multi-symplecticity, and there will be a field generalization of Noether’s theorem,which ensures conservation of momentum maps corresponding to symmetries.Numerically, instead of discretizing the equations (5), we prefer the approach of varia-tional integrators because they are intrinsically multi-symplectic and therefore structure-preserving [21, 22, 23, 20]. These integrators are obtained as follows: first discretize the7ction (9) using quadratures, then apply variational principle to the discrete action(which depends on finitely many arguments), and finally, solve the algebraic systemobtained from the variational principle, i.e., the discrete Euler-Lagrange equations.For an illustration, consider a nonlinear wave equation: u tt − u xx = V (cid:48) ( u ) (11)with periodic boundary condition u ( x + L, t ) = u ( x, t ) and compatible initial conditions u ( x, t = 0) = f ( x ) and u t ( x, t = 0) = g ( x ). Suppose we are interested in the solution ina domain [0 , L ] × [0 , T ].Rewrite the high order PDE as a system of first order PDEs (notice these covariantequations can be obtained through an intrinsic procedure, which works on manifolds aswell [5]): v t − w x = V (cid:48) ( u ) (12) u t = v (13) u x = w (14)The corresponding Lagrangian density is: L = 12 u t − u x + V ( u ) (15)Using a forward time forward space approximation, we obtain the following discreteLagrangian: L di,j (cid:44) h ij k ij (cid:34) (cid:18) u i,j +1 − u i,j h ij (cid:19) − (cid:18) u i +1 ,j − u i,j k ij (cid:19) + V ( u i,j ) (cid:35) (16) ≈ (cid:90) t j +1 = t j + h ij t j dt (cid:90) x i +1 = x i + k ij x i dx (cid:20) u t − u x + V ( u ) (cid:21) (17)where space step k ij and time step h ij define a rectangular grid of size k ij × h ij . Thesimplest single-scale choice would be k ij = k and h ij = h for some k and h .As a consequence, the continuous action S is approximated by a discrete action: S d = N (cid:88) α =1 M (cid:88) β =1 L dα,β ≈ S = (cid:90) (cid:90) L dt dx (18)and Hamilton’s principle of least action δ S d = 0 gives ∂∂u i,j N (cid:88) α =1 M (cid:88) β =1 L dα,β = 0 (19)for 1 ≤ i ≤ N and 1 ≤ j ≤ M , where N and M are such that (cid:80) Nα =1 k αβ = L for any β and (cid:80) Mβ =1 h αβ = T for any α . 8aking derivative with respect to u i,j , we obtain the following discrete Euler-Lagrangeequations: k ij u i,j − u i,j +1 h ij − u ij u i,j − u i +1 ,j k ij + h ij k ij V (cid:48) ( u i,j )+ k i,j − u i,j − u i,j − h i,j − − h i − ,j u i,j − u i − ,j k i − ,j = 0(20)The system of above equations is explicitly solvable when equipped with boundaryconditions and initial conditions; for instance, below is a consistent discretization of thecontinuous version: u i,j = u i + N,j , ∀ i, ju i, = f (cid:16)(cid:80) iα =1 k α (cid:17) , ∀ iu i, = u i, + h i g (cid:16)(cid:80) iα =1 k α (cid:17) , ∀ i (21)This numerical receipt is convergent. In fact, multi-symplectic integrators obtainedfrom variational principles can be viewed as special members of finite difference methods,whose error analysis is classical.We wish to point out that the above procedure works for any Hamiltonian PDEsof form (5). Also, notice that high-order derivatives are dealt with in an intrinsic wayregardless of whether mesh is uniform. Now consider a multiscale Hamiltonian PDE M (1 , (cid:15) − ) z t + K (1 , (cid:15) − ) z x = ∇ z H (1 , (cid:15) − , z ) (22)Any single-scale multi-symplectic integrator can be FLAVORized (to achieve com-putational acceleration) by using the following strategy: (i) Use the two-scale meshillustrated in Figure 1, and (ii) Turn off large coefficients when taking mesoscopic time-steps. Unlike FLAVORizing a general finite difference scheme, we FLAVORize the action S d instead of the PDE. Specifically, choose k ij = K, ∀ i, jh ij = h, ∀ i and odd jh ij = H − h, ∀ i and even j (23)and let (cid:15) − = 0 in L di,j for even j ’s and all i ’s, while the large value of (cid:15) − is kept in L di,j for odd j ’s and all i ’s. h and H correspond to a small and a mesoscopic time-step, and K corresponds to a mesoscopic space-step; the same rule of thumb for choosing them inSection 2 applies.After applying the discrete Hamilton’s principle, the resulting discrete Euler Lagrange-equations corresponding to a multi-symplectic integrator will still be (20), except thatstiffness is turned off in half of the grids. Multisymplecticity is automatically gained,because the updating equations originate from a discrete variational principle [21].9 .3 Example: multiscale Sine-Gordon wave equation Consider a specific nonlinear wave equation (11) in which V ( u ) = − cos( ωu ) − cos( u ).If ω = 0, this corresponds to Sine-Gordon equation, which has been studied extensivelydue to its soliton solutions and its relationships with quantum physics (for instance, asa nonlinear version of Klein-Gordon equation). We are interested in the case in which ω (identified with (cid:15) − ) is big, so that a separation of timescale exhibits.Arbitrarily choose L = 2 and use periodic boundary condition u ( x + L, t ) = u ( x, t ),and let initial condition be u ( x,
0) = sin(2 πx/L ) and u t ( x,
0) = 0. Denote total simula-tion time by T . Use the FLAVOR mesh (23). In order to obtain a stable and accuratenumerical solution, k and h have to be o (1 /ω ), and K and H need to be o (1). ω =20x 0 2 4 6 8 1000.511.52−1−0.500.51 tFLAVOR: ω =20x Figure 4:
Numerical solutions to multiscale Sine-Gordon equation by single-scale 1st-order multi-symplectic integrator (left) and its FLAVORization (right). For clarity, the surface plots (but notsimulations) use the same mesh size.
A comparison between the benchmark of the single-scale forward time forward spacemulti-symplectic integrator (Eq. 20 with h ij = h and k ij = k ) and its FLAVORization(Eq. 20 with mesh (23) and V (cid:48) ( u ) = ω sin( ωu ) + sin( u ) for odd j and V (cid:48) ( u ) = sin( u )for even j ) is presented in Figure 4. ω = 20, k = L/ /ω and h = k/
2, and K = L/ H = K/
2. It is intuitive to say that the slow process of wave propagation is well-approximated by FLAVOR, although the fast process of local fluctuation is not capturedin the strong sense. Error quantification is not done, because what the slow and fastprocesses are is not rigorously known here.
HK/ hk = 50-fold acceleration is obtainedby FLAVOR.Readers familiar with the splitting theory of ODEs [24] might question whetherFLAVORS are equivalent to an averaged stiffness of ˜ ω = ω hH (which corresponds ˜ ω = 2in the numerical experiment described above). The answer is no, because the equivalencygiven by the splitting theory is only local. In fact, the same single-scale forward timeforward space multi-symplectic integration of the case ω = 2 is shown in Figure 5, which10 ω =2x Figure 5:
Numerical solutions to multiscale Sine-Gordon equation with the ‘equivalent’ stiffness bysingle-scale 1st-order multi-symplectic integrator. For clarity, the surface plot (but not the simulation)uses the same mesh size (as in Figure 4). is clearly distinct from the FLAVOR result in Figure 4. Moreover, because of the e CωT error term, changing stiffness alone will not result in a converging method (and resultin a O (1) error on slow variables). Consider a PDE u t ( x, t ) = L u ( x, t ) (24)with periodic boundary condition u ( x, t ) = u ( x + L, t ) and initial condition u ( x,
0) = f ( x ), where L is a differential operator involving only spatial derivatives.The Fourier collocation method approximates the solutions by the truncated Fourierseries: u N ( x, t ) = (cid:88) | n |≤ N/ a n ( t ) e in πx/L (25)and solves for a n ( t )’s by requiring the PDE to hold at collocation points y j : ∂ t u N ( y j , t ) − L u N ( y j , t ) = 0 (26)This yields N ODEs, which can be integrated by any favorite ODE solver. Of course,specific choices of collocations points will affect the numerical approximation. Often-times, the simplest choice of y j = Lj/N, j = 0 , . . . , N − .2 FLAVORization of pseudospectral methods When the PDE is stiff (for instance, when L contains a large parameter (cid:15) − ), FLAVORScan be employed to integrate the stiff ODEs (which will still contain (cid:15) − ) resulting froma pseudospectral discretization.Similarly, for the FLAVORization of a pseudospectral method, it is sufficient tochoose N (cid:29) L instead of N (cid:29) (cid:15) − L , i.e., the space-step can be coarse ( K = o (1)).For time stepping, alternatively switching between h = o ( (cid:15) ) and H − h for a mesoscopic H = o (1) is again needed, and stiffness has to be turned off over the mesoscopic step of H − h . In a sense, we are still using the same FLAVOR ‘mesh’ (Figure 1), except thathere we do not discretize space, but instead truncate Fourier series to resolve the samespatial grid size. Consider the following system of PDEs u t + u x − q = 0 q t + q x − p = 0 p t + p x + ω q = 0 (27)with periodic boundary conditions u ( x, t ) = u ( x + L, t ), q ( x, t ) = q ( x + L, t ), and p ( x, t ) = p ( x + L, t ), and initial conditions u ( x,
0) = f u ( x ), q ( x,
0) = f q ( x ), and p ( x,
0) = f p ( x ).The integration domain is restricted to [0 , T ] × [0 , L ]. The stiffness (cid:15) − is identified with ω . We choose the initial condition of f u ( x ) = f q ( x ) = cos(2 πx/L ) and f p ( x ) = 0.In this system, q and p correspond to a fast process, which is a field theory versionof a harmonic oscillator with high frequency ω . u is a slow process, into which energy ispumped by the fast process in a non-trivial way.We have chosen to FLAVORize (27) because it does not fall into the (simpler) cat-egory of systems with fast processes converging towards Dirac (single point support)invariant distributions [15].We use the classical 4th order Runga-Kutta scheme (see, for instance, [16]) for the(single-step) time integration of the pseudospectrally discretized system of ODEs (26).Write φ ω h : ˜ a u,q,pn ( t ) (cid:55)→ ˜ a u,p,qn ( t + h ) its numerical flow over a microscopic time step h (consisting of four sub-steps), where ˜ a u,q,pn ( t ) are numerical approximations to theFourier coefficients in (25), for the unknowns u , q and p at an arbitrary time t . Then,the corresponding FLAVOR update over a mesoscopic time step H will be φ H − h ◦ φ ω h ,which consists of eight sub-steps.We present in Figure 6 and Figure 7 a comparison between the benchmark of single-scale pseudospectral simulation and its FLAVORization. It can be seen that the slowprocess of u is captured in strong (point-wise) sense, whereas the fast process of q is onlyapproximated in a weak sense (i.e. as a measure, in the case, wave shape and amplitudeare correct, but not the period). We choose L = 2, T = 10 and ω = 1000. The single-step integration uses N = 20 and h = 0 . /ω (notice that this is already beyond the12 u u Figure 6:
Single-scale (left) and multiscale pseudospectral (right) integrations of slow u in system (27).Plotting mesh for the single-scale simulation is coarser than its computation mesh. stability/accuracy region of a single-scale finite difference, since the space step does notdepend on 1 /ω ; the spectral method is more stable/accurate for a large space-step), andFLAVOR uses N = 20, h = 1 /ω and H = 0 . H/ h = 50-fold acceleration is achievedby FLAVOR. All FLow AVeraging integratORS described in previous sections are illustrations of thefollowing (semi-discrete) strategy: first, space is discretized or interpolated; next, spatialdifferential operators are approximated by algebraic functions of finitely many spatialvariables; finally, the resulting system of ODEs is numerically integrated by a corre-sponding ODE-FLAVOR [32]. In this section, we will use the semi-discrete ODE systemas an intermediate link to demonstrate that these PDE-FLAVORS are convergent tothe exact PDE solution under reasonable assumptions (in a strong sense with respectto (possibly hidden) slow variables and in the sense of measures with respect to fastvariables).More precisely, consider a spatial mesh (vector) M S = [ x , x , . . . ], a temporal mesh(vector) M T = [ t , t , . . . ], and a domain mesh (matrix) M = M S × M T . Examples ofthese meshes include the FLAVOR mesh M S = [ K, K, . . . , N K ] and M T = [ h, H, H + h, H, . . . , ( M − H, ( M − H + h, M H ], and a usual single-scale (step) integrationmesh M S = [ k, k, . . . , L ] and M T = [ h, h, . . . , T ] (recall the domain size is L = N K by T = M H ). We will use the FLAVOR mesh throughout this section. We will comparethe solution of the PDE (28) with the solution obtained with the FLAVOR strategy atthese discrete points. 13 q q Figure 7:
Single-scale (left) and multiscale pseudospectral (right) integrations of fast q in system (27).Plotting mesh for the single-scale simulation is coarser than its computation mesh. The same color doesnot indicate the same value in these two plots. For simplicity, assume the PDE of interest is 1st-order in time derivative: u t ( x, t ) = F (1 , (cid:15) − , x, t, u ( x, t ) , u x ( x, t ) , . . . ) (28)Observe that a PDE (1) with higher-order time derivatives can be written as a systemof 1st-order (in time derivatives) PDEs.Now consider a consistent discretization of PDE (28) with space step K and timestep h (we refer to Page 20 of [31] for a definition of the notion of consistency, whichintuitively means vanishing local truncation error). Letting h ↓ ˙ u ( t ) = f ( u , u , . . . , u N , (cid:15) − , t )˙ u ( t ) = f ( u , u , . . . , u N , (cid:15) − , t ) · · · ˙ u N ( t ) = f N ( u , u , . . . , u N , (cid:15) − , t ) (29)Assuming existence and uniqueness of an exact C strong solution u to the PDE (28),and writing u ( M Si , t ) its values at the spatial discretization points, we define for each i the following remainder: R i ( (cid:15) − , t ) (cid:44) ∂u∂t ( M Si , t ) − f i ( u ( M S , t ) , u ( M S , t ) , . . . , u ( M SN , t ) , (cid:15) − , t ) (30)which is a real function of t indexed by (cid:15) − .Then, u i ( t ) approximates the exact solution u ( M Si , t ) evaluated at grid points in thesense that these remainders vanish as (cid:15) − K ↓ K := M Si − M Si − ):14 emma 5.1. Assume that F in (28) satisfies | F (1 , (cid:15) − , x, t, u ( x, t ) , u x ( x, t ) , . . . ) | ≤ (1 + (cid:15) − ) | F (1 , , x, t, u ( x, t ) , u x ( x, t ) , . . . ) | (31) Assume that the f i in (29) satisfy similar inequalities. Then, there exists a constant C i independent from (cid:15) , h , H or K , such that for bounded t and u |R i ( (cid:15) − , t ) | ≤ (1 + (cid:15) − ) C i K (32) Remark 5.1. (31) is true, for instance, in cases where F (1 , (cid:15) − , x, t, u ( x, t ) , . . . ) = F ( x, t, u ( x, t ) , . . . ) + (cid:15) − F ( x, t, u ( x, t ) , . . . ) . (33) Proof.
The linear scaling with K in (32) immediately follows from the definition ofconsistency, and the parameter 1 + (cid:15) − in (32) has its origin (31). Remark 5.2.
The consistency of finite difference methods can be easily shown usingTaylor expansions. For instance, applying a Taylor expansion to the solution of u t − (cid:15) − u x = a ( u ) leads to u ( iK, ( j + 1) h ) = u ( iK, jh ) + h (cid:16) (cid:15) − (cid:0) u (( i + 1) K, jh ) − u ( iK, jh ) K + O ( K ) (cid:1) + a ( u ( iK, jh )) (cid:17) + O ( h ) (34)which implies ∂∂t u ( iK, t ) = (cid:15) − u (( i + 1) K, t ) − u ( iK, t ) K + a ( u ( iK, t )) + (cid:15) − O ( K ) (35)and naturally establishes the correspondence of f i ( u , . . . , u N , (cid:15) − , t ) = (cid:15) − u i +1 ( t ) − u i ( t ) K + a ( u i ( t )) and R i = (cid:15) − O ( K ) for a 1st-order finite difference scheme. Notice that theremainders are still stiff, but we will see later that this is not a problem, since theycan be handled by ODE-FLAVORs. The consistency of pseudospectral method can beshown similarly using Fourier analysis.With R i defined in (30), consider the following system of ODEs: ˙ u ( t ) = f ( u , u , . . . , u N , (cid:15) − , t ) + R ( (cid:15) − , t ) · · · ˙ u N ( t ) = f N ( u , u , . . . , u N , (cid:15) − , t ) + R N ( (cid:15) − , t ) (36)with initial condition u i (0) = u ( M Si , u i ( t )) ≤ i ≤ N is the exactPDE solution sampled at spatial grid points, i.e., u i ( t ) = u ( M Si , t ).We will now establish the accuracy of PDE-FLAVOR by showing that an ODE-FLAVOR integration of (36) leads to an accurate approximation of ( u i ( t )) ≤ i ≤ N . Since15pace (with fixed width L ) is discretized by N grid points, we use the following (nor-malized by N ) norm in our following discussion (suppose v i ( t ) = v ( M Si , t ) for a function v ): (cid:107) [ v ( t ) , v ( t ) , . . . , v N ( t )] (cid:107) (cid:44) N (cid:107) [ v ( t ) , v ( t ) , . . . , v N ( t )] (cid:107) (37)Observe that if v ( · , t ) is Riemann integrable, thenlim K ↓ (cid:13)(cid:13) [ v ( M S , t ) , v ( M S , t ) , . . . , v ( M SN , t )] (cid:13)(cid:13) → L (cid:107) v ( · , t ) (cid:107) L (recall L = N K is fixed) , (38)and hence the norm (37) does not blow up or vanish as N → ∞ . We will now prove the accuracy of PDE-FLAVORs under the assumption of existenceof (possibly hidden) slow and locally ergodic fast variables. The convergence of PDE-FLAVORs will be expressed using the notion of two-scale flow convergence introducedin [32] (corresponding to a strong convergence with respect to slow variables and weakone with respect to fast ones).
Condition 5.1.
Assume that the ODE system (36) satisfies the following conditions:1. (Existence of hidden slow and fast variables): There exists a (possibly time-dependent)diffeomorphism η t : [ u ( t ) , . . . , u N ( t )] (cid:55)→ [ x ( t ) , y ( t )] from R N onto R N − p × R p withuniformly bounded C , C derivatives with respect to u i ’s and t , and such that forall (cid:15) > , ( x ( t ) , y ( t )) satisfies (cid:40) ˙ x ( t ) = f ( x ( t ) , y ( t ) , t )˙ y ( t ) = (cid:15) − g ( x ( t ) , y ( t ) , t ) , (39) where f and g have bounded C derivatives with respect to x , y and t .2. (Local ergodicity of vast variables): There exists a family of probability measures µ t ( x, dy ) on R p indexed by x ∈ R N − p and t ∈ R , and a family of positive functions T (cid:55)→ E t ( T ) satisfying lim T →∞ E t ( T ) = 0 for all bounded t , such that for all x , y , t , T bounded and φ uniformly bounded and Lipschitz, the solution to ˙ Y t = g ( x , Y t , t ) Y = y (40) satisfies (cid:12)(cid:12)(cid:12) T (cid:90) T φ ( Y s ) ds − (cid:90) R p φ ( y ) µ t ( x , dy ) (cid:12)(cid:12)(cid:12) ≤ χ t (cid:0) (cid:107) ( x , y ) (cid:107) (cid:1) E t ( T )( (cid:107) φ (cid:107) L ∞ + (cid:107)∇ φ (cid:107) L ∞ )(41) where r (cid:55)→ χ t ( r ) is bounded on compact sets, and µ t has bounded derivative withrespect to t in total variation norm. αt,t + τ the numerical flow of a given (legacy) ODE integrator for (29):Φ αt,t + τ : [˜ u ( t ) , . . . , ˜ u N ( t )] (cid:55)→ [˜ u ( t + τ ) , . . . , ˜ u N ( t + τ )] , (42)where ˜ u i ( s ) approximates u i ( s ) for all s , τ is the integration time step, and α is acontrollable parameter that replaces the stiff parameter (cid:15) − in (29) and takes values of (cid:15) − (stiffness ‘on’) or 0 (stiffness ‘off’). Definition 5.1 (ODE-FLAVORS) . The FLow AVeraging integratOR associated withΦ is defined as the algorithm simulating the process:[¯ u ( t ) , . . . , ¯ u N ( t )]= (cid:0) Φ k − H + h,kH ◦ Φ (cid:15) ( k − H, ( k − H + h (cid:1) ◦ · · · ◦ (cid:0) Φ H + h, H ◦ Φ (cid:15) H,H + h (cid:1) ◦ (cid:0) Φ h,H ◦ Φ (cid:15) ,h (cid:1) ([ u (0) , . . . , u N (0)])(43)where (the number of steps) k is a piece-wise constant function of t satisfying kH ≤ t < ( k + 1) H , h is a microscopic time step resolving the fast timescale ( h (cid:28) (cid:15) ), H is amesoscopic time step independent of the fast timescale satisfying h (cid:28) (cid:15) (cid:28) H (cid:28) h(cid:15) ) (cid:28) H (cid:28) h(cid:15) (44) Condition 5.2.
Consider the legacy ODE integrator with one-step update map Φ αt,t + τ introduced in (42) . Suppose there exists constants C > and H > independent of N and α , such that for any τ ≤ H min(1 /α, and bounded vector [ u , . . . , u N ] , (cid:107) Φ αt,t + τ ( u , . . . , u N ) − [ u , . . . , u N ] − τ [ f ( u , . . . , u N , α, t ) , . . .. . . , f N ( u , . . . , u N , α, t )] (cid:107) ≤ Cτ (1 + α ) , (45)Condition 5.2 corresponds to the assumption that the integrator Φ αt,t + τ is consistentfor (29).Observe that we are integrating (29) but not (36), since the remainders R i ’s are a-priori unknown unless the exact PDE solution is known. However, the following lemmaimplies that the FLAVORization of this integration is in fact convergent to the solutionof (36), even though R i ’s are possibly stiff. Lemma 5.2.
Assume that Φ αt,t + τ , introduced in (42) , satisfies Condition 5.2. Let h and H be the time steps used in the FLAVORization 5.1. If h (cid:28) (cid:15) , H (cid:28) h/(cid:15) , and K = O ( H ) , then (cid:107) Φ αt,t + τ ( u , . . . , u N ) − [ u , . . . , u N ] − τ [ f ( u , . . . , u N , α, t ) + R ( α, t ) , . . .. . . , f N ( u , . . . , u N , α, t ) + R N ( α, t )] (cid:107) ≤ Cτ (1 + α ) (46) where τ = h when α = (cid:15) − and τ = H − h when α = 0 . roof. By Condition 5.2, we have (cid:107) Φ αt,t + τ ( u , . . . , u N ) − [ u , . . . , u N ] − τ [ f ( u , . . . , u N , α, t ) , . . .. . . , f N ( u , . . . , u N , α, t )] (cid:107) ≤ Cτ (1 + α ) (47)for any τ ≤ min(1 /α, H . In addition, Lemma 5.1 gives a bound on the remainders:when α = (cid:15) − , there exists a constant ˜ C > N and (cid:15) − , such that forall i , | τ R i ( (cid:15) − , t ) | ≤ τ ˜ CK(cid:15) − (48)Because we use τ = h in this case and K (cid:28) (cid:15) − τ , the above is bounded by τ ˜ C ( ˆ C(cid:15) − τ ) (cid:15) − ≤ Cτ (1 + α ) for some constants ˆ C (cid:28) C = ˜ C ˆ C . When α = 0 on the other hand,there exists a constant ˜ C > i | τ R i ( (cid:15) − , t ) | ≤ τ ˜ CK (49)Because K = O ( H ) and we use τ = H − h = O ( H ) in this case, the above is bounded by τ ˜ C ˆ Cτ ≤ Cτ (1 + α ) for some constants ˆ C and we let C = ˜ C ˆ C . Notice that the valueof K is fixed in both cases but τ has different values: the flow map used in FLAVORassociated with α = 0 is the one with mesoscopic step Φ t + h,t + H , i.e., τ = H − h ; when α = (cid:15) − on the other hand, the flow map is Φ (cid:15) − t,t + h and τ = h . Finally, the triangleinequality gives (cid:107) Φ αt,t + τ ( u , . . . , u N ) − [ u , . . . , u N ] − τ [ f ( u , . . . , u N , α, t ) + R ( α, t ) , . . .. . . , f N ( u , . . . , u N , α, t ) + R N ( α, t )] (cid:107) ≤ (cid:107) Φ αt,t + τ ( u , . . . , u N ) − [ u , . . . , u N ] − τ [ f ( u , . . . , u N , α, t ) , . . . , f N ( u , . . . , u N , α, t )] (cid:107) + 1 N N (cid:88) i =1 | τ R i ( α, t ) | ≤ Cτ (1 + α ) , (50)which finished the proof after absorbing the coefficient 2 into C .We also need the usual regularity and stability assumptions to prove the accuracy ofFLAVORS for (36). Condition 5.3.
Assume that1. f , f , . . . , f N are Lipschitz continuous.2. For all bounded initial condition [ u (0) , . . . , u N (0)] ’s, the exact trajectories ([ u ( t ) , . . . , u N ( t )]) ≤ t ≤ T (i.e., solution to (36) ) are uniformly bounded in (cid:15) .3. For all bounded initial condition [ u (0) , . . . , u N (0)] ’s, the numerical trajectories ([¯ u ( t ) , . . . , ¯ u N ( t )]) ≤ t ≤ T (defined by (43) ) are uniformly bounded in (cid:15) , < H ≤ H , h ≤ min( H (cid:15), H ) . x and in the sense of measures on fast ones y , see [32]) of FLAVORs under the aboveconditions. Theorem 5.1.
Consider FLAVOR trajectories in Definition 5.1. Under Conditions 5.1,5.2 and 5.3, there exist
C > , ˆ C > and H > independent from (cid:15) − and N , suchthat for K/ ˆ C < H < H , h < H (cid:15) and t > , (cid:107) x ( t ) − [ η t ] x (¯ u ( t ) , . . . , ¯ u N ( t )) (cid:107) ≤ Ce Ct χ ( u (0) , . . . , u N (0) , (cid:15), H, h ) (51) and for all bounded and uniformly Lipschitz continuous test functions ϕ : R N (cid:55)→ R , (cid:12)(cid:12)(cid:12)(cid:12) t (cid:90) t +∆ tt ϕ ([¯ u ( s ) , . . . , ¯ u N ( s )]) ds − (cid:90) R p ϕ ([ η t ] − ( x ( t ) , y )) µ t ( x ( t ) , dy ) (cid:12)(cid:12)(cid:12)(cid:12) ≤ χ ( u (0) , . . . , u N (0) , (cid:15), H, h, ∆ t, t )( (cid:107) ϕ (cid:107) L ∞ + (cid:107)∇ ϕ (cid:107) L ∞ ) (52) where χ and χ are bounded functions converging towards zero as (cid:15) ≤ H/ ( C ln H ) , h(cid:15) ↓ , (cid:15)h H ↓ and ( h(cid:15) ) H ↓ (and ∆ t ↓ for χ ). Recall notations:
N K = L is the fixed spatial width, [ η t ] x and [ η t ] − respectivelydenote the x (slow) component and the inverse of the diffeomorphism η t (defined inCondition 5.1), x ( t ) = [ η t ] x ( u ( t ) , . . . , u N ( t )) corresponds to the slow component of theexact PDE solution sampled at grid points. u i ( t ) and ¯ u i ( t ) represent the exact and theFLAVOR approximation of the solution to the semi-discrete system with the remainders(36). Proof.
The proof of Theorem 5.1 is analogous to that of Theorem 1.2 of [32] (which willnot be repeated here). The proof requires (46), which is guarantied from Condition 5.2by Lemma 5.2. It is easy to check that the slow dependence on time of f , g , η and µ does not affect the proof given in [32]. Remark 5.3.
Condition 5.2 implies that the constant C in Theorem 5.1 does not dependon N or K . This is important because although using a finer mesh leads to a smaller K and a larger N = L/K , Condition 5.2 (which is equivalent to the accuracy of the semi-discrete approximation of the PDE) ensures that, as long as K = O ( H ) and h (cid:29) (cid:15)H ,the constant C in the error bounds on the slow component (51) and the fast component(52) will not blow up. Remark 5.4.
Observe that the application of the FLAVOR strategy does not require theidentification of the diffeomorphism η (which may depend on the spatial discretization). The convergence result of the previous section is based on the semi-discretization ofthe original PDE. PDEs and ODEs are also naturally connected via the method of19haracteristics, and henceforth it is natural to wonder whether a numerical integrationof those characteristics by FLAVORs would lead to an accurate approximation of thesolution of the original PDE. The answer to this question will be illustrated by analyzingthe following (generic) PDE: (cid:40) F ( Du, u, q, (cid:15) − ) = 0 , q ∈ Uu ( q ) = γ ( q ) , q ∈ Γ (53)where U ⊂ R d is the domain in which solution is defined, Γ and γ define initial/boundaryconditions.The following condition corresponds to assuming that characteristics are well-posed. Condition 6.1.
Assume that1. The PDE F ( Du, u, q, (cid:15) − ) = 0 admits characteristics: ˙ q = f ( q, z, (cid:15) − ) (54)˙ z = g ( q, z ) (55) u ( q ( t )) = z ( t ) (56) where q ∈ U is a vector corresponding to coordinates of characteristics in thedomain of the PDE, and z corresponds to the unknown’s value along the charac-teristics.2. For arbitrary (cid:15) , any point in U is reachable from the initial condition via one andonly one characteristics. The following conditions correspond to the assumption of existence of (possibly hid-den) slow and locally ergodic fast variables for those characteristics.
Condition 6.2.
Consider ODE (54) . Assume that:1. There exists a z -dependent diffeomorphism η z : q (cid:55)→ [ x, y ] from R d onto R d − p × R p with uniformly bounded C , C derivatives with respect to both q and t , such that ( x, y ) satisfies (with z ( t ) given by (55) ) (cid:40) ˙ x = f ( x, y, z )˙ y = (cid:15) − f ( x, y, z ) (57) where f , f , and g have bounded C derivatives with respect to x , y and z , and u ([ η z ] − ( x, y )) has bounded C derivatives with respect to the (slow) variables x and z .2. There exists a family of probability measures µ z ( x, dy ) on R p indexed by x ∈ R d − p and z ∈ R , as well as a family of positive functions T (cid:55)→ E z ( T ) satisfying T →∞ E z ( T ) = 0 , such that for all x , y , z , T bounded and φ uniformly boundedand Lipschitz, the solution to ˙ Y t = f ( x , Y t , z ) Y = y (58) satisfies (cid:12)(cid:12)(cid:12) T (cid:90) T φ ( Y s ) ds − (cid:90) R p φ ( y ) µ z ( x , dy ) (cid:12)(cid:12)(cid:12) ≤ χ z (cid:0) (cid:107) ( x , y ) (cid:107) (cid:1) E z ( T )( (cid:107) φ (cid:107) L ∞ + (cid:107)∇ φ (cid:107) L ∞ )(59) where r (cid:55)→ χ z ( r ) is bounded on compact sets, and µ z has bounded derivative withrespect to z in total variation norm. The second item of Condition 6.2 corresponds to the assumption that the fast variable y is locally ergodic with respect to a family of measures µ drifted by the slow variables x and z .The following lemma shows that, under the above conditions, the solution of PDE(53) is nearly constant on the orbit of the fast components ( y ) of the characteristics. Lemma 6.1.
Under Conditions 6.1 and 6.2, for any fixed constant C (independent of (cid:15) − ), there exists a constant C independent of (cid:15) − , such that for any ≤ t ≤ C , ≤ t ≤ C and (fixed) x and z , (cid:12)(cid:12) u (cid:0) [ η z ] − ( x , Y ( t )) (cid:1) − u (cid:0) [ η z ] − ( x , Y ( t )) (cid:1)(cid:12)(cid:12) ≤ C (cid:15) (60) where Y ( t ) and Y ( t ) are two points on the orbit of ˙ Y ( t ) = f ( x , Y ( t ) , z ) .Proof. Under Conditions 6.1 and 6.2, it is known (we refer for instance to [29] or toTheorem 14, Section 3 of Chapter II of [30] or to [28]) that x and z converge as (cid:15) → x and ˜ z defined as the solution to the following ODEs with initial condition x and z (cid:40) ˙˜ x = (cid:82) f (˜ x, y, ˜ z ) µ ˜ z (˜ x, dy )˙˜ z = (cid:82) g ([ η ˜ z ] − (˜ x, y ) , ˜ z ) µ ˜ z (˜ x, dy ) (61)Therefore, writing y ( t ) the solution of ˙ y = (cid:15) − f (˜ x, y, ˜ z ), we have as (cid:15) → u ([ η ˜ z ( t ) ] − (˜ x ( t ) , y ( t ))) → ˜ z ( t ) (62)Now, taking the time derivative of ˆ u = u ◦ η − , we obtainˆ u x ˙˜ x + ˆ u y ˙ y + ˆ u z ˙˜ z = ˙˜ z + ˙ R ( (cid:15) ) (63)where R ( (cid:15) ) is a function of t that goes to 0 as (cid:15) → Y ( t ) = f ( x , Y ( t ) , z )= f (˜ x ( (cid:15)t ) , y ( (cid:15)t ) , ˜ z ( (cid:15)t )) + ∂f ∂ ˜ x (˜ x ( (cid:15)t ) − x ) + ∂f ∂ ˜ z (˜ z ( (cid:15)t ) − z ) + ∂f ∂y ( y ( (cid:15)t ) − Y ( t ))+ o ( (cid:15) ) + o ( y ( (cid:15)t ) − Y ( t )) 21y Taylor expansion, ˜ x ( (cid:15)t ) − x and ˜ z ( (cid:15)t ) − z are obviously O ( (cid:15) ). Applying Gronwall’slemma, we also obtain that y ( (cid:15)t ) − Y ( t ) = O ( (cid:15) ). Therefore,˙ Y ( t ) = f (˜ x ( (cid:15)t ) , y ( (cid:15)t ) , ˜ z ( (cid:15)t )) + O ( (cid:15) ) = (cid:15) ˙ y ( t ) + o ( (cid:15) ) (64)Combining Eq. 63 with Eq. 64, we obtain u (cid:0) η − ( x , Y ( t )) (cid:1) − u (cid:0) η − ( x , Y ( t )) (cid:1) = (cid:90) t t ˆ u y · ˙ Y ( t ) dt = (cid:15) (cid:90) t t ˆ u y · ˙ y dt + o ( (cid:15) )= (cid:15) (cid:18)(cid:90) t t ( ˙˜ z − ˆ u x ˙˜ x − ˆ u z ˙˜ z ) dt + R ( (cid:15) ) (cid:12)(cid:12)(cid:12) t t (cid:19) + o ( (cid:15) ) (65)Since ˆ u x , ˙˜ x , ˆ u t and ˙˜ z are bounded, and R ( (cid:15) ) is vanishing (and hence bounded), weconclude that the right hand side is O ( (cid:15) ). Condition 6.3.
Assume that the domain U is bounded (independently from (cid:15) − ). Lemma 6.2.
If Conditions 6.1, 6.2, and 6.3 hold, then every point in U is reachable bya characteristics from the initial condition in bounded time (independently from (cid:15) − ).Proof. From Condition 6.1, we already know that every point is reachable, and thereforeit suffices to show that hitting times do not blow up as (cid:15) →
0. Since x ( · ) converges to˜ x ( · ) (see proof of Lemma 6.1), by considering the x component of the characteristics(projected by η ), it becomes trivial to show that the hitting time converges to a fixedvalue (and hence, does not blow up). Using Condition 6.3, we conclude that that anypoint in U can be hit in (uniformly) bounded time from the initial condition.Analogously to the Integrator 5.1, a legacy integrator for (54) and (55) can be FLA-VORized, and shown to be convergent under regularity and stability conditions (anal-ogous to Condition 5.3) requiring f , f and g to be Lipschitz continuous and ˜ q ( t ) and˜ z ( t ) to be bounded. The convergence result is analogous to Theorem 5.1, modulo thefollowing change of notation: the slow index is now z instead of t , the original coor-dinates are q instead of u i , the vector field of the original coordinates is f instead of f i , and the dynamics of the slow index comes from the non-trivial drift of ˙ z = g ( q, z )instead of the trivial ˙ t = 1. We define ˜ u (˜ q ( t )) := ˜ z ( t ) for all t on each FLAVORizedcharacteristics [˜ q ( t ) , ˜ z ( t )]. Naturally, ˜ u is only defined at discrete points in the domain U . These discrete points, however, densely ‘fill’ the space in the sense that (as shown bythe proof of the following theorem) FLAVORied characteristics remain very close to ex-act characteristics ( x components are close in Euclidean distance, and y components areclose as well in terms of orbital distance induced by the infimum of point-wise Euclideandistances).By the two-scale convergence theorem, we can quantify: the strong convergence ofthe slow coordinate of the characteristics and the unknown’s value along the character-istics, and the weak convergence of fast coordinate of the characteristics. Finally, thesesingle characteristics’ ODE approximation error bounds can be transferred to the PDEapproximation error bounds by considering the entire family of characteristics startingfrom all points (in initial condition). 22 heorem 6.1. Write ˜ u (˜ q ) the solution obtained by FLAVORizing all characteristics.Under Conditions 6.1, 6.2, 6.3, the consistency and regularity and stability Conditionscorresponding to Conditions 5.2 and 5.3 (with the change of notation described above),there exist a constant C independent of (cid:15) − and q ∈ Γ , such that | ˜ u (˜ q ) − u (˜ q ) | ≤ Cχ ( q , γ ( q ) , (cid:15), δ, τ )(1 + χ ( q , γ ( q ) , (cid:15), δ, τ, T, t )) (66) for any ˜ q on any FLAVORized characteristics, where q ∈ Γ and γ ( q ) correspond to theinitial condition that leads to ˜ q via a FLAVORized characteristics, and χ and χ arevanishing error bound functions. Remark 6.1.
When Γ is compact (such as in the case of periodic boundary condition), χ and χ can be further chosen to be independent of q (hence ˜ q ) by taking a supremumover Γ. Proof.
By Condition 6.1, all q ∈ U can be traced back to q ∈ Γ through a charac-teristics. By Lemma 6.2, characteristics starting from q reach q in bounded time T .Using the two-scale convergence of the FLAVORization of these characteristics (a resultanalogous to Theorem 5.1), we deduce that the approximation error associated with ˜ z T (on FLAVORized characteristics) can be bounded Cχ (with respect to the true value u ( q ) = z T , the error Ce CT has been replaced by C because T is bounded).Now observe that ˜ q T (cid:54) = q T , where ˜ q T is the coordinate of the FLAVORized char-acteristics starting from q . As before, let [ x T , y T ] = η ( q T ) and [˜ x T , ˜ y T ] = η (˜ q T ). Theerror on the slow component is (cid:107) x T − ˜ x T (cid:107) ≤ Cχ . The possible large error on the fastcomponent is not a problem because we can look for a near-by point on the fast orbitwith introducing only an O ( (cid:15) ) error on the unknown’s value (Lemma 6.1): (cid:40) u ( η ( x T , y T )) = u ( η ( x T , y ∗ T )) + O ( (cid:15) ) y ∗ T = arg min Y t | ˙ Y t = f ( x T ,Y t ) (cid:107) ˜ y T − Y t (cid:107) (67)Since (cid:107) ˜ x T − x t (cid:107) is small, the local ergodic measures that represent the orbits given by˙ Y t = f ( x T , Y t ) and ˙ Y t = f (˜ x T , Y t ) will be small: (cid:107) µ ( x T , dy ) − µ (˜ x T , dy ) (cid:107) T.V. ≤ Cχ χ isby chain rule. Because ˜ y T is on the orbit of ˙ Y t = f (˜ x T , Y t ), we will have (cid:107) y ∗ T − η y (˜ q T ) (cid:107) ≤ Cχ χ .All together, we obtain | ˜ u (˜ q T ) − u (˜ q T ) | = | ˜ z T − u (˜ q T ) |≤ | ˜ z T − u ( q ) | + | u ( q T ) − u (˜ q T ) |≤ Cχ + C (cid:107)∇ ( u ◦ η ) (cid:107) ∞ ( (cid:107) x T − η x (˜ q T ) (cid:107) + (cid:107) y T − η y (˜ q T ) (cid:107) ) ≤ Cχ + C ( χ + χ χ ) = Cχ + Cχ χ (68) Remark 6.2.
To keep the presentation concise, we have written C all constants thatdo not depend on essential parameters. 23 emark 6.3. As shown above, u will be captured strongly. Du , on the other hand,depends on a derivative with respect to the fast variable, and therefore will only beconvergent in a weak sense. Relevance to an error analysis for PDE-FLAVORS
The above result guaran-tees the convergence of FLAVORized characteristics. It is also possible to establish anerror bound on the difference between a specific PDE-FLAVOR discretization and theapproximation given by the above FLAVORized characteristics (and hence prove theconvergence of this specific PDE-FLAVOR discretization). Such an error bound couldbe obtained by first transforming FLAVORized characteristics to PDE-FLAVOR gridpoints via interpolating functions, and then using the fact that coordinate transforma-tions do not affect the efficiency of FLAVORS. For the sake of conciseness, we did notelaborate on this point here.
This work is supported by NSF grant CMMI-092600. We thank Guo Luo for stimulatingdiscussions and Sydney Garstang for proofreading the manuscript.
References [1]
Z. Artstein, C. W. Gear, I. G. Kevrekidis, M. Slemrod, and E. S. Titi , Analysis and computation of a discrete kdv-burgers type equation with fast dispertionand slow diffusion . arXiv:0908.2752, 2009.[2]
L. Berlyand and H. Owhadi , Flux norm approach to finite dimensional homog-enization approximations with non-separated scales and high contrast , Archives forRational Mechanics and Analysis, 198 (2010), pp. 677–721.[3]
T. J. Bridges , A geometric formulation of the conservation of wave action andits implications for signature and the classification of instabilities , Proc. Roy. Soc.Lond. A, 453 (1997), pp. 1365–1395.[4]
T. J. Bridges , Multi-symplectic structures and wave propagation , Math. Proc.Camb. Phil. Soc., 121 (1997), pp. 147–190.[5]
T. J. Bridges , Canonical multi-symplectic structure on the total exterior algebrabundle , Proc. Roy. Soc. Lond. A, 462 (2006), pp. 1531–1551.[6]
T. J. Bridges and G. Derks , Unstable eigenvalues and the linearization aboutsolitary waves and fronts with symmetry , Proc. Roy. Soc. Lond. A, 455 (1999),pp. 2427–2469. 247]
T. J. Bridges and S. Reich , Multi-symplectic integrators: numerical schemes forHamiltonian PDEs that conserve symplecticity , Phys. Lett. A, 284 (2001), pp. 184–193.[8]
T. J. Bridges and S. Reich , Numerical methods for Hamiltonian PDEs , J. Phys.A, 39 (2006), pp. 5287–5320.[9]
J. G. Charney, R. Fj ¨Ortoft, and J. von Neumann , Numerical integration ofthe barotropic vorticity equation , Tellus, 2 (1950), pp. 237–254.[10]
J.-B. Chen , Symplectic and multisymplectic Fourier pseudospectral discretizationsfor the Klein–Gordon equation , Lett. Math. Phys., 75 (2006), pp. 293–305.[11]
R. Courant, K. Friedrichs, and H. Lewy , ¨Uber die partiellen Differenzengle-ichungen der mathematischen Physik , Math. Ann., 100 (1928), pp. 32–74.[12] W. E, B. Engquist, X. Li, W. Ren, and E. Vanden-Eijnden , Heterogeneousmultiscale methods: a review , Commun. Comput. Phys., 2 (2007), pp. 367–450.[13]
W. E, D. Liu, and E. Vanden-Eijnden , Nested stochastic simulation algorithmsfor chemical kinetic systems with multiple time scales , J. Comput. Phys., 221 (2007),pp. 158–180.[14]
B. Engquist and Y.-H. Tsai , Heterogeneous multiscale methods for stiff ordinarydifferential equations , Math. Comp., 74 (2005), pp. 1707–1742 (electronic).[15]
F. Filbet and S. Jin , A class of asymptotic preserving schemes for kinetic equa-tions and related problems with stiff sources , (2010). arXiv:0905.1378. Accepted byJ. Comput. Phys.[16]
E. Hairer, S. P. Nørsett, and G. Wanner , Solving ordinary differential equa-tions. I , vol. 8 of Springer Series in Computational Mathematics, Springer-Verlag,Berlin, second ed., 1993. Nonstiff problems.[17]
J. S. Hesthaven, S. Gottlieb, and D. Gottlieb , Spectral Methods for Time-Dependent Problems , vol. 21 of Cambridge Monographs on Applied and Computa-tional Mathematics, Cambridge University Press, United Kingdom, 2007.[18]
S. Jin, L. Pareschi, and G. Toscani , Uniformly accurate diffusive relaxationschemes for multiscale transport equations , SIAM J. Numer. Anal., 38 (2000),pp. 913–936 (electronic).[19]
P. Lax and R. Richtmyer , Survey of the stability of linear finite difference equa-tions , Comm. Pure Appl. Math., 9 (1956), pp. 267–293.[20]
A. Lew, J. Marsden, M. Ortiz, and M. West , Asynchronous variational in-tegrators , Arch. Ration. Mech. Anal., 167 (2003), pp. 85–146.2521]
J. E. Marsden, G. W. Patrick, and S. Shkoller , Multisymplectic geometry,variational integrators, and nonlinear PDEs , Commun. Math. Phys., 199 (1998),pp. 351–395.[22]
J. E. Marsden and S. Shkoller , Multisymplectic geometry, covariant hamilto-nians, and water waves , Math. Proc. Camb. Phil. Soc., 125 (1999), pp. 553–575.[23]
J. E. Marsden and M. West , Discrete mechanics and variational integrators ,Acta Numerica, (2001), pp. 357–514.[24]
R. McLachlan, G. Reinout, and W. Quispel , Splitting methods , Acta Numer-ica, (2002), pp. 341–434.[25]
G. Naldi and L. Pareschi , Numerical schemes for hyperbolic systems of con-servation laws with stiff diffusive relaxation , SIAM J. Numer. Anal., 37 (2000),pp. 1246–1270 (electronic).[26]
H. Owhadi and L. Zhang , Metric-based upscaling , Comm. Pure Appl. Math., 60(2007), pp. 675–723.[27] ,
Localized bases for finite dimensional homogenization approximations withnon-separated scales and high-contrast , (2010). arXiv:1011.0986.[28]
G. A. Pavliotis and A. M. Stuart , Multiscale methods , vol. 53 of Texts inApplied Mathematics, Springer, New York, 2008. Averaging and homogenization.[29]
J. A. Sanders and F. Verhulst , Averaging methods in nonlinear dynamicalsystems , vol. 59 of Applied Mathematical Sciences, Springer-Verlag, New York,1985.[30]
A. Skorokhod , Asymptotic methods in the theory of stochastic differential equa-tions , vol. 78 of Translations of Mathematical Monographs, American MathematicalSociety, Providence, RI, 1989. Translated from the Russian by H. H. McFaden.[31]
J. C. Strikwerda , Finite difference schemes and partial differential equations ,Society for Industrial and Applied Mathematics, 2nd ed., 2004.[32]