#  On a new class of BDF and IMEX schemes for parabolic type equations

Fukeng Huang† and Jie Shen‡

###### Abstract.

When applying the classical multistep schemes for solving differential
equations, one often faces the dilemma that smaller time steps are needed with
higher-order schemes, making it impractical to use high-order schemes for
stiff problems. We construct in this paper a new class of BDF and implicit-
explicit (IMEX) schemes for parabolic type equations based on the Taylor
expansions at time $t^{n+\beta}$ with $\beta>1$ being a tunable parameter.
These new schemes, with a suitable $\beta$, allow larger time steps at higher-
order for stiff problems than that is allowed with a usual higher-order
scheme. For parabolic type equations, we identify an explicit uniform
multiplier for the new second- to fourth-order schemes, and conduct rigorously
stability and error analysis by using the energy argument. We also present
ample numerical examples to validate our findings.

###### Key words and phrases:

stability; error analysis; implicit-explicit schemes; parabolic systems

###### 2000 Mathematics Subject Classification:

65M12; 76D05; 65M15

† Department of Mathematics, National University of Singapore, Singapore,
119076 (hfkeng@nus.edu.sg).
‡ School of Mathematical Science, Eastern Institute of Technology, Ningbo,
315200, China (jshen@eitech.edu.cn). This work is supported in part by NSFC
grant 11971407.

##  1\. Introduction

We consider in this paper numerical methods of a class of nonlinear ordinary
or partial differential equations in the form

$\begin{split}&u_{t}+\mathcal{L}u(t)+\mathcal{G}[u(t)]=f(t),\quad 0<t<T,\\\
&u(0)=u^{0},\end{split}$ (1.1)

where $\mathcal{L}$ is a linear (or possibly nonlinear) positive operator and
$\mathcal{G}$ is a nonlinear operator, whose exact descriptions can be found
in the next section.

Numerical approximation of ordinary differential equations (ODEs) is a very
mature field (see, for instance, [9, 10, 15, 18]), and the numerical methods
developed for ODEs have been playing important roles in solving partial
differential equations (PDEs) in the form of (1.1) through the method of lines
[27], or the so called method of lines transpose [20], i.e., discretizing
first in time followed by the discretization in space. In particular, the
backward difference formulae (BDF) and the implicit-explicit (IMEX) schemes
are frequently used to deal with (1.1) which exhibit stiff behaviors [8, 16,
21].

Two key issues of numerical methods for (1.1) are stability and accuracy. In
order to obtain highly accurate solution with less computational costs, it is
highly desirable to be able to use higher-order schemes with larger time
steps. However, as we increase the order of accuracy of BDF or IMEX type
schemes, their stability regions usually decrease, i.e., smaller time steps
need to be used with higher-order schemes, particularly for stiff problems,
making high-order schemes impractical for many complex nonlinear systems. A
natural question arises: is it possible to develop higher-order multi-step
schemes such that their stability regions are comparable or even larger than
lower-order classical BDF or IMEX schemes?

The main purposes of this paper are two-fold:

  * •

to construct a new class of BDF and IMEX schemes with a tunable parameter such
that larger time steps can be used in higher-order schemes;

  * •

to carry out a rigorous stability and error analysis for this new class of
IMEX schemes.

Furthermore, we provide convincing numerical evidences to validate our
theoretical findings.

We recall that the classical BDF and IMEX schemes for approximating solution
at time $t^{n+1}$ are usually constructed using the Taylor expansion formulae
at time $t^{n+\beta}$ with $\beta\in\\{0,1\\}$. In this paper, we shall
construct a new class of BDF and IMEX schemes based on the Taylor expansion
formulae at time $t^{n+\beta}$ with $\beta\geq 1$ being a tunable parameter.
The new schemes are a simple generalization of the classical BDF or IMEX
schemes with essentially the same computational efforts. However, they enjoy a
remarkable property that their stability regions increase as the parameter
$\beta$ increases, making it possible, by choosing a suitably large $\beta$,
to use high-order schemes with reasonably larger time steps. The price to pay
with a larger $\beta$ is increased truncation errors which can be more than
compensated with higher-order of accuracy.

On the other hand, it is well known that a rigorous stability and error
analysis  by using the energy technique of the classical BDF (and the related
IMEX) schemes of order up to five (cf. [5, 6, 14, 22, 24]) relies on a result
by Nevanlinna and Odeh [25] (see also [3] for the extension to the six-order
BDF scheme) in which the existence of suitable multiplier that can lead to
energy stability was established. It is therefore natural to ask whether such
a multiplier exists for the new class of BDF schemes. We shall construct
explicitly suitable multipliers in a more general form for the new class of
BDF schemes of orders two to four, and derive explicit telescoping formulae
associated with these multipliers. Furthermore, for nonlinear parabolic type
equations, we show rigorously that the stability condition of the new class of
IMEX schemes becomes less restrictive as $\beta$ increases, particularly
compared with the classical case of $\beta=1$.

The idea behind the new class of BDF and IMEX schemes is very simple but
original, and can be easily extended to other type numerical schemes. However,
our stability and error analysis rely on the explicit formulae for the uniform
multipliers and telescoping decomposition whose derivations are totally
nontrivial and original. On the other hand, the new schemes can be easily
implemented with a minimal effort by modifying the code based on the classical
BDF or IMEX schemes, and provide a much needed improvement on the stability of
higher-order schemes.

The rest of the paper is organized as follows. In Section 2, we describe the
abstract setting and construct the new class of BDF and IMEX methods based on
the Taylor expansion at time $t^{n+\beta}$ and investigate their stability
regions. In Section 3, we identify an explicit and uniform multiplier for the
new class of BDF and IMEX schemes, which plays an essential role in the
stability and error analysis. In Section 4, we establish the unconditional
stability for the linear parabolic equations and the stability, followed by
error analysis for the nonlinear parabolic equations in Section 5. In section
6, we discuss extension to the fifth-order scheme. In section 7, we provide
numerical examples to show the advantages of our new schemes, followed by some
concluding remarks in section 8.

##  2\. A new class of BDF and IMEX schemes

###  2.1. The abstract setting

We first describe the functional setting. For the sake of simplicity, we
consider a simpler setting than that used in [6], although our analysis would
also work for the more general setting there.

Let $V$ and $H$ be two real Hilbert spaces such that $V\subset
H=H^{\prime}\subset V^{\prime}$, with $V$ densely and continuously embedded in
$H$ and $V^{\prime}$ being the dual space of $V$. We consider (1.1) with
$\mathcal{L}$: $V\rightarrow V^{\prime}$ being a positive definite, self-
adjoint, linear operator, and $f$ in $V^{\prime}$ is a given source term. We
denote the inner product in $H$ by $(\cdot,\cdot)$, and the induced norm in
$H$ by $|\cdot|$. We also denote the norm in $V$ by $\|\cdot\|$ which is
defined as $\|u\|:=|\mathcal{L}^{1/2}u|=(\mathcal{L}u,u)^{1/2}$. The dual norm
in $V^{\prime}$ is defined by

$\|v\|_{\star}:=\mathop{\rm{sup}}\limits_{u\in
V\backslash\\{0\\}}\frac{|(v,u)|}{\|u\|},\quad\forall v\in V^{\prime}.$ (2.1)

We assume that the nonlinear operator $\mathcal{G}$ satisfies the following
local Lipschitz condition [6] in a ball  $\mathcal{B}_{u(t)}:=\\{v\in
V:\|v-u(t)\|\leq 1\\}$, centered at the exact solution $u(t)$,

$\|\mathcal{G}(v)-\mathcal{G}(\tilde{v})\|^{2}_{\star}\leq\gamma\|v-\tilde{v}\|^{2}+\mu|v-\tilde{v}|^{2},\quad\forall
v,\tilde{v}\in\mathcal{B}_{u(t)},\quad\forall t\in[0,T],$ (2.2)

with a non-negative constant $\gamma$ and an arbitrary constant $\mu$.

###  2.2. Construction of the new schemes

We shall first construct the new schemes for (1.1) based on the Taylor
expansion at time $t^{n+\beta}$. Given an integer $k\geq 2$, denoting
$t^{n}=n\Delta t$, it follows from the Taylor expansion at time $t^{n+\beta}$
that

$\phi(t^{n+1-i})=\sum_{m=0}^{k-1}[(1-i-\beta)\Delta
t]^{m}\frac{\phi^{(m)}(t^{n+\beta})}{m!}+\mathcal{O}(\Delta t^{k}),\quad\text{
for }k\geq i\geq 0.$ (2.3)

Then we can derive from the above an implicit difference formula to
approximate $\partial_{t}\phi(t^{n+\beta})$:

$\frac{1}{\Delta
t}\sum_{q=0}^{k}a_{k,q}(\beta)\phi(t^{n+1-k+q})=\partial_{t}\phi(t^{n+\beta})+\mathcal{O}(\Delta
t^{k}),$ (2.4)

where $a_{k,q}(\beta)$ can be uniquely determined by solving the following
linear system with a Vandermonde matrix:

${\left[\begin{array}[]{ccccc}1&1&...&...&1\\\
\beta-1&\beta&...&...&\beta+k-1\\\
(\beta-1)^{2}&\beta^{2}&...&...&(\beta+k-1)^{2}\\\
\vdots&\vdots&\vdots&\vdots&\vdots\\\
(\beta-1)^{k}&\beta^{k}&...&...&(\beta+k-1)^{k}\end{array}\right]}{\left[\begin{array}[]{c}a_{k,k}(\beta)\\\
a_{k,k-1}(\beta)\\\ a_{k,k-2}(\beta)\\\ \vdots\\\
a_{k,0}(\beta)\end{array}\right]}={\left[\begin{array}[]{c}0\\\ -1\\\ 0\\\
\vdots\\\ 0\end{array}\right]}.$ (2.5)

Similarly, we can derive an implicit difference formula to approximate
$\phi(t^{n+\beta})$:

$\sum_{q=0}^{k-1}b_{k,q}(\beta)\phi(t^{n+2-k+q})=\phi(t^{n+\beta})+\mathcal{O}(\Delta
t^{k}),$ (2.6)

with $b_{k,q}(\beta)$ being the unique solution of the following Vandermonde
system:

${\left[\begin{array}[]{ccccc}1&1&...&...&1\\\
\beta-1&\beta&...&...&\beta+k-2\\\ \vdots&\vdots&\vdots&\vdots&\vdots\\\
(\beta-1)^{k-1}&\beta^{k-1}&...&...&(\beta+k-2)^{k-1}\end{array}\right]}{\left[\begin{array}[]{c}b_{k,k-1}(\beta)\\\
b_{k,k-2}(\beta)\\\ \vdots\\\
b_{k,0}(\beta)\end{array}\right]}={\left[\begin{array}[]{c}1\\\ 0\\\ \vdots\\\
0\end{array}\right]}.$ (2.7)

To deal with the nonlinear term in (1.1), we also need the following explicit
difference formula to approximate $\phi(t^{n+\beta})$:

$\sum_{q=0}^{k-1}c_{k,q}(\beta)\phi(t^{n+1-k+q})=\phi(t^{n+\beta})+\mathcal{O}(\Delta
t^{k}),$ (2.8)

where $c_{k,q}(\beta)$ can be uniquely determined from:

${\left[\begin{array}[]{ccccc}1&1&...&...&1\\\
\beta&\beta+1&...&...&\beta+k-1\\\ \vdots&\vdots&\vdots&\vdots&\vdots\\\
\beta^{k-1}&(\beta+1)^{k-1}&...&...&(\beta+k-1)^{k-1}\end{array}\right]}{\left[\begin{array}[]{c}c_{k,k-1}(\beta)\\\
c_{k,k-2}(\beta)\\\ \vdots\\\
c_{k,0}(\beta)\end{array}\right]}={\left[\begin{array}[]{c}1\\\ 0\\\ \vdots\\\
0\end{array}\right]}.$ (2.9)

Then, a new class of BDF schemes for (1.1) with $\mathcal{G}=0$ is

$\frac{1}{\Delta
t}\sum_{q=0}^{k}a_{k,q}(\beta)\phi^{n+1-k+q}+\mathcal{L}(\sum_{q=0}^{k-1}b_{k,q}(\beta)\phi^{n+2-k+q})=f(t^{n+\beta}),\quad
k\geq 2,$ (2.10)

and a new class of IMEX schemes for (1.1) is

$\frac{1}{\Delta
t}\sum_{q=0}^{k}a_{k,q}(\beta)\phi^{n+1-k+q}+\mathcal{L}(\sum_{q=0}^{k-1}b_{k,q}(\beta)\phi^{n+2-k+q})+\mathcal{G}(\sum_{q=0}^{k-1}c_{k,q}(\beta)\phi^{n+1-k+q})=f(t^{n+\beta}),\quad
k\geq 2.$ (2.11)

######  Remark 1.

When $\beta=1$, (2.11) (resp. (2.10)) becomes the classical semi-implicit IMEX
(resp. BDF) schemes, and there have been extensive works regarding its
stability and error analysis [2, 4, 6, 22, 23] in the literature.  For
$\forall\beta>1$, (2.10) and (2.11) still involve values at the same
$k+1$-levels as the classical one (with $\beta=1$) on the left hand side while
they involve values at time $t^{n+\beta}$ on the right hand side.

For the reader’s convenience, we list below the coefficients in (2.11) for
$k=2,3,4$.

$k=2$:

$\displaystyle a_{2,2}(\beta)=\frac{2\beta+1}{2},\quad
a_{2,1}(\beta)=-2\beta,\quad a_{2,0}(\beta)=\frac{2\beta-1}{2},$ (2.12a)
$\displaystyle b_{2,1}(\beta)=\beta,\quad b_{2,0}(\beta)=-(\beta-1),$ (2.12b)
$\displaystyle c_{2,1}(\beta)=\beta+1,\,c_{2,0}(\beta)=-\beta.$ (2.12c)

$k=3$:

$\displaystyle
a_{3,3}(\beta)=\frac{3\beta^{2}+6\beta+2}{6},\,a_{3,2}(\beta)=\frac{-(9\beta^{2}+12\beta-3)}{6},\,a_{3,1}(\beta)=\frac{9\beta^{2}+6\beta-6}{6},\,a_{3,0}(\beta)=\frac{-(3\beta^{2}-1)}{6},$
(2.13a) $\displaystyle
b_{3,2}(\beta)=\frac{\beta^{2}+\beta}{2},\,b_{3,1}(\beta)=-(\beta^{2}-1),\,b_{3,0}(\beta)=\frac{\beta^{2}-\beta}{2},$
(2.13b) $\displaystyle
c_{3,2}(\beta)=\frac{\beta^{2}+3\beta+2}{2},\,c_{3,1}(\beta)=-(\beta^{2}+2\beta),\,c_{3,0}(\beta)=\frac{\beta^{2}+\beta}{2}.$
(2.13c)

$k=4$:

$\displaystyle
a_{4,4}(\beta)=\frac{2\beta^{3}+9\beta^{2}+11\beta+3}{12},\,a_{4,3}(\beta)=\frac{-8\beta^{3}-30\beta^{2}-20\beta+10}{12},\,a_{4,2}(\beta)=\frac{12\beta^{3}+36\beta^{2}+6\beta-18}{12}$
$\displaystyle
a_{4,1}(\beta)=\frac{-8\beta^{3}-18\beta^{2}+4\beta+6}{12},\,a_{4,0}(\beta)=\frac{2\beta^{3}+3\beta^{2}-\beta-1}{12},$
(2.14a) $\displaystyle
b_{4,3}(\beta)=\frac{\beta^{3}+3\beta^{2}+2\beta}{6},\,b_{4,2}(\beta)=\frac{-\beta^{3}-2\beta^{2}+\beta+2}{2},\,b_{4,1}(\beta)=\frac{\beta^{3}+\beta^{2}-2\beta}{2},\,b_{4,0}(\beta)=\frac{-\beta^{3}+\beta}{6},$
(2.14b) $\displaystyle
c_{4,3}(\beta)=\frac{\beta^{3}+6\beta^{2}+11\beta+6}{6},\,c_{4,2}(\beta)=\frac{-\beta^{3}-5\beta^{2}-6\beta}{2},\,c_{4,1}(\beta)=\frac{\beta^{3}+4\beta^{2}+3\beta}{2},$
$\displaystyle c_{4,0}(\beta)=\frac{-\beta^{3}-3\beta^{2}-2\beta}{6}.$ (2.14c)

