In this paper, we revisit two-step backward differentiation formula (BDF2) with variable time steps for solving the linear reaction-diffusion equation:

$ \begin{equation} \begin{aligned} & u_{t} = \Delta u +\kappa u +f(x, t), &&x\in\Omega, t\in(0, T], \\ &u(x, 0) = u_0(x), && x\in\bar{\Omega}, \\ &u(x, t) = 0, &&x\in\partial\Omega, t\in[0, T], \end{aligned} \end{equation} $

(1.1)

where the reaction coefficient $ \kappa \in \mathbb{R} $, and $ \Omega $ is a bounded domain.

Set the generally nonuniform time levels $ 0 = t_0<t_1<t_2<\cdots<t_N = T $ with the $ k $th time-step size $ \tau_k : = t_k-t_{k-1} $, the maximum step size $ \tau : = \max_{1\le k\le N}\tau_k $, and the adjacent time-step ratios

$ r_k = \frac{\tau_k}{\tau_{k-1}}, \quad \quad \quad \quad 2\leq k \leq N. $

The BDF1 and BDF2 formulas with variable time steps are respectively defined by

$ \begin{align*} \mathcal{D}_1 u^n& = \frac{1}{\tau_n} \nabla_\tau u^n, \quad \mathcal{D}_2 u^n = \frac{1+2r_n}{\tau_n(1+r_n)}\nabla_\tau u^n - \frac{r_n^2}{\tau_n(1+r_n)}\nabla_\tau u^{n-1}, \end{align*} $

where the difference operator $ \nabla_\tau u^n: = u^n-u^{n-1} $ for $ 1\leq n\leq N $.

By taking $ b_{0}^{(1)} : = 1/\tau_1 $, and for $ n>1 $

$ \begin{align} b_0^{(n)} = \frac{1+2r_n}{\tau_n(1+r_n)}, \quad b_1^{(n)} = - \frac{r_n^2}{\tau_n(1+r_n)}, \quad \text{and} \quad b_j^{(n)} = 0 \;\text{ (for } 2\leq j \leq n-1), \end{align} $

(1.2)

the BDF1 and BDF2 can be written as a unified discrete convolution form

$ \begin{equation} \mathcal{D}_2 u^n : = \sum\limits_{k = 1}^n b_{n-k}^{(n)} \nabla_\tau u^k, \quad \quad n\geq 1. \end{equation} $

(1.3)

For $ n = 1 $, we use BDF1 scheme to obtain the solution $ u^1 $, and for $ n>1 $, we use BDF2 scheme. Based on the unified notation (1.3), the BDF2 scheme with variable time steps is given as

$ \begin{align} \mathcal{D}_2 u^n = \Delta u^n +\kappa u^n +f^n, \quad \text{for} \quad 1\leq n\leq N. \end{align} $

(1.4)

The BDF2 with variable time steps is widely used to solve stiff or differential-algebraic problems [3, 4, 6, 15, 16, 17] as it has the nice property of the strong stability. One can refer to [1, 2, 3, 12, 16] for the details. While the practical use of BDF2 is well developed, the theoretical analysis seems to be difficult. Even so, many excellent mathematicians still make a big progress on the analysis of BDF2 scheme with variable time steps.

For the stability analysis of problem (1.1) with $ \kappa = 0 $, twenty years ago Becker [1] (one also refers to Thomée's classical book [16, Lemma 10.6]) presented the bound under the ratio condition $ 0< r_k \leq (2+\sqrt{13})/3\approx 1.868 $ that

$ \begin{equation} \left\|{u^n}\right\| \leq C \exp(C\Gamma_n)\big(\left\|{u_0}\right\| + \sum\limits_{j = 1}^n \tau_j \left\|{f^j}\right\| \big)\quad \text{for} \quad n\geq 1, \end{equation} $

(1.5)

where $ \Gamma_n : = \sum_{k = 2}^{n-2} \max\{0, r_k-r_{k+2}\} $ and $ \|\cdot\| $ denotes the $ L_2 $-norm. As pointed out in [16] and [2], the magnitudes of $ \Gamma_n $ can be zero, bounded [16, pp. 175] and unbounded [2, Remark 4.1] by selecting certain step-ratio sequence and vanishing step sizes. After that, Emmrich [3] improves the Becker's constrained condition to $ 0< r_k \leq 1.91 $, but still keeps the undesirable factor $ \exp(C\Gamma_n) $ in the $ L_2 $-norm stability. Recently, Chen et al. [2] present the energy stability for the Cahn-Hilliard equation under a new ratio condition $ 0\leq r_k \leq (3+\sqrt{17})/2 \approx 3.561 $, and then Liao and Zhang [11] propose the technique of the discrete orthogonal convolution (DOC) kernels and present the stability estimate for the linear problem (1.1) under the same ratio condition. The new ratio condition improves Grigorieff's stability condition $ 0\leq r_k < 1+\sqrt{2} $ given nearly forty years ago [5]. One also refers to [13] and [6, Section Ⅲ.5] a classical book by Hairer et al. In addition, Liao and Zhang [11] obtain the first-order convergence under $ 0\leq r_k \leq 3.561 $. If the second-order convergence is expected in [11], they need an extra restriction condition $ |\mathfrak{R}_p| \leq N_0 \ll N $ with the index set

$ \begin{equation} \mathfrak{R}_p = \left \{ k \;\left | \;1+\sqrt{2} \leq r_k \leq (3+\sqrt{17})/2 \right.\right\}. \end{equation} $

(1.6)

As pointed out in [11], the condition $ |\mathfrak{R}_p| \leq N_0 \ll N $ seems to be more theoretical rather than practical. For more practical applications, it is natural to ask if we can theoretically extend the restriction on adjacent time-step ratios and without any extra restriction condition like $ |\mathfrak{R}_p| \leq N_0 \ll N $, meanwhile, keep the sharp error estimate.

The aim of this paper is to extend $ 0\leq r_k \leq 3.561 $ to a new ratio condition $ 0< r_k <4.864{5} $ and prove that the second-order convergence of BDF2 scheme is sharp and robust. To do so, we use the concept of DOC kernels originally developed in [11], and also introduce the concept of discrete complementary convolution (DCC) kernels developed for solving fractional PDEs [7, 8, 9, 10]. Based on DOC and DCC kernels, we first present the corresponding energy ($ H^1 $-norm) stability estimate with a new adjacent time-step ratio condition

$ \begin{array}{l} {\bf{A1}}:0 < {r_k} \le {r_{\max }} = \frac{1}{6}\left( {\sqrt[3]{{1196 - 12\sqrt {177} }} + \sqrt[3]{{1196 + 12\sqrt {177} }}} \right) + \frac{4}{3} \approx 4.8645,\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\rm{for}}\;2 \le k \le N. \end{array}$

Here $ r_{\max} $ is the positive real solution of $ x^3 = (2x+1)^2 $, see the details in Lemma 2.1.

For the sharp and robust convergence, we further express the local truncated error by an error convolution structure (ECS) with the BDF2 kernel, see more details in Lemma 3.9. Using the definition of DOC kernels, the ECS can significantly circumvent the complex calculation of BDF2 and DOC kernels. Thus, we have the sharp and robust second-order convergence under the ratio condition A1 (i.e. $ 0< r_k \leq 4.86{45} $). The robustness means the error estimate only depends on the adjacent step ratio restriction A1, and does not suffer from other conditions on the mesh sizes, like the restricted condition $ |\mathfrak{R}_p| \leq N_0 \ll N $ in [11]. In this sense, the second-order convergence is robust for variable time step sizes. On the other hand, our analysis also shows that the first-order BDF1 for the first level solution $ u^1 $ is enough to have the sharp second-order convergence. Thus, our analysis removes the doubt of the choice of the first level solution $ u^1 $ computed by BDF1, which further improves the nice results in [11, 14].

The organization of the paper is given as follows. In Section 2, we present the semi-positive definition of BDF2 kernels under condition A1, and the properties of DOC and DCC kernels. The stability analysis and second-order convergence of the BDF2 scheme (1.3) are given in Section 3. Numerical examples are provided to demonstrate our theoretical analysis.

