In this paper, we revisit the two-step backward differentiation formula (BDF2) with variable time-steps for solving the molecular beam epitaxial (MBE) model [6, 9] without slope selection

$ \begin{equation} \begin{aligned} & u_{t} + \varepsilon\Delta^2 u + \nabla \cdot \mathit{\boldsymbol{f}}(\nabla u) =0, && {{\bf x}} \in\Omega, t\in(0, T], \\ &u({\mathit{\boldsymbol{x}}}, 0)=u_0(\mathit{\boldsymbol{x}}), && \mathit{\boldsymbol{x}}\in\bar{\Omega}, \end{aligned} \end{equation} $

(1.1)

with the periodic boundary conditions. Here the periodic solution $ u = u(\mathit{\boldsymbol{x}}, t) $ represents the scaled height function of a thin film in a co-moving frame, the fourth-order term models surface diffusion with a surface diffusion constant $ \varepsilon >0 $ and the nonlinear force vector $ \mathit{\boldsymbol{f}}({\bf v}) :={\bf v}/(1+|{\bf v}|^2) $ models the well-known Ehrlich-Schwoebel effect.

The MBE model (1.1) has been widely applied in various fields such as physics, biology, ecology and chemistry [7, 22], and can be derived from the gradient flow with the following energy functional in the $ L^2(\Omega) $ inner product

$ \begin{align} E(u) = \int_\Omega \left(\frac{\varepsilon}{2}|\Delta u|^2-\frac12\ln(1+|\nabla u|^2)\right) \, \mathrm{d} \mathit{\boldsymbol{x}}. \end{align} $

(1.2)

The logarithmic term $ -\frac12\ln(1+|\nabla u|^2) $ in the energy functional (1.2) can be bounded by zero but unbounded blew, which means the logarithmic term has no relative minima. The well-posedness of problem (1.1) is studied by Li and Liu [16] using the perturbation analysis and Galerkin spectral approximations.

Recently, to investigate the evolution process of thin-film epitaxial growth, various numerical schemes for MBE model (1.1) have been developed including the first and second order convex splitting schemes [2, 24], the nonlinear Crank-Nicolson type scheme [21], the stabilized semi-implicit scheme [28] and so on. However, the analysis in those mentioned literatures was based on uniform time steps.

A feature of phase field models is that the solutions admit multiple time scales, namely, the dynamics evolves on a fast time scale at the beginning and coarsening evolves slowly on a time later. In this situation, the coarse-grained and refined time steps are useful to capture the multi-scale dynamics according to the slow and fast change of the solution itself. Thus, the BDF2 scheme with variable time steps is a good choice due to its strong stability for solving stiff or differential-algebraic problems [5, 10, 11, 23, 25, 26].

The BDF2 scheme with variable time steps has been widely developed [1, 3, 5, 13, 20, 29] for the stability and convergence analysis, including linear diffusion problems [20, 29], semilinear parabolic problems [5, 13] and the Cahn-Hilliard (CH) equation [3]. More specifically for the stability analysis of linear diffusion equations, twenty years ago Becker [1] presented the bound under the adjacent time-step ratio condition $ 0< r_k:=\tau_k/\tau_{k-1} \leq (2+\sqrt{13})/3\approx 1.868 $ that

$ \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 \rm{for} \quad n\geq 1, $

where $ \Gamma_n := \sum{k=2}^{n-2} \max\{0, r_k-r_{k+2}\} $. The result is also given in Thomée's classical book [25, Lemma 10.6]. As shown in [25] and [3], the magnitudes of $ \Gamma_n $ can be bounded [25, pp. 175] and unbounded [3, Remark 4.1] by choosing certain step-ratio sequence and vanishing step sizes. Emmrich [5] extends the Becker's condition to $ 0< r_k \leq 1.91 $, but still keeps the undesirable factor $ \exp(C\Gamma_n) $. Recently, Liao and Zhang [20] introduce the technique of the discrete orthogonal convolution (DOC) kernels, and improve Grigorieff's stability condition [8] nearly forty years ago {(one also refers to [4] and [11, Section III.5] a classical book by Hairer et al.)} from $ 0\leq r_k \leq 1+\sqrt{2} $ to $ 0\leq r_k \leq (3+\sqrt{17})/2 \approx 3.561 $. However, the second-order convergence in [20] suffers from 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.3)

While the stability and convergence analysis of BDF2 with variable time steps has brought the great challenge for linear problems, the analysis for nonlinear problems is even hard and still has a great progress. For instance, Chen et al. [3] replace $ \exp(C\Gamma_n) $ in Becker's estimate with a bounded factor $ \exp(Ct_n) $ with $ 0< r_k \leq 1.53 $ for CH equations.

Liao et al. [19] consider MBE model (1.1) with variable-time-steps BDF2 scheme, and obtain the second-order convergence under the ratio condition $ 0<r_k<3.561 $, but they still require an additional condition $ |\mathfrak{R}_p| \leq N_0 \ll N $.

The aim of this paper is to achieve the robust and sharp second-order error estimate for the variable time-steps BDF2 scheme under a new ratio condition $ 0<r_k< r_{\max} \approx 4.8645 $. Under this new ratio condition, we first prove the BDF2 scheme with BDF1 as starting step to preserve a modified energy dissipation law. After that, we carefully analyze the positive definiteness of discrete convolution kernels [20], and then introduce the discrete complementary convolution (DCC) kernels (defined in (4.21)) and the error convolution structure (ECS) with the BDF2 kernels (see Lemma 5.3), and finally obtain the sharp second-order convergence given as

$ \begin{align} \|e_h^{n}\|\leq &2\exp(16\mathcal{Q}_\delta^2t_{n-1}/\varepsilon)(\|e_h^0\|+2C_ut_nh^2+\sum\limits_{k=1}^n \tau_k^2\int_{t_{k-1}}^{t_k}\| u_{ttt}\| \, \mathrm{d} t\\ &+2t_n\max\limits_{1\leq k\leq n }{\tau_k}\int_{t_{k-1}}^{t_k}\| u_{ttt}\| \, \mathrm{d} t+4\tau\int^{t_1}_{0} \| u_{tt}\| \, \mathrm{d} t). \end{align} $

(1.4)

For brevity, we list the adjacent time step ratio condition as

A1:

$0< r_k \leq r_{\max}-\delta \quad {\rm{for\ any\ small\ constant}} 0<\delta < r_{\max}{\rm{and}}2\leq k \leq N$,

where the maximum ratio $ r_{\max} = \frac{1}{6} \left(\sqrt[3]{1196\!-\!12\sqrt{177}}+\sqrt[3]{1196\!+\!12\sqrt{177}}\right)+\frac43 \approx 4.8645 $ is the root of the cubic equation

$ \begin{align} r_{\max}^3=(2r_{\max}+1)^2. \end{align} $

(1.5)

The $ \mathcal{Q}_\delta $ in (1.4) is a constant depending on the choice of adjacent time steps, and has a upper bound of the form $ \mathcal{O}(1/\delta) $, where the parameter $ \delta $ is given in A1 and generally taken as a given small constant, for example $ \delta = 0.1 $ or other any small constant.

Comparing with recent results in [19], our second-order convergence is sharp and robust with the new ratio condition A1. The robustness means the convergence does not suffer from other extra conditions on the time step sizes, like the constrained condition $ |\mathfrak{R}_p| \leq N_0 \ll N $ in [19], expect for A1. In addition, our analysis shows that the sharp second-order convergence is consistent to the first-order BDF1 scheme for the first step solution $ u^1 $. It is the first time to make clear that the BDF1 scheme as start step to compute $ u^1 $ is enough to guarantee the second-order convergence of BDF2 schemes with variable time steps for MBE models. Numerical examples are provided to demonstrate our theoretical analysis.

The remainder is organized as follows. In section 2, we present the fully discrete scheme with variable time steps by using the finite difference method in space and BDF2 scheme in time. The solvability of the BDF2 scheme and the energy stability are presented in section 3. In section 4, we introduce the concepts of DOC and DCC kernels, and also present the properties of DOC kernels and DCC kernels. In section 5, we give the stability and second-order convergence analysis. Numerical simulations are carried out in section 6. We end the paper with a conclusion.

We take the generally variable time grids $ 0=t_0<t_1<t_2<\cdots<t_N=T $ {and denote} the $ k $th time-step size by $ \tau_k := t_k-t_{k-1} $ and the maximum time step size by $ \tau :=\max_{1\le k\le N}\tau_k $. The adjacent time-step ratio is defined by

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

Set $ u^k=u(\cdot, t_k) $ and the difference operator $ \nabla_\tau u^k=u^k-u^{k-1} $ for $ 1\leq k\leq N $. The BDF1 and BDF2 formulas with variable time steps are defined respectively 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*} $

Set the discrete convolution kernels $ b_{n-k}^{(n)} $ as $ b_{0}^{(1)} : =1/\tau_1 $ and

$ \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 \rm{and} \quad b_j^{(n)} = 0 \rm{ for } n, j\geq 2 . \end{align} $

(2.6)

Thus, we may reformulate the BDF1 and BDF2 into 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} $

(2.7)

The spacial domain $ \Omega=(0, L)^2 $ considered here is approximated by a uniform grid $ h = L/M $ for a positive integer $ M $, and the discrete domains are denoted by

$ \begin{align*} \Omega_h:=\{{\mathit{\boldsymbol{x}}}_h: (ih, jh), 1\leq i, j\leq M-1\}, \quad \bar{\Omega}_h:=\{{\mathit{\boldsymbol{x}}}_h: (ih, jh), 0\leq i, j\leq M\}. \end{align*} $

The partial derivatives $ \partial_x w $ and $ \partial_{xx} w $ are respectively approximated by the following operators

$ \begin{align*} \Delta_x w_{i, j}:=(w_{i+1, j}-w_{i-1, j})/(2h), \quad {\delta_x^2:=(w_{i+1, j}-2w_{i, j}+w_{i-1, j})/h^2.} \end{align*} $

The operators $ \Delta_y w_{i, j} $ and $ \delta_y^2 w_{i, j} $ can be defined similarly. Moreover, the discrete gradient operator and the discrete Laplacian operator are accordingly defined by

$ \begin{align*} \nabla_h w_{i, j}:=(\Delta_x w_{i, j}, \Delta_y w_{i, j})^T, \quad \Delta_h w_{i, j}:=(\delta_x^2+\delta_y^2)w_{i, j}. \end{align*} $

For the vector $ \mathit{\boldsymbol{u}}_{i, j}=(u^1_{i, j}, u^2_{i, j})^T $, the discrete divergence is defined by

$ \begin{align*} \nabla_h\cdot \mathit{\boldsymbol{u}}_{i, j}:=\Delta_x u^1_{i, j}+\Delta_y u^2_{i, j}. \end{align*} $

By using the finite difference method in space and BDF2 scheme in time, we have the fully discrete scheme with variable time steps as

$ \begin{align} \mathcal{D}_2 u^n_h + \varepsilon\Delta_h^2 u_h^n +\nabla_h \cdot \mathit{\boldsymbol{f}}(\nabla_h u^n_h)=0, \quad \rm{for} \quad 1\leq n\leq N. \end{align} $

(2.8)

Define the space of L-periodic grid functions as

