Space-time FLAVORS: finite difference, multisymlectic, and pseudospectral integrators for multiscale PDEs
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  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 pre-existing 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) Structure-preserving: for stiff Hamiltonian PDEs (possibly on manifolds), they can be made to be multi-symplectic, symmetry-preserving (symmetries are group actions that leave the system invariant) in all variables and variational.
Multi-scale PDEs can be divided into two (possibly over-lapping) 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) . When fast variables converge toward Dirac (single point) distributions, asymptotic-preserving schemes  allow for simulations with large time steps. We also refer to [18, 25] for multi-scale transport equations and hyperbolic systems of conservation laws with stiff diffusive relaxation. Well-identified slow variables can be simulated with large time-steps using the two-scale structure of the original stiff PDEs (we refer to  and  for existing examples; slow variables satisfy a non-stiff PDE that can be identified in analogy to equations (A.9) and (A.13) of ; we also refer to  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  for stiff ODEs and SDEs) to stiff PDEs. Multi-scale integrators for stiff PDEs are obtained without the identification of slow variables by turning on and off stiff coefficients in single-step (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, multi-symplectic 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 semi-discretization in space, is analyzed in Subsection 5.1, where a non-asymptotic error bound indicates the two-scale convergence (, i.e., strong with respect to hidden slow variables and weak with respect to hidden fast variables) of PDE-FLAVORS. As illustrated by numerical (Figure 2) and theoretical results (Section 5), an explicit tuning () between microscopic and mesoscopic () time-steps 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 space-time FLAVOR mesh
2.1 Single-scale method and limitation
Consider a multiscale PDE:
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 single-scale 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  has to be met to ensure stability, which is also a neccessary condition for accuracy . Intuitively, the CFL condition guarantees that information does not propagate faster than what the numerical integrator can handle. The Von Neumann stability analysis  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  for additional discussions on single-scale 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 non-uniform mesh, grid sizes have to be taken into consideration when derivatives are approximated by finite differences. 1st-order 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 space-step , a microscopic time-step , and a mesoscopic time-step .
The intuition is as follows: adopt the point of view of semi-discrete 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 Ginzburg-Landau source
Consider a specific stiff PDE:
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 Lax-Friedrichs finite difference scheme:
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:
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 Lax-Friedrichs with different and values are computed by comparing the results to a benchmark Lax-Friedrichs integration with fine steps and . More precisely, we calculated the distance between two vectors respectively corresponding to FLAVOR and Lax-Friedrichs integrations, which contain ordered values on the intersection of FLAVOR and Lax-Friedrichs meshes (which is in fact the FLAVOR mesh as long as is a multiple of ). 1-norm 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 Single-scale method
where is a n-dimensional vector, and are arbitrary skew-symmetric matrices on , and is an arbitrary smooth function. The solution preserves the multi-symplectic structure in the following sense:
where and are differential 2-forms defined by
and and are two arbitrary solutions to the variational equation (the solution is identified with ):
Preservation of multi-symplecticity 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, fluid-structure interactions, etc. [4, 3, 6, 7]. We also refer to  and references therein for surveys on numerical recipes, and to  for an application to numerical nonlinear elastodynamics.
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 ) to the following action:
where the Lagrangian density is given by
This variational view of Hamiltonian PDEs will intrinsically guarantee the preservation 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 variational integrators because they are intrinsically multi-symplectic and therefore structure-preserving [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 Euler-Lagrange equations.
For an illustration, consider a nonlinear wave equation:
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 ):
The corresponding Lagrangian density is:
Using a forward time forward space approximation, we obtain the following discrete Lagrangian:
where space step and time step define a rectangular grid of size . The simplest single-scale choice would be and for some and .
As a consequence, the continuous action is approximated by a discrete action:
and Hamilton’s principle of least action gives
for and , where and are such that for any and for any .
Taking derivative with respect to , we obtain the following discrete Euler-Lagrange equations:
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:
This numerical receipt is convergent. In fact, multi-symplectic 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 high-order derivatives are dealt with in an intrinsic way regardless of whether mesh is uniform.
3.2 FLAVORization of multi-symplectic integrators
Now consider a multiscale Hamiltonian PDE
Any single-scale multi-symplectic integrator can be FLAVORized (to achieve computational acceleration) by using the following strategy: (i) Use the two-scale mesh illustrated 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 instead of the PDE. Specifically, choose
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 time-step, and corresponds to a mesoscopic space-step; the same rule of thumb for choosing them in Section 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 that stiffness is turned off in half of the grids. Multisymplecticity is automatically gained, because the updating equations originate from a discrete variational principle .
3.3 Example: multiscale Sine-Gordon wave equation
Consider a specific nonlinear wave equation (11) in which . If , this corresponds to Sine-Gordon equation, which has been studied extensively due to its soliton solutions and its relationships with quantum physics (for instance, as a nonlinear version of Klein-Gordon 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 single-scale forward time forward space multi-symplectic 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 well-approximated 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  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 single-scale forward time forward space multi-symplectic 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 Single-scale method
Consider a PDE
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:
and solves for ’s by requiring the PDE to hold at collocation points :
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  for additional details on Fourier collocation methods. It is worth mentioning that pseudospectral methods can also be multi-symplectic when applied to Hamiltonian PDEs .
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 space-step 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 non-Dirac fast process
Consider the following system of PDEs
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 non-trivial way.
We use the classical 4th order Runga-Kutta scheme (see, for instance, ) for the (single-step) time integration of the pseudospectrally discretized system of ODEs (26). Write its numerical flow over a microscopic time step (consisting of four sub-steps), 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 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 slow process of is captured in strong (point-wise) 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 single-step integration uses and (notice that this is already beyond the stability/accuracy region of a single-scale finite difference, since the space step does not depend on ; the spectral method is more stable/accurate for a large space-step), and FLAVOR uses , and . -fold acceleration is achieved by FLAVOR.
5 Convergence analysis
5.1 Semi-discrete system
All FLow AVeraging integratORS described in previous sections are illustrations of the following (semi-discrete) 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 ODE-FLAVOR . In this section, we will use the semi-discrete ODE system as an intermediate link to demonstrate that these PDE-FLAVORS 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 single-scale (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 1st-order in time derivative:
Observe that a PDE (1) with higher-order time derivatives can be written as a system of 1st-order (in time derivatives) PDEs.
Now consider a consistent discretization of PDE (28) with space step and time step (we refer to Page 20 of  for a definition of the notion of consistency, which intuitively means vanishing local truncation error). Letting in this discretization, we obtain a semi-discrete system (continuous in time and discrete in space). This semi-discrete system is denoted by the following system of ODEs, with approximated spatial derivatives:
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:
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 ):
(31) is true, for instance, in cases where
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
and naturally establishes the correspondence of and for a 1st-order 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 ODE-FLAVORs. The consistency of pseudospectral method can be shown similarly using Fourier analysis.
With defined in (30), consider the following system of ODEs:
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 PDE-FLAVOR by showing that an ODE-FLAVOR 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 ):
Observe that if is Riemann integrable, then
and hence the norm (37) does not blow up or vanish as .
5.2 Sufficient conditions for convergence, ODE-FLAVORS, and two-scale convergence of PDE-FLAVORS
We will now prove the accuracy of PDE-FLAVORs under the assumption of existence of (possibly hidden) slow and locally ergodic fast variables. The convergence of PDE-FLAVORs will be expressed using the notion of two-scale flow convergence introduced in  (corresponding to a strong convergence with respect to slow variables and weak one with respect to fast ones).
Assume that the ODE system (36) satisfies the following conditions:
(Existence of hidden slow and fast variables): There exists a (possibly time-dependent) diffeomorphism from onto with uniformly bounded derivatives with respect to ’s and , and such that for all , satisfies
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
where is bounded on compact sets, and has bounded derivative with respect to in total variation norm.
Write the numerical flow of a given (legacy) ODE integrator for (29):
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 (Ode-Flavors).
The FLow AVeraging integratOR associated with is defined as the algorithm simulating the process:
where (the number of steps) is a piece-wise 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
Consider the legacy ODE integrator with one-step update map introduced in (42). Suppose there exists constants and independent of and , such that for any and bounded vector ,
Observe that we are integrating (29) but not (36), since the remainders ’s are a-priori 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.
By Condition 5.2, we have
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 ,
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
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
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).
are Lipschitz continuous.
For all bounded initial condition ’s, the exact trajectories