In this section, we first consider the positive semi-definiteness of BDF2 convolution kernels and the properties of DOC and DCC kernels, which are useful for the analysis of stability and convergence of BDF2 scheme in section 3.

We first consider the positive semi-definiteness of BDF2 convolution kernels $ b^{(n)}_{n-k} $. It has been proven in [11] that the BDF2 convolution kernels $ b^{(n)}_{n-k} $ are positive semi-definite under the ratio condition $ 0<r_k\leq 3.561 $. A natural question is whether the ratio condition can be relaxed. In this subsection, we will prove the positive semi-definiteness of BDF2 convolution kernels $ b^{(n)}_{n-k} $ under a new ratio condition A1 (i.e., $ 0<r_k\leq 4.8645 $). To this end, we first present the following lemma which plays a key role in the proof of the positive semi-definiteness of $ b^{(n)}_{n-k} $.

Lemma 2.1  There exist typical values $ \epsilon_*>0 $ and $ x_{\max}>0 $ such that

$ \begin{align} \mathfrak{F}(x, y, \epsilon_*): = \frac{2\epsilon_*+4\epsilon_* x -\epsilon_*^2 x^2}{(1+x)} - \frac{y}{(1+y)}\geq 0, \quad \quad 0< x, y \leq x_{\max}, \end{align} $

(2.7)

where $ \epsilon_* = 1/\sqrt{x_{\max}} $ and $ x_{\max} $ is the positive root of the equation $ </p> <p>x^3 = (1+2x)^2. $

Proof  We now present the details how to find $ \epsilon_* $ and $ x_{\max} $. Set $ y = x $ in (2.7) and consider the quadratic function

$ \begin{equation} \mathfrak{H}(x) = -\epsilon^2 x^2 +(4\epsilon -1) x +2\epsilon. \end{equation} $

(2.8)

The positive root of $ \mathfrak{H}(x) = 0 $ is given by

$ \begin{equation} x = \frac{4\epsilon -1 +\sqrt{8\epsilon^3 + 16\epsilon^2 - 8\epsilon + 1} }{2\epsilon^2}. \end{equation} $

(2.9)

Noting $ x $ is a function of $ \epsilon $. We can produce its maximum by searching $ \epsilon_* $ such that $ x^\prime = 0 $. To do so, we take the derivative of $ x $ with respect to $ \epsilon $ as

$ \begin{align*} x^\prime & = \frac{-2 \epsilon^3 - 8\epsilon^2+6\epsilon -1 -(2\epsilon-1) {\sqrt{8\epsilon^3 + 16\epsilon^2 - 8\epsilon + 1}}}{\epsilon^3 {\sqrt{8\epsilon^3 + 16\epsilon^2 - 8\epsilon + 1}}}. \end{align*} $

To find an $ \epsilon_* $ such that $ x' = 0 $, we only need to consider

$ \begin{align*} &-2 \epsilon_*^3 - 8\epsilon_*^2+6\epsilon_* -1 -(2\epsilon_*-1) {\sqrt{8\epsilon_*^3 + 16\epsilon_*^2 - 8\epsilon_* + 1}} = 0, \\ \Rightarrow \quad & (2\epsilon_*-1)^2(8\epsilon_*^3+16\epsilon_*^2-8\epsilon_*+1) = (2\epsilon_*^3+8\epsilon_*^2-6\epsilon_*+1)^2, \\ \Rightarrow \quad &(\epsilon_*^3+2\epsilon_*-1)\epsilon_*^3 = 0, \\ (\epsilon_*\neq 0)\Rightarrow \quad & \epsilon_*^3+2\epsilon_*-1 = 0. \end{align*} $

Thus, we have the positive root of $ \epsilon_*^3+2\epsilon_*-1 = 0 $ given as

$ \epsilon_* = \frac{\sqrt[3]{12}}{6} \big( \sqrt[3]{\sqrt{177} + 9} - \sqrt[3]{\sqrt{177} - 9} \;\big) \approx 0.4534. $

Set $ g(\epsilon) = {\sqrt{8\epsilon^3 + 16\epsilon^2 - 8\epsilon + 1}} $. From $ x^\prime (\epsilon_* ) = 0 $, we have

$ \begin{align*} &x^\prime (\epsilon_* ) = \frac{-2 \epsilon_*^3 - 8\epsilon_*^2+6\epsilon_* -1 -(2\epsilon_*-1)g(\epsilon_*)}{\epsilon_*^3g(\epsilon_*)} = 0, \\ \Rightarrow \quad & g(\epsilon_*) = \frac{-2 \epsilon_*^3 - 8\epsilon_*^2+6\epsilon_* -1 }{2\epsilon_*-1} = \frac{ - 8\epsilon_*^2+10\epsilon_* -3 }{2\epsilon_*-1} , \end{align*} $

where we have used $ \epsilon_*^3 = 1-2\epsilon_* $ in the last identity.

From (2.9), we have the maximum value $ x_{\max} $ at $ \epsilon_* $ as

$ \begin{align} x_{\max} & = \frac{4\epsilon_* -1 +g(\epsilon_*) }{2\epsilon_*^2} = \frac{4\epsilon_* -1 +\frac{ - 8\epsilon_*^2+10\epsilon_* -3 }{2\epsilon_*-1} }{2\epsilon_*^2} = \frac{1}{\epsilon_*^2} \\ & = \frac{1}{6} \left(\sqrt[3]{1196\!-\!12\sqrt{177}}+\sqrt[3]{1196\!+\!12\sqrt{177}}\right)+\frac43 \approx 4.864{5}. \end{align} $

(2.10)

From $ \epsilon_*^3+2\epsilon_*-1 = 0 $, we have $ \epsilon_*^2 (\epsilon_*^2+2)^2 = 1 $ and substitute $ \epsilon_* = \frac{1}{\sqrt{x_{\max}}} $ into the resulting. Then, the direct calculation shows that $ x_{\max} $ satisfies the equation

$ x_{\max}^3 = (1+2x_{\max})^2. $

For the given value $ \epsilon_* $ and $ r_{\max} $, we now prove (2.7) holds. Considering the function

$ \begin{align*} \mathfrak{G}(x) & = (2\epsilon_*+4\epsilon_* x -\epsilon_*^2 x^2) (1+x)^{-1} , \end{align*} $

its derivative $ \mathfrak{G}^\prime (x) $ is given by

$ \mathfrak{G}^\prime (x) = -\frac{\big(x-(\sqrt{\epsilon_*(\epsilon_*+2)}-\epsilon_*)\big)\big(x+(\sqrt{\epsilon_*(\epsilon_*+2)}-\epsilon_*)\big)}{\epsilon_*^{2}(1+x)^{2} }. $

When $ 0<x \leq \sqrt{\epsilon_*(\epsilon_*+2)}-\epsilon_* $, we have $ \mathfrak{G}^\prime (x) \geq 0 $. Hence, it holds that

$ \mathfrak{G}(x)\geq\mathfrak{G}(0) = 2\epsilon_* (\approx 0.9068) > \frac{x_{\max}}{1+x_{\max}} (\approx 0.8295). $

When $ \sqrt{\epsilon_*(\epsilon_*+2)}-\epsilon_* < x \leq x_{\max} $, we have $ \mathfrak{G}^\prime (x) \leq 0 $. Hence, it holds that

$ \mathfrak{G}(x) \geq \mathfrak{G}(x_{\max}) = \frac{x_{\max}}{1+x_{\max}}. $

Thus, we prove $ \mathfrak{F}(x, x_{\max}, \epsilon_*) \geq 0 $ for $ 0<x\leq x_{\max} $. Noting $ \mathfrak{F}(x, y, \epsilon) $ is a decreasing function with respect to $ y $, we have $ \mathfrak{F}(x, y, \epsilon_*) \geq\mathfrak{F}(x, x_{\max}, \epsilon_*)\geq 0 $ for $ 0<x, y\leq x_{\max} $. The proof is complete.