$ \mathcal{V}_h = \{v_h| v_h \ \rm{ is\ L-periodic\ for }\ {{\bf x}}_h \in \bar{\Omega}_h \}. $

For any $ v, w\in \mathcal{V}_h $, the discrete $ L^2 $ inner product and norm are defined by

$ \langle{v, w}\rangle := h^2\sum\limits_{{{\bf x}}_h \in \Omega_h} v_hw_h, \quad \big\|{v}\big\|^2 := \langle{v, v}\rangle. $

The discrete norms $ \big\|{\nabla_h v}\big\| $ and $ \big\|{\Delta_h v}\big\| $ are defined by

$ \begin{align*} \|\nabla_h v\|:=\sqrt{h^2\sum\limits_{{\mathit{\boldsymbol{x}}}_h\in \Omega_h}|\nabla_h v_h|^2}, \quad \|\Delta_h v\|:=\sqrt{h^2\sum\limits_{{\mathit{\boldsymbol{x}}}_h\in \Omega_h}|\Delta_h v_h|^2}. \end{align*} $

For any $ v, w\in \mathcal{V}_h $, one has the discrete Green's formula with periodic boundary conditions

$ \begin{align} \langle -\nabla_h\cdot\nabla_h v, w\rangle=\langle \nabla_h v, \nabla_h w\rangle. \end{align} $

(3.9)

Lemma 3.1 (^[21])  For any grid function $ v\in \mathcal{V}_h $ and $ \epsilon>0 $, we have

$ \begin{align} \|\nabla_h v\|^2\leq \langle -\Delta_h v, v\rangle\leq \|\Delta_h v\|\|v\|\leq \epsilon\|\Delta_h v\|^2+\frac1{4\epsilon}\| v\|^2. \end{align} $

(3.10)

Theorem 3.2  If the time-step sizes $ \tau_n \leq 4\varepsilon $, the BDF2 time-stepping scheme (2.8) is convex and uniquely solvable.

The solvability of BDF2 scheme can be established by introducing a discrete energy functional $ G $ on the space $ \mathcal{V}_h $:

$ G[z] = \frac{1}{2} b_0^{(n)} \big\|{z-u^{n-1}_h}\big\|^2 + b_1^{(n)} \langle{\nabla_\tau u^{n-1}, z}\rangle \frac\varepsilon2 \big\|{\Delta_h z}\big\|^2 - \frac12 \langle \ln (1+|\nabla_h z|^2), 1\rangle. $

The detailed proof for the solvability of BDF2 scheme is given in [19] {and the key technique is referred to [27] }.

We now consider the energy stability of BDF2 scheme (2.8). To do so, we first present the positive definiteness of discrete convolution kernels $ b^{(n)}_{n-k} $ in the following lemma.

Lemma 3.2   Assume the time-step ratio $ r_k $ satisfies A1. For any real sequence $ \{w_k\}_{k=1}^n $, it holds for $ \epsilon_* = \sqrt[3]{12} \big(\sqrt[3]{\sqrt{177} + 9} - \sqrt[3]{\sqrt{177} - 9} \;\big)/6 \approx 0.4534 $ and { for any small constant $ 0<\delta < r_{\max} $ (see A1) } that

$ \begin{align} &2w_k \sum\limits_{j= 1}^k b^{(k)}_{k-j} w_j \geq \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{\delta w_k^2}{(1+r_{\max})^2\epsilon_*\tau_k}, \quad k \geq 2, \end{align} $

(3.11)

$ \begin{align} &2\sum\limits_{k=1}^nw_k \sum\limits_{j= 1}^k b^{(k)}_{k-j} w_j \geq \sum\limits_{k=1}^n\frac{\delta w_k^2}{(1+r_{\max})^2\epsilon_*\tau_k}\geq 0, \quad \rm{for } n \geq 1. \end{align} $

(3.12)

Proof Denote the multi-variable functions

$ \begin{align*} &\mathfrak{F}(x, y, \epsilon) = \frac{2\epsilon+4\epsilon x -\epsilon^2 x^2}{(1+x)} - \frac{y}{(1+y)}, \quad \rm{for } x, y, \epsilon \geq 0. \end{align*} $

It follows from [29, Lemma 2.1] and the proof of [29, Lemma 2.2] that

$ \begin{align} 2w_k \sum\limits_{j= 1}^k b^{(k)}_{k-j} w_j \geq& \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}, \quad k \geq 2, \end{align} $

(3.13)

and

$ \begin{align*} \frac{2\epsilon_*+4\epsilon_* r_k -\epsilon_*^2 r_k^2}{1+r_k}\geq \frac{r_{\max}}{1+r_{\max}}, \quad \forall \; 0<r_k\leq r_{\max}. \end{align*} $

Hence, for any $ 0<r_k\leq r_{\max}-\delta $, one has

$ \begin{align*} \mathfrak{F}(r_k, r_{k+1}, \epsilon_*)\geq \frac{r_{\max}}{1+r_{\max}}-\frac{r_{\max}-\delta}{1+r_{\max}-\delta}=\frac{\delta}{(1+r_{\max}-\delta)(1+r_{\max})}\geq \frac{\delta}{(1+r_{\max})^2}, \quad {k\geq 2} \end{align*} $

where the monotony of the function $ h(x) = x/(1+x) $ is used. Inserting above inequality to (3.13), one immediately has the inequality (3.11). Summing the inequality (3.11) from 1 to $ n $, one has

$ \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}}{(1+r_{n+1})}\frac{w_n^2}{\epsilon_*\tau_n} - \frac{r_2}{ (1+r_2)}\frac{w_{1}^2 }{\epsilon_*\tau_{1}}+\sum\limits_{k=2}^n\frac{\delta w_k^2}{(1+r_{\max})^2\epsilon_*\tau_k}\\ &\geq \left(2\epsilon_*-\frac{ r_{\max}-\delta}{ 1+r_{\max}-\delta}\right){\frac{w_{1}^2 }{\epsilon_*\tau_{1}}}+\sum\limits_{k=2}^n\frac{\delta w_k^2}{(1+r_{\max})^2\epsilon_*\tau_k}\\ &\geq \left(2-\frac{ r_{\max}}{ (1+r_{\max})\epsilon_*}\right){\frac{w_{1}^2 }{\tau_{1}}}+\sum\limits_{k=1}^n\frac{\delta w_k^2}{(1+r_{\max})^2\epsilon_*\tau_k}\\ &\geq \sum\limits_{k=1}^n\frac{\delta w_k^2}{(1+r_{\max})^2\epsilon_*\tau_k}, \quad n \geq 1. \end{align*} $

The proof is complete.

We now consider the energy stability of BDF2 scheme (2.8) by defining the modified discrete energy

$ \begin{align} E^n:=\frac{r_{n+1}}{2(1+r_{n+1})\epsilon_*\tau_n}\|\nabla_\tau u^n\|^2+\frac{\varepsilon}{2}\|\Delta_h u^n\|^2-\frac1{2}\langle\ln(1+|\nabla_h u^n|^2), 1\rangle, \quad n\geq 1, \end{align} $

(3.14)

and the initial energy $ E^0 :=\frac{\varepsilon}{2}\|\Delta_h u^0\|^2-\frac1{2}\langle\ln(1+|\nabla_h u^0|^2), 1\rangle $.

To establish the energy dissipation law, we need the time-step ratio $ r_k $ to hold A1 and the time step size $ \tau_n $ to satisfy

$ \begin{align} \tau_n\leq 4\varepsilon \min\{2-\frac{r_{\max}}{(1+r_{\max})\epsilon_*}, \frac{\delta}{(1+r_{\max})^2\epsilon_*}\}. \end{align} $

(3.15)

Theorem 3.2   Assume the time-step ratio condition A1 holds with the time-step condition (3.15), then the discrete energy $ E_n $ defined in (3.14) satisfies

$ \begin{align*} E^n\leq E^{n-1}\leq E^0, \quad n\geq 1. \end{align*} $

Proof   Taking the inner product on both sides of (2.8) with $ \nabla_\tau u^n $, one has

$ \begin{align} \langle\mathcal{D}_2 u^n, \nabla_\tau u^n\rangle+ \varepsilon\langle \Delta_h^2 u^n, \nabla_\tau u^n\rangle +\langle\nabla_h \cdot \mathit{\boldsymbol{f}}(\nabla_h u^n), \nabla_\tau u^n\rangle=0, \quad \rm{for} \quad 1\leq n\leq N. \end{align} $

(3.16)

Due to the periodic boundary conditions, the summation by parts argument holds, which implies

$ \begin{align} \varepsilon\langle \Delta_h^2 u^n, \nabla_\tau u^n\rangle=\varepsilon\langle \Delta_h u^n, \Delta_h\nabla_\tau u^n\rangle=\frac{\varepsilon}{2}(\|\Delta_h u^n\|^2-\|\Delta_h u^{n-1}\|^2+\|\Delta_h \nabla_\tau u^n\|^2). \end{align} $

(3.17)

By using the inequality $ \frac{x}{1+x}\leq \ln(1+x) $ with $ x=(\nabla_\tau|\nabla_h u^n|^2)/(1+|\nabla_h u^{n-1}|^2) $, one has

$ \frac{\nabla_\tau|\nabla_h u^n|^2}{1+|\nabla_h u^{n}|^2}\leq \ln(\frac{1+|\nabla_h u^n|^2}{1+|\nabla_h u^{n-1}|^2}), $

which together with the discrete Green's formula (3.9) and the inequality (3.10) with $ \epsilon=\varepsilon $ imply that

$ \begin{align} &\langle\nabla_h \cdot \mathit{\boldsymbol{f}}(\nabla_h u^n), \nabla_\tau u^n\rangle=-\langle\mathit{\boldsymbol{f}}(\nabla_h u^n), \nabla_h \nabla_\tau u^n\rangle\\ =&{-\langle\frac{\nabla_\tau|\nabla_h u^n|^2}{2(1+|\nabla_h u^n|^2)}, 1\rangle-\langle\frac{|\nabla_h\nabla_\tau u^n|^2}{2(1+|\nabla_h u^n|^2)}, 1\rangle}\\ \geq& -\frac12\langle\ln(1+|\nabla_h u^n|^2), 1\rangle+\frac12\langle\ln(1+|\nabla_h u^{n-1}|^2), 1\rangle-\frac12\|\nabla_h\nabla_\tau u^n\|^2\\ \geq& -\frac12\langle\ln(1+|\nabla_h u^n|^2), 1\rangle+\frac12\langle\ln(1+|\nabla_h u^{n-1}|^2), 1\rangle-\frac{\varepsilon}{2}\|\Delta_h\nabla_\tau u^n\|^2-\frac{1}{8\varepsilon}\|\nabla_\tau u^n\|^2. \end{align} $

(3.18)

For $ n\geq 2 $, it follows from Lemma 3.2 and the time-step condition (3.15) that

