Spacetime FLAVORS: finite difference, multisymlectic, and pseudospectral integrators for multiscale PDEs
Abstract
We present a new class of integrators for stiff PDEs. These integrators are generalizations of FLow AVeraging integratORS (FLAVORS) for stiff ODEs and SDEs introduced in [32] with the following properties: (i) Multiscale: they are based on flow averaging and have a computational cost determined by mesoscopic steps in 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 have hidden slow and fast processes). This bypasses the need for identifying explicitly (or numerically) the slow variables or reduced effective PDEs; (iii) Nonintrusive: A preexisting numerical scheme resolving the microscopic time scale can be used as a black box and easily turned into one of the integrators in this paper by turning the large coefficients on over a microscopic timescale and off during a mesoscopic timescale; (iv) Convergent over two scales: strongly over slow processes and in the sense of measures over fast ones; (v) Structurepreserving: for stiff Hamiltonian PDEs (possibly on manifolds), they can be made to be multisymplectic, symmetrypreserving (symmetries are group actions that leave the system invariant) in all variables and variational.
1 Introduction
Multiscale PDEs can be divided into two (possibly overlapping) categories: PDEs with highly oscillating or rough coefficients and PDEs with large (or stiff) coefficients. Classical numerical methods are usually: (i) stable but arbitrarily inaccurate for the former category (consider, for instance, a finite element method for the elliptic operator with a rapidly changing coefficient ), or (ii) unstable for the latter category. Accurate numerical methods for the former category, called numerical homogenization methods, are, in absence of local ergodicity or scale separation, based on the compactness of the solution space (we refer, for instance, to [26, 2, 27]). Numerical methods for the latter category are, in essence, based on the existence of slow and fast variables (or components) [14]. When fast variables converge toward Dirac (single point) distributions, asymptoticpreserving schemes [15] allow for simulations with large time steps. We also refer to [18, 25] for multiscale transport equations and hyperbolic systems of conservation laws with stiff diffusive relaxation. Wellidentified slow variables can be simulated with large timesteps using the twoscale structure of the original stiff PDEs (we refer to [1] and [12] for existing examples; slow variables satisfy a nonstiff PDE that can be identified in analogy to equations (A.9) and (A.13) of [32]; we also refer to [14] for a definition of slow variables).
In this paper, we consider the second category of PDEs and propose a generalization of FLow AVeraging integratORS (FLAVORS) (introduced in [32] for stiff ODEs and SDEs) to stiff PDEs. Multiscale integrators for stiff PDEs are obtained without the identification of slow variables by turning on and off stiff coefficients in singlestep (legacy) integrators (used as black boxes) and alternating microscopic and mesoscopic time steps (Subsection 2.2). We illustrate the generality of the proposed strategy by applying it to finite difference methods in Section 2, multisymplectic integrators in Section 3, 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 volume methods). The convergence of the proposed strategy, after semidiscretization in space, is analyzed in Subsection 5.1, where a nonasymptotic error bound indicates the twoscale convergence ([32], i.e., strong with respect to hidden slow variables and weak with respect to hidden fast variables) of PDEFLAVORS. As illustrated by numerical (Figure 2) and theoretical results (Section 5), an explicit tuning () between microscopic and mesoscopic () timesteps and the stiff parameter is necessary and sufficient for convergence. We also show in Section 6 that applying the FLAVOR strategy to characteristics leads to accurate approximations of solutions of stiff PDEs.
2 Finite difference and spacetime FLAVOR mesh
2.1 Singlescale method and limitation
Consider a multiscale PDE:
(1) 
where is a given function (possibly nonlinear), is a small positive real parameter and and are spatial and temporal coordinates.
To obtain a numerical solution of (1), the simplest singlescale finite difference approach employs a uniform rectangular mesh with time step length and space step length , and approximates the solution at its values at discrete grid points. Differential operators will be approximated by finite differences; for instance, according to forward space forward time rules: and , where is the numerical solution at discrete grid point with space index and time index . After this discretization, the original PDE is approximated by a finite dimensional algebraic system, which can be solved to yield the numerical solution.
Of course, a necessary condition for obtaining stability and accuracy in the numerical solution is that and have to be small enough. A quantitative statement on how small they need to be will depend on the specific PDE and discretization. For 1D linear advection equations and forward time forward space discretizations, the CFL condition [11] has to be met to ensure stability, which is also a neccessary condition for accuracy [19]. Intuitively, the CFL condition guarantees that information does not propagate faster than what the numerical integrator can handle. The Von Neumann stability analysis [9] helps determine analogous CFL conditions for linear equations with arbitrary discretizations. The stability of numerical schemes for general nonlinear equations remains a topic of study. We refer to [31] for additional discussions on singlescale finite difference schemes. In general, the presence of a stiff coefficient in equation (1) requires and to scale with in order to guarantee the stability of numerical integration schemes. This makes the numerical approximation of the solution of (1) computationally untractable when is close to 0.
2.2 Multiscale FLAVORization and general methodology
FLAVORs are multiscale in the sense that they accelerate computation by adopting both larger time and space steps. A finite difference scheme can be FLAVORized by employing two rules:
First, instead of a uniform mesh, use a mesh as depicted in Figure 1, in which a uniform spatial grid corresponds to a mesoscopic space step that does not scale with , and an alternating temporal grid corresponds to two time steps, microscopic (scaling with ) and mesoscopic ( independent from ). It is worth mentioning that when using this nonuniform mesh, grid sizes have to be taken into consideration when derivatives are approximated by finite differences. 1storder derivatives are straightforward to obtain, and we refer to Section 3 for approximations of higher order derivatives.
Second, the stiff parameter should be temporarily set to be 0 (i.e., turned off) when the current time step is the mesoscopic ; if the small time step is used instead, the large value of needs to be restored, or in other words, stiffness should be turned on again.
The rule of thumb is that and 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 stiffness turned off is stable and accurate, and and should be chosen to be an order of magnitude smaller than these values. FLAVORS does not require a microscopic , but only a mesoscopic spacestep , a microscopic timestep , and a mesoscopic timestep .
The intuition is as follows: adopt the point of view of semidiscrete approach for PDE integration, in which space is discretized first and the PDE is approximated by a system of ODEs. The integration (in the time) of the resulting finite dimensional ODE system 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 and alternating microscopic time steps (stiffness on) with mesoscopic time steps (stiffness on) preserves the symmetries of that scheme and at the same time induces an averaging of the dynamic of (possibly hidden) slow variables with respect to the fast ones. With this strategy, the FLAVORized scheme advances in mesoscopic time steps without losing stability. The (possibly hidden) slow dynamic is captured in a strong sense, while the fast one is captured only in the (weak) sense of measures. A rigorous proof of convergence of the proposed method relies on the assumption of existence of (possibly hidden) slow variables 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 identification of slow variables.
2.3 Example: conservation law with GinzburgLandau source
Consider a specific stiff PDE:
(2) 
in which and . Use the boundary condition of and the initial condition of . This system contains two scales: the fast process corresponds to quickly converging towards or , and the slow process corresponds to the front (with steep gradients) separating from propagating at an velocity.
We will FLAVORize the following LaxFriedrichs finite difference scheme:
(3) 
where and . If the domain of integration is restricted to , then , and . We use and for our purposes, both of which we found numerically at the order of the stability limit. In our experiment, we chose , and therefore and .
The FLAVORized version of this scheme is:
(4) 
where and . If the domain of integration is restricted to , then , and . We use the same as before, and choose and , which ensures that the stability of the integration remains independent of .
Errors of FLAVOR based on LaxFriedrichs with different and values are computed by comparing the results to a benchmark LaxFriedrichs integration with fine steps and . More precisely, we calculated the distance between two vectors respectively corresponding to FLAVOR and LaxFriedrichs integrations, which contain ordered values on the intersection of FLAVOR and LaxFriedrichs meshes (which is in fact the FLAVOR mesh as long as is a multiple of ). 1norm is used and normalized by the number of discrete points to mimic the norm for the continuous solution. Experimental settings are , and . As we can see in Figure 2, FLAVOR is indeed uniformly convergent in the sense that the error scales with , as long as takes an appropriate value. This is not surprising, because we have already proven in the ODE case that the error is bounded by a function of (uniformly in ) as long as , and this error can be made arbitrarily small as (notice can still be much larger than as ).
Also, a typical run of FLAVOR ( and ) in comparison to the benchmark ( and ) is shown in Figure 3. FLAVOR captured the slow process strongly in the sense that it obtained the correct speeds of both steep gradients’ propagations (up to arithmetic error and fringing). In this setting, FLAVOR achieves a fold acceleration. It is worth restating that both spatial and temporal step lengths of FLAVOR are mesocopic, whereas the counterparts in a single scale finite difference method have to be both microscopic for stability. The computational gain by FLAVOR will go to infinity as , and this statement will be true for all FLAVOR examples shown in this paper.
3 Multisymplectic integrator for Hamiltonian PDEs
3.1 Singlescale method
We refer to [7, 21, 22] for a discussion on the geometry of Hamiltonian PDEs (e.g., multisymplectic structure). We will now recall the Euclidean coordinate form of a Hamiltonian PDE:
(5) 
where is a ndimensional vector, and are arbitrary skewsymmetric matrices on , and is an arbitrary smooth function. The solution preserves the multisymplectic structure in the following sense:
(6) 
where and are differential 2forms defined by
(7) 
and and are two arbitrary solutions to the variational equation (the solution is identified with ):
(8) 
Preservation of multisymplecticity can be partially and intuitively interpreted as a conservation of infinitesimal volume in the jet bundle, which generalizes the conservation of phase space volume in Hamiltonian ODE settings to field theories.
A broad spectrum of PDEs fall in the class of Hamiltonian PDEs, including generalized KdV, nonlinear Schrödinger models, nonlinear wave equations, atmospheric flows, fluidstructure interactions, etc. [4, 3, 6, 7]. We also refer to [8] and references therein for surveys on numerical recipes, and to [20] for an application to numerical nonlinear elastodynamics.
Hamiltonian PDEs (5) can be viewed as EulerLagrange equations for field theories, which are obtained by applying Hamilton’s principle (i.e., a variational principle of ) to the following action:
(9) 
where the Lagrangian density is given by
(10) 
This variational view of Hamiltonian PDEs will intrinsically guarantee the preservation of multisymplecticity, 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 variational integrators because they are intrinsically multisymplectic and therefore structurepreserving [21, 22, 23, 20]. These integrators are obtained as follows: first discretize the action (9) using quadratures, then apply variational principle to the discrete action (which depends on finitely many arguments), and finally, solve the algebraic system obtained from the variational principle, i.e., the discrete EulerLagrange equations.
For an illustration, consider a nonlinear wave equation:
(11) 
with periodic boundary condition and compatible initial conditions and . Suppose we are interested in the solution in a domain .
Rewrite the high order PDE as a system of first order PDEs (notice these covariant equations can be obtained through an intrinsic procedure, which works on manifolds as well [5]):
(12)  
(13)  
(14) 
The corresponding Lagrangian density is:
(15) 
Using a forward time forward space approximation, we obtain the following discrete Lagrangian:
(16)  
(17) 
where space step and time step define a rectangular grid of size . The simplest singlescale choice would be and for some and .
As a consequence, the continuous action is approximated by a discrete action:
(18) 
and Hamilton’s principle of least action gives
(19) 
for and , where and are such that for any and for any .
Taking derivative with respect to , we obtain the following discrete EulerLagrange equations:
(20) 
The system of above equations is explicitly solvable when equipped with boundary conditions and initial conditions; for instance, below is a consistent discretization of the continuous version:
(21) 
This numerical receipt is convergent. In fact, multisymplectic integrators obtained from 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 PDEs of form (5). Also, notice that highorder derivatives are dealt with in an intrinsic way regardless of whether mesh is uniform.
3.2 FLAVORization of multisymplectic integrators
Now consider a multiscale Hamiltonian PDE
(22) 
Any singlescale multisymplectic integrator can be FLAVORized (to achieve computational acceleration) by using the following strategy: (i) Use the twoscale mesh illustrated in Figure 1, and (ii) Turn off large coefficients when taking mesoscopic timesteps. Unlike FLAVORizing a general finite difference scheme, we FLAVORize the action instead of the PDE. Specifically, choose
(23) 
and let in for even ’s and all ’s, while the large value of is kept in for odd ’s and all ’s. and correspond to a small and a mesoscopic timestep, and corresponds to a mesoscopic spacestep; the same rule of thumb for choosing them in Section 2 applies.
After applying the discrete Hamilton’s principle, the resulting discrete Euler Lagrangeequations corresponding to a multisymplectic integrator will still be (20), except that stiffness is turned off in half of the grids. Multisymplecticity is automatically gained, because the updating equations originate from a discrete variational principle [21].
3.3 Example: multiscale SineGordon wave equation
Consider a specific nonlinear wave equation (11) in which . If , this corresponds to SineGordon equation, which has been studied extensively due to its soliton solutions and its relationships with quantum physics (for instance, as a nonlinear version of KleinGordon equation). We are interested in the case in which (identified with ) is big, so that a separation of timescale exhibits.
Arbitrarily choose and use periodic boundary condition , and let initial condition be and . Denote total simulation time by . Use the FLAVOR mesh (23). In order to obtain a stable and accurate numerical solution, and have to be , and and need to be .
A comparison between the benchmark of the singlescale forward time forward space multisymplectic integrator (Eq. 20 with and ) and its FLAVORization (Eq. 20 with mesh (23) and for odd and for even ) is presented in Figure 4. , and , and and . It is intuitive to say that the slow process of wave propagation is wellapproximated by FLAVOR, although the fast process of local fluctuation is not captured in the strong sense. Error quantification is not done, because what the slow and fast processes are is not rigorously known here. fold acceleration is obtained by FLAVOR.
Readers familiar with the splitting theory of ODEs [24] might question whether FLAVORS are equivalent to an averaged stiffness of (which corresponds in the numerical experiment described above). The answer is no, because the equivalency given by the splitting theory is only local. In fact, the same singlescale forward time forward space multisymplectic integration of the case is shown in Figure 5, which is clearly distinct from the FLAVOR result in Figure 4. Moreover, because of the error term, changing stiffness alone will not result in a converging method (and result in a error on slow variables).
4 Pseudospectral methods
4.1 Singlescale method
Consider a PDE
(24) 
with periodic boundary condition and initial condition , where is a differential operator involving only spatial derivatives.
The Fourier collocation method approximates the solutions by the truncated Fourier series:
(25) 
and solves for ’s by requiring the PDE to hold at collocation points :
(26) 
This yields ODEs, which can be integrated by any favorite ODE solver. Of course, specific choices of collocations points will affect the numerical approximation. Oftentimes, the simplest choice of is used, and in this case, the method is also called a pseudospectral method. We refer to [17] for additional details on Fourier collocation methods. It is worth mentioning that pseudospectral methods can also be multisymplectic when applied to Hamiltonian PDEs [10].
4.2 FLAVORization of pseudospectral methods
When the PDE is stiff (for instance, when contains a large parameter ), FLAVORS can be employed to integrate the stiff ODEs (which will still contain ) resulting from a pseudospectral discretization.
Similarly, for the FLAVORization of a pseudospectral method, it is sufficient to choose instead of , i.e., the spacestep can be coarse (). For time stepping, alternatively switching between and for a mesoscopic is again needed, and stiffness has to be turned off over the mesoscopic step of . In a sense, we are still using the same FLAVOR ‘mesh’ (Figure 1), except that here we do not discretize space, but instead truncate Fourier series to resolve the same spatial grid size.
4.3 Example: a slow process driven by a nonDirac fast process
Consider the following system of PDEs
(27) 
with periodic boundary conditions , , and , and initial conditions , , and . The integration domain is restricted to . The stiffness is identified with . We choose the initial condition of and .
In this system, and correspond to a fast process, which is a field theory version of a harmonic oscillator with high frequency . is a slow process, into which energy is pumped by the fast process in a nontrivial way.
We have chosen to FLAVORize (27) because it does not fall into the (simpler) category of systems with fast processes converging towards Dirac (single point support) invariant distributions [15].
We use the classical 4th order RungaKutta scheme (see, for instance, [16]) for the (singlestep) time integration of the pseudospectrally discretized system of ODEs (26). Write its numerical flow over a microscopic time step (consisting of four substeps), where are numerical approximations to the Fourier coefficients in (25), for the unknowns , and at an arbitrary time . Then, the corresponding FLAVOR update over a mesoscopic time step will be , which consists of eight substeps.
We present in Figure 6 and Figure 7 a comparison between the benchmark of singlescale pseudospectral simulation and its FLAVORization. It can be seen that the slow process of is captured in strong (pointwise) sense, whereas the fast process of is only approximated in a weak sense (i.e. as a measure, in the case, wave shape and amplitude are correct, but not the period). We choose , and . The singlestep integration uses and (notice that this is already beyond the stability/accuracy region of a singlescale finite difference, since the space step does not depend on ; the spectral method is more stable/accurate for a large spacestep), and FLAVOR uses , and . fold acceleration is achieved by FLAVOR.
5 Convergence analysis
5.1 Semidiscrete system
All FLow AVeraging integratORS described in previous sections are illustrations of the following (semidiscrete) strategy: first, space is discretized or interpolated; next, spatial differential operators are approximated by algebraic functions of finitely many spatial variables; finally, the resulting system of ODEs is numerically integrated by a corresponding ODEFLAVOR [32]. In this section, we will use the semidiscrete ODE system as an intermediate link to demonstrate that these PDEFLAVORS are convergent to the exact PDE solution under reasonable assumptions (in a strong sense with respect to (possibly hidden) slow variables and in the sense of measures with respect to fast variables).
More precisely, consider a spatial mesh (vector) , a temporal mesh (vector) , and a domain mesh (matrix) . Examples of these meshes include the FLAVOR mesh and , and a usual singlescale (step) integration mesh and (recall the domain size is by ). We will use the FLAVOR mesh throughout this section. We will compare the solution of the PDE (28) with the solution obtained with the FLAVOR strategy at these discrete points.
For simplicity, assume the PDE of interest is 1storder in time derivative:
(28) 
Observe that a PDE (1) with higherorder time derivatives can be written as a system of 1storder (in time derivatives) PDEs.
Now consider a consistent discretization of PDE (28) with space step and time step (we refer to Page 20 of [31] for a definition of the notion of consistency, which intuitively means vanishing local truncation error). Letting in this discretization, we obtain a semidiscrete system (continuous in time and discrete in space). This semidiscrete system is denoted by the following system of ODEs, with approximated spatial derivatives:
(29) 
Assuming existence and uniqueness of an exact strong solution to the PDE (28), and writing its values at the spatial discretization points, we define for each the following remainder:
(30) 
which is a real function of indexed by .
Then, approximates the exact solution evaluated at grid points in the sense that these remainders vanish as (where ):
Lemma 5.1.
Remark 5.1.
(31) is true, for instance, in cases where
(33) 
Proof.
Remark 5.2.
The consistency of finite difference methods can be easily shown using Taylor expansions. For instance, applying a Taylor expansion to the solution of leads to
(34) 
which implies
(35) 
and naturally establishes the correspondence of and for a 1storder finite difference scheme. Notice that the remainders are still stiff, but we will see later that this is not a problem, since they can be handled by ODEFLAVORs. The consistency of pseudospectral method can be shown similarly using Fourier analysis.
With defined in (30), consider the following system of ODEs:
(36) 
with initial condition . Obviously, its solution is the exact PDE solution sampled at spatial grid points, i.e., .
We will now establish the accuracy of PDEFLAVOR by showing that an ODEFLAVOR integration of (36) leads to an accurate approximation of . Since space (with fixed width ) is discretized by grid points, we use the following (normalized by ) norm in our following discussion (suppose for a function ):
(37) 
Observe that if is Riemann integrable, then
(38) 
and hence the norm (37) does not blow up or vanish as .
5.2 Sufficient conditions for convergence, ODEFLAVORS, and twoscale convergence of PDEFLAVORS
We will now prove the accuracy of PDEFLAVORs under the assumption of existence of (possibly hidden) slow and locally ergodic fast variables. The convergence of PDEFLAVORs will be expressed using the notion of twoscale flow convergence introduced in [32] (corresponding to a strong convergence with respect to slow variables and weak one with respect to fast ones).
Condition 5.1.
Assume that the ODE system (36) satisfies the following conditions:

(Existence of hidden slow and fast variables): There exists a (possibly timedependent) diffeomorphism from onto with uniformly bounded derivatives with respect to ’s and , and such that for all , satisfies
(39) where and have bounded derivatives with respect to , and .

(Local ergodicity of vast variables): There exists a family of probability measures on indexed by and , and a family of positive functions satisfying for all bounded , such that for all bounded and uniformly bounded and Lipschitz, the solution to
(40) satisfies
(41) where is bounded on compact sets, and has bounded derivative with respect to in total variation norm.
Under Conditions 5.1, the computation of the solution of PDE (28) can be accelerated by applying the FLAVOR strategy to a singlescale time integration of its semidiscretization (29).
Write the numerical flow of a given (legacy) ODE integrator for (29):
(42) 
where approximates for all , is the integration time step, and is a controllable parameter that replaces the stiff parameter in (29) and takes values of (stiffness ‘on’) or (stiffness ‘off’).
Definition 5.1 (OdeFlavors).
The FLow AVeraging integratOR associated with is defined as the algorithm simulating the process:
(43) 
where (the number of steps) is a piecewise constant function of satisfying , is a microscopic time step resolving the fast timescale (), is a mesoscopic time step independent of the fast timescale satisfying and
(44) 
Condition 5.2.
Consider the legacy ODE integrator with onestep update map introduced in (42). Suppose there exists constants and independent of and , such that for any and bounded vector ,
(45) 
Observe that we are integrating (29) but not (36), since the remainders ’s are apriori unknown unless the exact PDE solution is known. However, the following lemma implies that the FLAVORization of this integration is in fact convergent to the solution of (36), even though ’s are possibly stiff.
Lemma 5.2.
Proof.
By Condition 5.2, we have
(47) 
for any . In addition, Lemma 5.1 gives a bound on the remainders: when , there exists a constant independent of and , such that for all ,
(48) 
Because we use in this case and , the above is bounded by for some constants and . When on the other hand, there exists a constant such that for all
(49) 
Because and we use in this case, the above is bounded by for some constants and we let . Notice that the value of is fixed in both cases but has different values: the flow map used in FLAVOR associated with is the one with mesoscopic step , i.e., ; when on the other hand, the flow map is and . Finally, the triangle inequality gives
(50) 
which finished the proof after absorbing the coefficient into . ∎
We also need the usual regularity and stability assumptions to prove the accuracy of FLAVORS for (36).
Condition 5.3.
Assume that

are Lipschitz continuous.

For all bounded initial condition ’s, the exact trajectories