######  Remark 2.

Instead of deriving (2.11) from Taylor expansions, one may also derive it by
following the standard construction of the usual multistep methods using
interpolation formulae (see, e.g., Section 2 in [19]). In fact, it can be
shown that the coefficients $a_{k,q}(\beta),b_{k,q}(\beta),c_{k,q}(\beta)$ can
be determined by the values at $t^{n+\beta}$ of the corresponding Lagrange
polynomials and their derivatives. For example,

$a_{k,q}(\beta)=\Delta tL_{q}^{\prime}(t^{n+\beta}),\quad q=0,...,k,$ (2.15)

where $L_{q}$ is the Lagrange polynomials associated with
$t^{n+1-k},...,t^{n+1}$.

###  2.3. Linear stability regions

In this subsection, we investigate the regions of linear stability of the new
schemes (2.10). For the test equation $\phi_{t}=\lambda\phi$, (2.10) reduces
to

$\frac{1}{\Delta
t}\sum_{q=0}^{k}a_{k,q}(\beta)\phi^{n+1-k+q}=\lambda\sum_{q=0}^{k-1}b_{k,q}(\beta)\phi^{n+2-k+q},\quad
k\geq 2.$ (2.16)

In order to study the stability regions for $\beta\neq 1$, we set
$\phi^{n}=w^{n}$ (here, $``n"$ is an upper index in $\phi^{n}$ and an exponent
in $w^{n}$) and $z=\lambda\Delta t$ in (2.20) to obtain its characteristic
equation, e.g., in the case of $k=2$, it takes the form:

$(2\beta+1-2\beta z)w^{2}+(2(\beta-1)z-4\beta)w+(2\beta-1)=0.$ (2.17)

Then the region of absolute stability of method (2.20) is the set of all
$z\in\mathbb{C}$ such that the characteristic polynomial satisfies the root
condition. We recall that the second order case was already considered in
[17], and it was shown that the second-order case of (2.20) is A-stable for
$\beta\geq 1$, and more importantly, the stability regions increase as we
increase $\beta$. Hence, for the test equation $\phi_{t}=\lambda\phi$, a more
general second-order IMEX method is

$\frac{(2\beta+1)\phi^{n+1}-4\beta\phi^{n}+(2\beta-1)\phi^{n-1}}{2\Delta
t}=\lambda\beta\phi^{n+1}-\lambda(\beta-1)\phi^{n}.$ (2.18)

Note that with $\beta=1$, it reduces to the usual second-order IMEX scheme. In
order to study the stability regions for $\beta\neq 1$, we set
$\phi^{n}=w^{n}$ and $z=\lambda\Delta t$ in (2.18) to obtain its
characteristic polynomial

$(2\beta+1-2\beta z)w^{2}+(2(\beta-1)z-4\beta)w+(2\beta-1)=0.$ (2.19)

Then the region of absolute stability of method (2.18) is the set of all
$z\in\mathbb{C}$ such that (2.19) holds for all $|w|<1$. In Fig. 1, we plot
the stability regions of the general IMEX2 type method (2.18) for
$\beta=1,3,5$. We observe that the general IMEX2 method are A-stable for
$\beta\geq 1$, and more importantly, the stability regions increases as we
increases $\beta$.

(a) second order, $\beta=1$ (b) second order, $\beta=3$ (c) second order,
$\beta=5$

(d) $\beta=1$, zoom in around the origin (e) $\beta=3$, zoom in around the
origin (f) $\beta=5$, zoom in around the origin

Figure 1.  The pink parts show the regions of absolute stability of the
general second order IMEX type method with Taylor expansion at
$n+\beta,\,\beta=1,3,5$.

Similarly, we can obtain the coefficients $a_{k,q},\,b_{k,q}$ for $k=3,4$ by
solving the linear systems (2.5) and (2.7) as

Then for the test equation $\phi_{t}=\lambda\phi$, by performing the Taylor
expansions at $t^{n+\beta}$, a more general third- and fourth- order IMEX
method can be written as

$\frac{1}{\Delta
t}\sum_{q=0}^{k}a_{k,q}(\beta)\phi^{n+1-k+q}=\lambda\sum_{q=0}^{k-1}b_{k,q}(\beta)\phi^{n+2-k+q},\quad
k=3,4.$ (2.20)

Same as the second order case, we can set $\phi^{n}=w^{n}$ and
$z=\lambda\Delta t$ in (2.20) to obtain the corresponding characteristic
polynomials for the third and fourth order scheme.

In Fig. 2 and Fig. 3 , we plot the stability regions of the general third- and
fourth-order BDF schemes for $\beta=1,3,5$. We observe that, the stability
regions increase as we increase $\beta$.

(a) third order, $\beta=1$

(b) third order, $\beta=3$

(c) third order, $\beta=5$

(d) $\beta=1$, zoom in around the origin

(e) $\beta=3$, zoom in around the origin

(f) $\beta=5$, zoom in around the origin

Figure 2.  The pink parts show the region of absolute stability of the general
third order BDF scheme with Taylor expansion at $n+\beta,\,\beta=1,3,5$.

(a) fourth order, $\beta=1$

(b) fourth order, $\beta=3$

(c) fourth order, $\beta=5$

(d) $\beta=1$, zoom in around the origin

(e) $\beta=3$, zoom in around the origin

(f) $\beta=5$, zoom in around the origin

Figure 3.  The pink parts show the region of absolute stability of the general
fourth-order BDF scheme with Taylor expansion at $n+\beta,\,\beta=1,3,5$.

In order to have a better sense on how the stability regions vary with
different $\beta$ and $k$, we plot in Table 1 a comparison of stability
regions in the same scale. We observe that (i) the stability regions increase
faster when $\beta$ is closer to 1; and (ii)  the area of the stability region
with $k=4$ and $\beta=3$ is already bigger than that of the classical second-
order BDF. Hence, we can expect that the general fourth-order scheme with
$\beta=3$ allows similar or larger time steps for nonlinear problems than the
classical second-order IMEX, avoiding the usual scenario that smaller time
step has to be used when increasing the accuracy order.

| $\beta=1$ | $\beta=3$ | $\beta=5$
---|---|---|---
second order |  |  |
third order |  |  |
fourth order |  |  |
Table 1. Comparison of stability regions for different $k$ and $\beta$ on the
same scale.

we can derive the following formulas for the third order scheme based on the
Taylor expansion at time $t^{n+\beta}$,

$\displaystyle\frac{(3\beta^{2}+6\beta+2)\phi(t^{n+1})-(9\beta^{2}+12\beta-3)\phi(t^{n})+(9\beta^{2}+6\beta-6)\phi(t^{n-1})-(3\beta^{2}-1)\phi(t^{n-2})}{6\Delta
t}$ (2.21a) $\displaystyle=\partial_{t}\phi(t^{n+\beta})+\mathcal{O}(\Delta
t^{3}),$ (2.21b)
$\displaystyle\frac{\beta^{2}+\beta}{2}\phi(t^{n+1})-(\beta^{2}-1)\phi(t^{n})+\frac{\beta^{2}-\beta}{2}\phi(t^{n-1})=\phi(t^{n+\beta})+\mathcal{O}(\Delta
t^{3}),$ (2.21c)

and for the fourth order scheme,

$\displaystyle\frac{(2\beta^{3}+9\beta^{2}+11\beta+3)\phi(t^{n+1})+(-8\beta^{3}-30\beta^{2}-20\beta+10)\phi(t^{n})+(12\beta^{3}+36\beta^{2}+6\beta-18)\phi(t^{n-1})}{12\Delta
t}$ (2.22a)
$\displaystyle+\frac{(-8\beta^{3}-18\beta^{2}+4\beta+6)\phi(t^{n-2})+(2\beta^{3}+3\beta^{2}-\beta-1)\phi(t^{n-3})}{12\Delta
t}=\partial_{t}\phi(t^{n+\beta})+\mathcal{O}(\delta t^{4}),$ (2.22b)
$\displaystyle\frac{\beta^{3}+3\beta^{2}+2\beta}{6}\phi(t^{n+1})+\frac{-\beta^{3}-2\beta^{2}+\beta+2}{2}\phi(t^{n})+\frac{\beta^{3}+\beta^{2}-2\beta}{2}\phi(t^{n-1})+\frac{-\beta^{3}+\beta}{6}\phi(t^{n-2})$
(2.22c) $\displaystyle=\phi(t^{n+\beta})+\mathcal{O}(\Delta t^{4}).$ (2.22d)

##  3\. Multipliers for the new BDF and IMEX schemes

In order to conduct the stability and error analysis for the BDF and IMEX
schemes by using energy techniques, a key step is to find a suitable
multiplier. A key result which allows one to prove energy stability of the
classical BDF schemes of order up to five is established in [25] where the
existence of such multiplier is shown, see [3] for extension of this result to
six-order BDF. In this section, we identify an explicit multiplier, and show
that it is suitable for the new BDF and IMEX schemes of second to fourth
order.

###  3.1. Notations and a key lemma

To simplify the presentations, we introduce the following notations:

$\quad
A_{k}^{\beta}(\phi^{i})=\sum_{q=0}^{k}a_{k,q}(\beta)\phi^{i-k+q},\,\quad
B_{k}^{\beta}(\phi^{i})=\sum_{q=0}^{k-1}b_{k,q}(\beta)\phi^{i-k+1+q},\,\quad
C_{k}^{\beta}(\phi^{i})=\sum_{q=0}^{k-1}c_{k,q}(\beta)\phi^{i-k+1+q}.$ (3.1)

with $a_{k,q},b_{k,q},c_{k,q}$ defined in (2.12), (2.13) and (2.14). We also
consider the characteristic polynomials of the new BDF and IMEX schemes (2.10)
and (2.11):

$\displaystyle\tilde{A}_{k}^{\beta}(\zeta)=\sum_{q=0}^{k}a_{k,q}(\beta)\zeta^{q},\quad
k=2,3,4;$ (3.2a)
$\displaystyle\tilde{C}_{k}^{\beta}(\zeta)=\sum_{q=0}^{k-1}c_{k,q}(\beta)\zeta^{q},\quad
k=2,3,4.$ (3.2b)

We first recall the following result from Dahlquist’s G-stability theory [13]
which plays a key role in establishing energy stability of multistep methods.

######  Lemma 1.

Let $\alpha(\zeta)=\alpha_{k}\zeta^{k}+...+\alpha_{0}$ and
$\mu(\zeta)=\mu_{k}\zeta^{k}+...+\mu_{0}$ be polynomials of degree at most $k$
(and at least one of them of degree $k$) that have no common divisors. Let
$(\cdot,\cdot)$ be an inner product with associated norm $|\cdot|$. If

$\text{Re}\frac{\alpha(\zeta)}{\mu(\zeta)}>0\quad\text{for}\,|\zeta|>1,$ (3.3)

then there exists a symmetric positive definite matrix
$G=(g_{ij})\in\mathbb{R}^{k\times k}$ and real $\delta_{0},...,\delta_{k}$
such that for $\upsilon^{0},...,\upsilon^{k}$ in the inner product space,

$\big{(}\sum_{i=0}^{k}\alpha_{i}\upsilon^{i},\sum_{j=0}^{k}\mu_{j}\upsilon^{j}\big{)}=\sum_{i,j=1}^{k}g_{ij}(\upsilon^{i},\upsilon^{j})-\sum_{i,j=1}^{k}g_{ij}(\upsilon^{i-1},\upsilon^{j-1})+\big{|}\sum_{i=0}^{k}\delta_{i}\upsilon^{i}\big{|}^{2}.$
(3.4)

It is clear from the above Lemma that the key for establishing the energy
stability of (2.11) is to find a suitable multiplier
$\mu(\zeta)=\mu_{k}\zeta^{k}+...+\mu_{0}$ such that (3.3) is satisfied with
$\alpha(\zeta)=\tilde{A}_{k}^{\beta}(\zeta)$. To this end, we first split
$B_{k}^{\beta}(\phi^{n+1})$ into two parts:

$B_{k}^{\beta}(\phi^{n+1})=\eta_{k}(\beta)C_{k}^{\beta}(\phi^{n+1})+D_{k}^{\beta}(\phi^{n+1}),\quad
k=2,3,4,$ (3.5)

with

$\eta_{2}(\beta)=\frac{\beta-1}{\beta},\quad\eta_{3}(\beta)=\frac{\beta-1}{\beta+1},\quad\eta_{4}(\beta)=\frac{\beta-1}{\beta+3},\quad\beta\geq
1,$ (3.6)

and $D_{k}^{\beta}$ can be written as

$D^{\beta}_{k}(\phi^{n+1})=\sum_{q=0}^{k-1}d_{k,q}(\beta)\phi^{n+2-k+q},\quad
k=2,3,4,$ (3.7)

with

$\displaystyle d_{2,1}(\beta)=\frac{1}{\beta},\,d_{2,0}(\beta)=0,$ (3.8a)
$\displaystyle
d_{3,2}(\beta)=1,\,d_{3,1}(\beta)=\frac{1-\beta}{1+\beta},\,d_{3,0}(\beta)=0,$
(3.8b) $\displaystyle
d_{4,3}(\beta)=\frac{\beta^{2}}{6}+\frac{\beta}{2}+\frac{1}{3},\,d_{4,2}(\beta)=-(\frac{\beta^{2}}{2}+\frac{\beta}{2}-1),\,d_{4,1}(\beta)=\frac{\beta(\beta-1)}{2},\,d_{4,0}(\beta)=-\frac{\beta(\beta^{2}-1)}{6(\beta+3)}.$
(3.8c) $\displaystyle D_{2}^{\beta}(\phi^{n+1})=\frac{1}{\beta}\phi^{n+1},$
(3.9a) $\displaystyle
D_{3}^{\beta}(\phi^{n+1})=\phi^{n+1}+\frac{1-\beta}{1+\beta}\phi^{n},$ (3.9b)
$\displaystyle
D_{4}^{\beta}(\phi^{n+1})=(\frac{\beta^{2}}{6}+\frac{\beta}{2}+\frac{1}{3})\phi^{n+1}-(\frac{\beta^{2}}{2}+\frac{\beta}{2}-1)\phi^{n}+\frac{\beta(\beta-1)}{2}\phi^{n-1}-\frac{\beta(\beta^{2}-1)}{6(\beta+3)}\phi^{n-2}.$
(3.9c)

Similar to $A_{i}^{\beta}$, $B_{i}^{\beta}$ and $C_{i}^{\beta}$, we introduce
$d_{k,q}$ to denote the coefficients appears in (3.9):

and hence (3.9) can also be written as

$D^{\beta}_{k}(\phi^{n+1})=\sum_{q=0}^{k-1}d_{k,q}(\beta)\phi^{n+2-k+q},\quad
k=2,3,4.$ (3.10)

We also define

$\displaystyle\tilde{D}_{k}^{\beta}(\zeta)=\sum_{q=0}^{k-1}d_{k,q}(\beta)\zeta^{q},\quad
k=2,3,4.$ (3.11)

######  Remark 3.

The choices of $\eta_{i}(\beta)$ are not unique. We choose
$\eta_{2}(\beta),\,\eta_{3}(\beta)$ defined in (3.6) to make
$D_{2}^{\beta},\,D_{3}^{\beta}$ as simple as possible and the choice of
$\eta_{4}(\beta)$ defined in (3.6) allows us to prove (3.15) in the next
subsection.