$ \begin{align} \langle\mathcal{D}_2 u^n, \nabla_\tau u^n\rangle&\geq \frac{r_{n+1}}{2(1+r_{n+1})}\frac{\|\nabla_\tau u^n\|^2}{\epsilon_*\tau_n} - \frac{r_n}{ 2(1+r_n)}\frac{\|\nabla_\tau u^{n-1}\|^2 }{\epsilon_*\tau_{n-1}}+\frac{\delta \|\nabla_\tau u^n\|^2}{2(1+r_{\max})^2\epsilon_*\tau_n}\\ &\geq \frac{r_{n+1}}{2(1+r_{n+1})}\frac{\|\nabla_\tau u^n\|^2}{\epsilon_*\tau_n} - \frac{r_n}{2 (1+r_n)}\frac{\|\nabla_\tau u^{n-1}\|^2 }{\epsilon_*\tau_{n-1}}+\frac1{8\varepsilon}\|\nabla_\tau u^n\|^2. \end{align} $

(3.19)

Hence, by inserting (3.17)-(3.19) into (3.16), we have

$ \begin{align*} E^n\leq E^{n-1}, \quad n\geq 2. \end{align*} $

For $ n=1 $, it follows from the condition A1 that

$ \frac{r_2}{(1+r_2)\epsilon_*}\leq \frac{r_{\max}}{(1+r_{\max})\epsilon_*} \leq 2, $

which together with the time step condition (3.15) imply that

$ \begin{align*} \langle\mathcal{D}_2 u^1, \nabla_\tau u^1\rangle=\frac1{\tau_1}\|\nabla_\tau u^1\|^2\geq \frac{r_2}{2(1+r_2)\epsilon_*\tau_1}\|\nabla_\tau u^1\|^2+\frac1{8\varepsilon}\|\nabla_\tau u^1\|^2. \end{align*} $

Thus, by inserting the above inequality and (3.17)-(3.18) into (3.16), we have

$ \begin{align*} E^1\leq E^{0}. \end{align*} $

The proof is complete.

To obtain stability analysis of {BDF2} scheme (2.8), we introduce the discrete complementary convolution (DCC) kernels $ p_{n-j}^{(n)} $ 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} $

(4.20)

As the identity (4.20) holds for all $ n\geq 1 $, it only requires

$ \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} $

(4.21)

The idea of DCC kernels has been successfully applied to the stability analysis for subdiffusion problems [14, 15, 17] and reaction-diffusion problem [29].

We now introduce the discrete orthogonal convolution (DOC) kernels by

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

(4.22)

where $ \delta_{nk} $ is the Kronecker delta symbol with $ \delta_{nk} = 1 $ if $ n=k $ and $ \delta_{nk} = 0 $ if $ n\neq k $. From (4.22), it holds

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

(4.23)

The DCC and DOC kernels have a close connection. In fact, summing from 1 to $ n $ with $ k $ on both sides of (4.23), and then exchanging the order of summation, we have

$ \begin{equation} \sum\limits_{j=1}^n\mathcal{D}_2 u^j\sum\limits_{k=j}^n\theta_{k-j}^{(k)} = u^n - u^0. \end{equation} $

(4.24)

One can compare the identity (4.24) with (4.20) to find that (also see [29])

$ \begin{align} &p_{n-j}^{(n)} = \sum\limits_{k=j}^n\theta_{k-j}^{(k)} \; \;\;(1\leq j\leq n). \end{align} $

(4.25)

From (4.25), the direct calculation leads to another relation between DOC and DCC kernels

$ \begin{align} &\theta_{0}^{(n)} = p_{0}^{(n)}, \quad \theta_{n-k}^{(n)} = p_{n-k}^{(n)} - p_{n-k-1}^{(n)} \; \;\;(1\leq k\leq n-1). \end{align} $

(4.26)

To establish the stability and error estimate, here we streamline the useful results in [20, 29].

Lemma 4.1 [20, Lemma 2.2]  Assume the BDF2 kernels $ b_{n-k}^{(n)} $ defined in (2.6) are positive definite. Then the DOC kernels $ \theta^{(n)}_{n-k} $ defined in (4.22) are also positive definite. This is, for any real sequence $ \{\omega_j\}_{j=1}^n $, it holds 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 4.2 [20, Corollary 2.1]  The DOC kernels $ \theta_{n-j}^{(n)} $ have the following properties:

$ \begin{align} &\theta_{n-j}^{(n)}>0, \quad \rm{for any}\quad 1\leq j\leq n, 1\leq n\leq N, \end{align} $

(4.27)

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

(4.28)

proposition 4.1  [29, Proposition 2.2]  Let $ \tau $ be the maximum time step size and $ r_{*} $ be any given positive constant. If the time step ratio satisfies $ 0<r_k\leq r_{*} $, then the DCC kernels $ p^{(n)}_{n-k} $ defined in (4.21) satisfy

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

(4.29)

$ \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} $

(4.30)

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

(4.31)

$ \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} $

(4.32)

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

The BDF2 kernels $ b^{(n)}_{n-k} $ and DOC kernels $ \theta_{n-k}^{(n)} $ defined in (2.6) and (4.22) respectively can be represented as the following matrix forms [18]

$ \begin{align*} B_2:=\left( \begin{array}{cccc} b^{(1)}_0&&&\\ b^{(2)}_1&b^{(2)}_0&&\\ &\ddots&\ddots&\\ &&b^{(n)}_1&b^{(n)}_0 \end{array} \right)_{n\times n}\ \rm{and}\quad \Theta_2:=\left( \begin{array}{cccc} \theta^{(1)}_0&&&\\ \theta^{(2)}_1&\theta^{(2)}_0&&\\ \vdots&\ddots&\ddots&\\ \theta^{(n)}_{n-1}&\cdots&\theta^{(n)}_1&\theta^{(n)}_0 \end{array} \right)_{n\times n} \; . \end{align*} $

It follows from the definition of DOC kernels $ \theta^{(n)}_{n-k} $ in (4.22) that

$ \begin{align} \Theta_2 = B_2^{-1}. \end{align} $

(4.33)

Assume A1 holds, then Lemmas 3.2 and 4.1 imply the real symmetric matrices

$ \begin{align*} B:=B_2+B_2^T\ \quad \rm{and}\ \quad \Theta:=\Theta_2+\Theta_2^T \end{align*} $

are positive definite.

Define the diagonal matrix $ \Lambda_\tau:=\rm{diag}(\sqrt{\tau_1}, \cdots, \sqrt{\tau_n}) $ and

$ \begin{align} \tilde{B}_2:=\Lambda_\tau B_2\Lambda_\tau =\left(\begin{array}{cccc} \tilde{b}^{(1)}_0&&&\\ \tilde{b}^{(2)}_1&\tilde{b}^{(2)}_0&&\\ &\ddots&\ddots&\\ &&\tilde{b}^{(n)}_1&\tilde{b}^{(n)}_0 \end{array}\right)_{n\times n, } \end{align} $

(4.34)

with

$ \begin{align*} \tilde{b}^{(1)}_0 = 1, \quad \tilde{b}^{(k)}_0 =\frac{1+2r_k}{1+r_k}, \quad \tilde{b}^{(k)}_1=-\frac{r_k^{3/2}}{1+r_k}, \quad 2\leq k\leq n. \end{align*} $

Moreover, we define the real symmetric matrix $ \tilde{B}:=\tilde{B}_2+\tilde{B}^T_2 $, which has the following properties.

Lemma 4.3  Assume A1 holds, the minimum eigenvalue of $ \tilde{B} $ can be bounded by

$ \begin{align*} \lambda_{\min}(\tilde{B})\geq \min\limits_{1\leq k\leq n} \tilde{R}(r_k, r_{k+1})\geq C_\delta, \end{align*} $

where

$ \begin{align} &\tilde{R}(x, y) = \frac{2+4x-x^{3/2}}{1+x}-\frac{y^{3/2}}{1+y}, \quad 0<x, y\leq r_{\max}, \end{align} $

(4.35)

$ \begin{align} &{C}_\delta = \min\{\tilde{R}(0, r_{\max}-\delta), \tilde{R}(r_{\max}-\delta, r_{\max}-\delta)\}. \end{align} $

(4.36)

Thus $ \tilde{B} $ is positive definite and there exists a non-singular upper triangular matrix $ L $ such that

$ \begin{align*} \tilde{B}=\Lambda_\tau B \Lambda_\tau=L^TL\quad\quad \rm{or}\quad\quad B =(L\Lambda_\tau^{-1})^TL\Lambda_\tau^{-1}. \end{align*} $

For brevity, we leave the detailed proof in the Appendix.

Remark 1  It follows from the A1 and the monotony of $ \tilde{R}(x, y) $ for $ x>1, y>0 $ that

$ \begin{align*} \tilde{R}(0, r_{\max}-\delta)\geq \tilde{R}(0&, r_{\max}) = 2-\frac{r_{\max}^{3/2}}{1+r_{\max}}>0\\ \tilde{R}(r_{\max}-\delta, r_{\max}-\delta)&\geq \tilde{R}(r_{\max}, r_{\max}-\delta)=\frac{r_{\max}^{3/2}}{1+r_{\max}}-\frac{(r_{\max}-\delta)^{3/2}}{1+r_{\max}-\delta}\\ & =\frac{\sqrt{r_{\max}}(r_{\max}+r_{\max}^2-\delta r_{\max})-\sqrt{r_{\max}-\delta}(r_{\max}+r_{\max}^2-\delta-\delta r_{\max})}{(1+r_{\max})(1+r_{\max}-\delta)}\\ & \geq \frac{\sqrt{r_{\max}}}{(1+r_{\max})^2}\delta, \quad \rm{for} \quad 0<\delta<r_{\max}-1, \end{align*} $

which implies the positive constant $ C_\delta $ depends on $ \delta $. Moreover, if $ 0<\delta\leq r_{\max}-4 $, one has

$ \begin{align*} \tilde{R}(0, r_{\max}-\delta)\geq \tilde{R}(r_{\max}-\delta, r_{\max}-\delta). \end{align*} $

Then, for any $ 0<\delta\leq r_{\max}-4 $, the constant $ C_\delta $ can be estimated by

$ \begin{align*} C_\delta=\tilde{R}(r_{\max}-\delta, r_{\max}-\delta)\geq \frac{\sqrt{r_{\max}}}{(1+r_{\max})^2}\delta. \end{align*} $

Thus, the lower bound of $ C_\delta $ is $ \mathcal{O}(\delta) $ for small $ \delta $.

The next Lemma gives an upper bound of the maximum singular value of $ \tilde{B}_2 $.

Lemma 4.4  If A1 holds, then the maximum eigenvalue of the real symmetric matrix $ \tilde{B}_2^T\tilde{B}_2 $ can be bounded by

$ \begin{align*} \lambda_{\max}(\tilde{B}_2^T\tilde{B_2})\leq \hat{R}(r_{\max}-\delta, r_{\max}-\delta)\leq \frac{4r_{\max}^3}{(1+r_{\max})^2}<14, \end{align*} $

where $ \tilde{B}_2 $ is defined by (4.34) and $ \hat{R}(x, y) $ is defined by

$ \begin{align*} \hat{R}(x, y):=\frac{(1+2x)(1+2x+x^{3/2})}{(1+x)^2}+\frac{y^{3/2}(1+2y+y^{3/2})}{(1+y)^2}. \end{align*} $

Again for brevity, we leave the detailed proof in the Appendix.

To deal with the nonlinear term $ \nabla_h \cdot \mathit{\boldsymbol{f}}(\nabla_h u^n) $, we now introduce the following matrices