Lemma 2.2  Assume the time step ratio $ r_k $ satisfy A1 (i.e., $ 0< r_k \leq 4.86{45} $). For any real sequence $ \{w_k\}_{k = 1}^n $, it holds that

$ \begin{align} &2w_k \sum\limits_{j = 1}^k b^{(k)}_{k-j} w_j \geq \frac{r_{k+1}\sqrt{r_{\max}}}{(1+r_{k+1})}\frac{w_k^2}{\tau_k} - \frac{r_k\sqrt{r_{\max}}}{ (1+r_k)}\frac{w_{k-1}^2 }{\tau_{k-1}}, \quad k \geq 2, \end{align} $

(2.11)

$ \begin{align} &2\sum\limits_{k = 1}^nw_k \sum\limits_{j = 1}^k b^{(k)}_{k-j} w_j \geq 0, \quad \text{for } n \geq 1. \end{align} $

(2.12)

Proof  Noting $ 2ab \leq \epsilon a^2 + b^2/ \epsilon \; (\forall \epsilon >0) $ and $ b_1^{(k)}<0 $, we have for $ k \geq 2 $ that

$ \begin{align*} 2w_k \sum\limits_{j = 1}^k b^{(k)}_{k-j} w_j & = 2b^{(k)}_0 w_k^2 + 2b_1^{(k)} w_k w_{k-1} \geq (2b^{(k)}_0 + \epsilon b_1^{(k)}) w_k^2 + \frac{b_1^{(k)}}{\epsilon} w_{k-1}^2\\ & = \frac{2\epsilon+4\epsilon r_k -\epsilon^2 r_k^2}{1+r_k}\frac{w_k^2}{\epsilon\tau_k} - \frac{r_k}{ 1+r_k}\frac{w_{k-1}^2 }{\epsilon\tau_{k-1}} \\ & = \frac{2\epsilon+4\epsilon r_k -\epsilon^2 r_k^2}{1+r_k}\frac{w_k^2}{\epsilon\tau_k} {- \frac{r_{k+1}}{1+r_{k+1}}\frac{w_k^2}{\epsilon\tau_k} +\frac{r_{k+1}}{1+r_{k+1}}\frac{w_k^2}{\epsilon\tau_k} }- \frac{r_k}{1+r_k}\frac{w_{k-1}^2 }{\epsilon\tau_{k-1}} \\ & = {\underbrace{\left(\frac{2\epsilon+4\epsilon r_k -\epsilon^2 r_k^2}{1+r_k}- \frac{r_{k+1}}{1+r_{k+1}}\right)}_{{\large = :\mathfrak{F}(r_k, r_{k+1}, \epsilon)}}\frac{w_k^2}{\epsilon\tau_k} } +\frac{r_{k+1}}{1+r_{k+1}}\frac{w_k^2}{\epsilon\tau_k} - \frac{r_k}{ 1+r_k}\frac{w_{k-1}^2 }{\epsilon\tau_{k-1}} \\ & = \frac{r_{k+1}}{1+r_{k+1}}\frac{w_k^2}{\epsilon\tau_k} - \frac{r_k}{ 1+r_k}\frac{w_{k-1}^2 }{\epsilon\tau_{k-1}} + {\mathfrak{F}(r_k, r_{k+1}, \epsilon)}\frac{w_k^2}{\epsilon\tau_k}. \end{align*} $

Set $ \epsilon = \epsilon_* = 1/\sqrt{r_{\max}} $, it follows from Lemma 2.1 that

$ \begin{equation} \mathfrak{F}(r_k, r_{k+1}, {\epsilon_*})\geq 0, \quad \forall 0< r_k, r_{k+1} \leq r_{\max}. \end{equation} $

(2.13)

Thus, the inequality (2.11) holds true.

From the inequality (2.11), the direct calculation yields

$ \begin{align*} 2\sum\limits_{k = 1}^nw_k \sum\limits_{j = 1}^k b^{(k)}_{k-j} w_j & \geq \frac{2}{\tau_1} w_1^2 + \frac{r_{n+1}\sqrt{r_{\max}}}{(1+r_{n+1})}\frac{w_n^2}{\tau_n} - \frac{r_2\sqrt{r_{\max}}}{ (1+r_2)}\frac{w_{1}^2 }{\tau_{1}}\\ & = \frac{r_{n+1}\sqrt{r_{\max}}}{(1+r_{n+1})}\frac{w_n^2}{\tau_n} +\frac{2 + (2 -\sqrt{r_{\max}}) r_2}{ (1+r_2)}\frac{w_{1}^2 }{\tau_{1}}\\ & \geq \frac{r_{n+1}\sqrt{r_{\max}}}{(1+r_{n+1})}\frac{w_n^2}{\tau_n} +\frac{2 + (2 -\sqrt{r_{\max}}) r_{\max}}{ (1+r_{\max})}\frac{w_{1}^2 }{\tau_{1}}\\ &\geq \frac{r_{n+1}\sqrt{r_{\max}}}{(1+r_{n+1})}\frac{w_n^2}{\tau_n} + \frac{w_{1}^2}{(1+r_{\max})\tau_{1}} \geq 0, \quad n \geq 1, \end{align*} $

where the monotonicity of function $ l(x) = \frac{x}{1+x} $ and the fact $ 2 + (2 -\sqrt{r_{\max}}) r_{\max} = 1 $ are used. The proof is complete.

The DCC kernels $ p_{n-j}^{(n)} $ are introduced in analogy of $ \int_0^t v^\prime (s) ds = v(t) - v(0) $ such that

$ \begin{equation} \sum\limits_{j = 1}^np_{n-j}^{(n)}\mathcal{D}_2 u^j = \sum\limits_{j = 1}^np_{n-j}^{(n)} \sum\limits_{l = 1}^j b_{j-l}^{(j)} \nabla_\tau u^l = \sum\limits_{l = 1}^n\nabla_\tau u^l \sum\limits_{j = l}^np_{n-j}^{(n)}b_{j-l}^{(j)} = u^n - u^0, \quad \forall n\geq 1. \end{equation} $

(2.14)

From the identity (2.14) holds for all $ n\geq 1 $, we define the DCC kernels by

$ \begin{equation} \sum\limits_{j = k}^n p_{n-j}^{(n)} b_{j-k}^{(j)} \equiv1, \quad \forall 1\leq k \leq n, \; 1 \leq n \leq N. \end{equation} $

(2.15)

From (2.15), the DCC kernels $ p_{n-j}^{(n)} $ can be explicitly expressed by the BDF2 kernels $ b_{j-k}^{(j)} $, namely,

$ \begin{align} { p_{0}^{(n)}} = 1/b_0^{(n)}, \; p_{n-j}^{(n)} = \frac{1}{b_0^{(j)}} \sum\limits_{k = j+1}^n (b^{(n)}_{k-j-1} - b^{(n)}_{k-j}) p^{(n)}_{n-k} {\quad (1\leq j\leq n-1)}. \end{align} $

(2.16)

The discrete orthogonal convolution (DOC) kernels are given in [11] as

$ \begin{align} & \sum\limits_{j = k}^n \theta_{n-j}^{(n)} b_{j-k}^{(j)} = \delta_{nk}, \quad \text{for all } 1\leq k \leq n, \end{align} $

(2.17)

where $ \delta_{nk} $ represents the Kronecker delta symbol with $ \delta_{nk} = 1 $ if $ n = k $ and $ \delta_{nk} = 0 $ if $ n\neq k $. From the DOC kernels (2.17), we have

$ \begin{equation} \sum\limits_{j = 1}^n\theta_{n-j}^{(n)}\mathcal{D}_2 u^j = \sum\limits_{l = 1}^n\nabla_\tau u^l \sum\limits_{j = l}^n\theta_{n-j}^{(n)}b_{j-l}^{(j)} = u^n - u^{n-1}, \quad 1\leq n \leq N. \end{equation} $

(2.18)

The two kernels have the following intimate relationship.

Proposition 2.1  The DCC kernels defined by (2.15) and DOC kernels defined in (2.17) have the following relationships