###  3.2. A uniform multiplier

Note that in [25], it was shown that there exists a multiplier in the form of
$\phi^{n+1}-\tilde{\eta}_{k}\phi^{n}$ with $\tilde{\eta}_{k}\geq 0$ for the
usual BDF schemes of order 2 to 5. Surprisingly, we can find a uniform
multiplier for the new BDF and IMEX schemes of order 2 to 4. More precisely,
we have the following results.

######  Theorem 1.

Given $\beta\geq 1$, then

$\rm{gcd}\big{(}\tilde{A}^{\beta}_{k}(\zeta),\zeta\tilde{C}^{\beta}_{k}(\zeta)\big{)}=\rm{gcd}\big{(}\tilde{D}^{\beta}_{k}(\zeta),\tilde{C}^{\beta}_{k}(\zeta)\big{)}=1,\,k=2,3,4,$
(3.12)

i.e. they have no common divisor, and

$\text{Re}\frac{\tilde{A}^{\beta}_{k}(\zeta)}{\zeta\tilde{C}^{\beta}_{k}(\zeta)}>0,\quad\text{for}\,|\zeta|>1,\,k=2,3,4.$
(3.13)

Moreover, we also have

$\text{Re}\frac{\tilde{D}_{k}^{\beta}(\zeta)}{\tilde{C}_{k}^{\beta}(\zeta)}>0,\quad\text{for}\,|\zeta|>1,\,k=2,3;$
(3.14)

and finally if $\beta\geq 2$, then we also have

$\text{Re}\frac{\tilde{D}^{\beta}_{4}(\zeta)}{\tilde{C}_{4}^{\beta}(\zeta)}>0,\quad\text{for}\,|\zeta|>1.$
(3.15)

###### Proof.

The proof follows the basic process in [3]. We will provide the proof for the
case $k=4$ in detail as it includes some technical estimations and then we
will point out the key steps for the cases $k=2,3$, which are easier to
handle. To simplify the notation, we often omit the dependence on $\beta$ for
the coefficients $a_{k,q}(\beta),c_{k,q}(\beta),d_{k,q}(\beta)$, i.e., we only
write them as $a_{k,q},c_{k,q},d_{k,q}$.

Case I: $k=4$. Firstly, we show
$\rm{gcd}\big{(}\tilde{A}_{4}^{\beta}(\zeta),\zeta\tilde{C}_{4}^{\beta}(\zeta)\big{)}=1$
by using the Sylvester Resultant [1] as follows. The Sylvester matrix [1] of
$\tilde{A}_{4}^{\beta}(\zeta)$ and $\tilde{C}_{4}^{\beta}(\zeta)$ is

$Sly(\tilde{A}_{4}^{\beta},\tilde{C}_{4}^{\beta})={\left(\begin{array}[]{ccccccc}a_{4,4}&a_{4,3}&a_{4,2}&a_{4,1}&a_{4,0}&0&0\\\
0&a_{4,4}&a_{4,3}&a_{4,2}&a_{4,1}&a_{4,0}&0\\\
0&0&a_{4,4}&a_{4,3}&a_{4,2}&a_{4,1}&a_{4,0}\\\
c_{4,3}&c_{4,2}&c_{4,1}&c_{4,0}&0&0&0\\\
0&c_{4,3}&c_{4,2}&c_{4,1}&c_{4,0}&0&0\\\
0&0&c_{4,3}&c_{4,2}&c_{4,1}&c_{4,0}&0\\\
0&0&0&c_{4,3}&c_{4,2}&c_{4,1}&c_{4,0}\end{array}\right)}.$ (3.16)

It is easy to verify that its determinant is

$\text{det}Sly(\tilde{A}^{\beta}_{4},\tilde{C}^{\beta}_{4})=-\frac{1}{5184}\big{(}18\beta^{6}+144\beta^{5}+426\beta^{4}+566\beta^{3}+321\beta^{2}+55\beta+3\big{)}\neq
0,\;\text{ for }\beta\geq 1,$ (3.17)

which implies that
$\rm{gcd}\big{(}\tilde{A}_{4}^{\beta}(\zeta),\tilde{C}_{4}^{\beta}(\zeta)\big{)}=1$.
Combined with $\tilde{A}_{4}^{\beta}(0)=a_{4,0}\neq 0$, it also implies that
$\tilde{A}_{4}^{\beta}(\zeta)$ and $\zeta\tilde{C}_{4}^{\beta}(\zeta)$ have no
common divisor.

Next, we show
$\frac{\tilde{A}_{4}^{\beta}(\zeta)}{\zeta\tilde{C}_{4}^{\beta}(\zeta)}$ is
holomorphic outside the unit disk in the complex plane. To this end, it
suffices to show that all three zeros of $\tilde{C}_{4}^{\beta}(\zeta)$ are
inside the unit disk. Note that

$\frac{d\tilde{C}_{4}^{\beta}}{dx}(x)=3c_{4,3}x^{2}+2c_{4,2}x+c_{4,1},\,$
(3.18)

with

$c_{4,3}=\frac{\beta^{3}+6\beta^{2}+11\beta+6}{6}>0,\,\Delta_{4}:=4c_{4,2}^{2}-12c_{4,3}c_{4,1}=-\beta(\beta+2)(\beta+3)^{2}<0,$
(3.19)

which means $\tilde{C}_{4}(x)$ is monotonically increasing in the real axis.
Note also that

$\tilde{C}_{4}^{\beta}(0)=c_{4,0}=-\frac{\beta^{3}+3\beta^{2}+2\beta}{6}<0,\quad\tilde{C}^{\beta}_{4}(1)=c_{4,3}+c_{4,2}+c_{4,1}+c_{4,0}=1.$
(3.20)

Therefore, $\tilde{C}_{4}^{\beta}(\zeta)=0$ has exactly one real root, denoted
as $x_{1}$, and two complex roots, denoted as $z_{2},\,z_{3}=\bar{z}_{2}$, in
the complex plane. Next, we denote

$x_{0}:=\frac{-c_{4,0}}{c_{4,3}-1}=\frac{\beta^{2}+3\beta+2}{\beta^{2}+6\beta+11}.$
(3.21)

Then we can find with $\beta\geq 1$,

$\tilde{C}_{4}^{\beta}(x_{0})=-\frac{2\beta^{6}+27\beta^{5}+141\beta^{4}+351\beta^{3}+405\beta^{2}+162\beta-8}{(\beta^{2}+6\beta+11)^{3}}<0.$
(3.22)

Combining (3.20) and (3.22), we have $x_{0}<x_{1}<1$. On the other hand, by
Vieta’s formulae, we have

$x_{1}z_{2}z_{3}=x_{1}|z_{2}|^{2}=-\frac{c_{4,0}}{c_{4,3}},\quad\text{then}\,\,|z_{2}|^{2}=\frac{1}{x_{1}}\frac{-c_{4,0}}{c_{4,3}}<\frac{1}{x_{0}}\frac{-c_{4,0}}{c_{4,3}}=\frac{c_{4,3}-1}{c_{4,3}}<1.$
(3.23)

As a result, we have $|x_{1}|,|z_{2}|,|z_{3}|<1$ and hence
$\frac{\tilde{A}_{4}^{\beta}(\zeta)}{\zeta\tilde{C}_{4}^{\beta}(\zeta)}$ and
$\frac{\tilde{D}_{4}^{\beta}(\zeta)}{\tilde{C}_{4}^{\beta}(\zeta)}$ are
holomorphic outside the unit disk.

On the other hand, we have

$\lim_{|\zeta|\rightarrow\infty}\frac{\tilde{A}_{4}^{\beta}(\zeta)}{\zeta\tilde{C}_{4}^{\beta}(\zeta)}=\frac{a_{4,4}}{c_{4,3}}=\frac{2\beta^{3}+9\beta^{2}+11\beta+3}{2(\beta^{3}+6\beta^{2}+11\beta+6)}>0.$
(3.24)

Therefore, it follows from the maximum principle for harmonic functions,
$\text{Re}\frac{\tilde{A}_{4}^{\beta}(\zeta)}{\zeta\tilde{C}_{4}^{\beta}(\zeta)}>0,\,\forall|\zeta|>1$
is equivalent to

$\text{Re}\frac{\tilde{A}_{4}^{\beta}(\zeta)}{\zeta\tilde{C}_{4}^{\beta}(\zeta)}\geq
0,\quad\forall\zeta\in\mathbb{S}^{1},$ (3.25)

with $\mathbb{S}^{1}$ being the unit circle in the complex plane, and which is
equivalent to

$Re[\tilde{A}_{4}^{\beta}(e^{i\theta})e^{-i\theta}\tilde{C}_{4}^{\beta}(e^{-i\theta})]\geq
0,\quad\theta\in[0,2\pi).$ (3.26)

Letting $y:=\cos(\theta)$ and using the trigonometric identities

$\cos(2\theta)=2y^{2}-1,\,\cos(3\theta)=4y^{3}-3y,\,\sin(2\theta)=2y\sin(\theta),\,\sin(3\theta)=(4y^{2}-1)\sin(\theta),$
(3.27)

we find

$\begin{split}\tilde{C}_{4}^{\beta}(e^{-i\theta})&=c_{4,3}\cos(3\theta)+c_{4,2}\cos(2\theta)+c_{4,1}\cos(\theta)+c_{4,0}-i[c_{4,3}\sin(3\theta)+c_{4,2}\sin(2\theta)+c_{4,1}\sin(\theta)]\\\
&=c_{4,3}(4y^{3}-3y)+c_{4,2}(2y^{2}-1)+c_{4,1}y+c_{4,0}-i[c_{4,3}(4y^{2}-1)+2c_{4,2}y+c_{4,1}]\sin(\theta),\end{split}$
(3.28)

and

$\begin{split}\tilde{A}_{4}^{\beta}(e^{i\theta})e^{-i\theta}&=a_{4,4}e^{3i\theta}+a_{4,3}e^{2i\theta}+a_{4,2}e^{i\theta}+a_{4,1}+a_{4,0}e^{-i\theta}\\\
&=a_{4,4}\cos(3\theta)+a_{4,3}\cos(2\theta)+a_{4,2}\cos(\theta)+a_{4,1}+a_{4,0}\cos(\theta)\\\
&+i[a_{4,4}\sin(3\theta)+a_{4,3}\sin(2\theta)+a_{4,2}\sin(\theta)-a_{4,0}\sin(\theta)]\\\
&=a_{4,4}(4y^{3}-3y)+a_{4,3}(2y^{2}-1)+(a_{4,2}+a_{4,0})y+a_{4,1}\\\
&+i[a_{4,4}(4y^{2}-1)+2a_{4,3}y+a_{4,2}-a_{4,0}]\sin(\theta).\end{split}$
(3.29)

It follows from (3.28), (3.29) and $\tilde{A}_{4}^{\beta}(1)=0$,
$\sin^{2}(\theta)=1-y^{2}$ that

$Re[\tilde{A}_{4}^{\beta}(e^{i\theta})e^{-i\theta}\tilde{C}_{4}^{\beta}(e^{-i\theta})]=\frac{1}{9}(1-y)(\omega_{3}(\beta)y^{3}+\omega_{2}(\beta)y^{2}+\omega_{1}(\beta)y+\omega_{0}(\beta))=:\frac{1}{9}(1-y)f_{4}(y),$
(3.30)

with

$\begin{split}&f_{4}(y)=\omega_{3}(\beta)y^{3}+\omega_{2}(\beta)y^{2}+\omega_{1}(\beta)y+\omega_{0}(\beta),\\\
&\omega_{0}(\beta)=2\beta^{6}+15\beta^{5}+39\beta^{4}+39\beta^{3}+10\beta^{2}+15,\\\
&\omega_{1}(\beta)=-6\beta^{6}-45\beta^{5}-117\beta^{4}-116\beta^{3}-21\beta^{2}+17\beta+9,\\\
&\omega_{2}(\beta)=6\beta^{6}+45\beta^{5}+117\beta^{4}+115\beta^{3}+12\beta^{2}-34\beta-12,\\\
&\omega_{3}(\beta)=-2\beta^{6}-15\beta^{5}-39\beta^{4}-38\beta^{3}-\beta^{2}+17\beta+6.\end{split}$
(3.31)

In the following, we omit the dependence on $\beta$ for
$\omega_{i},\,i=0,1,2,3$.

It is clear that (3.26) is equivalent to

$f_{4}(y)\geq 0,\quad\forall y\in[-1,1].$ (3.32)

With $\omega_{i}$ defined in (3.31) and $\beta\geq 1$, we have

$\begin{split}f_{4}(1)&=\omega_{0}+\omega_{1}+\omega_{2}+\omega_{3}=18>0,\\\
f_{4}(-1)&=\omega_{0}-\omega_{1}+\omega_{2}-\omega_{3}=16\beta^{6}+120\beta^{5}+312\beta^{4}+308\beta^{3}+44\beta^{2}-68\beta-12>0,\end{split}$
(3.33)

and

$f^{\prime}_{4}(y)=3\omega_{3}y^{2}+2\omega_{2}y+\omega_{1}.$ (3.34)

If $f_{4}^{\prime}(y)$ does not have zero in $[-1,1]$, then (3.33) implies
(3.32). Otherwise, supposing there exists $-1\leq y_{0}\leq 1$ such that
$f_{4}^{\prime}(y_{0})=0$, we only need to show $f_{4}(y_{0})\geq 0$. Indeed,
with $f_{4}^{\prime}(y_{0})=0$, we have

$3f_{4}(y_{0})=3f_{4}(y_{0})-y_{0}f_{4}^{\prime}(y_{0})=\omega_{2}y_{0}^{2}+2\omega_{1}y_{0}+3\omega_{0}.$
(3.35)

Denote

$g_{4}(y):=\omega_{2}y^{2}+2\omega_{1}y+3\omega_{0};$ (3.36)

then with $\beta\geq 1$, we have

$\begin{split}g_{4}(1)&=\omega_{2}+2\omega_{1}+3\omega_{0}=51>0,\\\
g_{4}(-1)&=\omega_{2}-2\omega_{1}+3\omega_{0}=24\beta^{6}+180\beta^{5}+468\beta^{4}+464\beta^{3}+84\beta^{2}-68\beta+15>0,\\\
\Delta_{g}&=4\omega_{1}^{2}-12\omega_{2}\omega_{0}=-1220\beta^{6}-9108\beta^{5}-23408\beta^{4}-22212\beta^{3}-1076\beta^{2}+7344\beta+2484<0,\end{split}$
(3.37)

which means $g_{4}(y)>0,\forall y\in[-1,1]$. In particular, we have
$f_{4}(y_{0})=\frac{1}{3}g_{4}(y_{0})>0$ which implies (3.32), which in turn
implies (3.26). Therefore, we proved (3.13) with $k=4$.

Next, we prove (3.15) with $\beta\geq 2$. The procedure is similar to the
proof of (3.13) above. First, the Sylvester matrix of
$\tilde{D}_{4}^{\beta}(\zeta)$ and $\tilde{C}_{4}^{\beta}(\zeta)$:

$Sly(\tilde{D}_{4}^{\beta},\tilde{C}_{4}^{\beta})={\left(\begin{array}[]{cccccc}d_{4,3}&d_{4,2}&d_{4,1}&d_{4,0}&0&0\\\
0&d_{4,3}&d_{4,2}&d_{4,1}&d_{4,0}&0\\\ 0&0&d_{4,3}&d_{4,2}&d_{4,1}&d_{4,0}\\\
c_{4,3}&c_{4,2}&c_{4,1}&c_{4,0}&0&0\\\ 0&c_{4,3}&c_{4,2}&c_{4,1}&c_{4,0}&0\\\
0&0&c_{4,3}&c_{4,2}&c_{4,1}&c_{4,0}\end{array}\right)},$ (3.38)