$ \begin{align*} \mathit{\boldsymbol{\widehat{B}}} := B\otimes I_2, \quad \mathit{\boldsymbol{\widehat{\Theta}}}:=\Theta\otimes I_2, \quad \mathit{\boldsymbol{\widetilde{B}_2}}:=\tilde{B}_2\otimes I_2, \quad \mathit{\boldsymbol{\widetilde{B}}}:=\tilde{B}\otimes I_2, \quad \mathit{\boldsymbol{\widehat{\Lambda}_\tau}} = \Lambda_\tau\otimes I_2, \quad \mathit{\boldsymbol{\widehat{L}}}:=L\otimes I_2. \end{align*} $

One can use Lemma 4.3 to derive that

$ \begin{align*} \|\widehat{\mathit{\boldsymbol{L}}}^{-1}\|_2^2 = \lambda_{\max}((\widehat{\mathit{\boldsymbol{L}}}^T\widehat{\mathit{\boldsymbol{L}}})^{-1})=\lambda^{-1}_{\min}((L^TL)\otimes I_2)=\lambda^{-1}_{\min}(\mathit{\boldsymbol{\widetilde{B}}}). \end{align*} $

For convenience, we define the vector norm $ \|\cdot\|_2 $ by $ \|\mathit{\boldsymbol{u}}\|_2=\sqrt{\mathit{\boldsymbol{u}}^T\mathit{\boldsymbol{u}}} $ and the associated matrix norm by $ \|A\|_2:=\sqrt{\lambda_{\max}(A^TA)} $. We also define

$ \begin{align} \mathcal{Q}_\delta:=\max\limits_{n\geq 1}\|\mathit{\boldsymbol{\widetilde{B_2}}}\|_2^2\|\mathit{\boldsymbol{\widehat{L}}}^{-1}\|_2^4=\max\limits_{n\geq 1}\frac{\lambda_{\max}(\mathit{\boldsymbol{\widetilde{B_2}}}^T\mathit{\boldsymbol{\widetilde{B_2}}})}{\lambda_{\min}^2(\mathit{\boldsymbol{\widetilde{B}}})}. \end{align} $

(4.37)

Note that under the codition A1, Lemmas 4.3 and 4.4 imply the positive constant $ \mathcal{Q}_\delta<\frac{14}{{C}_\delta} $, where $ {C}_\delta $ is defined by (4.36). By taking $ \delta = 1.303 $, one can obtain that $ \mathcal{Q_\delta}< 39 $, which is consistent with the one in [19].

We now present several important Lemmas, which plays a key role in dealing with the nonlinear term $ \nabla_h\cdot \mathit{\boldsymbol{f}}(\nabla_h u^n) $, and leave the proofs in the Appendix for brevity due to their similarity with [19].

Lemma 4.5  Assume A1 holds. For the positive definite matrix $ \mathit{\boldsymbol{\widehat{\Theta}}} = \mathit{\boldsymbol{\widehat{B}}}_2^{-1}\mathit{\boldsymbol{\widehat{B}}}(\mathit{\boldsymbol{\widehat{B}}}_2^{-1})^T $, and any vector sequences $ \{\mathit{\boldsymbol{v}}^k\}_{k=1}^n, \{\mathit{\boldsymbol{w}}^k\}_{k=1}^n, {\mathit{\boldsymbol{v}}}^k, {\mathit{\boldsymbol{w}}}^k\in \mathbb{R}^2 $, it holds

$ \begin{align} \sum\limits_{k=1}^n\sum\limits_{j=1}^k\theta^{(k)}_{k-j}(\mathit{\boldsymbol{v}}^k)^T\mathit{\boldsymbol{w}}^j\leq \frac{\epsilon}{2} {\mathit{\boldsymbol{v}}}^T\mathit{\boldsymbol{\widehat{\Theta}}} {\mathit{\boldsymbol{v}}}+ \frac1{2\epsilon} {\mathit{\boldsymbol{w}}}^T\mathit{\boldsymbol{\widehat{B}}}^{-1} {\mathit{\boldsymbol{w}}}, \quad \forall \epsilon>0, \end{align} $

(4.38)

where $ {\mathit{\boldsymbol{v}}} = (( {\mathit{\boldsymbol{v}}}^1)^T, \cdots, ( {\mathit{\boldsymbol{v}}}^n)^T)^T $ and $ {\mathit{\boldsymbol{w}}} = (( {\mathit{\boldsymbol{w}}}^1)^T, \cdots, ( {\mathit{\boldsymbol{w}}}^n)^T)^T $.

Lemma 4.6  ([12, Lemma 3.5]) For any $ \mathit{\boldsymbol{v}}, \mathit{\boldsymbol{w}}\in \mathbb{R}^2 $, there exists a symmetric matrix $ Q_{\mathit{\boldsymbol{f}}}\in \mathbb{R}^{2\times 2} $ such that

$ \begin{align*} \mathit{\boldsymbol{f}}(\mathit{\boldsymbol{v}})-\mathit{\boldsymbol{f}}(\mathit{\boldsymbol{w}}) = Q_{\mathit{\boldsymbol{f}}}(\mathit{\boldsymbol{v}}-\mathit{\boldsymbol{w}}). \end{align*} $

The eigenvalues $ \lambda_1, \lambda_2 $ of $ Q_{\mathit{\boldsymbol{f}}} $ satisfy $ -\frac18\leq \lambda_1, \lambda_2\leq 1 $. Consequently, it holds that

$ \begin{align} |\mathit{\boldsymbol{f}}(\mathit{\boldsymbol{v}})-\mathit{\boldsymbol{f}}(\mathit{\boldsymbol{w}})| \leq |\mathit{\boldsymbol{v}}-\mathit{\boldsymbol{w}}|. \end{align} $

(4.39)

We now give another important lemma as follows.

Lemma 4.7  Assume A1 holds, then for any vector sequence $ \mathit{\boldsymbol{v}}^k = (v^k_1, v^k_2)^T, \mathit{\boldsymbol{w}}^k = (w^k_1, w^k_2)^T, \mathit{\boldsymbol{z}}^k = (z^k_1, z^k_2)^T $ with $ 1\leq k\leq n $ and any $ \epsilon>0 $, it holds that

$ \begin{align*} \sum\limits_{k=1}^n\sum\limits_{j=1}^k \theta_{k-j}^{(k)}(\mathit{\boldsymbol{z}}^k)^T[ {\mathit{\boldsymbol{f}}}( {\mathit{\boldsymbol{v}}}^j+ {\mathit{\boldsymbol{w}}}^j)- {\mathit{\boldsymbol{f}}}( {\mathit{\boldsymbol{v}}}^j)]\leq \sum\limits_{k=1}^n\sum\limits_{j=1}^k \theta^{(k)}_{k-j} \left[\epsilon (\mathit{\boldsymbol{z}}^k)^T\mathit{\boldsymbol{z}}^j+\frac{\mathcal{Q}_\delta}{\epsilon}(\mathit{\boldsymbol{w}}^k)^T\mathit{\boldsymbol{w}}^j\right], \end{align*} $

where the positive constant $ \mathcal{Q}_\delta $ is defined in (4.37). Moreover,

$ \begin{align*} \sum\limits_{k=1}^n\sum\limits_{j=1}^k \theta_{k-j}^{(k)}(\mathit{\boldsymbol{z}}^k)^T[ {\mathit{\boldsymbol{f}}}( {\mathit{\boldsymbol{v}}}^j+\mathit{\boldsymbol{z}}^j)- {\mathit{\boldsymbol{f}}}( {\mathit{\boldsymbol{v}}}^j)]\leq 2\sqrt{\mathcal{Q}_\delta}\sum\limits_{k=1}^n\sum\limits_{j=1}^k \theta^{(k)}_{k-j} (\mathit{\boldsymbol{z}}^j)^T\mathit{\boldsymbol{z}}^k. \end{align*} $

In this section, we consider the stability and convergence analysis for the BDF2 scheme (2.8), which requires a discrete Grönwall inequality given as follows.

Lemma 5.1   Assume that $ \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 \rm{for} \quad 1\leq n \leq N, $

then it holds that

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

We now consider the $ L^2 $-norm stability analysis of BDF2 scheme (2.8) with variable time steps. It is noted that the inequality (3.10) plays an important role in [21] to obtain the $ L^2 $-norm error estimate on uniform time steps. But the method developed in [21] fails to prove the $ L^2 $-norm estimate if the DOC technique is used in this paper. In other words, the inequality (3.10) can not be used to have $ L^2 $-norm estimate since the DOC technique will lead to the cross inner product $ \langle\phi^k, \phi^j\rangle $ for different time levels. Alternatively, we here develop a new inequality to fill this gap as follows.

Lemma 5.2   Let the matrix $ \mathit{\boldsymbol{\beta}} := \mathit{\boldsymbol{\beta}}_2+\mathit{\boldsymbol{\beta}}_2^T $ be positive definite with $ \mathit{\boldsymbol{\beta}}_2:= (\beta_{k, j}), 1\leq j, k\leq n, \beta_{k, j}= 0 $ if $ j>k $. Assume $ \beta_{k, j}\geq0 $ for any $ 1\leq k, j\leq n $, then it holds

$ \begin{align} \sum\limits_{k=1}^n\sum\limits_{j=1}^k \beta_{k, j}\langle \nabla_h \phi_h^j, \nabla_h \phi_h^k\rangle \leq \sum\limits_{k=1}^n\sum\limits_{j=1}^k\beta_{k, j}\langle -\Delta_h \phi_h^j, \phi_h^k\rangle. \end{align} $

(5.40)

Proof  For notation convenience, we denote the operator $ \mathcal{L}_h := \nabla_h\cdot \nabla_h-\Delta_h $. It follows from the positive definiteness of the matrix $ \mathit{\boldsymbol{\beta}} $ and the standard Cholesky decomposition that there exists a non-singular upper triangular matrix $ A $ such that $ \mathit{\boldsymbol{\beta}} = A^TA $. Denote the matrix $ A = (\mathit{\boldsymbol{a}}_1, \cdots, \mathit{\boldsymbol{a}}_n) $, where $ \mathit{\boldsymbol{a}}_l = (a_{1, l}, \cdots, a_{n, l})^T, 1\leq l\leq n $ and $ a_{l, k}=0 $ if $ l>k $. It is easy to verify that $ \beta_{k, j} = \mathit{\boldsymbol{a}}_k^T\mathit{\boldsymbol{a}}_j = \sum_{l=1}^n a_{l, k}a_{l, j}, k\neq j $, $ \beta_{k, k} = \mathit{\boldsymbol{a}}_k^T\mathit{\boldsymbol{a}}_k/2= \sum_{l=1}^n a_{l, k}^2/2, 1\leq k\leq n $. The discrete Green's formula and the identity $ \langle \mathcal{L}_h v, w\rangle =\langle v, \mathcal{L}_hw\rangle $ yield

