1 Introduction
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 :
(1) 
Equation (1
) 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:
(2) 
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
(3) 
As is well known and immediate to see, a solution for the implicit scheme (3) can be found via the following optimization problem
(4) 
since (3) is the EulerLagrange equation for the optimization (4); here, . It follows that
(5) 
so that scheme (3) is unconditionally stable, provided that optimization problem (4) can be solved.
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 multistep methods and RungeKutta 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 blackbox 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 wellknown ordinary and partial differential equations that are gradient flows.
The code for section Section 5 is publicly available, and can be found at https://github.com/AZaitzeff/gradientflow.
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:
(6) 
where the intermediate stages , for , are given by
(7) 
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 twostage special case of scheme (6) & (7):
(8)  
(9) 
and impose the conditions
(10) 
on the parameters. Set . Note that (9) is equivalent to
(11) 
We have
establishing unconditional energy stability of scheme (8) & (9) under the condition (10) on its parameters. We will now extend this discussion to general, stage case of scheme (6) & (7):
Theorem 1
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:
Lemma 1
Proof
Note that
Then
Lemma 2
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 EulerLagrange equation:
(14) 
The consistency equations for the ’s are found by Taylor expanding around (or equivalently ). Set . We will calculate the onestep 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 :
(15) 
We now present the error at each stage of the multistage 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:
(16) 
where the coefficients obey the following recursive relation
(17) 
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:  
(18)  
Proof
We will now show by induction that the aforementioned consistency formulas, Eq. 16 and Eq. 17, hold.
Stage one:
(19) 
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
Noting that
completes stage one.
First we are going to solve for in Eq. 20:
(21) 
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 twostage method. However, it turns out that three stages are sufficient for unconditional stability:
(23) 
This choice of ’s that endows the threestage 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 :
(24) 
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:
(25) 
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 nonlinear ordinary and partial differential equations. The corresponding energies include convex and nonconvex forms. Careful numerical convergence studies are presented in each case to verify the anticipated convergence rates of previous sections.
Remark 1
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.
Number of  

time steps  
Error at  1.3e04  3.3e05  8.1e06  2.0e06  5.1e07  1.3e07 
Order    2.0  2.0  2.0  2.0  2.0 
Number of  

time steps  
Error at  1.1e06  1.3e07  1.7e08  2.1e09  2.6e10  3.2e11 
Order    3.0  3.0  3.0  3.0  3.0 
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.
Number of  

time steps  
Error at  5.2e04  1.3e04  3.3e05  8.2e06  2.0e06  5.1e07 
Order    2.0  2.0  2.0  2.0  2.0 
Number of  

time steps  
Error at  1.2e05  1.5e06  1.8e07  2.3e08  2.9e09  3.6e10 
Order    3.0  3.0  3.0  3.0  3.0 
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 (
Comments
There are no comments yet.