and its determinant is

$\text{det}\,Sly(\tilde{D}_{4}^{\beta},\tilde{C}_{4}^{\beta})=-\frac{\beta^{2}(\beta^{2}+3\beta+2)^{2}}{36}<0$
(3.39)

which implies $\tilde{D}_{4}^{\beta}(\zeta)$ and
$\tilde{C}_{4}^{\beta}(\zeta)$ have no common divisor. Since we have shown in
the above that
$\frac{\tilde{D}_{4}^{\beta}(\zeta)}{\tilde{C}_{4}^{\beta}(\zeta)}$ is
holomorphic outside the unit disk, following the same process as above, we
have that (3.15) is equivalent to:

$h_{4}(y)=\alpha_{3}(\beta)y^{3}+\alpha_{2}(\beta)y^{2}+\alpha_{1}(\beta)y+\alpha_{0}(\beta)\geq
0,\,\forall y\in[-1,1],$ (3.40)

with

$\begin{split}&\alpha_{0}(\beta)=\frac{1}{9(\beta+3)}(2\beta^{6}+15\beta^{5}+35\beta^{4}+15\beta^{3}-37\beta^{2}-39\beta+9),\\\
&\alpha_{1}(\beta)=-\frac{2\beta^{5}}{3}-3\beta^{4}-\frac{10\beta^{3}}{3}+\beta^{2}+2\beta+1,\\\
&\alpha_{2}(\beta)=\frac{1}{3}(\beta(2\beta^{4}+9\beta^{3}+12\beta^{2}+3\beta-2)),\\\
&\alpha_{3}(\beta)=-\frac{1}{9}(\beta(\beta+1)^{2}(2\beta^{2}+5\beta+2)).\end{split}$
(3.41)

In the following, we omit the dependence on $\beta$ for
$\alpha_{i},i=0,1,2,3$. Hence, we have

$\begin{split}&h_{4}(-1)=-\alpha_{3}+\alpha_{2}-\alpha_{1}+\alpha_{0}=\frac{2}{9(\beta+3)}(8\beta^{6}+60\beta^{5}+152\beta^{4}+132\beta^{3}-16\beta^{2}-57\beta-9)>0,\\\
&h_{4}(1)=\alpha_{3}+\alpha_{2}+\alpha_{1}+\alpha_{0}=\frac{4}{\beta+3}>0,\end{split}$
(3.42)

and

$h^{\prime}_{4}(y)=3\alpha_{3}y^{2}+2\alpha_{2}y+\alpha_{1}.$ (3.43)

Similarly as before, if $h^{\prime}_{4}(y)$ does not have zero in $[-1,1]$,
then (3.42) implies (3.40). Suppose $-1\leq y_{0}\leq 1$ such that
$h^{\prime}_{4}(y_{0})=0$, we only need to show $h_{4}(y_{0})\geq 0$.

With $h_{4}^{\prime}(y_{0})=0$ and $\alpha_{3}\neq 0$, we have

$\begin{split}3h_{4}(y_{0})=3h_{4}(y_{0})-y_{0}h_{4}^{\prime}(y_{0})&=\alpha_{2}y_{0}^{2}+2\alpha_{1}y_{0}+3\alpha_{0}\\\
&=\frac{\alpha_{2}}{3\alpha_{3}}h_{4}^{\prime}(y_{0})+(2\alpha_{1}-\frac{2\alpha_{2}^{2}}{3\alpha_{3}})y_{0}+3\alpha_{0}-\frac{\alpha_{1}\alpha_{2}}{3\alpha_{3}}\\\
&=0+(2\alpha_{1}-\frac{2\alpha_{2}^{2}}{3\alpha_{3}})y_{0}+3\alpha_{0}-\frac{\alpha_{1}\alpha_{2}}{3\alpha_{3}}.\end{split}$
(3.44)

We define

$p(y):=(2\alpha_{1}-\frac{2\alpha_{2}^{2}}{3\alpha_{3}})y+3\alpha_{0}-\frac{\alpha_{1}\alpha_{2}}{3\alpha_{3}}.$
(3.45)

Then $p(y^{*})=0$ if we define $y^{*}$ as

$y^{*}:=\frac{\frac{\alpha_{1}\alpha_{2}}{3\alpha_{3}}-3\alpha_{0}}{2\alpha_{1}-\frac{2\alpha_{2}^{2}}{3\alpha_{3}}}=\frac{4\beta^{4}+30\beta^{3}+35\beta^{2}+3\beta}{4\beta^{4}+30\beta^{3}+71\beta^{2}+54\beta+9},$
(3.46)

and we also have

$2\alpha_{1}-\frac{2\alpha_{2}^{2}}{3\alpha_{3}}=\frac{8\beta^{2}+28\beta+6}{6\beta+3}>0.$
(3.47)

Therefore, to prove (3.40), it suffices to show $y_{0}\geq y^{*}$. However,
this is more complicated as (3.46) implies that $y^{*}$ can be arbitrarily
close to 1 by increasing $\beta$, and meanwhile, there indeed exists
$y^{*}<y_{0}<1$ such that $h_{4}^{\prime}(y_{0})=0$.

If follows from (3.43) that

$y_{0}=\frac{-2\alpha_{2}\pm\sqrt{\Delta_{h}}}{6\alpha_{3}},$ (3.48)

with

$\Delta_{h}:=4\alpha_{2}^{2}-12\alpha_{1}\alpha_{3}=\frac{4\beta(\beta+1)^{2}(4\beta^{3}+22\beta^{2}+31\beta+6)}{9}>0.$
(3.49)

We can estimate $\Delta_{h}$ as follows

$\Delta_{h}<\Delta_{h}+\frac{4(\beta+1)^{2}(2\beta^{3}+5\beta^{2}-6\beta)}{9}=\frac{4}{9}(\beta+1)^{2}(2\beta^{2}+6\beta)^{2}=:\Delta_{h}^{*}.$
(3.50)

To show $y_{0}\geq y^{*}$, we only consider the smallest root of
$h_{4}^{\prime}(y)=0$. Since we have $\alpha_{2}>0$ and $\alpha_{3}<0$, the
smallest root is

$y_{0}=\frac{-2\alpha_{2}+\sqrt{\Delta_{h}}}{6\alpha_{3}}>\frac{-2\alpha_{2}+\sqrt{\Delta_{h}^{*}}}{6\alpha_{3}}=\frac{2\beta^{3}+7\beta^{2}+3\beta-8}{2\beta^{3}+7\beta^{2}+7\beta+2}.$
(3.51)

Finally, we can prove $y_{0}\geq y^{*}$ as follows. It follows from (3.46) and
(3.51) that

$\begin{split}y_{0}-y^{*}&>\frac{2\beta^{3}+7\beta^{2}+3\beta-8}{2\beta^{3}+7\beta^{2}+7\beta+2}-\frac{4\beta^{4}+30\beta^{3}+35\beta^{2}+3\beta}{4\beta^{4}+30\beta^{3}+71\beta^{2}+54\beta+9}\\\
&=\frac{56\beta^{4}+138\beta^{3}-95\beta^{2}-339\beta-72}{8\beta^{6}+80\beta^{5}+300\beta^{4}+523\beta^{3}+430\beta^{2}+153\beta+18},\end{split}$
(3.52)

and given $\beta\geq 2$,

$\begin{split}56\beta^{4}+138\beta^{3}-95\beta^{2}-339\beta-72&\geq 56\times
2^{3}\beta+138\times 2\beta^{2}-95\beta^{2}-339\beta-72\\\
&=109\beta+181\beta^{2}-72>0.\end{split}$ (3.53)

Therefore, we have $y_{0}\geq y^{*}$. Hence (3.15) is proved for $\beta\geq
2$.

For the case $k=2$ and $3$, we can prove (3.13) and (3.14) by the same process
as above, so we only point out some related facts below, which are sufficient
to complete the proof.

Case II: $k=2$.

  * •

det$\,Sly(\tilde{A}^{\beta}_{2},\tilde{C}^{\beta}_{2})=-\frac{1}{2}\neq 0$,
det$\,Sly(\tilde{D}^{\beta}_{2},\tilde{C}^{\beta}_{2})=-1\neq
0,\,\tilde{A}^{\beta}_{2}(0)\neq 0$.

  * •

The only zero of $\tilde{C}_{2}^{\beta}(\zeta)$ is $\frac{\beta}{1+\beta}<1$,
which means
$\frac{\tilde{A}^{\beta}_{2}(\zeta)}{\zeta\tilde{C}^{\beta}_{2}(\zeta)}$ and
$\frac{\tilde{D}^{\beta}_{2}(\zeta)}{\tilde{C}^{\beta}_{2}(\zeta)}$ are
holomorphic outside the unit disk.

  * •

For $k=2$, (3.13) is equivalent to

$f_{2}(y)=(-2\beta^{2}-\beta+1)y+2\beta^{2}+\beta+1\geq 0,\quad\forall
y\in[-1,1],$ (3.54)

which is true since $f_{2}(y)$ is monotonically decreasing and $f_{2}(1)=2$.

  * •

For $k=2$, (3.14) is equivalent to

$h_{3}(y)=-y+1+\frac{1}{\beta}\geq 0,\quad\forall y\in[-1,1],$ (3.55)

which is obviously true.

Case III: $k=3$.

  * •

det$\,Sly(\tilde{A}_{3}^{\beta},\tilde{C}_{3}^{\beta})=\frac{\beta^{2}}{8}+\frac{5\beta}{24}+\frac{1}{36}\neq
0$,
det$\,Sly(\tilde{D}_{3}^{\beta},\tilde{C}_{3}^{\beta})=\frac{\beta(\beta+1)}{2}\neq
0,\,\tilde{A}^{\beta}_{3}(0)\neq 0,\,\forall\beta\geq 1.$

  * •

$\tilde{C}_{3}^{\beta}(\zeta)$ has two complex zeros $z_{1}$ and $z_{2}$ such
that $|z_{1}|^{2}=|z_{2}|^{2}=\frac{\beta}{\beta+2}<1$, which means
$\frac{\tilde{A}_{3}^{\beta}(\zeta)}{\zeta\tilde{C}_{3}^{\beta}(\zeta)}$ and
$\frac{\tilde{D}_{3}^{\beta}(\zeta)}{\tilde{C}_{3}^{\beta}(\zeta)}$ are
holomorphic outside the unit disk.

  * •

For $k=3$, (3.13) is equivalent to

$f_{3}(y)=\sigma_{2}(\beta)y^{2}+\sigma_{1}(\beta)y+\sigma_{0}\geq
0,\quad\forall y\in[-1,1],$ (3.56)

with

$\displaystyle\sigma_{2}(\beta)=3\beta^{4}+9\beta^{3}+5\beta^{2}-3\beta-2,$
(3.57a)
$\displaystyle\sigma_{1}(\beta)=-6\beta^{4}-18\beta^{3}-13\beta^{2}+\beta+4,$
(3.57b)
$\displaystyle\sigma_{0}(\beta)=3\beta^{4}+9\beta^{3}+8\beta^{2}+2\beta+4.$
(3.57c)

(3.56) is true since $\sigma_{2}(\beta)>0$ for $\beta\geq 1$ and

$\displaystyle f_{3}(-1)=12\beta^{4}+36\beta^{3}+26\beta^{2}-2\beta-2>0,$
(3.58a) $\displaystyle f_{3}(1)=6,$ (3.58b)
$\displaystyle\Delta_{3}:=\sigma_{1}^{2}-4\sigma_{0}\sigma_{2}=-63\beta^{4}-186\beta^{3}-95\beta^{2}+72\beta+48<0.$
(3.58c)

  * •

For $k=3$, (3.14) is equivalent to

$h_{3}(y)=\mu_{2}(\beta)y^{2}+\mu_{1}(\beta)y+\mu_{0}(\beta)\geq
0,\quad\forall y\in[-1,1],$ (3.59)

with

$\displaystyle\mu_{2}(\beta)=\beta(\beta+1),$ (3.60a)
$\displaystyle\mu_{1}(\beta)=-2\beta^{2}-2\beta+1,$ (3.60b)
$\displaystyle\mu_{0}(\beta)=\frac{\beta^{3}+2\beta^{2}+1}{\beta+1}.$ (3.60c)

(3.59) is true since $\mu_{2}(\beta)>0$ for $\beta\geq 1$ and

$\displaystyle h_{3}(-1)=\frac{2\beta(2\beta^{2}+4\beta+1)}{\beta+1}>0,$
(3.61a) $\displaystyle h_{3}(1)=\frac{2}{\beta+1}>0,$ (3.61b)
$\displaystyle\Delta_{3}^{*}:=\mu_{1}^{2}-4\mu_{0}\mu_{2}=1-8\beta<0.$ (3.61c)

The proof for all the cases is completed. ∎

######  Remark 4.

The restriction $\beta\geq 2$ is a sufficient condition for (3.15), which
comes from (3.53). One can easily show that (3.39) and (3.53) are true
whenever $\beta>1.6$. On the other hand, (3.40) is not true when $\beta=1$ as
$h_{4}(0.2)=-0.312<0$ with $h_{4}$ defined in (3.40).

###  3.3. Explicit telescoping formulae for the second and third order schemes

Note that Lemma 1 only provides the existence of a symmetric positive definite
matrix $G$ without giving the exact value of $g_{ij}$. In the following, we
provide explicit formulae for $g_{ij}$ in the second and third order cases.

######  Proposition 3.1.

For the second-order version of (2.11), we have

$\big{(}D_{2}^{\beta}(\phi^{n+1}),C_{2}^{\beta}(\phi^{n+1})\big{)}=\frac{1}{\beta}|\phi^{n+1}|^{2}+\frac{1}{2}|\phi^{n+1}|^{2}-\frac{1}{2}|\phi^{n}|^{2}+\frac{1}{2}|\phi^{n+1}-\phi^{n}|^{2},$
(3.62)

and

$\begin{split}\big{(}A_{2}^{\beta}(\phi^{n+1}),C_{2}^{\beta}(\phi^{n+1})\big{)}&=a_{2}|\phi^{n+1}|^{2}-a_{2}|\phi^{n}|^{2}+|b_{2}\phi^{n+1}+c_{2}\phi^{n}|^{2}-|b_{2}\phi^{n}+c_{2}\phi^{n-1}|^{2}\\\
&+|d_{2}\phi^{n+1}+e_{2}\phi^{n}+f_{2}\phi^{n-1}|^{2},\end{split}$ (3.63)

where the coefficients are given by

$\begin{split}&e_{2}=-\sqrt{2\beta(2\beta+1)},\,\Delta_{2}=2\beta(2\beta+1),\,c_{2}=f_{2}=\frac{-\sqrt{2}+\sqrt{\Delta_{2}}}{2},\\\
&d_{2}=\sqrt{2}+f_{2},\,E_{2}=-\beta(2\beta-1),\,b_{2}=\frac{E_{2}-2e_{2}f_{2}}{-2c},\,a_{2}=\frac{3\beta+1-2\sqrt{\beta(2\beta+1)}}{2(\beta+1)^{2}}.\end{split}$

Moreover, we have $a_{2}>0$ for all $\beta\geq 1$.

######  Proposition 3.2.

For the third-order version of (2.11), we have

$\begin{split}\big{(}D_{3}^{\beta}(\phi^{n+1}),C_{3}^{\beta}(\phi^{n+1})\big{)}&=\hat{a}_{3}|\phi^{n+1}|^{2}-\hat{a}_{3}|\phi^{n}|^{2}+|\hat{b}_{3}\phi^{n+1}+\hat{c}_{3}\phi^{n}|^{2}-|\hat{b}_{3}\phi^{n}+\hat{c}_{3}\phi^{n-1}|^{2}\\\
&+|\hat{d}_{3}\phi^{n+1}+\hat{e}_{3}\phi^{n}+\hat{f}_{3}\phi^{n-1}|^{2},\end{split}$
(3.64)