$ \begin{align*} \frac12\sum\limits_{l=1}^n \langle \mathcal{L}_h(\sum\limits_{k=1}^n a_{l, k}\phi_h^k), \sum\limits_{k=1}^n a_{l, k}\phi_h^k\rangle&=\frac12\sum\limits_{k=1}^n (\sum\limits_{l=1}^n a_{l, k}^2) \langle \mathcal{L}_h\phi_h^k, \phi_h^k\rangle +\sum\limits_{k=2}^n \sum\limits_{j=1}^{k-1}(\sum\limits_{l=1}^{n} a_{l, k}a_{l, j}) \langle \mathcal{L}_h\phi_h^k, \phi_h^j\rangle \\ &=\sum\limits_{k=1}^n\beta_{k, k} \langle \mathcal{L}_h\phi_h^k, \phi_h^k\rangle +\sum\limits_{k=2}^n \sum\limits_{j=1}^{k-1}\beta_{k, j}\langle \mathcal{L}_h\phi_h^k, \phi_h^j\rangle\\ &=\sum\limits_{k=1}^n\sum\limits_{j=1}^k \beta_{k, j}\langle \mathcal{L}_h \phi_h^k, \phi_h^j\rangle. \end{align*} $

It follows from the discrete Green's formula and the inequality (3.10) that $ \langle \mathcal{L}_h v, v\rangle\geq 0 $ for any $ v\in \mathbb{V}_h $, which implies

$ \begin{align*} \sum\limits_{k=1}^n\sum\limits_{j=1}^k \beta_{k, j}\langle \mathcal{L}_h \phi_h^k, \phi_h^j\rangle\geq 0. \end{align*} $

The proof is complete.

We now consider the $ L^2 $-norm stability of BDF2 with variable time steps.

Theorem 5.1   Assume A1 holds and the maximum time step $ \tau\leq 16\mathcal{Q}_\delta^2/{\varepsilon} $. Then the BDF2 scheme (2.8) is stable in $ L^2 $-norm, and has the estimate

$ \begin{align} \|\hat{u}_h^n-u_h^n\|\leq 2\exp(16\mathcal{Q}_\delta t_{n-1}/\varepsilon)\|\hat{u}_h^0-u_h^0\|, \end{align} $

(5.41)

where $ u_h^n $ and $ \hat{u}_h^n $ are the solutions to (2.8) with the initial values $ u_h^0 $ and $ \hat{u}_h^0 $, respectively.

Proof   Let $ \phi_h^n:=\hat{u}_h^n-{u}_h^n $ ($ 0\leq n\leq N $) be the solution perturbation for $ \mathit{\boldsymbol{x}}_h\in \bar{\Omega}_h $. The BDF2 scheme (2.8) implies that $ \phi_h^n $ solves

$ \begin{align*} \mathcal{D}_2\phi_h^j+\varepsilon \Delta_h^2\phi_h^j+\nabla_h\cdot ( {\mathit{\boldsymbol{f}}}(\nabla_h\hat{u}_h^j)- {\mathit{\boldsymbol{f}}}(\nabla_h{u}_h^j))=0. \end{align*} $

Multiplying both sides of the above identity by $ \theta_{k-j}^{(k)} $ and summing $ j $ from 1 to $ k $, one can use the identity (4.23) to obtain

$ \begin{align*} \nabla_\tau\phi_h^k+\varepsilon \sum\limits_{j=1}^k\theta_{k-j}^{(k)}\Delta_h^2\phi_h^j+\sum\limits_{j=1}^k\theta_{k-j}^{(k)}\nabla_h\cdot ( {\mathit{\boldsymbol{f}}}(\nabla_h\hat{u}_h^j)- {\mathit{\boldsymbol{f}}}(\nabla_h{u}_h^j))=0. \end{align*} $

Taking the inner product of the above identity with $ 2\phi^k $ and summing $ k $ from 1 to $ n $, one has

$ \begin{align} \|\phi_h^n\|^2-\|\phi_h^0\|^2 &+2\varepsilon\sum\limits_{k=1}^n\sum\limits_{j=1}^k\theta_{k-j}^{(k)}\langle\Delta_h\phi_h^j, \Delta_h\phi_h^k\rangle \\ &\leq \quad 2\sum\limits_{k=1}^n\sum\limits_{j=1}^k\theta_{k-j}^{(k)}\langle {\mathit{\boldsymbol{f}}}(\nabla_h\hat{u}_h^j)- {\mathit{\boldsymbol{f}}}(\nabla_h{u}_h^j), \nabla_h\phi_h^k\rangle, \end{align} $

(5.42)

where the discrete Green's formula (3.9) and $ 2a^2-2ab\geq a^2-b^2 $ are used. We now deal with the rightmost term of (5.42). In view of $ \nabla_h\hat{u}_h^j=\nabla_h{u}_h^j+ \nabla_h\phi_h^j $, one can apply Lemmas 5.2 and 4.7 to obtain

$ \begin{align} 2\sum\limits_{k=1}^n\sum\limits_{j=1}^k&\theta_{k-j}^{(k)}\langle ( {\mathit{\boldsymbol{f}}}(\nabla_h\hat{u}_h^j)- {\mathit{\boldsymbol{f}}}(\nabla_h{u}_h^j)), \nabla_h\phi_h^k\rangle\\ &\leq 4\sqrt{\mathcal{Q}_\delta}\sum\limits_{k=1}^n\sum\limits_{j=1}^k \theta^{(k)}_{k-j} \langle \nabla_h \phi_h^j, \nabla_h \phi_h^k\rangle\\ &\leq 4\sqrt{\mathcal{Q}_\delta}\sum\limits_{k=1}^n\sum\limits_{j=1}^k \theta^{(k)}_{k-j} \langle -\Delta_h \phi_h^j, \phi_h^k\rangle \end{align} $

(5.43)

$ \begin{align} &\leq 2\varepsilon \sum\limits_{k=1}^n\sum\limits_{j=1}^k\theta_{k-j}^{(k)}\langle \Delta_h \phi_h^j, \Delta_h \phi_h^k\rangle+\frac{8\mathcal{Q}_\delta^2}{\varepsilon}\sum\limits_{k=1}^n\sum\limits_{j=1}^k\theta_{k-j}^{(k)}\langle \phi_h^j, \phi_h^k\rangle. \end{align} $

(5.44)

The last inequality holds by taking $ \epsilon=\varepsilon/(2\sqrt{\mathcal{Q}_\delta}) $ in Lemma 4.7 since Lemma 4.7 still holds for the linear function $ \mathit{\boldsymbol{f}}(\mathit{\boldsymbol{v}}):=\mathit{\boldsymbol{v}} $. Inserting above inequality to (5.42), one can use the discrete Cauchy-Schwartz inequality to obtain

$ \begin{align*} \|\phi_h^n\|^2-\|\phi_h^0\|^2\leq \frac{8\mathcal{Q}_\delta^2}{\varepsilon}\sum\limits_{k=1}^n\sum\limits_{j=1}^k\theta_{k-j}^{(k)}\langle \phi_h^j, \phi_h^k\rangle\leq \frac{8\mathcal{Q}_\delta^2}{\varepsilon}\sum\limits_{k=1}^n\sum\limits_{j=1}^k\theta_{k-j}^{(k)}\|\phi_h^j\|\phi_h^k\|. \end{align*} $

Select a $ n_0 $ such that $ \|\phi_h^{n_0}\|=\max_{0\leq k\leq n} \|\phi_h^k\| $. Then for the { time level $ n_0 $}, we have

$ \begin{align*} \|\phi_h^{n_0}\|^2\leq \|\phi_h^0\|\|\phi_h^{n_0}\|+ \frac{8\mathcal{Q}_\delta^2}{\varepsilon}\|\phi_h^{n_0}\|\sum\limits_{k=1}^{n_0}\|\phi_h^k\|\sum\limits_{j=1}^k\theta_{k-j}^{(k)}, \end{align*} $

which together with Lemma 4.2 imply that

$ \begin{align*} \|\phi_h^{n}\|\leq \|\phi_h^{n_0}\|\leq \|\phi_h^0\|+ \frac{8\mathcal{Q}_\delta^2}{\varepsilon}\|\sum\limits_{k=1}^{n_0}\tau_k\|\phi_h^k\|\leq\|\phi_h^0\|+ \frac{8\mathcal{Q}_\delta^2}{\varepsilon}\|\sum\limits_{k=1}^{n}\tau_k\|\phi_h^k\|. \end{align*} $

With the help of the maximum time step $ \tau\leq 16\mathcal{Q}_\delta^2/{\varepsilon} $, we arrive at

$ \begin{align*} \|\phi_h^{n}\|\|\leq 2\|\phi_h^0\|+ \frac{16\mathcal{Q}_\delta^2}{\varepsilon}\sum\limits_{k=1}^{n}\tau_k\|\phi_h^k\|. \end{align*} $

Thus, the proof is completed by using the Grönwall inequality in Lemma 5.1.

Let $ e^n_h:= u(t_n, \mathit{\boldsymbol{x}}_h)-u^n_h\; (n\geq 0) $ be the error between the numerical solution $ u^n_h $ and exact solution $ u({\mathit{\boldsymbol{x}}}_h, t_n) $ for $ {\mathit{\boldsymbol{x}}}_h\in \bar{\Omega}_h $. Then $ e_h^n $ solves the error function

$ \begin{align} \mathcal{D}_2 e_h^j +\varepsilon \Delta_h^2 e_h^j+\nabla_h\cdot (\mathit{\boldsymbol{f}}(\nabla_h (e^j_h+u_h^j))-\mathit{\boldsymbol{f}}(\nabla_h u_h^j))= \xi_h^j+\eta_h^j, \quad \rm{for} \quad 1\leq j\leq N, \end{align} $

(5.45)

where $ \eta_h^n:=\mathcal{D}_2 u({\mathit{\boldsymbol{x}}}_h, t_j)- u_{t}({\mathit{\boldsymbol{x}}}_h, t_j) \; (1\leq j\leq N) $ denotes the temporal truncation error and $ \xi_h^j:=C_uh^2 $ denotes the spatial truncation error. Here $ C_u $ is a constant depending on the exact solution $ u $, but independent of the mesh sizes $ h $ and $ \tau_j $.

The following Lemma implies that the temporal truncation error $ \eta_h^j $ can be divided into a {convolution part and a rest part}, which plays a key role in simplifying the complicate calculation for the estimate of convergence order when using the techniques of DOC and DCC kernels.

Lemma 5.3   ([29, Lemma 3.2])Denote the temporal truncation error by $ \eta_h^j:=\mathcal{D}_2 u(t_j)- u_{t}(t_j) \;(1\leq j\leq N) $ and set

$ \begin{equation} \begin{aligned} G^l &:= -\frac12\int_{t_{l-1}}^{t_l} (t-t_{l-1})^2 u_{ttt} \, \mathrm{d} t, \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 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} t u_{tt} \, \mathrm{d} t. \end{aligned} \end{equation} $

(5.46)

Then it holds that

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

(5.47)

Theorem 5.2  Let $ u(t, {\mathit{\boldsymbol{x}}}) $ and $ u^n_h $ be the solutions to problem (1.1) and discrete problem (2.8), respectively. Assume A1 holds and the maximum time step $ \tau\leq \varepsilon/16\mathcal{Q}_\delta $, then the discrete solution $ u_h^n $ of BDF2 scheme (2.8) is convergent in $ L^2 $-norm in the sense that