$ \begin{align} p^{(n)}_{n-j}& = \sum\limits_{l = j}^n \theta^{(l)}_{l-j}, \quad \forall 1\leq j\leq n, \end{align} $

(2.19)

$ \begin{align} \theta^{(n)}_{n-j} & = p^{(n)}_{n-j}-p^{(n-1)}_{n-1-j}, \quad \forall 1\leq j\leq n, \end{align} $

(2.20)

where $ p^{(n)}_{-1}: = 0 \; (\forall n\geq 0) $ are defined.

Proof  Set $ q^{(n)}_{n-j} = \sum_{l = j}^n\theta^{(l)}_{l-j}\; (\forall 1\leq j\leq n). $ Then from the definition (2.17), we have

$ \sum\limits_{j = k}^n q^{(n)}_{n-j} b^{(j)}_{j-k} = \sum\limits_{j = k}^n \sum\limits_{l = j}^n \theta^{(l)}_{l-j}b^{(j)}_{j-k} = \sum\limits_{l = k}^n \sum\limits_{j = k}^l \theta^{(l)}_{l-j}b^{(j)}_{j-k} = \sum\limits_{l = k}^n \delta_{lk} = 1. $

Hence, $ q^{(n)}_{n-j} = \sum_{l = j}^n\theta^{(l)}_{l-j} \;( 1\leq j\leq n) $ are solutions to (2.15). Noting the DCC kernels uniquely exist due to the explicit expression (2.16). Thus, we have $ p^{(n)}_{n-j} = q^{(n)}_{n-j} = \sum_{l = j}^n\theta^{(l)}_{l-j} $. The equality (2.20) can be directly yielded by (2.19) and the proof is complete.

We now rewrite the definitions of DCC kernel (2.15) and DOC kernel (2.17) by

$ \begin{align} &{\bf B}{\bf P} = {\bf I}, \quad\quad \quad {\bf B}{\bf \Theta} = {\bf I_0}, \end{align} $

(2.21)

where we denote

$ \begin{array}{l} {\bf B} = \left[ \begin{array}{ccccccc} b_{0}^{(n)} & 0&0& \cdots & 0& 0& 0\\ b_{1}^{(n)} & b_{0}^{(n-1)} &0& \cdots & 0& 0& 0\\ 0 & b_{1}^{(n-1)} &b_{0}^{(n-2)}& \cdots & 0& 0& 0\\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\ 0&0& 0& \cdots & b_{1}^{(3)} & b_{0}^{(2)}& 0\\ 0 &0& 0& \cdots & 0 & b_{1}^{(2)} & b_{0}^{(1)} \end{array} \right], \; \end{array} \begin{array}{l} {\bf P} = \left[ \begin{array}{c} p_0^{(n)} \\ p_1^{(n)} \\ p_2^{(n)} \\ \vdots \\ \vdots\\ {p_{n-2}^{(n)}} \\ p^{(n)}_{n-1} \end{array} \right], \; {\bf \Theta} = \left[ \begin{array}{c} \theta_0^{(n)} \\ \theta_1^{(n)} \\ \theta_2^{(n)} \\ \vdots \\ \vdots\\ {\theta_{n-2}^{(n)}} \\ \theta^{(n)}_{n-1} \end{array} \right], \; \end{array} $

and $ {\bf I} = [1, 1, 1, \cdots, 1, 1]^{{T}} $ and $ {\bf I}_0 = [1, 0, 0, \cdots, 0, 0]^{{T}} $. It is easy to verify that each component of $ {\bf P} $ and $ {\bf \Theta} $ is positive by using mathematical induction.

The positive semi-definitiveness of the DOC kernel $ \theta^{(n)}_{n-k} $ can be derived in [11] by the positive semi-definitiveness of $ b^{(n)}_{n-k} $.

Lemma 2.3 ([11])  If the BDF2 kernels $ b_{n-k}^{(n)} $ defined in (1.2) are positive semi-definite, then the DOC kernels $ \theta^{(n)}_{n-k} $ defined in (2.17) are also positive semi-definite. This is, it holds for any real sequence $ \{\omega_j\}_{j = 1}^n $ that

$ \sum\limits_{k = 1}^n\omega_k\sum\limits_{j = 1}^k\theta_{k-j}^{(k)}\omega_j \geq 0, \quad \forall n\geq 1. $

Lemma 2.4 ([11])  The DOC kernels $ \theta_{n-j}^{(n)} $ have the following properties:

$ \begin{align} \sum\limits_{j = 1}^{n} \theta_{n-j}^{(n)} & = \tau_n, \quad \quad \text{for } n\geq 1, \end{align} $

(2.22)

$ \begin{align} \sum\limits_{k = 1}^{n} \sum\limits_{j = 1}^{k}\theta_{k-j}^{(k)} & = t_n , \quad \quad \text{for } n\geq 1. \end{align} $

(2.23)

Lemma 2.5 ([11])  The DOC kernel $ \theta_{n-j}^{(n)} $ can be explicitly represented by