where the coefficients are given by

$\begin{split}&\hat{M}=\frac{2\beta^{3}+4\beta^{2}+\beta+1}{\beta+1},\,\hat{N}=\frac{(2\beta^{2}+2\beta-1)^{2}}{4},\,\hat{\Delta}_{3}=\hat{M}^{2}-4\hat{N}=\frac{4\beta(2\beta^{2}+4\beta+1)}{(\beta+1)^{2}}\\\
&\hat{e}_{3}=-\sqrt{\frac{\hat{M}-\sqrt{\hat{\Delta}_{3}}}{2}},\,\hat{P}=\frac{\beta^{3}+2\beta^{2}+1}{\beta+1}-\hat{e}_{3}^{2},\,\hat{Q}=\frac{2\beta^{3}+4\beta^{2}+\beta+1}{\beta+1}-\hat{e}_{3}^{2},\\\
&\hat{f}_{3}=\frac{-\sqrt{\hat{P}}+\sqrt{\hat{Q}}}{2},\,\hat{c}_{3}=\hat{f}_{3},\,\hat{d}_{3}=\sqrt{\hat{P}}+\hat{f}_{3},\,\hat{b}_{3}=\frac{\beta(\beta-1)+4\hat{e}_{3}\hat{f}_{3}}{4\hat{c}_{3}},\,\hat{a}_{3}=\frac{\beta^{2}}{2}+\frac{3\beta}{2}+1-\hat{b}_{3}^{2}-\hat{d}_{3}^{2};\end{split}$
(3.65)

and

$\begin{split}&\big{(}A_{3}^{\beta}(\phi^{n+1}),C_{3}^{\beta}(\phi^{n+1})\big{)}=a_{3}|\phi^{n+1}|^{2}-a_{3}|\phi^{n}|^{2}+|b_{3}\phi^{n+1}+c_{3}\phi^{n}|^{2}-|b_{3}\phi^{n}+c_{3}\phi^{n-1}|^{2}\\\
&+|d_{3}\phi^{n+1}+e_{3}\phi^{n}+f_{3}\phi^{n-1}|^{2}-|d_{3}\phi^{n}+e_{3}\phi^{n-1}+f_{3}\phi^{n-2}|^{2}+|g_{3}\phi^{n+1}+h_{3}\phi^{n}+i_{3}\phi^{n-1}+j_{3}\phi^{n-2}|^{2},\end{split}$
(3.66)

where the coefficients are given by

$\begin{split}&M=2\beta^{4}+6\beta^{3}+\frac{13\beta^{2}}{3}-\frac{\beta}{3}-\frac{1}{3},\,N=-(\frac{\beta^{2}}{2}-\frac{1}{6})(\frac{\beta^{2}}{2}+\frac{3\beta}{2}+1),\,P=\frac{\sqrt{M}+1}{2},\\\
&Q=-\frac{1}{2}(\beta(\frac{\beta^{2}}{2}-\frac{1}{6})(\beta+1)),\,R=\beta^{4}+\frac{7\beta^{3}}{2}+\frac{19\beta^{2}}{6}-\frac{\beta}{3}-1,\,S=\frac{7\beta^{4}}{4}+\frac{25\beta^{3}}{4}+\frac{17\beta^{2}}{3}+\frac{\beta}{2}+\frac{1}{3},\\\
&W=(\frac{\beta^{2}}{2}+\frac{3\beta}{2}+1)(\frac{\beta^{2}}{2}+\beta+\frac{1}{3}),\,U=\frac{1}{2}-\frac{79\beta^{2}}{12}-\frac{21\beta^{3}}{4}-\frac{5\beta^{4}}{4}-\frac{23\beta}{12},\\\
&f_{3}=\frac{\sqrt{P^{2}+2N}+P}{2},\,j_{3}=f_{3},\,g_{3}=f_{3}-P,\,i_{3}=-\sqrt{M}-g_{3},\,h_{3}=\sqrt{M}-f_{3},\,e_{3}=\frac{2i_{3}j_{3}-Q}{2f_{3}},\\\
&d_{3}=\frac{R-2g_{3}i_{3}}{2f_{3}},\,c_{3}=\sqrt{S-e_{3}^{2}-g_{3}^{2}-h_{3}^{2}},\,b_{3}=\frac{U-2d_{3}e_{3}-2g_{3}h_{3}}{2c_{3}},\,a_{3}=W-g_{3}^{2}-d_{3}^{2}-b_{3}^{2}.\end{split}$
(3.67)

Moreover, it is numerically verified that all variables appearing in (3.65)
and (3.67) are real and bounded, and $\hat{a}_{3},\,a_{3}>0$ for
$1\leq\beta\leq 100$ (cf. Fig. 4).

The proof of the above two propositions is based on the method of undetermined
coefficients, more precisely, we assume a desired form and use the method of
undetermined coefficients to find the suitable coefficients. The detail of the
proof is tedious but straightforward so we leave it to the interested readers.

Figure 4.  Values of $\hat{a}_{3}$ and $a_{3}$ with different $\beta$.

##  4\. Stability of (2.10) for linear parabolic type equations

We consider in this section the new BDF schemes for the linear case (2.10),
which can be written as

$\frac{A_{k}^{\beta}(\phi^{n+1})}{\Delta
t}+\mathcal{L}B_{k}^{\beta}(\phi^{n+1})=f^{n+\beta},\quad k=2,3,4,$ (4.1)

and establish a stability result based on Theorem 1.

######  Theorem 2.

Assuming $\|f(t)\|_{\star}^{2}\leq C_{f},\,\forall t\leq T$, $\beta>1$ for
$k=2,3$, and $\beta\geq 2$ for $k=4$, then the scheme (4.1) is stable in the
sense that

$g_{k}|\phi^{n+1}|^{2}+\frac{1}{2}\Delta
t\eta_{k}(\beta)\sum_{q=k}^{n+1}\|C_{k}^{\beta}(\phi^{q})\|^{2}\leq
C\sum_{q=0}^{k-1}\big{(}|\phi^{q}|^{2}+\Delta
t\|\phi^{q}\|^{2}\big{)}+\frac{TC_{f}}{2\eta_{k}(\beta)},\quad\forall k\leq
n+1\leq\frac{T}{\Delta t},$ (4.2)

with $g_{k}$ a positive constant depending only on $k$, $C$ a constant
independent of $\Delta t$ and $\eta_{k}(\beta)$ is defined in (3.6).

###### Proof.

We denote $f^{i}=f(t^{i}),\,\forall i\leq\frac{T}{\Delta t}$. Taking the inner
product of (4.1) with $\Delta tC_{k}^{\beta}(\phi^{n+1})$ and splitting
$B_{k}^{\beta}(\phi^{n+1})$ as in (3.5), we obtain

$\big{(}A_{k}^{\beta}(\phi^{n+1}),C_{k}^{\beta}(\phi^{n+1})\big{)}+\Delta
t\eta_{k}(\beta)\|C_{k}^{\beta}(\phi^{n+1})\|^{2}+\Delta
t\big{(}\mathcal{L}D_{k}^{\beta}(\phi^{n+1}),C_{k}^{\beta}(\phi^{n+1})\big{)}=\Delta
t\big{(}f^{n+\beta},C_{k}^{\beta}(\phi^{n+1})\big{)},$ (4.3)

where we used
$\big{(}\mathcal{L}C_{k}^{\beta}(\phi^{n+1}),C_{k}^{\beta}(\phi^{n+1})\big{)}=\|C_{k}^{\beta}(\phi^{n+1})\|^{2}$.
We estimate the terms in (4.3) as follows.

It follows from (2.1) and the assumption on $f$ that

$\begin{split}\big{(}f^{n+\beta},C_{k}^{\beta}(\phi^{n+1})\big{)}&\leq\|f^{n+\beta}\|_{\star}\|C_{k}^{\beta}(\phi^{n+1})\|\\\
&\leq\frac{1}{2\eta_{k}(\beta)}\|f^{n+\beta}\|_{\star}^{2}+\frac{\eta_{k}(\beta)}{2}\|C_{k}^{\beta}(\phi^{n+1})\|^{2}\\\
&\leq\frac{C_{f}}{2\eta_{k}(\beta)}+\frac{\eta_{k}(\beta)}{2}\|C_{k}^{\beta}(\phi^{n+1})\|^{2}.\\\
\end{split}$ (4.4)

Denote $\Phi^{n+1}_{k}:=(\phi^{n-k+1},...,\phi^{n+1})^{T}$. It follows from
Lemma 1 and Theorem 1 that there exist symmetric positive definite matrices
$G=(g_{ij})\in\mathbb{R}^{k\times k}$ and
$H=(h_{ij})\in\mathbb{R}^{(k-1)\times(k-1)}$ such that

$\begin{split}\big{(}A_{k}^{\beta}(\phi^{n+1}),C_{k}^{\beta}(\phi^{n+1})\big{)}&\geq\sum_{i,j=1}^{k}g_{ij}(\phi^{n+1+i-k},\phi^{n+1+j-k})-\sum_{i,j=1}^{k}g_{ij}(\phi^{n+i-k},\phi^{n+j-k})\\\
&=:|\Phi^{n+1}_{k}|^{2}_{G}-|\Phi^{n}_{k}|^{2}_{G},\end{split}$ (4.5)

and

$\begin{split}\big{(}\mathcal{L}D_{k}^{\beta}(\phi^{n+1}),C_{k}^{\beta}(\phi^{n+1})\big{)}&\geq\sum_{i,j=1}^{k-1}h_{ij}(\mathcal{L}\phi^{n+2+i-k},\phi^{n+2+j-k})-\sum_{i,j=1}^{k-1}h_{ij}(\mathcal{L}\phi^{n+1+i-k},\phi^{n+1+j-k})\\\
&=:\|\Phi^{n+1}_{k}\|^{2}_{H}-\|\Phi^{n}_{k}\|^{2}_{H}.\end{split}$ (4.6)

Now, combining (4.3)-(4.6), we obtain

$|\Phi^{n+1}_{k}|^{2}_{G}-|\Phi^{n}_{k}|^{2}_{G}+\Delta
t(\|\Phi^{n+1}_{k}\|^{2}_{H}-\|\Phi^{n}_{k}\|^{2}_{H})+\frac{1}{2}\Delta
t\eta_{k}(\beta)\|C_{k}^{\beta}(\phi^{n+1})\|^{2}\leq\frac{C_{f}\Delta
t}{2\eta_{k}(\beta)}.$ (4.7)

Summing up (4.7) from $n=k-1$ to $n=m$, we obtain

$|\Phi^{m+1}_{k}|^{2}_{G}+\Delta t\|\Phi^{m+1}_{k}\|^{2}_{H}+\frac{1}{2}\Delta
t\eta_{k}(\beta)\sum_{q=k-1}^{m}\|C_{k}^{\beta}(\phi^{q+1})\|^{2}\leq|\Phi^{k-1}_{k}|^{2}_{G}+\Delta
t\|\Phi^{k-1}_{k}\|^{2}_{H}+\frac{TC_{f}}{2\eta_{k}(\beta)}.$ (4.8)

Let $g_{k}$ be the smallest eigenvalue of the matrix $G\in\mathbb{R}^{k,k}$,
then we have

$|\Phi^{m+1}_{k}|^{2}_{G}\geq g_{k}|\phi^{m+1}|^{2},$ (4.9)

and we can choose a constant $C$ large enough such that

$\displaystyle|\Phi^{k-1}_{k}|^{2}_{G}\leq C\sum_{i=0}^{k-1}|\phi^{i}|^{2},$
(4.10a) $\displaystyle\Delta t\|\Phi^{k-1}_{k}\|^{2}_{H}\leq C\Delta
t\sum_{i=0}^{k-1}\|\phi^{i}\|^{2}.$ (4.10b)

Finally, combining (4.8) and (4.10) leads to

$g_{k}|\phi^{m+1}|^{2}+\frac{1}{2}\Delta
t\eta_{k}(\beta)\sum_{q=k-1}^{m}\|C_{k}^{\beta}(\phi^{q+1})\|^{2}\leq
C\sum_{i=0}^{k-1}(|\phi^{i}|^{2}+\Delta
t\|\phi^{i}\|^{2})+\frac{TC_{f}}{2\eta_{k}(\beta)},$ (4.11)

which implies (4.2). ∎

######  Remark 5.

Note that in order to obtain (4.6), the linear operator $\mathcal{L}$ is
required to be self-adjoint while using the Nevanlinna-Odeh approach in [25]
can also deal with $\mathcal{L}$ which is not self-adjoint.

##  5\. Stability and error analysis of (2.11) for nonlinear parabolic type
equations

In this section, we use the stability result established in the last section
to carry out a stability and error analysis of (2.11) for nonlinear parabolic
equations.

###  5.1. Stability

Under the local Lipschitz condition (2.2) on the nonlinear operator
$\mathcal{G}$, we can derive a local stability result for (2.11) similarly as
in the proof of the linear case (cf. Theorem 2) if we further assume

$C_{k}^{\beta}(\phi^{n})\in\mathcal{B}_{\phi(t^{n+\beta})},$ (5.1)

with $\beta>1$ for $k=2,3$, and $\beta\geq 2$ for $k=4$. Note that formally
(5.1) must be true when $\Delta t$ small enough since
$C_{k}^{\beta}(\phi^{n})$ is a $k$-th order approximation to
$\phi(t^{n+\beta})$. We shall defer the rigorous proof of (5.1) to subsection
5.3 by induction together with the error analysis.

######  Theorem 3.

Assume (2.2) and (5.1), then for $k=2,3,4$, under the stability condition

$\eta_{k}(\beta)-\sqrt{\gamma}\geq\rho>0,$ (5.2)

the general implicit-explicit IMEX method (2.11) is locally stable in the
sense that, with $\vartheta^{m}:=\phi^{m}-\psi^{m}$, for $\Delta t$
sufficiently small, we have

$\delta_{k}|\vartheta^{m+1}|^{2}+\rho\Delta
t\sum_{q=k-1}^{m}\|C_{k}(\vartheta^{q+1})\|^{2}\leq
C_{0}\exp(C_{0}T)\sum_{i=0}^{k-1}(|\vartheta^{i}|^{2}+\Delta
t\|\vartheta^{i}\|^{2}),\quad k\leq m+1\leq\frac{T}{\Delta t},$ (5.3)

with $\delta_{k}$ a positive constant depending only on $k$, $C_{0}$ a
constant independent of $\Delta t$ and $C_{k}(\cdot)$ is defined in section
LABEL:notation.

###### Proof.

Denoting $g^{m}:=\mathcal{G}[C_{k}(\phi^{m})]-\mathcal{G}[C_{k}(\psi^{m})]$
and subtracting (LABEL:Gscheme2) from (2.11), we obtain

$\frac{A_{k}(\vartheta^{n+1})}{\Delta
t}+\mathcal{L}B_{k}(\vartheta^{n+1})+g^{n}=0.$ (5.4)

Now, we take in (5.4) the inner product with $\Delta tC_{k}(\vartheta^{n+1})$
and split $B_{k}(\vartheta^{n+1})$ as in (3.5) to obtain

$\big{(}A_{k}(\vartheta^{n+1}),C_{k}(\vartheta^{n+1})\big{)}+\Delta
t\eta_{k}(\beta)\|C_{k}(\vartheta^{n+1})\|^{2}+\Delta
t\big{(}\mathcal{L}D_{k}(\vartheta^{n+1}),C_{k}(\vartheta^{n+1})\big{)}=-\Delta
t\big{(}g^{n},C_{k}(\vartheta^{n+1})\big{)},$ (5.5)

where we use
$\big{(}\mathcal{L}C_{k}(\vartheta^{n+1}),C_{k}(\vartheta^{n+1})\big{)}=\|C_{k}(\vartheta^{n+1})\|^{2}$.
We estimate the terms in (5.5) as follows. Firstly, it follows from (2.1) and
(2.2) that