$ \begin{align*} \|e_h^{n}\|\leq &2\exp(16\mathcal{Q}_\delta^2t_{n-1}/\varepsilon)(\|e_h^0\|+2C_ut_nh^2+\sum\limits_{k=1}^n \tau_k^2\int_{t_{k-1}}^{t_k}\| u_{ttt}\| \, \mathrm{d} t\nonumber\\ &+2t_n\max\limits_{1\leq k\leq n }{\tau_k}\int_{t_{k-1}}^{t_k}\| u_{ttt}\| \, \mathrm{d} t+4\tau\int^{t_1}_{0} \| u_{tt}\| \, \mathrm{d} t). \end{align*} $

Proof  Multiplying both sides of (5.45) by $ \theta_{k-j}^{(k)} $, and summing $ j $ from 1 to $ k $, then taking the inner product with $ e_h^k $ on both sides and summing $ k $ from 1 to $ n $, one has

$ \begin{align} \|e_h^n\|^2-&\|e_h^0\|^2+2\varepsilon\sum\limits_{k=1}^n\sum\limits_{j=1}^k\theta_{k-j}^{(k)}\langle\Delta_h e_h^j, \Delta_h e_h^k\rangle\\ &\leq 2\sum\limits_{k=1}^n\sum\limits_{j=1}^k\theta_{k-j}^{(k)}\langle {\mathit{\boldsymbol{f}}}(\nabla_h(e_h^j+{u}_h^j))- {\mathit{\boldsymbol{f}}}(\nabla_h{u}_h^j), \nabla_he_h^k\rangle+2\sum\limits_{k=1}^n\sum\limits_{j=1}^k\theta_{k-j}^{(k)}\langle \xi_h^j+\eta_h^j, e_h^k\rangle \\ &:=I_1+I_2+I_3, \end{align} $

(5.48)

where the discrete Green's formula (3.9) and $ 2a^2-2ab\geq a^2-b^2 $ are used. The first term $ I_1 $ has been estimated by (5.44). One can use Lemma 4.2 and the discrete Cauchy–Schwartz inequality to obtain

$ \begin{align} I_2\leq 2t_n\max\limits_{1\leq k\leq n}\|\xi_h^k\|\|e_h^k\|\leq 2C_ut_nh^2\max\limits_{1\leq k\leq n}\|e_h^k\|. \end{align} $

(5.49)

We now consider $ I_3 $. By using Lemma 5.3 and exchanging the order of summation, we have

$ \begin{align*} I_3&=2\sum\limits_{k=1}^n\langle\sum\limits_{j=1}^k\theta_{k-j}^{(k)}\sum\limits_{l=1}^{j}b^{(j)}_{j-l}G^l, e_h^k\rangle +2\sum\limits_{k=1}^n\langle \sum\limits_{j=1}^k\theta_{k-j}^{(k)}R^j, e_h^k\rangle\\ &=2\sum\limits_{k=1}^n\langle\sum\limits_{l=1}^kG^l\sum\limits_{j=l}^{k}\theta_{k-j}^{(k)}b^{(j)}_{j-l}, e_h^k\rangle +2\sum\limits_{k=1}^n\langle\sum\limits_{j=1}^k\theta_{k-j}^{(k)}R^j, e_h^k\rangle\\ &=2\sum\limits_{k=1}^n\langle G^k, e_h^k\rangle+2\sum\limits_{k=1}^n\langle\sum\limits_{j=1}^k\theta_{k-j}^{(k)}R^j, e_h^k\rangle, \\ &\leq 2\sum\limits_{k=1}^n \|G^k\|\|e_h^k\|+2\sum\limits_{k=1}^n\sum\limits_{j=1}^k\theta_{k-j}^{(k)}\|R^j\|\|e_h^k\|. \end{align*} $

Inserting (5.44), (5.49) and above inequality into (5.48), one has

$ \begin{align*} \|e_h^n\|^2\leq& \|e_h^0\|^2+\frac{8\mathcal{Q}_\delta^2}{\varepsilon}\sum\limits_{k=1}^n\sum\limits_{j=1}^k\theta_{k-j}^{(k)}\|e_h^j\|\|e_h^k\|\\ &+2C_ut_nh^2\max\limits_{1\leq k\leq n}\|e_h^k\|+2\sum\limits_{k=1}^n \|G^k\|\|e_h^k\|+2\sum\limits_{k=1}^n\sum\limits_{j=1}^k\theta_{k-j}^{(k)}\|R^j\|\|e_h^k\|. \end{align*} $

Choose $ n_0 $ such that $ \|e_h^{n_0}\|=\max_{1\leq k\leq n} \|e_h^n\| $. Then, we have

$ \begin{align*} \|e_h^{n_0}\|^2\leq& \|e_h^0\|\|e_h^{n_0}\|+\frac{8\mathcal{Q}_\delta^2}{\varepsilon}\|e_h^{n_0}\|\sum\limits_{k=1}^{n_0}\sum\limits_{j=1}^k\theta_{k-j}^{(k)}\|e_h^k\|\\ &+2C_ut_nh^2\|e_h^{n_0}\|+2\|e_h^{n_0}\|\sum\limits_{k=1}^{n_0} \|G^k\|+2\|e_h^{n_0}\|\sum\limits_{k=1}^{n_0}\sum\limits_{j=1}^k\theta_{k-j}^{(k)}\|R^j\|. \end{align*} $

Thus, by exchanging the order of summation, we arrive at

$ \begin{align*} \|e_h^{n}\|&\leq \|e_h^{n_0}\|\leq \|e_h^0\|+\frac{8\mathcal{Q}_\delta^2}{\varepsilon}\sum\limits_{k=1}^{n}\tau_k\|e_h^k\|\\ &+2C_ut_nh^2+2\sum\limits_{k=1}^n \|G^k\|+2\sum\limits_{j=1}^np_{n-j}^{(k)}\|R^j\|. \end{align*} $

One may use Lemma 4.2 and the assumption of the maximum time step $ \tau\leq \varepsilon/16\mathcal{Q}_\delta $ to obtain

$ \begin{align*} \|e_h^{n}\|&\leq 2\|e_h^0\|+\frac{16\mathcal{Q}_\delta^2}{\varepsilon}\sum\limits_{k=1}^{n-1}\tau_k\|e_h^k\|\\ &+4C_ut_nh^2+4\sum\limits_{k=1}^n \|G^k\|+4\sum\limits_{j=1}^np_{n-j}^{(k)}\|R^j\|, \end{align*} $

which together with Lemma 4.2 and the discrete Grönwall inequality in Lemma 5.1 derives

$ \begin{align} \|e_h^{n}\|\leq 2\exp(16\mathcal{Q}_\delta^2t_{n-1}/\varepsilon)(\|e_h^0\|+2C_ut_nh^2+2\sum\limits_{k=1}^n \|G^k\|+2\sum\limits_{k=1}^np^{(n)}_{n-k}\|R^k\|). \end{align} $

(5.50)

Note that $ G^l, R^j $ defined in Lemma 5.3 satisfy

$ \begin{align*} \|G^l\| &\leq \frac{\tau_l^2}{2}\int_{t_{l-1}}^{t_l}\| u_{ttt}\| \, \mathrm{d} t, \\ \|R^j\|&\leq \frac{r_j}{2+2r_j}(2\tau_j+\tau_{j-1})\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*} $

Then it follows from Proposition 4.1 that

$ \begin{aligned}\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}\|\xi_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{aligned} $

Thus, we derive that

$ \begin{align} \|e_h^{n}\|\leq &2\exp(16\mathcal{Q}_\delta^2t_{n-1}/\varepsilon)(\|e_h^0\|+2C_ut_nh^2+\sum\limits_{k=1}^n \tau_k^2\int_{t_{k-1}}^{t_k}\| u_{ttt}\| \, \mathrm{d} t\\ &+2t_n\max\limits_{1\leq k\leq n}{\tau_k}\int_{t_{k-1}}^{t_k}\| u_{ttt}\| \, \mathrm{d} t+4\tau\int^{t_1}_{0} \| u_{tt}\| \, \mathrm{d} t). \end{align} $

(5.51)

The proof is complete.

Remark 2  Recently, [19] gives a result for the error estimate of the molecular beam epitaxial model (1.1) in form of

$ \begin{align} \|e_h^{n}\|\leq &C_u\exp(16\mathcal{M}_r^2t_{n-1}/\varepsilon)(\|e_h^0\|+t_n \tau_1+t_n(\tau^2+h^2)). \end{align} $

(5.52)

with $ 0\leq r_k \leq 3.561 $. One can see that the right-hand-side second term is the first-order convergence for large $ t_n $. To obtain the second-order convergence in [19], it requires another restriction condition $ |\mathfrak{R}_p| \leq N_0 \ll N $, where the index set $ |\mathfrak{R}_p| $ defined by (1.3). Our result in Theorem 5.2 shows the robust second-order (optimal) convergence remains valid to a new ratio condition $ 0<r_k\leq r_{\max}\approx 4.8645-\delta $, where $ \delta>0 $ is a arbitrarily small constant. The robustness means that the convergence does not need other conditions on the time step like the constrained condition $ |\mathfrak{R}_p| \leq N_0 \ll N $.

We now present numerical example to investigate the convergence order of BDF2 with variable time steps, which is tested on random time meshes. To do so, we set the computational domain $ \Omega = (0, S)^2 $ with $ S $ a positive constant, the final time $ T=1 $, and consider the following exterior-forced MBE model

$ \begin{align} u_t = -\varepsilon \Delta^2 u-\nabla\cdot \mathit{\boldsymbol{f}}(\nabla u)+g(\mathit{\boldsymbol{x}}, t). \end{align} $

(6.53)

By choosing a suitable function $ g $, the exact solution to (6.53) is constructed as follows

$ \begin{align*} u = \cos(t)\sin(\frac{2\pi x}{S})\sin(\frac{2\pi y}{S}). \end{align*} $

The random time meshes are given by $ \tau_k=T\chi_k/C $, where $ C = \sum_{k=1}^N\chi_k $ with $ \chi_k $ randomly drawn from the uniform distribution on $ (0, 1) $ and the number of spatial meshes is chosen by $ M=N $. In each run, the error $ e(N)= \|u(T)-u^N\| $ and the numerical rate of convergence at the final time $ T=1 $ are recorded in Tables 1 and 2, in which the maximum time step $ \tau $ and maximum adjacent time step ratio are also listed, where the convergence rate is calculated by

$ \begin{align*} \rm{Order}= \log_2(e(N)/e(2N)). \end{align*} $

As shown in Tables 1 and 2, even though the time step is randomly chosen beyond the constrained condition A1, the BDF2 scheme with variable time steps is robustly stable and convergent in the second order. In addition for the simulations, the first-step BDF1 does not bring the loss of accuracy and is consistent with our theoretical analysis.

The BDF2 scheme with variable time steps is considered to solve the MBE model without slope selection. Our proposed BDF2 scheme with BDF1 for the first step is proved to preserve the discrete energy dissipation law under a new adjacent time-step ratio condition: $ r_k:=\tau_k/\tau_{k-1}\leq4.8645-\delta $, where $ \delta>0 $ is a arbitrarily small constant. By using the techniques of DOC and DCC kernels, we achieve the robust and sharp second-order error estimate under the condition A1.

Our second-order convergence analysis shows two folds: (i) the variable time steps under condition A1 will not influence on the convergence order; (ii) the BDF1 scheme for the first step is consistent to the globally optimal convergence order. This conclusion removes the doubt of the classical choice of the first level solution with first-order consistent BDF1 scheme for the sharp second-order convergence. Numerical results demonstrate the theoretical analysis.

