We are concerned with numerical schemes for evolution equations that arise as gradient flow (steepest descent) for an energy , where is a Hilbert space with inner product :
) may represent a (scalar or vectorial) ordinary or partial differential equation. Our focus is on unconditionally energy stable, high order in time discretizations. To be precise, by energy stable we mean the following dissipation property:
where denotes the approximation to the solution at the -th time step. Thus, in the context of PDEs, where is infinite dimensional, we are concerned with discrete in time, continuum in space schemes.
The backward Euler method for the abstract equation (1), with time step size , reads
As is well known and immediate to see, a solution for the implicit scheme (3) can be found via the following optimization problem
Energetic formulation (4) of the backward Euler scheme (3) is often referred to as minimizing movements. It enables extending numerical schemes for the stationary optimization problem to the dynamic, evolutionary problem (1) provided an additional, typically quadratic term in the cost function can be accommodated. The quadratic term in (4) is often referred to as the movement limiter, as it opposes deviation from the current configuration . It encodes the inner product with respect to which the gradient flow is being generated. Beyond numerical analysis and computation, minimizing movements approximation of gradient flows have been instrumental in the analysis of evolution equations of the form (1), e.g. in defining and finding weak solutions beyond the formation of singularities when classical notions of solution cease to exist.
There are many general purpose numerical methods that can certainly be used for solving Eq. 1. For instance, among high order schemes, some of the most well known are linear multi-step methods and Runge-Kutta methods, which can be regarded as special cases of the wider class of generalized linear methods [butcher2016numerical]. However, the energy stability of such general purpose numerical schemes is not immediate, and needs to be studied on a case by case basis. Among methods that focus on unconditional energy stability are convexity splitting [eyre1998unconditionally] and the more recent scalar auxiliary variable approach [shen2018scalar]. The following combination of desirable properties distinguish the new schemes introduced in this paper:
Complete generality. There is no assumption (e.g. convexity) on the energy in (1) beyond sufficient differentiability.
Unconditional energy stability.
High (at least up to third) order accuracy.
Each time step requires a few standard minimizing movements solves, equivalent to backward Euler substeps, or optimization of the associated energy plus a quadratic term.
Property 4 is perhaps the most unique and appealing aspect of the new framework: There are many existing schemes that can be recognized as some form of minimizing movements, sometimes relying on efficient optimization algorithms to solve (3) via (4). Our contribution shows how to painlessly jack up the order of accuracy of these schemes while preserving unconditional stability, relying only on a black-box implementation of the standard backward Euler scheme. In that sense, our new schemes can be understood as a variational analogue of Richardson extrapolation on Eq. 3, which in its standard form lacks the stability guarantees of our new schemes.
The rest of the paper is organized as follows:
Section 2 presents the general framework for the new scheme, focusing on unconditional energy stability.
Section 3 focuses on consistency, showing how to attain 2nd and 3rd order accuracy.
Section 4 gives 2nd and 3rd order examples of the new schemes.
Section 5 presents numerical convergence studies on a number of well-known ordinary and partial differential equations that are gradient flows.
2 The New Schemes: Stability
In this section, we formulate a wide class of numerical schemes that are energy stable by construction. We thus place stability front and center, leaving consistency to be dealt with subsequently. It is therefore important to allow many degrees of freedom in the scheme at this stage, in the form of a large number of coefficients, that will eventually be chosen, in the next section, to attain consistency at a high order of accuracy.
Our method is a linear -stage scheme of the following form:
where the intermediate stages , for , are given by
with the proviso . Notice that the proposed scheme Eq. 6 & Eq. 7, as promised, merely requires the solution of exactly the same type of problem at every time step as the standard backward Euler scheme: minimization of the associated energy plus a quadratic term.
At this point, it is not clear why a scheme such as Eq. 6 & (7) should dissipate energy at every iteration as in Eq. 2. However, in this section we establish quite broad conditions on the coefficients that ensure energy dissipation Eq. 2; this is the essential observation at the heart of the present paper. To demonstrate the idea, consider the following two-stage special case of scheme (6) & (7):
and impose the conditions
on the parameters. Set . Note that (9) is equivalent to
As we will see in Section 3, the conditions on the parameters of scheme (6) & (7) imposed in Claim 1 are loose enough to enable meeting consistency conditions to high order. We will establish Claim 1 with the help of the following couple of lemmas:
Let the auxiliary quantities , and be defined as in Theorem 1. We have
Let the auxiliary quantities , be given in Theorem 1 and let for . Then
3 The New Schemes: Consistency
We now turn to the question of whether the coefficients in scheme Eq. 6 and Eq. 7 can be chosen to ensure its high order consistency with the abstract evolution law Eq. 1. From Eq. 7, each stage satisfies the Euler-Lagrange equation:
The consistency equations for the ’s are found by Taylor expanding around (or equivalently ). Set . We will calculate the one-step error. For , let denote the multilinear form given by
so that the linear functional may be identified with an element of , and so on.
We begin with the Taylor expansion of the exact solution around :
We now present the error at each stage of the multi-stage algorithm, Eq. 6 and Eq. 7, and the conditions required to achieve various orders of accuracy: Let and be given in Eq. 6 and Eq. 7. The Taylor expansion of at each stage has the same form as Eq. 15, namely:
where the coefficients obey the following recursive relation
with . Furthermore, the following conditions for in scheme Eq. 6 are necessary and sufficient for various orders of accuracy:
|First Order:||Second Order:||Third Order:|
We first will Taylor expand around in Eq. 19:
Now we plug in an ansatz for the expansion on around , , and solve for , and :
Matching terms of the same order we get
completes stage one.
First we are going to solve for in Eq. 20:
Now Taylor expand around in Eq. 21:
Plug in the ansatz for and equation Eq. 16 for , and retaining up to terms of third order, we have that
Solving for , , by matching terms of the same order in Eq. 22, we arrive at:
completing the induction step.
4 The New Schemes: Examples
In this section, we exhibit second order and third order examples of scheme Eq. 6 that satisfy concurrently the hypothesis guaranteeing unconditional energy stability (Theorem 1) and the consistency equations (Section 3) up to second and third order. We found the ’s numerically and then sought a nearby algebraic solution to the consistency equations that still satisfied the conditions in Theorem 1.
4.1 Second order examples
It can be shown that there is no unconditionally energy stable second order two-stage method. However, it turns out that three stages are sufficient for unconditional stability:
This choice of ’s that endows the three-stage method Eq. 6 and Eq. 7 with unconditional stability and second order accuracy is by no means unique; indeed, here is another that has the additional benefit of having each one of its stages depend only on the previous one and :
4.2 Third order examples
We now exhibit a six stage version of scheme Eq. 6 and Eq. 7 that concurrently satisfies the conditions for unconditional energy stability (Theorem 1) as well the consistency equations (Section 3) up to third order:
The exact values of ’s above are given in the appendix (Section 7); they are all rational numbers but with long fractional representations. Again, we cannot rule out other solutions for , possibly with fewer stages.
5 The New Schemes: Numerical Tests
In this section, we will apply the second order Eq. 23 and third Eq. 25 order accurate unconditionally stable schemes to a variety of gradient flows. We found Eq. 23 before Eq. 24 and therefore ran all our numerical tests with the former. The gradient flows considered span linear and non-linear ordinary and partial differential equations. The corresponding energies include convex and non-convex forms. Careful numerical convergence studies are presented in each case to verify the anticipated convergence rates of previous sections.
5.1 Ordinary Differential Equations
Our first test is on the simple equation that corresponds to gradient flow for the scalar energy with respect to the standard inner product on . We take the initial condition in our numerical tests, so that the exact solution is . Table 1 and Table 2 show the error in the solution at time computed by the second order scheme (6), (7) & (23) and the third order scheme (6), (7) & (25), respectively, at various choices of the time step size. The anticipated order or convergence is clearly observed for both schemes.
Next, for a less trivial example, we turn to the ODE with the corresponding energy . With initial condition , the exact solution is . The errors for the two new schemes are tabulated in Table 3 and Table 4, and once again bear out the anticipated convergence rates.
For PDEs, we start with a preliminary test on the one dimensional heat equation on subject to periodic boundary conditions with initial data . This is gradient flow with respect to the inner product for the energy . The exact solution is . The spatial domain is discretized into a uniform grid of 2048 points, and a high order accurate discretization for the Laplacian is chosen so that the contribution to the error from spatial discretization is negligible. Table 5 and Table 6 show the error in the approximate solution at , computed via the second order accurate scheme (6), (7) & (23), and the third order accurate scheme (