$\begin{split}-\big{(}g^{n},C_{k}(\vartheta^{n+1})\big{)}&\leq\|g^{n}\|_{\star}\|C_{k}(\vartheta^{n+1})\|\\\
&\leq\frac{\varepsilon}{2}\|g^{n}\|_{\star}^{2}+\frac{1}{2\varepsilon}\|C_{k}(\vartheta^{n+1})\|^{2}\\\
&\leq\frac{\varepsilon}{2}(\gamma\|C_{k}(\vartheta^{n})\|^{2}+\mu|C_{k}(\vartheta^{n})|^{2})+\frac{1}{2\varepsilon}\|C_{k}(\vartheta^{n+1})\|^{2},\\\
\end{split}$ (5.6)

with $\varepsilon>0$ can be arbitrarily small.

On the other hand, we use the notation
$\Theta^{n+1}_{k}:=(\vartheta^{n-k+1},...,\vartheta^{n+1})^{T}$ and two norms
$|\Theta^{n+1}_{k}|_{G}$, $\|\Theta^{n+1}_{k}\|_{H}$ defined below. It follows
from Lemma 1 and Theorem 1 that there exists symmetric positive definite
matrices $G=(g_{ij})\in\mathbb{R}^{k\times k}$ and
$H=(h_{ij})\in\mathbb{R}^{(k-1)\times(k-1)}$ such that

$\begin{split}\big{(}A_{k}(\vartheta^{n+1}),C_{k}(\vartheta^{n+1})\big{)}&\geq\sum_{i,j=1}^{k}g_{ij}(\vartheta^{n+1+i-k},\vartheta^{n+1+j-k})-\sum_{i,j=1}^{k}g_{ij}(\vartheta^{n+i-k},\vartheta^{n+j-k})\\\
&:=|\Theta^{n+1}_{k}|^{2}_{G}-|\Theta^{n}_{k}|^{2}_{G}\end{split}$ (5.7)

and

$\begin{split}\big{(}\mathcal{L}D_{k}(\vartheta^{n+1}),C_{k}(\vartheta^{n+1})\big{)}&\geq\sum_{i,j=1}^{k-1}h_{ij}(\mathcal{L}\vartheta^{n+2+i-k},\vartheta^{n+2+j-k})-\sum_{i,j=1}^{k-1}h_{ij}(\mathcal{L}\vartheta^{n+1+i-k},\vartheta^{n+1+j-k})\\\
&:=\|\Theta^{n+1}_{k}\|^{2}_{H}-\|\Theta^{n}_{k}\|^{2}_{H}\end{split}$ (5.8)

Now, combining (5.5)-(5.8) and choose $\varepsilon=\frac{1}{\sqrt{\gamma}}$,
we obtain

$\begin{split}|\Theta^{n+1}_{k}|^{2}_{G}-|\Theta^{n}_{k}|^{2}_{G}&+\Delta
t(\|\Theta^{n+1}_{k}\|^{2}_{H}-\|\Theta^{n}_{k}\|^{2}_{H})+\Delta
t(\eta_{k}(\beta)-\frac{\sqrt{\gamma}}{2})\|C_{k}(\vartheta^{n+1})\|^{2}\\\
&\leq\Delta
t\big{(}\frac{\sqrt{\gamma}}{2}\|C_{k}(\vartheta^{n})\|^{2}+\frac{\mu}{2\sqrt{\gamma}}|C_{k}(\vartheta^{n})|^{2}\big{)}\end{split}$
(5.9)

Summing (5.9) from $n=k-1$ to $n=m$, we obtain

$\begin{split}&|\Theta^{m+1}_{k}|^{2}_{G}+\Delta
t\|\Theta^{m+1}_{k}\|^{2}_{H}+\Delta
t(\eta_{k}(\beta)-\sqrt{\gamma})\sum_{q=k-1}^{m}\|C_{k}(\vartheta^{q+1})\|^{2}\\\
&\leq|\Theta^{k-1}_{k}|^{2}_{G}+\Delta
t\|\Theta^{k-1}_{k}\|^{2}_{H}+\frac{\Delta
t\sqrt{\gamma}}{2}\|C_{k}(\vartheta^{k-1})\|^{2}+\frac{\mu\Delta
t}{2\sqrt{\gamma}}\sum_{q=k-1}^{m}|C_{k}(\vartheta^{q})|^{2}.\end{split}$
(5.10)

Suppose $\delta_{k}$ is the smallest eigenvalue of the matrix
$G\in\mathbb{R}^{k,k}$, then we have

$|\Theta^{m+1}_{k}|^{2}_{G}\geq\delta_{k}|\vartheta^{m+1}|^{2},$ (5.11)

and we can choose a constant $C_{0}$ large enough such that

$\displaystyle|\Theta^{k-1}_{k}|^{2}_{G}\leq
C_{0}\sum_{i=0}^{k-1}|\vartheta^{i}|^{2},$ (5.12a) $\displaystyle\Delta
t\|\Theta^{k-1}_{k}\|^{2}_{H}+\frac{\Delta
t\sqrt{\gamma}}{2}\|C_{k}(\vartheta^{k-1})\|^{2}\leq C_{0}\Delta
t\sum_{i=0}^{k-1}\|\vartheta^{i}\|^{2},$ (5.12b) $\displaystyle\frac{\mu\Delta
t}{2\sqrt{\gamma}}\sum_{q=k-1}^{m}|C_{k}(\vartheta^{q})|^{2}\leq C_{0}\Delta
t\sum_{q=0}^{m}|\vartheta^{q}|^{2},$ (5.12c)

then (5.10) and (5.19) together imply

$\delta_{k}|\vartheta^{m+1}|^{2}+\rho\Delta
t\sum_{q=k-1}^{m}\|C_{k}(\vartheta^{q+1})\|^{2}\leq
C_{0}\sum_{i=0}^{k-1}(|\vartheta^{i}|^{2}+\Delta
t\|\vartheta^{i}\|^{2})+C_{0}\Delta t\sum_{q=0}^{m}|\vartheta^{q}|^{2}.$
(5.13)

Finally, we get (5.3) by applying the discrete Gronwall lemma LABEL:Gron2 to
(5.13). ∎

###  5.2. Truncation errors

Using the notations introduced in previous sections, we define the truncation
errors for $k=2,3,4$ as

$\displaystyle E_{k}^{n+1}:=\Delta
t\phi_{t}(t^{n+\beta})-A_{k}^{\beta}(\phi(t^{n+1})),$ (5.14a) $\displaystyle
R_{k}^{n+1}:=\phi(t^{n+\beta})-B_{k}^{\beta}(\phi(t^{n+1})),$ (5.14b)
$\displaystyle P_{k}^{n}:=\phi(t^{n+\beta})-C_{k}^{\beta}(\phi(t^{n})).$
(5.14c)

It follows from (2.4), (2.6) and (2.8) that

$E_{k}^{n+1}=\mathcal{O}(\Delta t^{k+1}),\quad R_{k}^{n+1}=\mathcal{O}(\Delta
t^{k}),\quad\quad P_{k}^{n}=\mathcal{O}(\Delta t^{k}).\\\ $ (5.15)

More precisely, one can verify

$\displaystyle
E_{k}^{n+1}=\frac{1}{k!}\sum_{q=0}^{k}a_{k,q}(\beta)\int_{t^{n+1+q-k}}^{t^{n+\beta}}(t^{n+1+q-k}-s)^{k}\phi^{(k+1)}(s)ds,$
(5.16a) $\displaystyle
R_{k}^{n+1}=\frac{1}{(k-1)!}\sum_{q=0}^{k-1}b_{k,q}(\beta)\int_{t^{n+2+q-k}}^{t^{n+\beta}}(t^{n+2+q-k}-s)^{k-1}\phi^{(k)}(s)ds,$
(5.16b) $\displaystyle
P_{k}^{n}=\frac{1}{(k-1)!}\sum_{q=0}^{k-1}c_{k,q}(\beta)\int_{t^{n+1+q-k}}^{t^{n+\beta}}(t^{n+1+q-k}-s)^{k-1}\phi^{(k)}(s)ds.$
(5.16c)

Therefore, under suitable regularity requirements, we have

$|E_{k}^{n+1}|^{2}\leq C(\Delta t)^{2k+2},\quad\|R_{k}^{n+1}\|^{2}\leq
C(\Delta t)^{2k},\quad\|P_{k}^{n}\|^{2}\leq C(\Delta t)^{2k},\quad\forall
n+1\leq\frac{T}{\Delta t}.$ (5.17)

###  5.3. Error estimate

We denote $e^{m}:=\phi^{m}-\phi(t^{m})$, where $\phi(t^{m})$ is the exact
solution of (1.1) at time $t^{m}$, i.e.,

$\phi_{t}(t^{m})+\mathcal{L}\phi({t^{m}})+\mathcal{G}[\phi(t^{m})]=f(t^{m}).$
(5.18)

We will use the following discrete version of the Gronwall lemma [26].

######  Lemma 2.

Let $y^{k},\,h^{k},\,g^{k},\,f^{k}$ be four nonnegative sequences satisfying

$y^{n}+\Delta t\sum_{k=0}^{n}h^{k}\leq B+\Delta
t\sum_{k=0}^{n}(g^{k}y^{k}+f^{k})\;\text{ with }\;\Delta t\sum_{k=0}^{T/\Delta
t}g^{k}\leq M,\,\forall\,0\leq n\leq T/\Delta t.$

We assume $\Delta t\,g^{k}<1$ for all $k$, and let $\sigma=\max_{0\leq k\leq
T/\Delta t}(1-\Delta tg^{k})^{-1}$. Then

$y^{n}+\Delta t\sum_{k=1}^{n}h^{k}\leq\exp(\sigma M)(B+\Delta
t\sum_{k=0}^{n}f^{k}),\,\,\forall\,n\leq T/\Delta t.$

######  Theorem 4.

Assume (2.2) and the solution of (1.1) is sufficiently smooth such that (5.17)
is true, and the following stability condition

$\eta_{k}(\beta)-\sqrt{\gamma}\geq\rho>0$ (5.19)

is satisfied. Given $\phi^{0}=\phi(0)\in V$, we assume $\beta>1$ for $k=2,3$,
and $\beta\geq 2$ for $k=4$, and that $\phi^{i},i=1,...,k-1,$ are computed
with a proper initialization procedure such that

$|\phi^{i}-\phi(t^{i})|^{2},\,\|\phi^{i}-\phi(t^{i})\|^{2}\leq C(\Delta
t)^{2k},i=1,...,k-1,\,\text{and}\,\,C_{k}^{\beta}(\phi^{k-1})\in\mathcal{B}_{\phi(t^{k-1+\beta})};$
(5.20)

then for $\Delta t$ sufficiently small, we have

$C_{k}^{\beta}(\phi^{n+1})\in\mathcal{B}_{\phi(t^{n+1+\beta})},\quad\forall
n+1\leq\frac{T}{\Delta t},$ (5.21)

and

$g_{k}|e^{n+1}|^{2}+\frac{\rho}{2}\Delta
t\sum_{q=k-1}^{n+1}\|C_{k}(e^{q})\|^{2}\leq C\exp\big{(}(1-C\Delta
t)^{-1}T\big{)}(\Delta t)^{2k},\quad\forall n+1\leq\frac{T}{\Delta t},$ (5.22)

where $g_{k}$ is a positive constant depending only on $k$, $C$ is a constant
independent of $\Delta t$.

###### Proof.

We shall prove (5.21) and (5.22) by induction. Suppose we already have

$C_{k}^{\beta}(\phi^{n})\in\mathcal{B}_{\phi(t^{n+\beta})},\quad\forall n\leq
m,$ (5.23)

and (5.22) is satisfied with all $n\leq m-1$, we need to prove

$C_{k}^{\beta}(\phi^{m+1})\in\mathcal{B}_{\phi(t^{m+1+\beta})},$ (5.24)

and (5.22) is satisfied with all $n\leq m$.

Subtracting (5.18) with $m=n+\beta$ from (2.11) and multiplying by $\Delta t$,
we obtain

$A_{k}^{\beta}(e^{n+1})+\Delta t\mathcal{L}B_{k}^{\beta}(e^{n+1})=-\Delta
t\big{(}\mathcal{G}[C_{k}^{\beta}(\phi^{n})]-\mathcal{G}[\phi(t^{n+\beta})]\big{)}+E_{k}^{n+1}+\Delta
t\mathcal{L}R_{k}^{n+1}$ (5.25)

where $E_{k}^{n+1},\,R_{k}^{n+1}$ are given in (5.14). We split
$\mathcal{G}[C_{k}^{\beta}(\phi^{n})]-\mathcal{G}[\phi(t^{n+\beta})]$ as

$\begin{split}\mathcal{G}[C_{k}^{\beta}(\phi^{n})]-\mathcal{G}[\phi(t^{n+\beta})]&=\big{(}\mathcal{G}[C_{k}^{\beta}(\phi^{n})]-\mathcal{G}[C_{k}^{\beta}(\phi(t^{n}))]\big{)}+\big{(}\mathcal{G}[C_{k}^{\beta}(\phi(t^{n}))]-\mathcal{G}[\phi(t^{n+\beta})]\big{)}\\\
&=:T_{1}^{n}+T_{2}^{n}.\end{split}$ (5.26)

Taking the inner product of (5.25) with $C_{k}^{\beta}(e^{n+1})$, and
splitting $B_{k}^{\beta}(e^{n+1})$ as in (3.5), we obtain

$\begin{split}&\big{(}A_{k}^{\beta}(e^{n+1}),C_{k}^{\beta}(e^{n+1})\big{)}+\Delta
t\eta_{k}(\beta)\|C_{k}^{\beta}(e^{n+1})\|^{2}+\Delta
t\big{(}\mathcal{L}D_{k}^{\beta}(e^{n+1}),C_{k}^{\beta}(e^{n+1})\big{)}\\\
&=-\Delta t\big{(}T_{1}^{n},C_{k}^{\beta}(e^{n+1})\big{)}-\Delta
t\big{(}T_{2}^{n},C_{k}^{\beta}(e^{n+1})\big{)}+\big{(}E_{k}^{n+1},C_{k}^{\beta}(e^{n+1})\big{)}+\Delta
t\big{(}\mathcal{L}R_{k}^{n+1},C_{k}^{\beta}(e^{n+1})\big{)}.\end{split}$
(5.27)

Next, we bound the right hand side of (5.27) with the help of the consistency
estimate. First, it follows from (2.8) that with $\Delta t$ sufficiently
small, we have $C_{k}^{\beta}(\phi(t^{n}))\in\mathcal{B}_{\phi(t^{n+\beta})}$,
then for the terms with $T_{1}^{n}$ and $T_{2}^{n}$, it follows from (2.2) and
(5.23) that for any given $\varepsilon>0$,

$\big{|}\big{(}T_{1}^{n},C_{k}^{\beta}(e^{n+1})\big{)}\big{|}\leq\|T_{1}^{n}\|_{\star}\|C_{k}^{\beta}(e^{n+1})\|\leq\frac{\varepsilon}{2}(\gamma\|C_{k}^{\beta}(e^{n})\|^{2}+\mu|C_{k}^{\beta}(e^{n})|^{2})+\frac{1}{2\varepsilon}\|C_{k}^{\beta}(e^{n+1})\|^{2},$
(5.28)

With $P_{k}^{n}$ defined in (5.14), we have