The Proof of Lemma 4.3

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

$ \begin{align*} 2w_k \sum\limits_{j= 1}^k \tilde{b}^{(k)}_{k-j} w_j &= 2\tilde{b}^{(k)}_0 w_k^2 + 2\tilde{b}_1^{(k)} w_k w_{k-1} \geq (2\tilde{b}^{(k)}_0 + \tilde{b}_1^{(k)}) w_k^2 +\tilde{b}_1^{(k)} w_{k-1}^2\\ & = \frac{2+4r_k - r_k^{3/2}}{(1+r_k)}w_k^2 - \frac{r_k^{3/2}}{(1+r_k)}w_{k-1}^2 \\ &= \frac{2+4 r_k - r_k^{3/2}}{(1+r_k)}w_k^2 - \frac{r_k^{3/2}}{ (1+r_k)}w_{k-1}^2 \\ & = \frac{r_{k+1}^{3/2}}{(1+r_{k+1})}w_k^2- \frac{r_k^{3/2}}{ (1+r_k)}w_{k-1}^2 +\tilde{R}(r_k, r_{k+1})w_k^2, \quad k \geq 2. \end{align*} $

Note that $ \partial_x \tilde{R}(x, y)=\frac12 (1+x)^{-2}(1-\sqrt{x})(x+\sqrt{x}+4). $ Hence, $ \tilde{R}(x, y) $ is increasing in $ (0, 1) $ and decreasing in $ (1, r_{\max}) $ with respect to $ x $. And it is easy to verify that $ \tilde{R}(x, y) $ is decreasing in $ (0, r_{\max}) $ with respect to $ y $. Thus, $ \tilde{R}(x, y) $ attains its minimum at $ (x, y)=(0, r_{\max}-\delta) $ or $ (x, y) =(r_{\max}, r_{\max}-\delta) $, namely,

$ \begin{align*} \min\limits_{0<x, y\leq r_{\max}-\delta} \tilde{R}(x, y) &= \min\{\tilde{R}(0, r_{\max}-\delta), \tilde{R}(r_{\max}-\delta, r_{\max}-\delta)\}= C_\delta>0. \end{align*} $

Due to the symmetry of matrix $ \tilde{B} $, for any $ \mathit{\boldsymbol{v}}\in \mathbb{R}^n $, we have

$ \begin{align*} \lambda_{\min}(\tilde{B}) = \min\limits_{\mathit{\boldsymbol{v}}\in \mathbb{R}^n}\frac{\mathit{\boldsymbol{v}}^T\tilde{B}\mathit{\boldsymbol{v}}}{\mathit{\boldsymbol{v}}^T\mathit{\boldsymbol{v}}}= \min\limits_{\mathit{\boldsymbol{v}}\in \mathbb{R}^n}\frac2{\mathit{\boldsymbol{v}}^T\mathit{\boldsymbol{v}}}\sum\limits_{k=1}^nw_k \sum\limits_{j= 1}^k \tilde{b}^{(k)}_{k-j} w_j\geq \min\limits_{1\leq k\leq n} \tilde{R}(r_k, r_{k+1})\geq C_\delta, \end{align*} $

which implies that the real symmetric matrix $ \tilde{B} $ is positive definite. The last claim holds by applying the standard Cholesky decomposition to $ \tilde{B} $. The proof is complete.

The proof of Lemma 4.4:

The direct calculation from the definition of $ \tilde{B}_2 $ (4.34) produces that

$ \begin{align*} \tilde{B}_2^T\tilde{B}_2=\left(\begin{array}{ccccc} d_0^{(1)}&d_1^{(2)}&&&\\ d_1^{(2)}&d_0^{(2)}&d_1^{(3)}\\ &\ddots&\ddots&\ddots&\\ && d_1^{(n-1)}&d_0^{(n-1)}&d_1^{(n)}\\ &&&d_1^{(n)}&d_0^{(n)} \end{array}\right)_{n\times n} \; , \end{align*} $

where the elements $ d_0^{(k)} $ and $ d_1^{(k)} $ are given by (set $ r_1\equiv 0 $)

$ \begin{align*} d_0^{(n)}&=(\frac{1+2r_n}{1+r_n})^2, \quad d_0^{(k)}=(\frac{1+2r_k}{1+r_k})^2+\frac{r_{k+1}^3}{(1+r_{k+1})^2}, \quad 1\leq k\leq n-1\\ d_1^{(k)}&=-\frac{r_k^{3/2}(1+2r_k)}{(1+r_k)^2}, \quad 1\leq k\leq n. \end{align*} $

Using the Gerschgorin circle theorem, one can obtain the upper bound of the maximum eigenvalue of matrix $ \tilde{B}_2^T\tilde{B}_2 $ as

$ \begin{align} \lambda_{\max}(\tilde{B}_2^T\tilde{B}_2)&\leq \max\limits_{1\leq k\leq n-1}\{d_0^{(k)}-d_1^{(k)}-d_1^{(k+1)}, \; d_0^{(n)}-d_1^{(n)}\}\\ &= \max\limits_{1\leq k\leq n-1}\{\hat{R}(r_k, r_{k+1}), \hat{R}(r_n, 0)\}\leq \max\limits_{1\leq k\leq n}\hat{R}(r_k, r_{k+1}). \end{align} $

(7.54)

It is easy to verify that $ \hat{R}(x, y) $ is increasing with respect to $ x $ and $ y $. Hence, by the time-step ratio condition A1, we have

$ \begin{align*} \lambda_{\max}(\tilde{B}_2^T\tilde{B}_2)&\leq \max\limits_{1\leq k\leq n}\hat{R}(r_k, r_{k+1})\leq \hat{R}(r_{\max}-\delta, r_{\max}-\delta)\\ &\leq \hat{R}(r_{\max}, r_{\max})= \frac{4r_{\max}^3}{(1+r_{\max})^2}\approx 13.3880<14. \end{align*} $

The proof is complete

The Proof of Lemma 4.5:

The first claim holds by a simple calculation that

$ \begin{align*} \mathit{\boldsymbol{\widehat{B}}}_2^{-1}\mathit{\boldsymbol{\widehat{B}}}(\mathit{\boldsymbol{\widehat{B}}}_2^{-1})^T=\mathit{\boldsymbol{\widehat{B}}}_2^{-1}(\mathit{\boldsymbol{\widehat{B}}}_2^T +\mathit{\boldsymbol{\widehat{B}}}_2)(\mathit{\boldsymbol{\widehat{B}}}_2^{-1})^T= (\mathit{\boldsymbol{\widehat{B}}}_2^{-1})^T+(\mathit{\boldsymbol{\widehat{B}}}_2^{-1}) = \mathit{\boldsymbol{\widehat{\Theta}}}. \end{align*} $

We now prove the second claim. It follows from Lemma 4.3 that

$ \begin{align*} \mathit{\boldsymbol{\widehat{B}}} = (\mathit{\boldsymbol{\widehat{L}}}\mathit{\boldsymbol{\widehat{\Lambda}}}_\tau^{-1})^T\mathit{\boldsymbol{\widehat{L}}}\mathit{\boldsymbol{\widehat{\Lambda}}}_\tau^{-1}, \end{align*} $

together with the Young's inequality, one has

$ \begin{align*} \sum\limits_{k=1}^n\sum\limits_{j=1}^k\theta^{(k)}_{k-j}(\mathit{\boldsymbol{v}}^k)^T\mathit{\boldsymbol{w}}^j&= {\mathit{\boldsymbol{w}}}^T\mathit{\boldsymbol{\widehat{\Theta}}}_2^T {\mathit{\boldsymbol{v}}} = {\mathit{\boldsymbol{w}}}^T(\mathit{\boldsymbol{\widehat{L}}}\mathit{\boldsymbol{\widehat{\Lambda}}}_\tau^{-1})^{-1}\mathit{\boldsymbol{\widehat{L}}}\mathit{\boldsymbol{\widehat{\Lambda}}}_\tau^{-1} \mathit{\boldsymbol{\widehat{\Theta}}}_2^T {\mathit{\boldsymbol{v}}}\\ &\leq \frac{\epsilon}{2} {\mathit{\boldsymbol{v}}}^T(\mathit{\boldsymbol{\widehat{L}}}\mathit{\boldsymbol{\widehat{\Lambda}}}_\tau^{-1} \mathit{\boldsymbol{\widehat{\Theta}}}_2^T)^T\mathit{\boldsymbol{\widehat{L}}}\mathit{\boldsymbol{\widehat{\Lambda}}}_\tau^{-1} \mathit{\boldsymbol{\widehat{\Theta}^T}}_2 {\mathit{\boldsymbol{v}}}+\frac1{2\epsilon} {\mathit{\boldsymbol{w}}}^T(\mathit{\boldsymbol{\widehat{L}}}\mathit{\boldsymbol{\widehat{\Lambda}}}_\tau^{-1})^{-1}(\mathit{\boldsymbol{\widehat{L}}} \mathit{\boldsymbol{\widehat{\Lambda}}}_\tau^{-1})^{-T} {\mathit{\boldsymbol{w}}}\\ &\leq \frac{\epsilon}{2} {\mathit{\boldsymbol{v}}}^T\mathit{\boldsymbol{\widehat{\Theta}}} {\mathit{\boldsymbol{v}}}+\frac1{2\epsilon} {\mathit{\boldsymbol{w}}}^T\mathit{\boldsymbol{\widehat{B}}}^{-1} {\mathit{\boldsymbol{w}}}. \end{align*} $

The proof is complete.

The Proof of Lemma 4.6:

Consider the gradient matrix of $ \mathit{\boldsymbol{f}} $ with respect to $ \mathit{\boldsymbol{v}}=(v_1, v_2) $, namely

$ \begin{align*} \nabla \mathit{\boldsymbol{f}}( {\mathit{\boldsymbol{v}}})=\frac1{1+|\mathit{\boldsymbol{v}}|^2}\left( \begin{array}{cc} 1-v_1^2+v_2^2&-2v_1v_2\\ -2v_1v_2&1-v_1^2+v_2^2 \end{array}\right). \end{align*} $

By calculating the eigenvalues of $ \nabla {\mathit{\boldsymbol{f}}} $, one has

$ \begin{align*} \mu_1(\mathit{\boldsymbol{v}})=\frac{1-|\mathit{\boldsymbol{v}}|^2}{(1+|\mathit{\boldsymbol{v}}|^2)^2}, \quad \mu_2(\mathit{\boldsymbol{v}})=\frac{1}{1+|\mathit{\boldsymbol{v}}|^2}, \end{align*} $

with $ -\frac18\leq \mu_1(\mathit{\boldsymbol{v}})\leq 1 $ and $ 0< \mu_2(\mathit{\boldsymbol{v}})\leq 1 $. The application of Taylor expansion yields

$ \begin{align*} {\mathit{\boldsymbol{f}}}(\mathit{\boldsymbol{v}})= {\mathit{\boldsymbol{f}}}(\mathit{\boldsymbol{w}})+\int_0^1\nabla {\mathit{\boldsymbol{f}}}(\theta \mathit{\boldsymbol{v}}+(1-\theta) \mathit{\boldsymbol{w}}) \, \mathrm{d} \theta(\mathit{\boldsymbol{v}}-\mathit{\boldsymbol{w}}). \end{align*} $

