 # Variational Extrapolation of Implicit Schemes for General Gradient Flows

We introduce a class of unconditionally energy stable, high order accurate schemes for gradient flows in a very general setting. The new schemes are a high order analogue of the minimizing movements approach for generating a time discrete approximation to a gradient flow by solving a sequence of optimization problems. In particular, each step entails minimizing the associated energy of the gradient flow plus a movement limiter term that is, in the classical context of steepest descent with respect to an inner product, simply quadratic. A variety of existing unconditionally stable numerical methods can be recognized as (typically just first order accurate in time) minimizing movement schemes for their associated evolution equations, already requiring the optimization of the energy plus a quadratic term at every time step. Therefore, our approach gives a painless way to extend these to high order accurate in time schemes while maintaining their unconditional stability. In this sense, it can be viewed as a variational analogue of Richardson extrapolation.

## Authors

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 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 :

 u′=−∇HE(u). (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:

 E(un+1)≤E(un) (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

 un+1−unk=−∇HE(un+1). (3)

As is well known and immediate to see, a solution for the implicit scheme (3) can be found via the following optimization problem

 un+1=argminuE(u)+12k||u−un||2 (4)

since (3) is the Euler-Lagrange equation for the optimization (4); here, . It follows that

 E(un+1)≤E(un+1)+12k||un+1−un||2≤E(un)+12k||un−un||2=E(un) (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 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:

1. Complete generality. There is no assumption (e.g. convexity) on the energy in (1) beyond sufficient differentiability.

2. Unconditional energy stability.

3. High (at least up to third) order accuracy.

4. 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.

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:

 un+1=UM=argminuE(u)+M−1∑i=0γM,i2k||u−Ui||2 (6)

where the intermediate stages , for , are given by

 Um=argminuE(u)+m−1∑i=0γm,i2k||u−Ui||2. (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 two-stage special case of scheme (6) & (7):

 U1=argminuE(u)+γ1,02k||u−un||2 (8) un+1=argminuE(u)+γ2,02k||u−un||2+γ2,12k||u−U1||2 (9)

and impose the conditions

 γ1,0−γ22,0γ2,0+γ2,1≥0 % and γ2,0+γ2,1≥0 (10)

on the parameters. Set . Note that (9) is equivalent to

 un+1=argminuE(u)+γ2,0+γ2,12k∥∥u−(θun+(1−θ)U1)∥∥2. (11)

We have

 E(un+1) ≤E(un+1)+γ2,0+γ2,12k∥∥un+1−(θun+(1−θ)U1)∥∥2 (by (10)) ≤E(U1)+γ2,0+γ2,12k∥∥U1−(θun+(1−θ)U1)∥∥2 (by (11)) =E(U1)+γ22,0(γ2,1+γ2,0)2k∥U1−un∥2 ≤E(U1)+γ1,02k∥U1−un∥2 (by (10)) ≤E(un) (by (8))

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

Define the following auxiliary quantities in terms of the coefficients of scheme (6) & (7):

 ~γm,i=γm,i−M∑j=m+1~γj,i~Sj,m~Sj,j (12) ~Sj,m=m−1∑i=0~γj,i (13)

If for , then scheme (6) & (7) satisfies the energy stability condition (2): For every we have .

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

Let the auxiliary quantities , and be defined as in Theorem 1. We have

 argminE(u)+m−1∑i=0γm,i2k||u−Ui||2 = argminE(u)+12kM∑j=m~S2j,m~Sj,j||u−m−1∑i=0~γj,i~Sj,mUi||2

###### Proof
 m−1∑i=0γm,i2k||u−Ui||2 = ||u||22km−1∑i=0γm,i−1k⟨u,m−1∑i=0γm,iUi⟩+terms that do not % depend on u.

Note that

 M∑j=m~S2j,m~Sj,j =~Sm,m+M∑j=m+1~S2j,m~Sj,j =m−1∑i=0~γm,i+M∑j=m+1~S2j,m~Sj,j =m−1∑i=0[γm,i−M∑j=m+1~γj,i~Sj,m~Sj,j]+M∑j=m+1~S2j,m~Sj,j =m−1∑i=0γm,i−M∑j=m+1[m−1∑i=0~γj,i]~Sj,m~Sj,j+M∑j=m+1~S2j,m~Sj,j =m−1∑i=0γm,i

Then

 12kM∑j=m~S2j,m~Sj,j||u−m−1∑i=0~γj,i~Sj,mUi||2 = ||u||22kM∑j=m~S2j,m~Sj,j−1k⟨u,m−1∑i=0M∑j=m~γj,i~Sj,m~Sj,jUi⟩+terms % that do not depend on u = ||u||22km−1∑i=0γm,i−1k⟨u,m−1∑i=0γm,iUi⟩+terms that do not % depend on u.

###### Lemma 2

Let the auxiliary quantities , be given in Theorem 1 and let for . Then

 E(Um)+12kM∑j=m~S2j,m~Sj,j||Um−m−1∑i=0~γj,i~Sj,mUi||2 ≤ E(Um−1)+12kM∑j=m−1~S2j,m−1~Sj,j||Um−1−m−2∑i=0~γj,i~Sj,m−1Ui||2

###### Proof
 E(Um)+12kM∑j=m~S2j,m~Sj,j||Um−m−1∑i=0~γj,i~Sj,mUi||2 ≤ E(Um−1)+12kM∑j=m~S2j,m~Sj,j||Um−1−m−1∑i=0~γj,i~Sj,mUi||2 (by Eq. 7 & Lemma 1) = E(Um−1)+12kM∑j=m~S2j,m~Sj,j||Um−1(1−~γj,m−1~Sj,m)−m−2∑i=0~γj,i~Sj,mUi||2 = E(Um−1)+12kM∑j=m~S2j,m~Sj,j||Um−1(~Sj,m−1~Sj,m)−m−2∑i=0~γj,i~Sj,mUi||2 (by Eq. 13) = E(Um−1)+12kM∑j=m~S2j,m−1~Sj,j||Um−1−m−2∑i=0~γj,i~Sj,m−1Ui||2 ≤ E(Um−1)+12kM∑j=m−1~S2j,m−1~Sj,j||Um−1−m−2∑i=0~γj,i~Sj,m−1Ui||2 (since ~Sm−1,m−1>0)

###### Proof

(of theorem)

 E(un+1) ≤E(un+1)+12k~S2M,M~SM,M||UM−M−1∑i=0~γM,i~SM,MUi||2 (since ~SM,M>0) ≤E(U1)+12kM∑j=1~S2j,1~Sj,j||U1−~γj,0~Sj,1U0||2 (by Lemma 2 repeatedly) ≤E(U0)+12kM∑j=1~S2j,1~Sj,j||U0−U0||2=E(un) (by Eq. 7 & Lemma 1)

## 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:

 [m−1∑i=0γm,i]Um+k∇HE(Um)=m−1∑i=0γm,iUi. (14)

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

 DnE(u)(v1,…,vn)=∂n∂s1⋯∂snE(u+s1v1+s2v2+⋯+snvn)∣∣∣s1=s2=⋯=sn=0

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 :

 u(k+t0)=u(t0)+kut(t0)+12k2utt(t0)+16k3uttt(t0)+% h.o.t.=U0−kDE(U0)+12k2D2E(U0)DE(U0)−16k3[D2E(U0)(D2E(U0)(DE(U0)))+D3E(U0)(DE(U0),DE(U0))]+h.o.t. (15)

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:

 Ui=U0−β1,ikDE(U0)+β2,ik2D2E(U0)DE(U0)−k3[β3,iD2E(U0)(D2E(U0)(DE(U0)))+β4,iD3E(U0)(DE(U0),DE(U0))]+h.o.t. (16)

where the coefficients obey the following recursive relation

 β1,0=β2,0=β3,0=β4,0=0β1,m=1Sm[1+m−1∑i=1γm,iβ1,i]β2,m=1Sm[β1,m+m−1∑i=1γm,iβ2,i]β3,m=1Sm[β2,m+m−1∑i=1γm,iβ3,i]β4,m=1Sm[β21,m2+m−1∑i=1γm,iβ4,i] (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: β1,M=1 β1,M=1 β1,M=1 β2,M=1/2 β2,M=1/2 (18) β3,M=1/6 β4,M=1/6
###### Proof

We will now show by induction that the aforementioned consistency formulas, Eq. 16 and Eq. 17, hold.

Stage one:

 γ1,0U1+kDE(U1)=γ1,0U0 (19)

We first will Taylor expand around in Eq. 19:

 γ1,0U1+kDE(U0)+kD2E(U0)(U1−U0)+12kD3E(U0)(U1−U0,U1−U0)+h.o.t.=γ1,0U0.

Now we plug in an ansatz for the expansion on around , , and solve for , and :

 γ1,0(A1k+A2k2+A3k3)+kDE(U0)+k2D2E(U0)(A1)+k3D2E(U0)(A2)+12k3D3E(U0)(A1,A1)+h.o.t.=0.

Matching terms of the same order we get

 A1 =−1γ1,0DE(U0) A2 =−1γ1,0D2E(U0)(A1)=1γ21,0D2E(U0)(DE(U0)) A3 =−1γ1,0D2E(U0)(A2)−121γ1,0D3E(U0)(A1,A1)

Noting that completes stage one.

Stage m:

 [m−1∑i=0γm,i]Um+kDE(Um)=m−1∑i=0γm,iUi. (20)

and assume Eq. 16 and Eq. 17 up to .

First we are going to solve for in Eq. 20:

 Um−U0= −kSmDE(Um)+1Smm−1∑i=0γm,iUi−U0. (21)

Now Taylor expand around in Eq. 21:

 Um−U0= −kSm[DE(U0)+D2E(U0)(Um−U0)+12D3E(U0)(Um−U0,Um−U0)] +1Smm−1∑i=0γm,iUi−U0+h.o% .t.

Plug in the ansatz for and equation Eq. 16 for , and retaining up to terms of third order, we have that

 (22)

Solving for , , by matching terms of the same order in Eq. 22, we arrive at:

completing the induction step.

Matching the consistency equations, Eq. 16 and Eq. 17, at with the one step error Eq. 15 gives the conditions on for various orders of accuracy Section 3, completing the proof.

In the next section, we give examples of ’s that satisfy the consistency equations (Section 3) as well as the hypothesis of Theorem 1 concurrently.

## 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:

 γ=⎛⎜⎝γ1,000γ2,0γ2,10γ3,0γ3,1γ3,2⎞⎟⎠=⎛⎜ ⎜⎝500−260−2314447⎞⎟ ⎟⎠. (23)

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 :

 γ=⎛⎜ ⎜ ⎜⎝9200−1164470−2875911483060944163148306⎞⎟ ⎟ ⎟⎠. (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:

 γ≈⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝67600000−15213670000−2120−1945874200095121−47669500315−436−43138242210−17675162.4577−11.55176.680111.9455⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠ (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 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.

###### Remark 1

We note that equations Eq. 6 and Eq. 7 can be rewritten using only one quadratic movement limiter term, so a black box implementation for backward Euler Eq. 3, or equivalently Eq. 4, is all that is needed for our method, and is called once per stage.

### 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 (