$\begin{split}\big{|}\big{(}T_{2}^{n},C_{k}^{\beta}(e^{n+1})\big{)}\big{|}\leq\|T_{2}^{n}\|_{\star}\|C_{k}^{\beta}(e^{n+1})\|&\leq\frac{1}{\rho}(\gamma\|P_{k}^{n}\|^{2}+\mu|P_{k}^{n}|^{2})+\frac{\rho}{4}\|C_{k}^{\beta}(e^{n+1})\|^{2},\\\
&\leq C(\Delta
t)^{2k}+\frac{\rho}{4}\|C_{k}^{\beta}(e^{n+1})\|^{2}.\end{split}$ (5.29)

Similarly,

$\big{(}E_{k}^{n+1},C_{k}^{\beta}(e^{n+1})\big{)}\leq\frac{1}{2\Delta
t}|E_{k}^{n+1}|^{2}+\frac{\Delta t}{2}|C_{k}^{\beta}(e^{n+1})|^{2}\leq
C(\Delta t)^{2k+1}+\frac{\Delta t}{2}|C_{k}^{\beta}(e^{n+1})|^{2},$ (5.30)

and

$\big{(}\mathcal{L}R_{k}^{n+1},C_{k}^{\beta}(e^{n+1})\big{)}\leq\frac{1}{\rho}\|R_{k}^{n+1}\|^{2}+\frac{\rho}{4}\|C_{k}^{\beta}(e^{n+1})\|^{2}\leq
C(\Delta t)^{2k}+\frac{\rho}{4}\|C_{k}^{\beta}(e^{n+1})\|^{2}.$ (5.31)

Now, under the stability condition (5.19), combining the assumption on the
initial steps (5.20) and estimations in (5.28)-(5.31), taking
$\varepsilon=\frac{1}{\sqrt{\gamma}}$ in (5.28), and following the same
process as in the proof of Theorem 2 to handle the terms on the left hand side
of (5.27), we can obtain the following from (5.27):

$g_{k}|e^{n+1}|^{2}+\frac{\rho}{2}\Delta
t\sum_{q=k-1}^{n+1}\|C_{k}^{\beta}(e^{q})\|^{2}\leq C\Delta
t\sum_{q=0}^{n+1}|e^{q}|^{2}+C(\Delta t)^{2k},\quad\forall\,n\leq m.$ (5.32)

Therefore, by applying the discrete Gronwall lemma 2 to (5.32), we can obtain

$g_{k}|e^{m+1}|^{2}+\frac{\rho}{2}\Delta
t\sum_{q=k-1}^{m+1}\|C_{k}^{\beta}(e^{q})\|^{2}\leq C\exp\big{(}(1-C\Delta
t)^{-1}T\big{)}(\Delta t)^{2k},\quad\forall m+1\leq\frac{T}{\Delta t},$ (5.33)

with $C$ a constant independent of $\Delta t$ which implies (5.22). Finally,
it follows from (5.33) and (2.8) that

$\begin{split}\|C_{k}^{\beta}(\phi^{m+1})-\phi(t^{m+1+\beta})\|^{2}&\leq
2\|C_{k}^{\beta}(\phi^{m+1})-C_{k}^{\beta}(\phi(t^{m+1}))\|^{2}+2\|C_{k}^{\beta}(\phi(t^{m+1}))-\phi(t^{m+1+\beta})\|^{2}\\\
&\leq 2\|C_{k}^{\beta}(e^{m+1})\|^{2}+\mathcal{O}(\Delta t^{2k})\\\
&\leq\bar{C}\Delta t^{2k-1},\end{split}$ (5.34)

with $\bar{C}$ a constant independent of $\Delta t$, which implies (5.24) for
$\Delta t$ sufficiently small. Thus, the proof is complete with the induction.
∎

######  Remark 6.

Note that $\eta_{k}(\beta)$ in (3.6) monotonically increases as $\beta$
increases. On the other hand, for many applications, given $\delta>0$, one can
choose $\gamma<\delta$ with a suitable $\mu$ such that (2.2) is satisfied [6].
Hence, the stability condition (5.19) can always be satisfied with these
applications.

######  Remark 7.

The analysis in Theorem 2 and Theorem 3 can not be directly extended to the
standard BDF methods (with $\beta=1$) since $\eta_{k}(1)=0$.

###  5.4. Comparison to the classical BDF and IMEX schemes

In this subsection, we compare the stability condition (5.19) to that of the
classical BDF and IMEX methods (with Taylor expansion at time $t^{n+1}$) for
which the stability condition (5.19) does not apply. So we shall derive below
a corresponding stability condition for the classical BDF and IMEX methods. To
simplify the presentation, we assume $\mu=0$ in (2.2) since the general case
can be handled by applying the discrete Gronwall lemma as in Theorem 4.

The stability condition (5.19) in Theorem 4 is derived from

$\big{(}\mathcal{L}B_{k}^{\beta}(e^{n}),C_{k}^{\beta}(e^{n})\big{)}=\big{(}\eta_{k}(\beta)C_{k}^{\beta}(e^{n})+D_{k}^{\beta}(e^{n}),C_{k}^{\beta}(e^{n})\big{)}=\eta_{k}(\beta)\|C_{k}^{\beta}(e^{n})\|^{2}+\big{(}D_{k}^{\beta}(e^{n}),C_{k}^{\beta}(e^{n})\big{)},$
(5.35)

and

$\begin{split}\big{(}\mathcal{G}[C_{k}^{\beta}(\phi^{n})]-\mathcal{G}[C_{k}^{\beta}(\phi(t^{n}))],C_{k}^{\beta}(e^{n})\big{)}&\leq\mathop{\rm{min}}\limits_{\varepsilon>0}\big{(}\frac{\varepsilon}{2}\|\mathcal{G}[C_{k}^{\beta}(\phi^{n})]-\mathcal{G}[C_{k}^{\beta}(\phi(t^{n}))]\|_{\star}^{2}+\frac{1}{2\varepsilon}\|C_{k}(e^{n})\|^{2}\big{)}\\\
&\leq\mathop{\rm{min}}\limits_{\varepsilon>0}\big{(}\frac{\varepsilon\gamma}{2}\|C_{k}^{\beta}(e^{n})\|^{2}+\frac{1}{2\varepsilon}\|C_{k}^{\beta}(e^{n})\|^{2}\big{)}\\\
&\overset{\varepsilon=\frac{1}{\sqrt{\gamma}}}{=}\sqrt{\gamma}\|C_{k}^{\beta}(e^{n})\|^{2}.\end{split}$
(5.36)

As a result, the stability condition (5.19) is derived by requiring
$\eta_{k}(\beta)>\sqrt{\gamma}$ since the term
$\big{(}D_{k}^{\beta}(e^{n}),C_{k}^{\beta}(e^{n})\big{)}$ can be handled by
Lemma 1 and Theorem 1.

On the other hand, for the classical IMEX$k$ $(k=2,3,4)$ schemes, i.e., (2.11)
with $\beta=1$, the suitable multipliers are given as
$e^{n}-\tilde{\eta}_{k}e^{n-1}$ [25] and the smallest possible values of
$\tilde{\eta}_{k}$ are

$\tilde{\eta}_{2}=0,\quad\tilde{\eta}_{3}=0.0836,\quad\tilde{\eta}_{4}=0.2878.$
(5.37)

Hence, the corresponding versions of (5.36) and (5.35) become

$\begin{split}&\big{(}\mathcal{G}[\sum_{q=0}^{k-1}c_{k,q}(1)\phi^{n-k+1+q}]-\mathcal{G}[\sum_{q=0}^{k-1}c_{k,q}(1)\phi(t^{n-k+1+q})],e^{n}-\tilde{\eta}_{k}e^{n-1}\big{)}\\\
&\leq\mathop{\rm{min}}\limits_{\varepsilon>0}\big{(}\frac{\varepsilon}{2}\|\mathcal{G}[\sum_{q=0}^{k-1}c_{k,q}(1)\phi^{n-k+1+q}]-\mathcal{G}[\sum_{q=0}^{k-1}c_{k,q}(1)\phi(t^{n-k+1+q})]\|_{\star}^{2}+\frac{1}{2\varepsilon}\|e^{n}-\tilde{\eta}_{k}e^{n-1}\|^{2}\big{)}\\\
&\leq\mathop{\rm{min}}\limits_{\varepsilon>0}\big{(}\frac{\varepsilon\gamma}{2}\sum_{q=0}^{k-1}|c_{k,q}(1)|\|e^{n-k+1+q}\|^{2}+\frac{1}{2\varepsilon}\|e^{n}-\tilde{\eta}_{k}e^{n-1}\|^{2}\big{)}\\\
&\leq\mathop{\rm{min}}\limits_{\varepsilon>0}\big{(}\frac{\varepsilon\gamma}{2}\sum_{q=0}^{k-1}|c_{k,q}(1)|\|e^{n-k+1+q}\|^{2}+\frac{1}{2\varepsilon}(\|e^{n}\|^{2}+\tilde{\eta}_{k}^{2}\|e^{n-1}\|^{2})\big{)},\end{split}$
(5.38)

where $c_{k,q}(1)$ are defined in (2.12)-(2.14) with $\beta=1$, and

$\big{(}\mathcal{L}e^{n},e^{n}-\tilde{\eta}_{k}e^{n-1}\big{)}=\|e^{n}\|^{2}-\tilde{\eta}_{k}\big{(}\mathcal{L}e^{n},e^{n-1}\big{)}\geq\|e^{n}\|^{2}-\frac{\tilde{\eta}_{k}}{2}(\|e^{n}\|^{2}+\|e^{n-1}\|^{2}).$
(5.39)

Combining (5.38) and (5.39), we obtain the following stability condition for
the classical IMEX type scheme with multiplier
$e^{n}-\tilde{\eta}_{k}e^{n-1}$,

$1-\tilde{\eta}_{k}>\mathop{\rm{min}}\limits_{\varepsilon>0}\big{(}\frac{\varepsilon\gamma}{2}\sum_{q=0}^{k-1}|c_{k,q}(1)|+\frac{1}{2\varepsilon}(1+\tilde{\eta}_{k}^{2})\big{)}\geq\sqrt{\tilde{c}_{k}\gamma(1+\tilde{\eta}_{k}^{2})},$
(5.40)

with $\tilde{c}_{k}=\sum_{q=0}^{k-1}|c_{k,q}(1)|$. Comparing (5.19) with
(5.40), we have two remarks:

  * •

From (5.37) and (5.40), we observe that for the classical IMEX schemes,
higher-order (i.e., larger $k$) requires stronger stability condition on the
parameter $\gamma$ appearing in (2.2). It is this requirement on the time step
that limits the use of high order scheme in practice.

  * •

On the other hand, for the new class of IMEX schemes, we observe from (3.6)
and (5.19) that the stability condition on $\gamma$ becomes weaker as we
increase $\beta$. In particular, the new higher-order schemes with a suitable
$\beta$ can be stable with a larger time step than that is allowed with a
classical IMEX scheme of the same-order. For example, we have from (3.6) that
$\eta_{2}(2)=\eta_{3}(3)=\eta_{4}(5)=1/2$ which indicates that the stability
condition (5.19) of the new fourth-order scheme with $\beta=5$ and third-order
scheme with $\beta=3$ is the same as that of the second-order classical
scheme. Our numerical results in Example 3 below indicate that we can use the
maximum allowable time step of the second-order classical scheme in our new
third- and fourth-order schemes to obtain more accurate results.

######  Remark 8.

Note that a new multiplier $e^{n}-\frac{2}{169}e^{n-1}-\frac{11}{169}e^{n-2}$
for the classical BDF3 scheme is reported in [4] and since
$\hat{\eta}_{3}:=\frac{2}{169}+\frac{11}{169}<\tilde{\eta}_{3}=0.0836$, one
can obtain milder conditions on $\gamma$ compared to adopting the Nevanlinna-
Odeh multipliers. Nevertheless, we can derive even milder conditions on
$\gamma$ by choosing larger $\beta$ in our new methods.

##  6\. Extension to fifth-order

In Theorem 1, we found suitable multipliers for the second- and third-order
scheme with $\beta\geq 1$ and for the fourth-order scheme with $\beta\geq 2$.
In this section, we would like to show numerically that the multiplier we
found in section 3 also works for the fifth-order scheme.

Following the same notations as before, we can obtain the coefficients
$a_{5,q}(\beta),b_{5,q}(\beta),c_{5,q}(\beta)$ by solving the linear systems
(2.5), (2.7) and (2.9) with $k=5$, respectively. Then we can define
$A_{5}^{\beta}(\phi^{i}),B_{5}^{\beta}(\phi^{i}),C_{5}^{\beta}(\phi^{i})$ as
in (3.1). Next, we split $B_{5}^{\beta}(\phi^{n+1})$ as

$B_{5}^{\beta}(\phi^{n+1})=\eta_{5}(\beta)C_{5}^{\beta}(\phi^{n+1})+D_{5}^{\beta}(\phi^{n+1}),\quad\rm{with}\quad\eta_{5}(\beta)=\frac{\beta-1}{\beta+15},$
(6.1)

and define
$\tilde{A}_{5}^{\beta}(\zeta),\tilde{C}_{5}^{\beta}(\zeta),\tilde{D}_{5}^{\beta}(\zeta)$
as in (3.2). Following the key steps in the proof of Theorem 1, we present a
sequence of numerical results to show that $C_{5}^{\beta}(\phi^{n+1})$ is a
suitable multiplier for the fifth-order scheme with $6.5\leq\beta\leq 100$.

  * •

We have
$\rm{gcd}\big{(}\tilde{A}^{\beta}_{5}(\zeta),\zeta\tilde{C}^{\beta}_{5}(\zeta)\big{)}=\rm{gcd}\big{(}\tilde{D}^{\beta}_{5}(\zeta),\tilde{C}^{\beta}_{5}(\zeta)\big{)}=1$
since $\tilde{A}_{5}^{\beta}(0)=a_{5,0}\neq 0$ and

$\begin{split}&\rm{det}\,Sly(\tilde{A}^{\beta}_{5},\tilde{C}^{\beta}_{5})=\frac{\beta^{12}}{221184}+\frac{11\beta^{11}}{110592}+\frac{635\beta^{10}}{663552}+\frac{78937\beta^{9}}{14929920}+\frac{552809\beta^{8}}{29859840}+\frac{638383\beta^{7}}{14929920}\\\
&+\frac{9801769\beta^{6}}{149299200}+\frac{4912619\beta^{5}}{74649600}+\frac{765683\beta^{4}}{18662400}+\frac{225157\beta^{3}}{15552000}+\frac{6143\beta^{2}}{2488320}+\frac{2071\beta}{10368000}+\frac{1}{160000}>0,\end{split}$
(6.2)

and

$\rm{det}\,Sly(\tilde{D}^{\beta}_{5},\tilde{C}^{\beta}_{5})=\frac{\beta^{3}(\beta^{3}+6\beta^{2}+11\beta+6)^{3}}{13824}>0.$
(6.3)

  * •

Let $r_{1},r_{2},...,r_{5}$ be the five roots of
$\tilde{C}^{\beta}_{5}(\zeta)=0$, and denote $r_{\rm
max}=\mathop{\rm{max}}\limits_{1\leq i\leq 5}|r_{i}|$. In Fig. 5, we plot the
numerical values of $r_{\rm max}$ for $0\leq\beta\leq 100$. We observe that
$r_{\rm max}<1$ for $0\leq\beta\leq 100$, which implies
$\tilde{C}^{\beta}_{5}(\zeta)$ is holomorphic outside the unit disk in the
complex plane.

Figure 5. $r_{\rm max}$ with different $\beta$.

  * •

Following the same process as in the proof of Theorem 1, we can derive that
$\text{Re}\frac{\tilde{A}_{5}^{\beta}(\zeta)}{\zeta\tilde{C}_{5}^{\beta}(\zeta)}>0\,\text{for}\,|\zeta|>1$
is equivalent to

$\frac{1}{180}(1-y)f_{5}(y)\geq 0,\quad\forall y\in[-1,1],$

where