From the symmetry of $ \nabla {\mathit{\boldsymbol{f}}}(\mathit{\boldsymbol{v}}) $, there exists an orthogonal matrix $ \left( \begin{array}{cc} a_\theta&b_\theta\\ c_\theta&d_\theta \end{array}\right) $ such that

$ \begin{align*} \left( \begin{array}{cc} \lambda_1&\\&\lambda_2 \end{array}\right)=\int_0^1 \left( \begin{array}{cc} a_\theta&b_\theta\\ c_\theta&d_\theta \end{array}\right)\left( \begin{array}{cc} \mu_1(\mathit{\boldsymbol{\xi}}_\theta)&\\&\mu_2(\mathit{\boldsymbol{\xi}}_\theta) \end{array}\right)\left( \begin{array}{cc} a_\theta&c_\theta\\ b_\theta&d_\theta \end{array}\right) \, \mathrm{d} \theta, \end{align*} $

where $ \mathit{\boldsymbol{\xi}}_\theta = \theta {\mathit{\boldsymbol{v}}}+(1-\theta) {\mathit{\boldsymbol{w}}} $, and $ \lambda_1, \lambda_2 $ are the eigenvalues of matrix $ \int_0^1\nabla {\mathit{\boldsymbol{f}}}(\mathit{\boldsymbol{\xi}}_\theta) \, \mathrm{d} \theta $, and $ \mu_1(\mathit{\boldsymbol{\xi}}_\theta), \mu_2(\mathit{\boldsymbol{\xi}}_\theta) $ are the eigenvalues of matrix $ \nabla {\mathit{\boldsymbol{f}}}(\mathit{\boldsymbol{\xi}}_\theta) $. Without loss of generality, here we assume $ \mu_1(\mathit{\boldsymbol{\xi}}_\theta)\leq \mu_2(\mathit{\boldsymbol{\xi}}_\theta) $. It is easy to verify

$ \begin{align*} \lambda_1 = \int_0^1 (a_\theta^2\mu_1(\xi_\theta)+b_\theta^2 \mu_2(\xi_\theta)) \, \mathrm{d} \theta, \quad \lambda_2 = \int_0^1 (c_\theta^2\mu_1(\xi_\theta)+d_\theta^2 \mu_2(\xi_\theta)) \, \mathrm{d} \theta. \end{align*} $

The orthogonality of matrix $ \left( \begin{array}{cc} a_\theta&b_\theta\\ c_\theta&d_\theta \end{array}\right) $ yields $ a_\theta^2+b_\theta^2=1, c_\theta^2+d_\theta^2=1 $, which implies $ \mu_1(\xi_\theta)\leq \lambda_1, \lambda_2\leq \mu_2(\xi_\theta) $. Hence, we have $ Q_ {\mathit{\boldsymbol{f}}}=\int_0^1\nabla {\mathit{\boldsymbol{f}}}(\theta \mathit{\boldsymbol{v}}+(1-\theta) \mathit{\boldsymbol{w}}) \, \mathrm{d} \theta $ and its eigenvalues $ \lambda_1, \lambda_2 $ satisfy $ -\frac18\leq \lambda_1, \lambda_2\leq 1 $, which implies $ \|Q_ {\mathit{\boldsymbol{f}}}\|_2\leq 1 $. By using the property of matrix norm that $ |A\mathit{\boldsymbol{x}}|\leq \|A\|_2|\mathit{\boldsymbol{x}}| $, one has

$ \begin{align*} | {\mathit{\boldsymbol{f}}}(\mathit{\boldsymbol{v}})- {\mathit{\boldsymbol{f}}}(\mathit{\boldsymbol{w}})|\leq\|\nabla {\mathit{\boldsymbol{f}}}(\theta \mathit{\boldsymbol{v}}+(1-\theta) \mathit{\boldsymbol{w}})\|_2|\mathit{\boldsymbol{v}}-\mathit{\boldsymbol{w}}|\leq |\mathit{\boldsymbol{v}}-\mathit{\boldsymbol{w}}|. \end{align*} $

The proof is complete.

The Proof of Lemma 4.7:

From Lemma 4.6, there exist symmetric matrix sequences $ Q_f^j \in \mathbb{R}^{2\times 2} $ such that

$ \begin{align*} {\mathit{\boldsymbol{f}}}( {\mathit{\boldsymbol{v}}}^j+ {\mathit{\boldsymbol{w}}}^j)- {\mathit{\boldsymbol{f}}}( {\mathit{\boldsymbol{v}}}^j)=Q_f^j {\mathit{\boldsymbol{w}}}^j, \quad 1\leq j\leq n, \end{align*} $

where the eigenvalues $ \mu^j_1, \mu_2^j $ of $ Q_f^j $ satisfy $ -\frac18\leq \mu^j_1, \mu_2^j\leq 1 $.

Define the symmetric matrix

$ \begin{align*} Q:=\rm{diag}(Q^1_f, \cdots, Q^n_f)\in \mathbb{R}^{2n\times 2n}. \end{align*} $

The spectral radius of the symmetric matrix $ Q $ satisfies $ \rho(Q)\leq 1 $, which implies that $ \|Q\|_2\leq1. $ According to Lemma 4.5, one has

$ \begin{align} \sum\limits_{k=1}^n\sum\limits_{j=1}^k &\theta_{k-j}^{(k)}(\mathit{\boldsymbol{z}}^k)^T[ {\mathit{\boldsymbol{f}}}( {\mathit{\boldsymbol{v}}}^j+ {\mathit{\boldsymbol{w}}}^j)- {\mathit{\boldsymbol{f}}}( {\mathit{\boldsymbol{v}}}^j)]=\sum\limits_{k=1}^n\sum\limits_{j=1}^k \theta_{k-j}^{(k)}(\mathit{\boldsymbol{z}}^k)^T Q_f^j {\mathit{\boldsymbol{w}}}^j \\ &\leq \frac{\epsilon}{2}\mathit{\boldsymbol{z}}^T\mathit{\boldsymbol{\widehat{\Theta}}}\mathit{\boldsymbol{z}}+\frac1{2\epsilon} {\mathit{\boldsymbol{w}}}^TQ^T\mathit{\boldsymbol{\widehat{B}}}^{-1}Q {\mathit{\boldsymbol{w}}} \\ &=\epsilon \sum\limits_{k=1}^n\sum\limits_{j=1}^k \theta_{k-j}^{(k)}(\mathit{\boldsymbol{z}}^k)^T\mathit{\boldsymbol{z}}^j+\frac1{2\epsilon} {\mathit{\boldsymbol{w}}}^TQ^T\mathit{\boldsymbol{\widehat{B}}}^{-1}Q {\mathit{\boldsymbol{w}}}, \end{align} $

(7.55)

where $ \mathit{\boldsymbol{z}} = ((\mathit{\boldsymbol{z}}^1)^T, \cdots, (\mathit{\boldsymbol{z}}^n)^T)^T $ and $ \mathit{\boldsymbol{w}} = ((\mathit{\boldsymbol{w}}^1)^T, \cdots, (\mathit{\boldsymbol{w}}^n)^T)^T $. It follows from Lemmas 4.3 and 4.5 that

$ \begin{align*} \mathit{\boldsymbol{\widehat{\Theta}}} =(\mathit{\boldsymbol{\widehat{B}}}_2^{-1})^T\mathit{\boldsymbol{\widehat{B}}}\mathit{\boldsymbol{\widehat{B}}}_2^{-1}=(\mathit{\boldsymbol{\widehat{L}}}\mathit{\boldsymbol{\widehat{\Lambda}}}_\tau^{-1} \mathit{\boldsymbol{\widehat{B}}}_2^{-1})^T\mathit{\boldsymbol{\widehat{L}}}\mathit{\boldsymbol{\widehat{\Lambda}}}_\tau^{-1}\mathit{\boldsymbol{\widehat{B}}}_2^{-1}, \end{align*} $

which implies

$ \begin{align*} {\mathit{\boldsymbol{w}}}^T\mathit{\boldsymbol{\widehat{\Theta}}} {\mathit{\boldsymbol{w}}}=\|\mathit{\boldsymbol{\widehat{L}}}\mathit{\boldsymbol{\widehat{\Lambda}}}_\tau^{-1}\mathit{\boldsymbol{\widehat{B}}}_2^{-1} {\mathit{\boldsymbol{w}}}\|_2^2. \end{align*} $

Then the right term of the last identity in (7.55) can be estimated by

$ \begin{align*} \mathit{\boldsymbol{w}}^T\mathit{\boldsymbol{\widehat{Q}}}^T\mathit{\boldsymbol{\widehat{B}}}^{-1}Q\mathit{\boldsymbol{w}}&=\|(\mathit{\boldsymbol{\widehat{L}}}^{-1})^T\mathit{\boldsymbol{\widehat{\Lambda}}}_\tau \mathit{\boldsymbol{\widehat{Q}}}\mathit{\boldsymbol{w}}\|_2^2\\ &=\|(\mathit{\boldsymbol{\widehat{L}}}^{-1})^T\mathit{\boldsymbol{\widehat{\Lambda}}}_\tau Q\mathit{\boldsymbol{\widehat{B}}}_2\mathit{\boldsymbol{\widehat{\Lambda}}}_\tau \mathit{\boldsymbol{\widehat{L}}}^{-1}\mathit{\boldsymbol{\widehat{L}}}\mathit{\boldsymbol{\widehat{\Lambda}}}_\tau^{-1}\mathit{\boldsymbol{\widehat{B}}}_2^{-1}\mathit{\boldsymbol{w}}\|_2^2\\ &\leq \|(\mathit{\boldsymbol{\widehat{L}}}^{-1})^T\mathit{\boldsymbol{\widehat{\Lambda}}}_\tau Q\mathit{\boldsymbol{\widehat{B}}}_2\mathit{\boldsymbol{\widehat{\Lambda}}}_\tau \mathit{\boldsymbol{\widehat{L}}}^{-1}\|^2_2 {\mathit{\boldsymbol{w}}}^T\mathit{\boldsymbol{\widehat{\Theta}}} {\mathit{\boldsymbol{w}}}\\ &\leq \| Q\|_2^2\|\mathit{\boldsymbol{\widetilde{B}}}_2\|_2^2\|\mathit{\boldsymbol{\widehat{L}}}^{-1}\|_2^4 {\mathit{\boldsymbol{w}}}^T\mathit{\boldsymbol{\widehat{\Theta}}} {\mathit{\boldsymbol{w}}}\leq 2\mathcal{Q}_\delta \sum\limits_{k=1}^n\sum\limits_{j=1}^k \theta_{k-j}^{(k)}(\mathit{\boldsymbol{w}}^j)^T\mathit{\boldsymbol{w}}^k, \end{align*} $

where the commutativity $ \mathit{\boldsymbol{\widehat{\Lambda}}}_\tau Q=Q \mathit{\boldsymbol{\widehat{\Lambda}}}_\tau $ is used and $ \mathcal{Q}_\delta $ is defined in (4.37). Inserting the above estimate to (7.55), one can obtain the first claim. The second claim can be proved by taking $ \epsilon=\sqrt{\mathcal{Q}_\delta} $ in the first claim. The proof is complete.