$ \begin{equation} \begin{aligned} \theta_{n-k}^{(n)} & = \left \{\begin{array}{l} \frac{\tau_n (1+r_k)}{1+2r_k} \prod\limits_{i = k+1}^n \frac{r_i}{1+2r_i}, \quad\quad \text{ for } 2\leq k\leq n. \\ \tau_n \prod\limits_{i = k+1}^n \frac{r_i}{1+2r_i}, \quad\quad \text{ for } k = 1. \end{array} \right. \end{aligned} \end{equation} $

(2.24)

We point out the results in Lemma 2.5 play an important role for the following convergence analysis and the bound of DCC kernels. We now consider the properties of DCC kernels.

Proposition 2.2  Let $ \tau $ be the maximum time-step size and the time-step ratios satisfy $ 0<r_k\leq r_{*} $, where $ r_{*} $ is any given positive constant. The DCC kernels $ p^{(n)}_{n-k} $ defined in (2.15) satisfy

$ \begin{align} &p^{(n)}_{n-j} = {\frac{1+r_j}{1+2r_j}}\sum\limits_{k = j}^n\tau_k\prod\limits_{i = j+1}^k \frac{r_i}{1+2r_i}, \quad 2\leq j\leq n, \end{align} $

(2.25)

$ \begin{align} &p^{(n)}_{n-1} = \sum\limits_{k = 1}^n\tau_k\prod\limits_{i = 2}^k \frac{r_i}{1+2r_i}, \end{align} $

(2.26)

$ \begin{align} &\sum\limits_{j = 1}^n p^{(n)}_{n-j} = t_n, \end{align} $

(2.27)

$ \begin{align} &p^{(n)}_{n-j}\leq \sum\limits_{k = j}^n\tau_k\left(\frac{r_{*}}{1+2r_{*}}\right)^{k-j} \leq \sum\limits_{k = j}^n\frac{\tau_k}{2^{k-j}}\leq 2\tau, \end{align} $

(2.28)

where $ \prod_{i = j+1}^k = 1 $ for $ j\geq k $ is defined.

Proof  It follows from (2.19) in Proposition 2.1 and Lemma 2.5 that, for $ 2\leq j\leq n $,

$ p^{(n)}_{n-j} = \sum\limits_{k = j}^n \theta_{k-j}^{(k)} = {\frac{1+r_j}{1+2r_j}}\sum\limits_{k = j}^n\tau_k\prod\limits_{i = j+1}^k \frac{r_i}{1+2r_i}\leq \sum\limits_{k = j}^n\tau_k \left(\frac{r_{*}}{1+2r_{*}}\right)^{k-j}, $

and for $ j = 1 $,

$ p^{(n)}_{n-1} = \sum\limits_{k = j}^n \theta_{k-1}^{(k)} = \sum\limits_{k = 1}^n\tau_k\prod\limits_{i = 2}^k \frac{r_i}{1+2r_i}\leq \sum\limits_{k = j}^n\tau_k \left(\frac{r_{*}}{1+2r_{*}}\right)^{k-j}, $

where the monotonicity of function $ h(x) = \frac{x}{1+2x} $ is used. The application of $ \frac{r_{*}}{1+2r_{*}} \leq \frac12 $ for any $ r_*\geq 0 $ yields the last inequality in (2.28). The equality (2.27) can be derived directly by Proposition 2.1 and Lemma 2.4 since

$ \sum\limits_{j = 1}^n p^{(n)}_{n-j} = \sum\limits_{j = 1}^n \sum\limits_{k = j}^n \theta^{(k)}_{k-j} = \sum\limits_{k = 1}^n\sum\limits_{j = 1}^k\theta^{(k)}_{k-j} = \sum\limits_{k = 1}^n \tau_k = t_n. $

The proof is complete.

It is known that problem (1.1) with $ \kappa \leq 0 $ has the property of energy dissipation. We now present the corresponding energy stability for BDF2 scheme (1.4). To the end, we define a (modified) discrete energy $ E^k $ by

$ \begin{equation} E^k: = \frac{r_{k+1}\sqrt{r_{\max}}\tau_k}{1+r_{k+1}}\|\partial_\tau u^k\|^2+|u^k|_1^2-\kappa \| u^k\|^2, \quad \text{for}\ \kappa\leq 0 \ \text{and} \ k\geq 1, \end{equation} $

(3.29)

where the initial energy $ E^0: = |u^0|_1^2-\kappa\left\|{ u^0}\right\|^2 $ and $ \partial_\tau u^k = \nabla_\tau u^k /\tau_k $. Here we remark that our discrete energy $ E^k $ defined by (3.29) differs from the one in [11] due to the different modified formula, i.e., the first term in (3.29).

Here and below $ \langle{\cdot, \cdot}\rangle $ and $ \left\|{\cdot}\right\| $ represent the inner product and norm in $ L^2(\Omega) $ space.

Theorem 3.6  Assume the condition $ {\bf{A{1}}} $ holds and $ \kappa\leq 0 $, then the discrete solution $ u^n $ to the BDF2 scheme (1.4) with variable time steps satisfies

$ \begin{equation} \partial_\tau E^k\leq 2\langle f^k, \partial_\tau u^k\rangle, \quad \forall k\geq 1. \end{equation} $

(3.30)

Furthermore, the discrete energy has the following estimate

$ \begin{equation} \sqrt{E^n}\leq \sqrt{E^0}+ 4C_\Omega( \sum\limits_{k = 1}^{n}\|\nabla_\tau f^{k}\|+\|f^0\|), \quad \forall n\geq 1. \end{equation} $

(3.31)

Proof  For $ k\geq 2 $, the weak form of (1.4) is given as

$ \begin{align} \langle\mathcal{D}_2 u^k, v\rangle = -\langle\nabla u^k, \nabla v\rangle +\kappa \langle u^k, v\rangle +\langle f^k, v\rangle, \quad \forall v\in H^1_0(\Omega), \; 2\leq k\leq N. \end{align} $

(3.32)

Setting $ \upsilon = 2\nabla_\tau u^k $ in the weak form (3.32), we have

$ \begin{equation} 2\langle \mathcal{D}_2 u^k, \nabla_\tau u^k\rangle+2\langle\nabla u^k, \nabla_\tau \nabla u^k\rangle+2\langle -\kappa u^k, \nabla_\tau u^k\rangle = 2\langle f^k, \nabla_\tau u^k\rangle, \quad k\geq 2. \end{equation} $

(3.33)

It follows from Lemma 2.2 that

$ \begin{align} 2\langle \mathcal{D}_2 u^k, \nabla_\tau u^k\rangle\geq \frac{r_{k+1}\sqrt{r_{\max}}\tau_k}{1+r_{k+1}}\|\partial_\tau u^k\|^2-\frac{r_{k}\sqrt{r_{\max}}\tau_{k-1}}{1+r_{k}}\|\partial_\tau u^{k-1}\|^2. \end{align} $

(3.34)

Applying the inequality $ 2ab\leq a^2+b^2 $, one has

$ \begin{equation} \begin{aligned} & 2\langle {\nabla} u^k, \nabla_\tau \nabla u^k\rangle = 2\|\nabla u^k\|^2-2\|\nabla u^k\|\|\nabla u^{k-1}\| \geq \|\nabla u^k\|^2-\|\nabla u^{k-1}\|^2, \\ &2\langle -\kappa u^k, \nabla_\tau u^k\rangle = -\kappa(2\|u^k\|^2-2\| u^k\|\| u^{k-1}\|) \geq -\kappa(\|u^k\|^2-\| u^{k-1}\|^2). \end{aligned} \end{equation} $

(3.35)

Inserting (3.34) and (3.35) into (3.33) and using the definition (3.29), we arrive at

$ \nabla_\tau E^k\leq 2\langle f^k, \nabla_\tau u^k\rangle, \quad \text{for} \ k\geq 2. $

We now consider $ k = 1 $. For $ k = 1 $, the direct calculation produces

$ \begin{equation} 2\langle \mathcal{D}_2 u^1, \nabla_\tau u^1\rangle = 2\tau_1\|\partial_\tau u^1\|^2{\geq \frac{r_{\max}^{3/2}}{1+r_{\max}}\tau_1\|\partial_\tau u^1\|^2\geq \frac{r_2\sqrt{r_{\max}}}{1+r_2}\tau_1\|\partial_\tau u^1\|^2.} \end{equation} $

(3.36)

Noting the inequalities (3.35) also hold for $ k = 1 $, together with (3.36), we have

$ {\frac{r_2\sqrt{r_{\max}}}{1+r_2}\tau_1\|\partial_\tau u^1\|^2}+\|\nabla u^1\|^2-\kappa\|u^1\|^2\leq \|\nabla u^{0}\|^2-\kappa\| u^{0}\|^2+2\langle f^1, \nabla_\tau u^1\rangle, $

which implies

$ \nabla_\tau E^1\leq 2\langle f^1, \nabla_\tau u^1\rangle. $

Thus, we prove the inequality (3.30).

Taking summation from 1 to $ n $ for (3.30), we have

$ \begin{align} E^n &\leq E^0+2\sum\limits_{k = 1}^n \langle f^k, \nabla_\tau u^k\rangle \\ & = E^0+2\langle f^n, u^n\rangle-2\sum\limits_{k = 2}^{n}\langle \nabla_\tau f^{k}, u^{k-1}\rangle-2\langle f^{1}, u^0\rangle\\ & \leq E^0+ 2\|f^n\|\|u^n\|+2\sum\limits_{k = 2}^{n}\|u^{k-1}\|\|\nabla_\tau f^{k}\|+2\|f^1\|\|u^0\|, \end{align} $

(3.37)

where the Cauchy-Schwartz inequality is used to the last term in (3.37). On the other hand, noting the Poincaré inequality produces $ \|u^n\|\leq C_\Omega |u^n|_1 $, we have $ \|u^n\|\leq C_\Omega \sqrt{E^n} $. Thus, from (3.37), we arrive at

$ E^n\leq E^0+ 2C_\Omega( \|f^n\|\sqrt{E^n}+\sum\limits_{k = 2}^{n}\sqrt{E^{k-1}}\|\nabla_\tau f^{k}\|+\|f^1\|\sqrt{E^0}). $

Choose an integer $ n_0 $ $ (0\leq n_0\leq n) $ to satisfy $ E^{n_0} = \max_{0\leq k\leq n} E^k $. Then

$ \begin{aligned}E^{n_0}&\leq \sqrt{E^0}\sqrt{E^{n_0}}+ 2C_\Omega\sqrt{E^{n_0}}( \|f^{n_0}\|+\sum\limits_{k = 2}^{n_0}\|\nabla_\tau f^{k}\|+\|f^1\|)\\ & \leq \sqrt{E^0}\sqrt{E^{n_0}}+ 4C_\Omega\sqrt{E^{n_0}}( \sum\limits_{k = 2}^{n_0}\|\nabla_\tau f^{k}\|+\|f^1\|), \end{aligned} $

where the last inequality has used the fact that $ f^{n_0} = f^1+\sum_{k = 2}^{n_0} \nabla_\tau f^k $. Noting that $ \|f^1\|\leq \|\nabla_\tau f^1\|+\|f^0\| $, we have

$ \begin{aligned} \sqrt{E^n}&\leq \sqrt{E^{n_0}}\leq \sqrt{E^0}+ 4C_\Omega( \sum\limits_{k = 1}^{n_0}\|\nabla_\tau f^{k}\|+\|f^0\|)\\ &\leq \sqrt{E^0}+ 4C_\Omega( \sum\limits_{k = 1}^{n}\|\nabla_\tau f^{k}\|+\|f^0\|). \end{aligned} $

The proof is complete.

We first introduce a discrete Grönwall inequality.

Lemma 3.7  Assume $ \lambda>0 $ and the sequences $ \{v_j\}_{j = 1}^N $ and $ {\{\eta_j\}_{j = 0}^N} $ are nonnegative. If

$ v_n \leq \lambda \sum\limits_{j = 1}^{n-1} \tau_j v_j + \sum\limits_{j = 0}^n \eta_j, \quad \text{for} \quad 1\leq n \leq N, $

then it holds

$ v_n \leq \exp\big(\lambda t_{n-1}\big) \sum\limits_{j = 0}^n \eta^j, \quad \text{for} \quad 1\leq n \leq N. $

The standard induction hypothesis can give the proof of lemma 3.7, which is omitted here.

Theorem 3.8  If the condition A1 holds, the solution $ u^n $ of BDF2 scheme (1.4) is unconditionally stable in the $ L^2 $-norm. If $ \kappa>0 $ and the maximum time-step size $ \tau \leq \frac{1}{4\kappa} $, it holds

$ \begin{align} \left\|{u^n}\right\| &{\leq 2\exp\big(4\kappa t_{n-1}\big) \Big(\left\|{u^0}\right\| + 2\sum\limits_{k = 1}^{n}\|\sum\limits_{j = 1}^k\theta_{k-j}^{(k)} f^j\|\Big)}\\ &\leq 2\exp\big(4\kappa t_{n-1}\big) \Big(\left\|{u^0}\right\| + 2\sum\limits_{j = 1}^{n} p_{n-j}^{(n)} \left\|{f^j}\right\|\Big), \quad \text{for } 1\leq n\leq N. \end{align} $

(3.38)

If $ \kappa\leq 0 $, it holds

$ \begin{align} \left\|{u^n}\right\| &\leq \left\|{u^0}\right\|+2\sum\limits_{j = 1}^n\|\sum\limits_{j = 1}^k\theta_{k-j}^{(k)} f^j\| \leq \left\|{u^0}\right\| + 2\sum\limits_{j = 1}^{n}p_{n-j}^{(n)} \left\|{f^j}\right\|, \quad \text{for } 1\leq n\leq N. \end{align} $

(3.39)

Proof  Applying the property (2.18) of DOC kernels to scheme (1.4), we have

$ \begin{equation} \nabla_\tau u^k = \sum\limits_{j = 1}^k\theta_{k-j}^{(k)} (\Delta u^j +\kappa u^j) + \sum\limits_{j = 1}^k\theta_{k-j}^{(k)} f^j, \quad \text{for} \quad k\geq 1. \end{equation} $

(3.40)

Noting the positive semi-definiteness of the DOC kernels in Lemma 2.3, we have

$ \begin{equation} \sum\limits_{k = 1}^{n} \sum\limits_{j = 1}^{k} \langle u^k, \theta_{k-j}^{(k)} \Delta u^j\rangle = -\sum\limits_{k = 1}^{n} \sum\limits_{j = 1}^{k} \langle \nabla u^k, \theta_{k-j}^{(k)} \nabla u^j \rangle\leq 0. \end{equation} $

(3.41)

Taking the inner product with $ u^k $ on both sides of (3.40), summing the resulting from $ 1 $ to $ n $ and using (3.41), we have

$ \begin{align} \sum\limits_{k = 1}^n\langle u^k, \nabla_\tau u^k \rangle &\leq \sum\limits_{k = 1}^n\langle u^k, \sum\limits_{j = 1}^k\theta_{k-j}^{(k)} \kappa u^j\rangle + \sum\limits_{k = 1}^n\langle u^k, \sum\limits_{j = 1}^k\theta_{k-j}^{(k)} f^j \rangle, \quad \text{for } 1\leq n\leq N . \end{align} $

(3.42)

If $ \kappa\leq 0 $ in (3.42), we use $ \|u^k\|^2 - \|u^{k-1}\|^2\leq 2 \langle u^k, \nabla_\tau u^k \rangle $, the Cauchy-Schwarz inequality and Lemma 2.3 to have

$ \begin{align} \|u^n\|^2&\leq \|u^0\|^2+2\sum\limits_{k = 1}^n\|u^k\|{\| \sum\limits_{j = 1}^k\theta_{k-j}^{(k)}f^j\|}, \quad \text{for } 1\leq n\leq N. \end{align} $

(3.43)

Selecting an integer $ n_0 $ $ (0 \leq n_0 \leq n) $ such that $ \left\|{u^{n_0}}\right\| = \max_{0\leq k \leq n} \left\|{u^k}\right\| $. From (3.43), we have

$ \begin{align} \left\|{u^{n_0}}\right\|^2 &\leq \left\|{u^0}\right\| \left\|{u^{n_0}}\right\| + 2\|u^{n_0}\| \sum\limits_{k = 1}^{n_0} \|\sum\limits_{j = 1}^k\theta_{k-j}^{(k)} f^j\|. \end{align} $

(3.44)

Eliminating a $ \left\|{u^{n_0}}\right\| $ for both sides of (3.44) and noting $ n_0\leq n $, we have

$ \begin{align} \|u^n\|\leq \|u^{n_0}\|&\leq \left\|{u^0}\right\|+2\sum\limits_{k = 1}^{n_0}\|\sum\limits_{j = 1}^k\theta_{k-j}^{(k)} f^j\|\\ &\leq \left\|{u^0}\right\|+2\sum\limits_{k = 1}^n\|\sum\limits_{j = 1}^k\theta_{k-j}^{(k)} f^j\|\\ &\leq \left\|{u^0}\right\|+2\sum\limits_{j = 1}^n\| f^j\| \sum\limits_{k = j}^n\theta_{k-j}^{(k)} \\ &{ = \left\|{u^0}\right\| + 2\sum\limits_{j = 1}^n p_{n-j}^{(n)}\| f^j\|}, \end{align} $

where we have used the Cauchy-Schwarz inequality, exchanged the order of summation and used the property (2.19).

If $ \kappa>0 $ in (3.42), we apply the Cauchy-Schwarz inequality to have

$ \begin{align} \|u^n\|^2&\leq \|u^0\|^2+2\kappa\sum\limits_{k = 1}^n\|u^k\| \sum\limits_{j = 1}^k\theta_{k-j}^{(k)}\| u^j\|+2\sum\limits_{k = 1}^n\|u^k\| { \|\sum\limits_{j = 1}^k\theta_{k-j}^{(k)} f^j\|}, \quad \text{for } 1\leq n\leq N. \end{align} $

(3.45)

Similar to (3.44) by selecting $ n_0 $ $ (0 \leq n_0 \leq n) $ such that $ \left\|{u^{n_0}}\right\| = \max_{0\leq k \leq n} \left\|{u^k}\right\| $, one has

$ \begin{align} \left\|{u^{n_0}}\right\|^2 &\leq \left\|{u^0}\right\| \left\|{u^{n_0}}\right\| +2\kappa\|u^{n_0}\|\sum\limits_{k = 1}^{n_0}\| u^k\| \sum\limits_{j = 1}^k\theta_{k-j}^{(k)}+ 2\|u^{n_0}\| \sum\limits_{k = 1}^{n_0}{\| \sum\limits_{j = 1}^k\theta_{k-j}^{(k)}f^j\|}. \end{align} $

(3.46)

Eliminating a $ \left\|{u^{n_0}}\right\| $ for both sides of (3.46), we further have

$ \begin{align*} \left\|{u^{n}}\right\| &\leq \left\|{u^{n_0}}\right\| \leq \left\|{u^0}\right\| + 2\kappa \tau_k \sum\limits_{k = 1}^{n_0} \left\|{u^k}\right\| + 2\sum\limits_{k = 1}^{n_0}{\| \sum\limits_{j = 1}^k\theta_{k-j}^{(k)}f^j\|}\\ &\leq \left\|{u^0}\right\| + 2\kappa \tau_k \sum\limits_{k = 1}^{n} \left\|{u^k}\right\| + 2\sum\limits_{k = 1}^{n}{\|\sum\limits_{j = 1}^k\theta_{k-j}^{(k)} f^j\|}, \end{align*} $

where we use the facts $ n_0\leq n $ and $ \sum\limits_{j = 1}^k\theta_{k-j}^{(k)} = \tau_k $. Taking the maximum time-step size $ \tau\leq \frac{1}{4\kappa} $ in the above inequality, we finally arrive at

$ \begin{align*} \left\|{u^{n}}\right\| &\leq 2\left\|{u^0}\right\| + 4\kappa \tau_k \sum\limits_{k = 1}^{n-1} \left\|{u^k}\right\| + 4 \sum\limits_{k = 1}^{n}{\|\sum\limits_{j = 1}^k\theta_{k-j}^{(k)} f^j\|}\\ &\leq 2\left\|{u^0}\right\| + 4\kappa \tau_k \sum\limits_{k = 1}^{n-1} \left\|{u^k}\right\| + 4 \sum\limits_{j = 1}^{n}\| f^j\| \sum\limits_{k = j}^n\theta_{k-j}^{(k)}\\ &{ = 2\left\|{u^0}\right\| + 4\kappa \tau_k \sum\limits_{k = 1}^{n-1} \left\|{u^k}\right\| + 4 \sum\limits_{j = 1}^{n} p_{n-j}^{(n)}\| f^j\|.} \end{align*} $

The Grönwall inequality in Lemma 3.7 directly produces the result (3.38). The proof is complete.

Set $ e^n: = u(t_n, x)-u^n(x) \; (n\geq 1) $. From (1.4), the error function is governed by

$ \begin{align} \mathcal{D}_2 e^n = \Delta e^n +\kappa e^n +\eta^n, \quad \text{for} \quad 1\leq n\leq N, \end{align} $

(3.47)

where $ \eta^n: = \mathcal{D}_2 u(t_n)- u_{t} (t_n) {(1\leq n\leq N)} $ denotes the truncation error.

Lemma 3.9  Denote

$ \begin{equation} \begin{aligned} G^l &: = -\frac12\int_{t_{l-1}}^{t_l} (t-t_{l-1})^2 u_{ttt} \, \mathrm{d} t, \quad \quad 1\leq l \leq N, \\ R^j &: = -\frac1{2}b^{(j)}_1\tau_{j-1}\int_{t_{j-1}}^{t_j}(2(t-t_{j-1})+\tau_{j-1}) u_{ttt} \, \mathrm{d} t, \quad\quad 2 \leq j\leq {N}, \\ R^1&: = \frac1{2\tau_1}\int_0^{t_1} t^2 u_{ttt} \, \mathrm{d} t{-\frac1{\tau_1}\int_0^{t_1} tu_{tt} \, \mathrm{d} t}. \end{aligned} \end{equation} $

(3.48)

The truncation error $ \eta^j: = \mathcal{D}_2 u(t_j)- u_{t} (t_j) \;(1\leq j\leq N) $ can be expressed by the following form

$ \begin{equation} \eta^j = \sum\limits_{l = 1}^{j}b^{(j)}_{j-l}G^l+R^j, \quad 1\leq j\leq {N}. \end{equation} $

(3.49)

Moreover, we have the follwing estimate

$ \begin{align} &2\sum\limits_{k = 1}^{n}\|\sum\limits_{j = 1}^k\theta_{k-j}^{(k)} \eta^j\|\leq 2\sum\limits_{k = 1}^{n}\|G^k\|+2\sum\limits_{k = 1}^np_{n-k}^{(k)} \|R^k\|\\ \leq& 4\tau\int^{t_1}_{0} \| u_{tt} \| \, \mathrm{d} t+ \sum\limits_{k = 1}^{n} \tau_k^2\int_{t_{k-1}}^{t_k}\| u_{ttt} \| \, \mathrm{d} t + {2} t_n \max\limits_{{1}\leq k\leq n}\tau_k \int_{t_{k-1}}^{t_k} \| u_{ttt} \| \, \mathrm{d} t . \end{align} $

(3.50)

Proof  By using the Taylor's expansion (see [16]), one has

$ \begin{aligned} \eta^j & = \frac1{2}b^{(j)}_0 G^j+\frac1{2}b^{(j)}_1G^{j-1}-\frac1{2}b^{(j)}_1\tau_{j-1}\int_{t_{j-1}}^{t_j}(2(t-t_{j-1})+\tau_{j-1}) u_{ttt} \, \mathrm{d} t\\ & = \sum\limits_{l = 1}^jb^{(j)}_{j-l}G^l-\frac1{2}b^{(j)}_1\tau_{j-1}\int_{t_{j-1}}^{t_j}(2(t-t_{j-1})+\tau_{j-1}) u_{ttt} \, \mathrm{d} t \\ & = \sum\limits_{l = 1}^jb^{(j)}_{j-l}G^l+R^j, \quad 2\leq j\leq N, \end{aligned} $

where the property of BDF2 kernels (1.2) that $ b^{(j)}_k = 0 $ for $ k\geq 2 $ is used. For $ j = 1 $, using the Taylor's expansion again, one has

$ \begin{aligned} \eta^1 & = \frac{u(t_1)-u(0)}{\tau_1}-u_t(t_1) = {-b^{(1)}_0\int_0^{t_1} t u_{tt} \, \mathrm{d} t}\\ & = b^{(1)}_0G^1+\frac1{2\tau_1}\int_0^{t_1} t^2 u_{ttt} \, \mathrm{d} t{-\frac1{\tau_1}\int_0^{t_1} t u_{tt} \, \mathrm{d} t}\\ & = b^{(1)}_0G^1+R^1. \end{aligned} $

Hence, the equality (3.49) holds.

From (3.49), we now estimate

$ \begin{align} 2\sum\limits_{k = 1}^{n}\|\sum\limits_{j = 1}^k\theta_{k-j}^{(k)} \eta^j\|&\leq 2\sum\limits_{k = 1}^{n}\|G^k\|+2\sum\limits_{k = 1}^{n}\|\sum\limits_{j = 1}^k\theta_{k-j}^{(k)} R^j\|\\ &\leq 2\sum\limits_{k = 1}^{n}\|G^k\|+2\sum\limits_{k = 1}^np_{n-k}^{(k)} \|R^k\|, \end{align} $

(3.51)

where the last inequality uses (2.19). Note that $ G^l, R^j $ can be bounded by

$ \begin{align} \|G^l\|&\leq \frac{\tau_l^2}{2}\int_{t_{l-1}}^{t_l}\| u_{ttt} \| \, \mathrm{d} t, \quad 1\leq l\leq n \;, \\ \|R^j\|&\leq \frac{2r_j+1}{2r_j+2}\tau_j\int_{t_{j-1}}^{t_j}\| u_{ttt} \| \, \mathrm{d} t\\ &\leq {\tau_j} \int_{t_{j-1}}^{t_j}\| u_{ttt} \| \, \mathrm{d} t, \quad 2 \leq j\leq n. \end{align} $

(3.52)

Noting $ \sum_{k = 1}^np^{(n)}_{n-k} = t_n $ in (2.27) and $ p^{(n)}_{n-1}\leq 2\tau $ in (2.28), we have

$ \begin{align} \sum\limits_{k = 1}^{n} p_{n-k}^{(n)} \|R^k\|& = \sum\limits_{k = 2}^{n} p_{n-k}^{(n)} \|R^k\|+p^{(n)}_{n-1}\|R^1\|\\ &\leq \sum\limits_{k = 1}^{n} p_{n-k}^{(n)} {\tau_j} \int_{t_{j-1}}^{t_j}\| u_{ttt} \| \, \mathrm{d} t+p^{(n)}_{n-1}\int_0^{t_1} \|u_{tt}\| \, \mathrm{d} t\\ &\leq { t_n} \max\limits_{{1}\leq k\leq n}\tau_k \int_{t_{k-1}}^{t_k} \| u_{ttt} \| \, \mathrm{d} t+2\tau\int^{t_1}_{0} \| u_{tt} \| \, \mathrm{d} t. \end{align} $

(3.53)

Inserting (3.52) and (3.53) into (3.51), we have the inequality (3.54). The proof is complete.

Theorem 3.10  Let $ u(t, x) $ be the exact solution to problem (1.1). If the condition A1 holds, then the discrete solution $ u^n $ to BDF2 scheme (1.4) has the second-order convergence in the $ L^2 $-norm. If $ \kappa>0 $ and the maximum time-step size $ \tau<1/(4\kappa) $, it holds

$ \begin{align} \left\|{u(t_n)-u^n}\right\| \leq& 2\exp\big(4\kappa t_{n-1}\big) \left(\left\|{u(0)-u^0}\right\| +4\tau\int^{t_1}_{0} \| u_{tt} \| \, \mathrm{d} t+ \sum\limits_{k = 1}^{n} \tau_k^2\int_{t_{k-1}}^{t_k}\| u_{ttt} \| \, \mathrm{d} t\right. \\ &\left.+ {2} t_n \max\limits_{{1}\leq k\leq n}\tau_k \int_{t_{k-1}}^{t_k} \| u_{ttt} \| \, \mathrm{d} t\right), \quad \text{for } 1\leq n\leq N. \end{align} $

(3.54)

If $ \kappa\leq 0 $, it holds

$ \begin{align} \left\|{u(t_n)-u^n}\right\| \leq& \left\|{u(0)-u^0}\right\| +4\tau\int^{t_1}_{0} \| u_{tt} \| \, \mathrm{d} t + \sum\limits_{k = 1}^{n} \tau_k^2\int_{t_{k-1}}^{t_k}\| u_{ttt} \| \, \mathrm{d} t \\ &+ {2}t_n \max\limits_{{1}\leq k\leq n}\tau_k \int_{t_{k-1}}^{t_k} \| u_{ttt} \| \, \mathrm{d} t, \quad \text{for } 1\leq n\leq N. \end{align} $

(3.55)

Proof  From Theorem 3.8 that if $ \kappa>0 $ and the maximum time step $ \tau \leq \frac{1}{4\kappa} $, it holds

$ \begin{align} \left\|{u(t_n)-u^n}\right\| &{\leq 2\exp\big(4\kappa t_{n-1}\big) \Big(\left\|{u(0)-u^0}\right\| + 2\sum\limits_{k = 1}^{n}\|\sum\limits_{j = 1}^k\theta_{k-j}^{(k)} f^j\|\Big)}\quad \text{for } 1\leq n\leq N. \end{align} $

(3.56)

If $ \kappa\leq 0 $, it holds

$ \begin{align} \left\|{u(t_n)-u^n}\right\|\leq \left\|{u(0)-u^0}\right\|+2\sum\limits_{j = 1}^n\|\sum\limits_{j = 1}^k\theta_{k-j}^{(k)} f^j\|. \quad \text{for } 1\leq n\leq N. \end{align} $

(3.57)

The direct application of Lemma 3.9 to (3.56) and (3.57) produces the error estimates (3.54) and (3.55), respectively. The proof is complete.

Remark 1  For the error estimate of problem (1.1) with $ \kappa = 0 $, Becker [1] gives an estimate for $ 0<r_k<\frac{2+\sqrt{13}}{3}\approx1.868 $, which is improved to $ 0<r_k<1.91 $ in [3] later. By choosing different $ r_k $, the fact $ \Gamma_n $ in (1.5) can be bounded [16, pp. 175], unbounded [2, Remark 4.1] or zero. Recently, Liao and Zhang [11] give an improved estimate

$ \begin{align} \left\|{u(t_n)-u^n}\right\| \!\leq \! 2\exp\big(4\kappa t_{n-1}\big) \big(\left\|{u_0\!-\!u^0}\right\| \!+\!2t_n\int^{t_1}_{0} \| u_{tt} \| \, \mathrm{d} t \!+\! 3 t_n \!\max\limits_{1\leq k\leq n}\tau_k \int_{t_{k-1}}^{t_k} \| u_{ttt} \| \, \mathrm{d} t\big) \end{align} $

(3.58)

with $ 0\leq r_k \leq 3.561 $. One can see that the right-hand-side second term in (3.58) has the first-order convergence when $ t_n $ is large. If they expect to have the second-order convergence, they need an extra restriction condition $ |\mathfrak{R}_p| \leq N_0 \ll N $ with the index set defined by (1.6). A similar error estimate is given in [18] with $ 0<r_k<\sqrt{2}+1 $.

Our result in Theorem 3.10 shows the sharp second-order convergence under $ 0<r_k\leq 4.8645 $. And the second-order convergence is robust, which means the convergence order remains valid for any time step satisfying the ratio $ 0<r_k\leq r_{\max}\approx 4.864 $. As far as we know, it is a pioneer paper to clarify the robust and sharp second-order convergence under the new ratio $ 0<r_k\leq 4.864 $.

We now report two examples to investigate the convergence order of BDF2 scheme (1.4) with variable time-steps. In the simulations, we set the computational domain $ \Omega = (0, 2)^2 $, final time $ T = 1 $, the number of spatial mesh $ M $ chosen by $ M = N $. By taking

$ f = (\pi^2/2-\kappa-1)\exp(-t)\sin(\pi x/2)\sin(\pi y/2), $

we can construct an exact solution to problem (1.1) as a benchmark solution in the form of

$ u = \exp(-t)\sin(\pi x/2)\sin(\pi y/2). $

The time meshes are constructed by the random time-steps $ \tau_k = T\chi_k/C $, where $ C = \sum_{k = 1}^N \chi_k $ and $ \chi_k $ is randomly drawn from the uniform distribution on $ (0, 1) $. In each run, the discrete $ L^2 $-norm at the final time $ T = 1 $

$ e(N) = {h}\sqrt{\sum\limits_{1\leq i, j\leq M} (u(x_i, y_j, T)-u^N_h(x_i, y_j))^2} $

is recorded in Tables 1 and 2, in which we also list the maximum time-step $ \tau $ and maximum adjacent time-step ratio. The numerical rate of convergence is calculated by

$ \text{Order} = \log_2(e(N)/e(2N)). $

From the current data and more tests not listed here, we see that the BDF2 scheme is robustly stable and convergent in the second order, which is consistent with our theoretical analysis. Due to the time step is randomly chosen without any constrain condition, one can see the first-step BDF1 does not bring the loss of accuracy, which again implies the effectiveness of our analysis.

With the applications of DCC and DOC kernels, we present the stability and convergence analysis of BDF2 scheme with variable time-steps under condition A1. We extend the adjacent time step to a new ratio $ r_k: = \tau_k/\tau_{k-1} \leq r_{\max} = 4.86{45} $, and obtain the robust and sharp second-order convergence without the extra constrained condition on ratios in [11]. Our convergence results shows that the BDF1 scheme for first step is enough to have the globally optimal second-order convergence. This conclusion removes the doubt of the choice of the first level solution with first-order accuracy. Numerical results are provided to demonstrate the theoretical analysis.

The technique developed in this work can be extended to a family of multi-step schemes with variable time-steps for the stability and convergence analysis. The main challenge is how to explore the useful properties of DOC and DCC kernels, and multi-step schemes's kernels.

The authors would like to thank professor Tao Tang, professor Zhimin Zhang and Dr. Honglin Liao for their valuable discussions on this topic.