 # Stability and error estimates for the variable step-size BDF2 method for linear and semilinear parabolic equations

In this paper stability and error estimates for time discretizations of linear and semilinear parabolic equations by the two-step backward differentiation formula (BDF2) method with variable step-sizes are derived. An affirmative answer is provided to the question: whether the upper bound of step-size ratios for the l^∞(0,T;H)-stability of the BDF2 method for linear and semilinear parabolic equations is identical with the upper bound for the zero-stability. The l^∞(0,T;V)-stability of the variable step-size BDF2 method is also established under more relaxed condition on the ratios of consecutive step-sizes. Based on these stability results, error estimates in several different norms are derived. To utilize the BDF method the trapezoidal method and the backward Euler scheme are employed to compute the starting value. For the latter choice, order reduction phenomenon of the constant step-size BDF2 method is observed theoretically and numerically in several norms. Numerical results also illustrate the effectiveness of the proposed method for linear and semilinear parabolic equations.

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

In this paper we shall study stability and error estimates for time discretizations by the two-step backward differentiation formula (BDF2) with variable step-sizes for linear parabolic partial differential equations (PDEs)

 {u′(t)+Au(t)+Bu(t)=f(t),t∈J:=(0,T],u(0)=u0, (1)

and its semilinear extension, where : is a positive definite, self-adjoint, linear operator on a Hilbert space with domain dense in , the linear operator satisfies some structural assumptions; here the forcing term , and initial value .

Let the time interval for given , , be partitioned via . Let , be the time step-sizes which in general will be variable, and . We set

 rn=knkn−1,  n=2,3,…,N;kmax=maxn=2,3,…,Nkn;rmax=maxn=2,…,Nrn.

Assuming we are given starting approximations and , which is computed by the trapezoidal method or the backward Euler scheme, we discretize (1.1) in time by the variable step-size BDF2, i.e., we define nodal approximations to the values of the exact solution to (1.1) as follows:

 ¯∂2BUn+AUn+BUn=fn,n=2,3,…,N, (2)

where and

 ¯∂2BUn := (1+sn)¯∂1BUn−sn¯∂1BUn−1 = 1kn((1+sn)Un−(1+rn)Un−1+rnsnUn−2).

Here and are defined by

 ¯∂1BUn=¯∂Un:=Un−Un−1kn,and  sn:=rn1+rn=knkn+kn−1,

respectively. For an equidistant partition with , we have and the well-known formula

 ¯∂2BUn=1k(32Un−2Un−1+12Un−2).

The BDF2 method is one of the most popular time-stepping methods and many studies have been conducted on the stability and error estimates for it. Because of its good stability property (the scheme is -stable), the BDF2 method with constant step-size has been dealt with for various equations as, e.g., linear parabolic equations [24, 4], integro-differential equations , jump-diffusion model in finance , the Navier-Stokes equations [15, 20, 14]. When the solutions of time dependent differential equations have different time scales, i.e., solutions rapidly varying in some regions of time while slowly changing in other regions, variable step-sizes are often essential to obtain computationally efficient, accurate results. Owing to these prominent advantages, the variable step-size BDF2 method has been successfully applied to partial integro-differential equations  and Cahn-Hilliard equation  recently. An important result that the variable step-size BDF2 method is zero-stable if the step-size ratios are less than has been independently proved by several authors [26, 16, 10]. And the value of cannot be improved when dealing with arbitrary variable step-sizes (see, for example, [16, 6, 7]).

For the variable step-size BDF2 method applied to linear parabolic equations, Le Roux  derived stability and error bounds in the norm by using spectral techniques under the step-size conditions

 ckmax≤kj≤kmax,|rn−1|≤Ckn−11+|logkn−1|, (3)

with constants and . Palencia and García-Archilla  studied linear parabolic equations in a Banach space setting and obtained that the ratios should be close to such as, e.g., in (3) for the stability factor to be moderate. Grigorieff [17, 18] showed stability and optimal error estimates with smooth or non-smooth data under the assumption that the step-size ratios are less than in a Banach space setting. Becker improved the bound up to in Hilbert space  (see, also, ) by testing (2) with for a specified constant . Based on the same technique, i.e., testing (2) with , Emmrich  extended the results to semilinear parabolic problems and further improved the bound to using a more general identity for . Emmrich  also studied the stability and convergence of the variable step-size BDF2 method for nonlinear evolution problems governed by a monotone potential operator.

It is natural to ask what the upper bound of step-size ratios is and whether it is identical with the upper bound for the zero-stability. In this paper we will address this question and give an affirmative answer. We explore a new technique, which is very different from the one used by Becker , Thomée  and Emmrich [12, 13]. We first test (2) with to obtain -stability and -stability (their definitions will be introduced in Section 2). Then after we test (2) with , using -stability estimate, we obtain the usual stability result in the norm under the sharp zero-stability condition on the ratios of consecutive step-sizes. Following the approach of Chen et. al. , the and -stabilities of the variable step-size BDF2 method are also established under a more relaxed assumption that the step-size ratios are less than .

It is well known that the method (2) yields second order approximations to (in the norm) when the backward Euler method is used to compute the starting value , since it is applied only once. This choice for is quite popular in the multistep methods for computations of parabolic equations. However, the error bounds derived in this paper based on the obtained stability results suggest that this is not the best choice for the constant step-size BDF2 method, since it will cause the reduction of the convergence order in the and norms. This will be discussed in detail in Section 4.

The rest of this paper is organized as follows. We start in Section 2 by introducing the necessary notation and recalling a lemma which will be used in the following analysis. The stability of the method in several norms under the condition that the step-size ratios are less than (or ) is proved in Section 3. Error estimates in different norms are derived in Section 4. Since our error estimates will depend on the first step error , the error produced by the trapezoidal method or the backward Euler scheme will be analyzed in this section too. Section 5 will extend the analysis to the semilinear case

 {u′(t)+Au(t)=f(t,u),t∈J:=(0,T],u(0)=u0, (4)

with some assumptions on the nonlinear operator . Section 6 is devoted to numerical experiments, which confirm our theoretical results and illustrate the effectiveness of the proposed method for linear and semilinear parabolic equations. Section 7 contains a few concluding remarks.

## 2 Variable step-size BDF2 method for linear parabolic equations

Now we consider the variable two-step BDF method for solving (1.1). To do this, we first make some assumptions and introduce the necessary notation.

### 2.1 Linear parabolic equations

Let and denote the norms in and by and , , respectively. Let be the dual of (), and denote by the dual norm on , . We denote by the duality pairing between and . We define a bilinear form via . For the linear operator , we assume that

 |Bu|≤γ(t)∥u∥,∀u∈V,t∈J, (5)

with a smooth nonnegative function . Let . We may write the parabolic problem in variational form as

 (ut,v)+a(u,v)+(Bu,v)=(f,v),∀v∈V,  t∈J;u(0)=u0. (6)

Emmrich in  has shown that for given and , problem (6) admits a unique solution with .

Standard example. Let and be defined by

 Au:=−d∑i,j=1∂∂xi(αij∂u∂xj),Bu:=d∑j=1αj∂u∂xj+α0u,

respectively, where are sufficiently smooth functions in with being a bounded domain in with sufficiently smooth boundary . Let and be the usual Sobolev and Lebesgue space, respectively. Assume that is symmetric and uniformly positive definite. Then the operators is a positive definite, self-adjoint, linear operator, and satisfies the condition (5).

### 2.2 Variable step-size BDF2 method for linear parabolic equations

For the method (2) we need the starting values and . We set and perform an initial trapezoidal approximation to get

 ¯∂U1+AU12+BU12=f12 (7)

with . Note that with , the two-step BDF formally degenerates to a backward Euler step. It is also easy to see that

 ¯∂2BUn=kn¯∂2Un+1kn−1+kn(Un−Un−2)=kn¯∂2Un+sn¯∂Un+kn−1kn−1+kn¯∂Un−1=knsn¯∂2Un+¯∂Un, (8)

where

 ¯∂2vn:=¯∂¯∂vn=1kn[¯∂vn−¯∂vn−1].

With respect to the solvability of the time discrete problem, Emmrich has shown in  that for given and , the problem

 (¯∂2BUn,v)+a(Un,v)+(BUn,v)=(fn,v),∀v∈V (9)

admits a unique solution. For the obtained solution sequence , we define the , and norms as

 |U|l∞(tn1,tn2;H):=maxn1≤j≤n2|Uj|,∥U∥l∞(tn1,tn2;V):=maxn1≤j≤n2∥Uj∥, (10)

and

 ∥U∥l2(tn1,tn2;V):=(n2∑j=n1kj∥Uj∥2)1/2, (11)

respectively. It is well known that they are the discrete counterparts of the , and norms, respectively.

Remark. [The choice for ]. The starting value can be also obtained by the backward Euler method

 ¯∂U1+AU1+BU1=f1 (12)

with . It is well known that the constant step-size BDF2 method (2) with this initial approximation also yields second order approximations to in norm; see Corollary 4.6 in Section 4, or, [5, 24]. However, we find that order reduction will be caused in the and norms when the starting value is obtained by the backward Euler method (12). Here the norm of a solution sequence with is defined by

 ∥U∥l2(tn1,tn2;H,H):=(n2−1∑j=n1kjsj∣∣¯∂Uj∣∣2)1/2, (13)

which is the discrete counterpart of with nonuniform grid weights .

Because of the different choices for and , we pay special attention to and set .

In subsequent sections, by convention, we set and if . We will use the identity

 2(κa−εb)a=(2κ−ε)a2−εb2+ε(a−b)2,κ,ε≥0. (14)

We also need the following discrete Gronwall lemma proved in .

[Discrete Gronwall lemma ] Let , and with being monotonically increasing. Then

 an+bn≤n−1∑j=ϖλjaj+λan+gn,n=ϖ,ϖ+1,… (15)

implies for

 an+bn≤gn1−λn−1∏j=ϖ(1+λj1−λ)≤gn1−λexp(11−λn−1∑j=ϖλj).

## 3 Stability analysis

In this section we shall show stability of the variable step-size BDF2 method with applied to linear parabolic equations (1.1). As mentioned in Introduction, for the variable step-size BDF2 method applied to parabolic equations, the best known result is that it is stable in the norm when the step-size ratios are less than . To improve the bound to , the upper bound for the zero-stability, we first need the following stability results in the norm. Additionally, since the norm is an energy norm, from a physical point of view, stability is of utmost important.

[ and stability under ] Let with . If there exists a constant such that satisfies

 (4+2√2)γ2kmax≤c1<1, (16)

then we have, for any ,

 kn∣∣¯∂Un∣∣2+∥U∥2l2(t2,tn;H,H)+∥U∥2l∞(t2,tn;V)≤C1(n∑j=2kj∣∣fj∣∣2+k2s2∣∣¯∂U1∣∣2+∥U1∥2), (17)

where

 C1=4+2√21−c1exp(4+2√21−c1n∑j=2γ2(tj)kj)≤4+2√21−c1exp(4+2√21−c1γ2tn).

We first have from (2) and (8)

 knsn¯∂2Un+¯∂Un+AUn+BUn=fn. (18)

Taking in (18) the inner product with , using the relation (14), we obtain, for any ,

 knsn(∣∣¯∂Un∣∣2−∣∣¯∂Un−1∣∣2+∣∣¯∂Un−¯∂Un−1∣∣2)+2kn∣∣¯∂Un∣∣2+∥Un∥2−∥Un−1∥2+∥Un−Un−1∥2≤2γ(tn)kn∥Un∥∣∣¯∂Un∣∣+2kn|fn|∣∣¯∂Un∣∣≤1ϵγ2(tn)kn∥Un∥2+ϵkn∣∣¯∂Un∣∣2+1ϵkn|fn|2+ϵkn∣∣¯∂Un∣∣2. (19)

Summing up gives for ,

 [sn+2(1−ϵ)]kn∣∣¯∂Un∣∣2+n−1∑j=2(kjsj+2(1−ϵ)kj−kj+1sj+1)∣∣¯∂Uj∣∣2+∥Un∥2 (20) ≤ 1ϵn∑j=2γ2(tj)kj∥Uj∥2+1ϵn∑j=2kj∣∣fj∣∣2+k2s2∣∣¯∂U1∣∣2+∥U1∥2.

Now take such that

 sup0≤x≤R0x21+x≤2(1−ϵ).

Then we get

 (sn+1)kn∣∣¯∂Un∣∣2+n−1∑j=2kjsj∣∣¯∂Uj∣∣2+∥Un∥2≤(4+2√2)n∑j=2γ2(tj)kj∥Uj∥2+(4+2√2)n∑j=2kj∣∣fj∣∣2+k2s2∣∣¯∂U1∣∣2+∥U1∥2. (21)

An application of Lemma 2.2 leads to the desired result.

Recently, the stability of a linearly implicit stabilization BDF2 method with variable step-sizes has been established under the condition for the Cahn-Hilliard equation in . Following their approach, we can also improve the bound to for the stability of the variable step-size BDF2 method for the problem (1.1).

[ and stability under ] Let with , and let . If there exists a constant such that satisfies

 cRγ2kmax≤c1<1, (22)

then we have, for any ,

 R1+Rkn∣∣¯∂Un∣∣2+∥Un∥2≤C2(R1+Rk1∣∣¯∂U1∣∣2+∥U1∥2+cRn∑j=2kj∣∣fj∣∣2), (23)

and

 n∑j=2kj∥Uj∥2≤C2tn(k1∣∣¯∂U1∣∣2+∥U1∥2+cRn∑j=2kj∣∣fj∣∣2), (24)

where

 C2=cR1−c1exp(cR1−c1n∑j=2γ2(tj)kj)≤cR1−c1exp(cR1−c1γ2tn).

Taking in (18) the inner product with , we obtain, for any ,

 kn(2+4rn−r2n)1+rn|¯∂Un|2−kn1+rn|¯∂Un−1|2+rnsn|Un−2Un−1+Un−2|2kn+∥Un∥2−∥Un−1∥2+∥Un−Un−1∥2≤2γ(tn)kn∥Un∥∣∣¯∂Un∣∣+2kn|fn|∣∣¯∂Un∣∣≤1ϵγ2(tn)kn∥Un∥2+ϵkn∣∣¯∂Un∣∣2+1ϵkn|fn|2+ϵkn∣∣¯∂Un∣∣2. (25)

Ignoring some of the positive terms on the left-hand side, we have

 g(z,ϵ)kn∣∣¯∂Un∣∣2+∥Un∥2≤snkn−1∣∣¯∂Un−1∣∣2+∥Un−1∥2+1ϵγ2(tn)kn∥Un∥2+1ϵkn|fn|2, (26)

where . In the case , noting that , we can take such that

 g(rn,ϵ)>R1+R. (27)

In the case , since is a decreasing function, we take such that (27) holds. Thus in both cases, we have

 R1+Rkn∣∣¯∂Un∣∣2+∥Un∥2≤R1+Rkn−1∣∣¯∂Un−1∣∣2+∥Un−1∥2+1ϵγ2(tn)kn∥Un∥2+1ϵkn|fn|2≤R1+Rk1∣∣¯∂U1∣∣2+∥U1∥2+1ϵn∑j=2γ2(tj)kj∥Uj∥2+1ϵn∑j=2kj∣∣fj∣∣2, (28)

Apply Lemma 2.2 to (28) to obtain the desired inequality (23). The estimate (24) is a direct result of (23). This completes the proof.

Now we use the stability estimate (17) in Theorem 3.1 to show the stability of the variable step-size BDF2 method in the and norms.

[ and stability under ] Let with . If there exist constants and such that satisfies (16) and

 2c3γ2kmax≤c2<1, (29)

where , then the following estimate holds for :

 |U|2l∞(t2,tn;H)+∥U∥2l2(t2,tn;V)≤CfU, (30)

with

 fU=|U1|2+kmax∥U1∥2+(k21+kmaxk2)|¯∂U1|2+n∑j=2kj(∥fj∥2∗+kmax∣∣fj∣∣2).

Here, depends on , , , and , , with being defined by

 Φn:=n−2∑j=2[rj−rj+2]+,[x]+:=|x|+x2.

Taking in (2) the inner product with yields

 2kn1+rn(¯∂2BUn,Un)+2kn1+rna(Un,Un)+2kn1+rn(BUn,Un)=2kn1+rn(fn,Un),  n=2,3,…,N. (31)

By simple calculations, the first term on the left-hand side becomes

 2kn1+rn(¯∂2BUn,Un)=kn1+rn¯∂2B|Un|2+1+2rn(1+rn)2|Un−Un−1|2−s2n|Un−1−Un−2|2−2s2n(Un−Un−1,Un−1−Un−2)=kn1+rn¯∂2B|Un|2+1+2rn−r2n(1+rn)2|Un−Un−1|2−2s2n|Un−1−Un−2|2+s2n|Un−2Un−1+Un−2|2. (32)

Now

 2|(fn,Un