$f_{5}(y)=\sigma_{4}(\beta)y^{4}+\sigma_{3}(\beta)y^{3}+\sigma_{2}(\beta)y^{2}+\sigma_{1}(\beta)y+\sigma_{0}\geq
0,\quad\forall y\in[-1,1],$ (6.4)

with

$\displaystyle\sigma_{4}(\beta)=5\beta^{8}+70\beta^{7}+390\beta^{6}+1090\beta^{5}+1539\beta^{4}+820\beta^{3}-350\beta^{2}-540\beta-144,$
(6.5a)
$\displaystyle\sigma_{3}(\beta)=-20\beta^{8}-280\beta^{7}-1550\beta^{6}-4260\beta^{5}-5836\beta^{4}-3024\beta^{3}+950\beta^{2}+1396\beta+336,$
(6.5b)
$\displaystyle\sigma_{2}(\beta)=30\beta^{8}+420\beta^{7}+2310\beta^{6}+6240\beta^{5}+8244\beta^{4}+3932\beta^{3}-1260\beta^{2}-1340\beta-204,$
(6.5c)
$\displaystyle\sigma_{1}(\beta)=-20\beta^{8}-280\beta^{7}-1530\beta^{6}-4060\beta^{5}-5136\beta^{4}-2072\beta^{3}+1070\beta^{2}+652\beta+36,$
(6.5d)
$\displaystyle\sigma_{0}(\beta)=5\beta^{8}+70\beta^{7}+380\beta^{6}+990\beta^{5}+1189\beta^{4}+344\beta^{3}-410\beta^{2}-168\beta+336.$
(6.5e)

  * •

On the other hand, we can also show that
$\text{Re}\frac{\tilde{D}_{5}^{\beta}(\zeta)}{\tilde{C}_{5}^{\beta}(\zeta)}>0\,\text{for}\,|\zeta|>1$
is equivalent to

$h_{5}(y)=\mu_{4}(\beta)y^{4}+\mu_{3}(\beta)y^{3}+\mu_{2}(\beta)y^{2}+\mu_{1}(\beta)y+\mu_{0}(\beta)\geq
0,\quad\forall y\in[-1,1],$ (6.6)

with

$\displaystyle\mu_{4}(\beta)=\frac{\beta(\beta^{2}+3\beta+2)^{2}(6\beta^{3}+37\beta^{2}+48\beta-27)}{18(\beta+15)},$
(6.7a)
$\displaystyle\mu_{3}(\beta)=-\frac{\beta(24\beta^{7}+292\beta^{6}+1366\beta^{5}+3013\beta^{4}+2881\beta^{3}+193\beta^{2}-1391\beta-618)}{18(\beta+15)},$
(6.7b)
$\displaystyle\mu_{2}(\beta)=\frac{\beta(12\beta^{7}+146\beta^{6}+670\beta^{5}+1385\beta^{4}+1021\beta^{3}-553\beta^{2}-1127\beta-402)}{6(\beta+15)},$
(6.7c)
$\displaystyle\mu_{1}(\beta)=-\frac{(24\beta^{8}+292\beta^{7}+1314\beta^{6}+2527\beta^{5}+1203\beta^{4}-2405\beta^{3}-3117\beta^{2}-1008\beta-270)}{18(\beta+15)},$
(6.7d)
$\displaystyle\mu_{0}(\beta)=\frac{6\beta^{8}+73\beta^{7}+322\beta^{6}+571\beta^{5}+91\beta^{4}-926\beta^{3}-995\beta^{2}-312\beta+18}{18(\beta+15)}.$
(6.7e)

In Fig. 6, we plot the minimum values of $f_{5}(y)$ and $h_{5}(y)$ in $[-1,1]$
with $1\leq\beta\leq 100$, which show (6.4) is true for $1\leq\beta\leq 100$
and (6.6) is true for $6.5\leq\beta\leq 100$. Therefore, we have numerically
verified that Theorem 1 is also true for (2.11) with $k=5$ and
$6.5\leq\beta\leq 100$.

######  Remark 9.

The choice of $\eta_{5}(\beta)$ in (6.1) is not unique, and the range
$6.5\leq\beta\leq 100$ is not necessarily the largest possible. But our
numerical results indicate (6.4) and (6.6) do not hold for some $\beta>100$.

For the sixth-order scheme, our numerical results show there exists
$|r_{6}|>1$, which is one root of $\tilde{C}^{\beta}_{6}(\zeta)=0$ and this
implies that it is not holomporphic outside the unit disk. Hence, the proof in
Theorem 1 can not be extended to the sixth-order.

Figure 6. Minimum value of $f_{5}$ and $h_{5}$ in $[-1,1]$ with different
$\beta$.

##  7\. Numerical examples

In this section, we provide some numerical approximation of the Allen-Cahn [7]
and Cahn-Hilliard [11] equations to validate our theoretical results, and to
show the advantages of the new IMEX schemes (2.11).

Given a free energy

$\mathcal{E}[\phi]=\int\frac{1}{2}|\nabla\phi|^{2}+\frac{1}{4\varepsilon^{2}}(1-\phi^{2})^{2}d\bm{x}.$
(7.1)

We consider the $H^{-\alpha}$ gradient flow,

$\frac{\partial\phi}{\partial
t}=-m(-\Delta)^{\alpha}\big{(}-\Delta\phi-\frac{1}{\varepsilon^{2}}\phi(1-\phi^{2})\big{)}+f(t),\quad\alpha=0\quad\rm{or}\quad
1,$ (7.2)

where $f$ is the given source term. When $\alpha=0$, (7.2) is the standard
Allen-Cahn equation; when $\alpha=1$, it becomes the standard Cahn-Hilliard
equation.

Example 1. In the first example, we validate the convergence order of the new
schemes. Considering a two-dimensional domain  $(0,2)^{2}$ with periodic
boundary conditions, let $\alpha=0$, $m=\varepsilon=0.2$ in (7.2) and $f$ is
chosen such that the exact solution of (7.2) is

$\phi(x,y,t)=e^{\sin(\pi x)\sin(\pi y)}\sin(t).$ (7.3)

We use the Fourier Galerkin method with $Nx=Ny=40$ in space so that the
spatial discretization error is negligible compared to the time discretization
error. In Fig. 7, we plot the convergence rate of the $L^{2}$ error at $T=1$
by using the second- to fourth- order schemes (2.11). We observe the expected
convergence order for all the cases with different $\beta$. We also observe
that for the same order, the error increases slightly with larger $\beta$.

Figure 7.  Convergence test for the general IMEX type methods. From left to
right: second order, third order and fourth order schemes with different
$\beta$.

Example 2. In the second example, we solve a benchmark problem for the Allen-
Cahn equation [12]. Consider a two-dimensional domain $(-128,128)^{2}$ with a
circle of radius $R_{0}=100$. In other words, the initial condition is given
as

$\phi(x,y,0)=\left\\{\begin{aligned} &1,\quad x^{2}+y^{2}<100^{2},\\\
&-1,\quad x^{2}+y^{2}\geq 100^{2}.\end{aligned}\right.$ (7.4)

By mapping the domain to $(-1,1)^{2}$, the parameters in (7.2) are given by
$m=6.10351\times 10^{-5}$, $\varepsilon=0.0078$, $\alpha=0$ and $f=0$. In the
sharp interface limit, the radius at time $t$ is given by

$R=\sqrt{R_{0}^{2}-2t}.$ (7.5)

We use the Fourier Galerkin method with $Nx=Ny=512$ in space. Then we fix
$\Delta t=0.75$, which is the maximum time step we can use for the classical
second-order scheme to get acceptable numerical results, and use (2.11) with
different orders and different $\beta$. We plot the computed radius $R(t)$ in
Fig. 8, which shows that we can use higher-order schemes with the same large
time step as the second-order schemes by choosing $\beta>1$. More importantly,
we can get much more accurate results with higher-order schemes. Here,
$k=1,\beta=1$ represents the usual first-order scheme.

Figure 8.  The evolution of radius $R$ with $\Delta t=0.75$ under different
schemes.

Example 3. In the third example, we consider the Cahn-Hilliard equation in a
two-dimensional domain $(0,1)^{2}$ with periodic boundary condition and let
$\alpha=1$, $m=1$, $\varepsilon=0.02$ in (7.2). The initial condition is given
as $\phi(0)=0.2+r$ and $r$ is a random perturbation variable with uniform
distribution in $[-0.02,0.02]$. We use the Fourier Galerkin method with
$Nx=Ny=128$ in space. In Fig. 9, we compare the first- to the fourth-order
schemes with different $\beta$, the reference solution is generated by using
the classical fourth-order scheme with sufficiently small time step $\Delta
t=5\times 10^{-9}$.

Several observations are in order:

  * •

1\. We take $\Delta t=7.5\times 10^{-8}$ which is the maximum allowable time
step for the classical second-order scheme, and observe in Fig. 9(a) that we
can use the same time step for the higher-order schemes by choosing a suitable
$\beta>1$, and obtain more accurate results.

  * •

2\. We observe in Fig. 9(b) that the usual third- and fourth-order schemes
with $\beta=1$ are unstable, but we can get correct solutions with the third-
and fourth-order schemes by choosing a suitable $\beta>1$.

  * •

3\. We also observe in Fig. 9(b) that $\beta$ too large may lead to inaccurate
results due to larger truncation errors.

(a) Comparisons of different order schemes

(b) Effect of $\beta$ in high order schemes

Figure 9.  Comparisons of different order schemes with different $\beta$ for
the Cahn-Hilliard equation

Change the Table below to a figure.

| $t=2.5\times 10^{-4}$ | $t=5\times 10^{-4}$ | $t=2\times 10^{-3}$
---|---|---|---
Reference |  |  |
IMEX2, $\beta=2$, $\Delta t=10^{-7}$ |  |  |
IMEX3, $\beta=3$, $\Delta t=1.5\times 10^{-7}$ |  |  |
IMEX4, $\beta=2.3$, $\Delta t=7.5\times 10^{-8}$ |  |  |
Table 2. Example 3. Snapshots of $\phi$ at different time by different
schemes.

##  8\. Concluding remarks

We presented in this paper a new class of BDF and IMEX schemes for parabolic
type equations based on the Taylor expansion at time $t^{n+\beta}$ with
$\beta>1$ being a tunable parameter. The new schemes are a simple
generalization of the classical BDF or IMEX schemes with essentially the same
computational efforts. However, they enjoy a remarkable property that their
stability regions increase as the parameter $\beta$ increases, making it
possible, by choosing a suitably large $\beta$, to use high-order schemes with
larger time steps that are only allowed with lower-order classical schemes. We
also identified an explicit uniform multiplier for the new schemes of second-
to fourth-order, and carried out a rigorous stability and error analysis by
using the energy argument. We also presented numerical examples to show the
benefit of using higher-order schemes with a suitable $\beta>1$.

This class of new BDF and IMEX schemes makes it possible to use higher-order
schemes for highly stiff systems with reasonably large time steps, and can be
easily implemented with a minimal effort by modifying the code based on the
classical BDF or IMEX schemes. Thus, it provides a much needed improvement on
the stability of higher-order schemes. The idea behind the new class of BDF
and IMEX schemes is very simple but original, and can be extended to other
type of numerical schemes.

## References

  * [1] A. G. Akritas.  Sylvester’s forgotten form of the resultant.  Fibonacci Quart, 31(4):325–332, 1993.
  * [2] G. Akrivis.  Stability of implicit-explicit backward difference formulas for nonlinear parabolic equations.  SIAM J. Numer. Anal., 53(1):464–484, 2015.
  * [3] G. Akrivis, M. Chen, F. Yu, and Z. Zhou.  The energy technique for the six-step BDF method.  SIAM J. Numer. Anal., 59(5):2449–2472, 2021.
  * [4] G. Akrivis and E. Katsoprinakis.  Backward difference formulae: new multipliers and stability properties for parabolic equations.  Math. Comput., 85(301):2195–2216, 2016.
  * [5] G. Akrivis, B. Li, and C. Lubich.  Combining maximal regularity and energy estimates for time discretizations of quasilinear parabolic equations.  Math. Comput., 86(306):1527–1552, 2017.
  * [6] G. Akrivis and C. Lubich.  Fully implicit, linearly implicit and implicit–explicit backward difference formulae for quasi-linear parabolic equations.  Numer. Math., 131:713–735, 2015.
  * [7] S. M. Allen and J. W. Cahn.  A microscopic theory for antiphase boundary motion and its application to antiphase domain coarsening.  Acta metallurgica, 27(6):1085–1095, 1979.
  * [8] U. M. Ascher, S. J. Ruuth, and B. T. Wetton.  Implicit-explicit methods for time-dependent partial differential equations.  SIAM J. Numer. Anal., 32(3):797–823, 1995.
  * [9] J. C. Butcher.  The numerical analysis of ordinary differential equations: Runge-Kutta and general linear methods.  Wiley-Interscience, 1987.
  * [10] J. C. Butcher.  Numerical methods for ordinary differential equations.  John Wiley & Sons, 2016.
  * [11] J. W. Cahn and J. E. Hilliard.  Free energy of a nonuniform system. i. interfacial free energy.  J. Chem. Phys, 28(2):258–267, 1958.
  * [12] L. Chen and J. Shen.  Applications of semi-implicit fourier-spectral method to phase field equations.  Comput. Phys. Commun., 108(2-3):147–158, 1998.
  * [13] G. Dahlquist.  G-stability is equivalent to A-stability.  BIT Numer. Math., 18:384–401, 1978.
  * [14] M. Gunzburger, X. He, and B. Li.  On Stokes-Ritz projection and multistep backward differentiation schemes in decoupling the Stokes-Darcy model.  SIAM J. Numer. Anal., 56(1):397–427, 2018.
  * [15] E. Hairer, M. Hochbruck, A. Iserles, and C. Lubich.  Geometric numerical integration.  Oberwolfach Rep., 3(1):805–882, 2006.
  * [16] D. J. Higham and L. N. Trefethen.  Stiffness of ODEs.  BIT Numer. Math., 33:285–303, 1993.
  * [17] F. Huang and J. Shen.  Stability and error analysis of a second order consistent splitting scheme for the Navier-Stokes equation.  SIAM J. Numer. Anal., 61(5):2408–2433, 2023.
  * [18] W. H. Hundsdorfer and J. G. Verwer.  Numerical solution of time-dependent advection-diffusion-reaction equations, volume 33.  Springer, 2003.
  * [19] A. Iserles.  A first course in the numerical analysis of differential equations.  Number 44. Cambridge university press, 2009.
  * [20] J. Jia and J. Huang.  Krylov deferred correction accelerated method of lines transpose for parabolic problems.  J. Comput. Phys., 227(3):1739–1753, 2008.
  * [21] A. K. Kassam and L. N. Trefethen.  Fourth-order time-stepping for stiff PDEs.  SIAM J. Sci. Comput., 26(4):1214–1233, 2005.
  * [22] B. Li, K. Wang, and Z. Zhou.  Long-time accurate symmetrized implicit-explicit BDF methods for a class of parabolic equations with non-self-adjoint operators.  SIAM J. Numer. Anal., 58(1):189–210, 2020.
  * [23] C. Lubich.  On the convergence of multistep methods for nonlinear stiff differential equations.  Numer. Math., 58:839–853, 1990.
  * [24] C. Lubich, D. Mansour, and C. Venkataraman.  Backward difference time discretization of parabolic differential equations on evolving surfaces.  IMA J. Numer. Anal., 33(4):1365–1385, 2013.
  * [25] O. Nevanlinna and F. Odeh.  Multiplier techniques for linear multistep methods.  Numer. Funct. Anal. Optim., 3(4):377–423, 1981.
  * [26] A. Quarteroni and A. Valli.  Numerical approximation of partial differential equations, volume 23.  Springer Science & Business Media, 2008.
  * [27] W. E. Schiesser.  The numerical method of lines: integration of partial differential equations.  Elsevier, 2012.