Average energy dissipation rates of additive implicit-explicit Runge-Kutta methods for gradient flow problems

Hong-lin Liao Xuping Wang Cao Wen ORCID 0000-0003-0777-6832. School of Mathematics, Nanjing University of Aeronautics and Astronautics, Nanjing 211106, China; Key Laboratory of Mathematical Modeling and High Performance Computing of Air Vehicles (NUAA), MIIT, Nanjing 211106, China. Emails: [email protected] and [email protected]. This author’s work is supported by NSF of China under grant number 12071216.School of Mathematics, Nanjing University of Aeronautics and Astronautics, Nanjing 211106, China. Email: [email protected] of Mathematics, Nanjing University of Aeronautics and Astronautics, Nanjing 211106, China. Email: [email protected].

Abstract

A unified theoretical framework is suggested to examine the energy dissipation properties at all stages of additive implicit-explicit Runge-Kutta (IERK) methods up to fourth-order accuracy for gradient flow problems. We construct some parameterized IERK methods by applying the so-called first same as last method, that is, the diagonally implicit Runge-Kutta method with the explicit first stage and stiffly-accurate assumption for the linear stiff term, and applying the explicit Runge-Kutta method for the nonlinear term. The main part of the novel framework is to construct the differential forms and the associated differentiation matrices of IERK methods by using the difference coefficients of method and the so-called discrete orthogonal convolution kernels. As the main result, we prove that an IERK method can preserve the original energy dissipation law unconditionally if the associated differentiation matrix is positive semi-definite. The recent indicator, namely average energy dissipation rate, is also adopted for these multi-stage methods to evaluate the overall energy dissipation rate of an IERK method such that one can choose proper parameters in some parameterized IERK methods. It is found that the selection of method parameters in the IERK methods is at least as important as the selection of different IERK methods. Extensive numerical experiments are also included to support our theory.
Keywords: gradient flow problem, implicit-explicit Runge-Kutta methods, orthogonal convolution kernels, energy dissipation law, average dissipation rate
AMS subject classifications: 35K58, 65L20, 65M06, 65M12

1 Introduction

We propose a unified theoretical framework to examine the energy dissipation properties at all stages of additive implicit-explicit Runge-Kutta (in short, IERK) methods for solving the following semi-discrete semilinear parabolic problem, cf. [10, 11, 24, 32],

\displaystyle u_{h}^{\prime}(t)=M_{h}\left[L_{h}u_{h}(t)-g_{h}(u_{h}(t))\right],\quad u_{h}(t_{0})=u_{h}^{0},

(1.1)

where $L_{h}$ is a symmetric, positive definite matrix resulting from certain spatial discretization of stiff term, typically the Laplacian operator $-\Delta$ with periodic boundary conditions, $M_{h}$ is a symmetric, negative definite matrix (having the same size as $L_{h}$ ) resulting from certain spatial discretization of the mobility operator, and $g_{h}$ represents a nonlinear but non-stiff term. Without losing the generality, the Fourier pseudo-spectral method is assumed to approximate spatial operators and we define the discrete $L^{2}$ inner product $\left\langle u,v\right\rangle:=v^{T}u$ and the $L^{2}$ norm $\left\|v\right\|:=\sqrt{\left\langle v,v\right\rangle}$ . Assume that there exists a non-negative Lyapunov function $G_{h}$ such that $g_{h}(v)=-\frac{\delta}{\delta v}G_{h}(v).$ Then the problem (1.1) can be formulated into the following gradient flow system, see [1, 8, 12, 29, 30, 31],

\displaystyle\frac{\,\mathrm{d}u_{h}}{\,\mathrm{d}t}=M_{h}\frac{\delta E}{\delta u_{h}}\quad\text{with}\quad E[v_{h}]:=\frac{1}{2}\langle v_{h},L_{h}v_{h}\rangle+\langle G_{h}(v_{h}),1\rangle.

(1.2)

The dynamics approaching the steady-state solution $u_{h}^{*}$ , that is $L_{h}u_{h}^{*}=g_{h}(u_{h}^{*})$ , of this dissipative system (1.1) satisfies the following original energy dissipation law

\displaystyle\frac{\,\mathrm{d}E}{\,\mathrm{d}t}=\Big{\langle}\frac{\delta E}{\delta u_{h}},\frac{\,\mathrm{d}u_{h}}{\,\mathrm{d}t}\Big{\rangle}=\Big{\langle}M_{h}^{-1}\frac{\,\mathrm{d}u_{h}}{\,\mathrm{d}t},\frac{\,\mathrm{d}u_{h}}{\,\mathrm{d}t}\Big{\rangle}\leq 0.

(1.3)

Our recent interests in the Runge-Kutta (RK) methods are firstly to develop some high-order starting procedures for the backward differentiation formulas [22, 25], which are proven to possess certain discrete gradient structures, cf. [32, Section 5.6], and then preserve the energy dissipation law of phase field models with the original energy added with a small positive term in their original forms. Thus, it would be expected that the RK methods also have some discrete gradient structures and preserve certain discrete energy dissipation laws for gradient flow problems. Combining the scalar auxiliary variable (SAV) method [29] and the Gaussian-type RK methods, arbitrarily high-order energy-stable schemes were constructed [1, 12]. However, these SAV schemes only satisfy a modified energy law, which may not necessarily ensure the energy stability in the original form [17, 11]. Recently, Fu et al. [10, 11] proved that some exponential time difference RK methods maintain the decaying of original energy. For some parameterized explicit exponential RK methods in [15], the discrete version of (1.3) was established in [24] for certain range of the method parameters. Since the functions of matrix exponential are always involved, these explicit RK methods would be also computationally intensive. Actually, the efficient algorithm to accurately compute the matrix exponential is still limited [27, 9].

In this study, we will focus on whether and to what extent the IERK methods preserve the original energy dissipation law (1.3). Due to the relative ease of implementation, the diagonally implicit Runge-Kutta (DIRK) methods are possibly the most widely used implicit Runge-Kutta method in practical applications, see the recent review [21] and abundant references therein, especially involving stiff differential equations. The DIRK method has greatly reduced the computational complexity of the fully implicit Runge-Kutta method, but it still requires an iterative method to solve the nonlinear equations at each stage. In order to improve the computational efficiency, the IERK methods have attracted much attention and are widely used, see [2, 4, 3, 5, 6, 7, 16, 13, 18, 20, 19, 26, 31, 30]. Cooper and Sayfy [6] presented some IERK methods up to fourth-order accuracy with the implicit part is A-stable. Kennedy and Carpenter [20] constructed high-order IERK methods from third- to fifth-order to simulate convection-diffusion-reaction equations. The widespread ARS-type IERK methods by Ascher, Ruuth and Spiteri were developed in [2] to solve the convection-diffusion problems. These ARS-type methods had better stability regions than implicit-explicit multi-step schemes over a wide parameter range. Some ARS-type IERK methods from second- to fourth-order were also proposed in [26] for nonlinear differential equations with constraints, such as the Navier-Stokes equation. Cardone et.al. [5] proposed a class of ARS-type IERK methods up to fourth-order based on extrapolation of the stage values at the current step by stage values in the previous step. Izzo and Jackiewicz [16] constructed some parameterized IERK methods up to fourth-order by choosing the method parameters to maximize the regions of absolute stability for the explicit part in the method with the assumption of A-stable implicit part.

We will consider another class of parameterized IERK methods and choose the method parameters to maintain the original energy dissipation law (1.3) as much as possible. Let $u_{h}^{k}$ be the numerical approximation of $u_{h}(t_{k})$ at the mesh point $t_{k}$ for $0\leq k\leq N$ . For a $s$ -stage Runge-Kutta method, let $U^{n,i}$ be the approximation of $u_{h}(t_{n-1}+c_{i}\tau)$ at the abscissas $c_{1}:=0$ , $c_{i}>0$ for $2\leq i\leq s-1$ , and $c_{s}:=1$ . To integrate the semilinear parabolic problem (1.1) from the time $t_{n-1}$ ( $n\geq 1$ ) to the next grid point $t_{n}=t_{n-1}+\tau$ (typically, $\tau$ also represents a variable-step size), one can construct the following $s$ -stage IERK methods with explicit first stage, see [20, 26, 33]


	$\displaystyle~{}U^{n,i}:=U^{n,1}+\tau\sum_{j=1}^{i}a_{i,j}M_{h}L_{h}U^{n,j}-\tau\sum_{j=1}^{i-1}\hat{a}_{i,j}M_{h}g_{h}(U^{n,j}),\quad\text{$2\leq i\leq s-1$,}$		(1.4a)
	$\displaystyle~{}U^{n,s}:=U^{n,1}+\tau\sum_{j=1}^{s}b_{j}M_{h}L_{h}U^{n,j}-\tau\sum_{j=1}^{s-1}\hat{b}_{j}M_{h}g_{h}(U^{n,j}),$		(1.4b)

where $U^{n,1}:=u_{h}^{n-1}$ and $u_{h}^{n}:=U^{n,s}$ . The implicit part of IERK methods (1.4) approximating the stiff linear term $L_{h}$ can be represented by the abscissa vector $\mathbf{c}$ , the coefficient matrix $A$ and the vector of weights $\mathbf{b}$ ,

\begin{array}[]{c|c}\mathbf{c}&A\\ \hline\cr\\[-8.0pt] &\mathbf{b}^{T}\end{array}=\begin{array}[]{c|ccccc}c_{1}=0&0&&&&\\ c_{2}&a_{21}&a_{22}&&&\\ c_{3}&a_{31}&a_{32}&a_{33}&&\\ \vdots&\vdots&\vdots&\ddots&\ddots&\\[2.0pt] c_{s}=1&a_{s,1}&a_{s,2}&\cdots&a_{s,s-1}&a_{s,s}\\[2.0pt] \hline\cr&a_{s,1}&a_{s,2}&\cdots&a_{s,s-1}&a_{s,s}\end{array},

where the abscissa $c_{i}=\sum_{j=1}^{i}a_{i,j}$ for $1\leq i\leq s$ , and we use the so-called first same as last methods [28, 21], that is, the DIRK methods with the explicit first stage and the stiffly-accurate assumption: $b_{j}=a_{s,j}$ for $1\leq j\leq s$ such that

c_{s}=\sum_{j=1}^{s}a_{s,j}=\sum_{j=1}^{s}b_{j}=1.

That is, we consider Lobatto-type DIRK methods, with the number of implicit stage $s_{\mathrm{I}}:=s-1,$ which differs from the widespread Radau-type DIRK methods with $a_{11}\neq 0$ or $a_{i1}=0$ for $2\leq i\leq s$ .

The explicit part of IERK methods (1.4), namely explicit Runge-Kutta methods, for the nonlinear term $g_{h}$ can be represented by the abscissa vector $\hat{\mathbf{c}}$ , the strictly lower triangular coefficient matrix $\widehat{A}$ and the vector of weights $\hat{\mathbf{b}}$ ,

\begin{array}[]{c|c}\hat{\mathbf{c}}&\widehat{A}\\ \hline\cr\\[-8.0pt] &\hat{\mathbf{b}}^{T}\end{array}=\begin{array}[]{c|ccccc}\hat{c}_{1}=0&0&&&&\\ \hat{c}_{2}&\hat{a}_{21}&0&&&\\ \hat{c}_{3}&\hat{a}_{31}&\hat{a}_{32}&0&&\\ \vdots&\vdots&\vdots&\ddots&\ddots&\\[2.0pt] \hat{c}_{s}=1&\hat{a}_{s,1}&\hat{a}_{s,2}&\cdots&\hat{a}_{s,s-1}&0\\[1.0pt] \hline\cr\\[-8.0pt] &\hat{a}_{s,1}&\hat{a}_{s,2}&\cdots&\hat{a}_{s,s-1}&0\end{array},

where the abscissa $\hat{c}_{i}=\sum_{j=1}^{i-1}\hat{a}_{i,j}$ for $1\leq i\leq s$ , and introduce, to simplify our notations,

\displaystyle\hat{b}_{j}:=\hat{a}_{s,j},\quad 1\leq j\leq s-1\quad\text{such that}\quad\hat{c}_{s}=\sum_{j=1}^{s-1}\hat{a}_{s,j}=1.

(1.5)

Always, we assume that $\hat{a}_{k+1,k}\neq 0$ for any $1\leq k\leq s-1$ . The IERK methods (1.4) become

\displaystyle~{}U^{n,i}=U^{n,1}+\tau\sum_{j=1}^{i}a_{i,j}M_{h}L_{h}U^{n,j}-\tau\sum_{j=1}^{i-1}\hat{a}_{i,j}M_{h}g_{h}(U^{n,j})\quad\text{for $2\leq i\leq s$.}

(1.6)

Since the Lobatto-type DIRK methods are used in the implicit part, we call (1.6) as Lobatto-type IERK methods. They are quite different from the Radau-type IERK method containing the Radau-type DIRK methods in the implicit part, which are also known as ARS-type IERK methods [2, 11, 16, 31, 13, 26].

In general, the consistency of IERK methods (1.6) requires the canopy node condition

\displaystyle c_{i}=\hat{c}_{i}\quad\text{or}\quad\sum_{j=1}^{i}a_{i,j}=\sum_{j=1}^{i-1}\hat{a}_{i,j}\quad\text{for $1\leq i\leq s$},

(1.7)

which makes the form (1.6) invariant under the transformation of IERK method to the non-autonomous system. Under the canopy node condition $\hat{\mathbf{c}}={\mathbf{c}}$ , Table 1 lists the order conditions for the coefficient matrices and the weights vector to make the IERK method (1.6) accurate up to fourth-order. A detailed description of these order conditions can also be found in [2, 6, 16, 31].

Table 1: Order conditions for IERK methods up to fourth-order accuracy.

Order	Stand-alone conditions		Coupling condition
Order	Implicit part	Explicit part
1	$\mathbf{b}^{T}\mathbf{1}=1$	$\hat{\mathbf{b}}^{T}\mathbf{1}=1$	-
2	$\mathbf{b}^{T}\mathbf{c}=\tfrac{1}{2}$	$\hat{\mathbf{b}}^{T}\mathbf{c}=\tfrac{1}{2}$	-
3	$\mathbf{b}^{T}\mathbf{c}^{.2}=\tfrac{1}{3}$	$\hat{\mathbf{b}}^{T}\mathbf{c}^{.2}=\tfrac{1}{3}$
	$\mathbf{b}^{T}A\mathbf{c}=\tfrac{1}{6}$	$\hat{\mathbf{b}}^{T}\widehat{A}\mathbf{c}=\tfrac{1}{6}$	$\mathbf{b}^{T}\widehat{A}\mathbf{c}=\tfrac{1}{6},\;\hat{\mathbf{b}}^{T}A\mathbf{c}=\tfrac{1}{6}$
	$\mathbf{b}^{T}\mathbf{c}^{.3}=\frac{1}{4}$	$\hat{\mathbf{b}}^{T}\mathbf{c}^{.3}=\frac{1}{4}$
	$\mathbf{b}^{T}[\mathbf{c}\odot(A\mathbf{c})]=\frac{1}{8}$	$\hat{\mathbf{b}}^{T}[\mathbf{c}\odot(\widehat{A}\mathbf{c})]=\frac{1}{8}$	$\mathbf{b}^{T}[\mathbf{c}\odot(\widehat{A}\mathbf{c})]=\frac{1}{8},\;\hat{\mathbf{b}}^{T}[\mathbf{c}\odot(A\mathbf{c})]=\frac{1}{8}$
4	$\mathbf{b}^{T}A\mathbf{c}^{.2}=\frac{1}{12}$	$\hat{\mathbf{b}}^{T}\widehat{A}\mathbf{c}^{.2}=\frac{1}{12}$	$\mathbf{b}^{T}\widehat{A}\mathbf{c}^{.2}=\frac{1}{12},\;\hat{\mathbf{b}}^{T}A\mathbf{c}^{.2}=\frac{1}{12}$
	$\mathbf{b}^{T}A^{2}\mathbf{c}=\frac{1}{24}$	$\hat{\mathbf{b}}^{T}\widehat{A}^{2}\mathbf{c}=\frac{1}{24}$	$\mathbf{b}^{T}\widehat{A}^{2}\mathbf{c}=\frac{1}{24},\;\hat{\mathbf{b}}^{T}A\widehat{A}\mathbf{c}=\frac{1}{24},\;\hat{\mathbf{b}}^{T}\widehat{A}A\mathbf{c}=\frac{1}{24}$
			$\mathbf{b}^{T}A\widehat{A}\mathbf{c}=\frac{1}{24},\;\mathbf{b}^{T}\widehat{A}A\mathbf{c}=\frac{1}{24},\;\hat{\mathbf{b}}^{T}A^{2}\mathbf{c}=\frac{1}{24}$

^{* For the vectors $\mathbf{x}=(x_{1},x_{2},\cdots,x_{s})^{T}$ and $\mathbf{y}=(y_{1},y_{2},\cdots,y_{s})^{T}$ , $\mathbf{x}\odot\mathbf{y}:=(x_{1}y_{1},x_{2}y_{2},\cdots,x_{s}y_{s})^{T}$ and $\mathbf{x}^{.m}:=\mathbf{x}\odot\mathbf{x}^{.(m-1)}$ for $m>1$ .}

In simulating the semilinear parabolic problems (1.1) and related gradient flow problems (1.2), IERK methods turned out to be very competitive. Shin et al. [30] observed that the Radau-type IERK methods combined with convex splitting technique (CSRK) exhibit the original energy stability in numerical experiments. Subsequently, they proved [31] that the proposed CSRK schemes can unconditionally maintain the original energy stability by assuming $\hat{a}_{ij}=a_{i,j+1}$ in (1.4). The proposed second-order CSRK scheme requires at least four stages, while the third-order CSRK scheme requires at least seven stages. The proposed CSRK schemes in [31] are nonlinear and would be computationally expensive due to the required inner iterations at each time stage. Very recently, Fu et al. [11] derived some sufficient conditions of Radau-type IERK methods to maintain the decay of original energy for the gradient flow problems and presented some concrete schemes up to third-order accuracy in time. It is noted that, the Radau-type IERK methods are only focused in the mentioned works, while we will construct some Lobatto-type IERK methods and present a unified theoretical framework to examine the energy dissipation properties of IERK methods up to fourth-order accurate for gradient flow problems (1.2).

The unified theoretical framework, presented in the next section, is inspired by two aspects: one is the idea by treating the stiffly-accurate (but not necessarily algebraically stable) DIRK method as the composite linear multistep method; the other is the recent discrete energy technique for the stability and convergence of multistep backward differentiation formulas [22, 25]. We will construct the differential forms and the associated differentiation matrices of IERK methods by using the difference coefficients of method and the so-called discrete orthogonal convolution kernels. It is proven that an IERK method can preserve the original energy dissipation law unconditionally if the associated differentiation matrix is positive semi-definite, see Theorem 2.1. The recent indicator in [24], namely average energy dissipation rate, is also defined for the IERK methods to evaluate the overall energy dissipation rate when applied to the gradient flow problems (1.2) such that one can choose proper parameters in some parameterized IERK methods, see Lemma 2.2. As described later, this framework will be fit for both Lobatto-type and Radau-type IERK methods, but is quite different from those in previous studies [31, 30, 11].

Section 3 addresses the details in constructing three parameterized second-order IERK (IERK2) methods and establishes the original energy dissipation laws for the resulting IERK2 methods. Also, we choose the optimal parameter in these IERK2 methods by using the concept of average dissipation rates, and present extensive tests to support our theoretical predictions, see Table 2, in which the parameter choices for the energy stability of six IERK2 methods are summarized. Third-order parameterized IERK (IERK3) methods are discussed and tested in Section 4 and, at the end of section, Table 3 collects some parameter choices for the energy stability of eight IERK3 methods. To show the existence of energy stable IERK methods with fourth-order time accuracy, two approximately fourth-order IERK (IERK4-A) methods are presented and tested in Section 5. In the last section, we summarize this article and present some further issues to be studied.

2 Stage energy laws of IERK methods

Requiring $U^{n,i}=u_{h}^{n}=u_{h}^{*}$ for all $i$ and $n\geq 0$ immediately shows that the canopy node condition (1.7) makes the IERK method (1.6) preserve naturally the equilibria $u_{h}^{*}$ of the gradient flow problem (1.2), that is, $L_{h}u_{h}^{*}=g_{h}(u_{h}^{*})$ . So one can reformulate the form (1.6) into the following steady-state preserving form

	$\displaystyle U^{n,i}=$	$\displaystyle\,U^{n,1}+\tau\sum_{j=1}^{i}a_{i,j}M_{h}\big{(}L_{h}U^{n,j}-L_{h}U^{n,1}\big{)}-\tau\sum_{j=1}^{i-1}\hat{a}_{i,j}M_{h}\big{[}g_{h}(U^{n,j})-L_{h}U^{n,1}\big{]}$
	$\displaystyle=$	$\displaystyle\,U^{n,1}+\tau\sum_{j=2}^{i}a_{i,j}M_{h}\big{(}L_{h}U^{n,j}-L_{h}U^{n,1}\big{)}-\tau\sum_{j=1}^{i-1}\hat{a}_{i,j}M_{h}\big{[}g_{h}(U^{n,j})-L_{h}U^{n,1}\big{]}$		(2.1)

for $2\leq i\leq s$ , in which we drop the terms with the coefficients $a_{i1}$ for $2\leq i\leq s$ . In this sense, we define the lower triangular coefficient matrices for the implicit and explicit parts, respectively,

A_{\mathrm{I}}:=\big{(}a_{i+1,j+1}\big{)}_{i,j=1}^{s_{\mathrm{I}}}\quad\text{and}\quad A_{\mathrm{E}}:=\big{(}\hat{a}_{i+1,j}\big{)}_{i,j=1}^{s_{\mathrm{I}}}\,.

The two matrices are always required in our theory on the energy dissipation property, while the coefficient vector $\mathbf{a}_{1}:=(a_{21},a_{31},\cdots,a_{s1})^{T}$ would be not involved directly. Due to the addition of $s_{\mathrm{I}}$ coefficients (degree of freedom) in the vector $\mathbf{a}_{1}$ , as shown later, it would be useful in designing some computationally effective Lobatto-type IERK methods.

2.1 Our theoretical framework

Motivated by the stabilization idea from Du et al. [8] and Fu et al. [10, 11], we introduce the following stabilized operators with a parameter $\kappa\geq 0$ ,

\displaystyle L_{\kappa}:=L_{h}+\kappa I\quad\text{and}\quad g_{\kappa}(u):=g_{h}(u)+\kappa u,

(2.2)

such that the problem (1.1) becomes the stabilized version

\displaystyle u_{h}^{\prime}(t)=M_{h}\left[L_{\kappa}u_{h}(t)-g_{\kappa}(u_{h})\right],\quad u_{h}(t_{0})=u_{h}^{0}.

(2.3)

Thus, applying (2) to (2.3), we have the following IERK method

\displaystyle U^{n,i+1}=U^{n,1}+\tau\sum_{j=1}^{i}a_{i+1,j+1}\sum_{\ell=1}^{j}M_{h}L_{\kappa}\delta_{\tau}U^{n,\ell+1}-\tau\sum_{j=1}^{i}\hat{a}_{i+1,j}M_{h}\big{[}g_{\kappa}(U^{n,j})-L_{\kappa}U^{n,1}\big{]}

(2.4)

for $1\leq i\leq s_{\mathrm{I}}$ , where the (stage) time difference $\delta_{\tau}U^{n,\ell+1}:=U^{n,\ell+1}-U^{n,\ell}$ for $\ell\geq 1$ .

To make our idea more concise, we assume further that the nonlinear function $g_{h}$ is Lipschitz continuous with a constant $\ell_{g}>0$ , cf. [32] or the recent discussions in [11]. In theoretical manner, the stabilization parameter $\kappa$ in (2.2) is chosen properly large (determining the minimum stabilized parameter is out of our current scope although it is also practically useful). In this sense, if an IERK method is proven to maintain the original energy dissipation law (1.3) unconditionally, we mean that this IERK method can be stabilized by setting a properly large parameter $\kappa$ . To derive the energy dissipation law of the general IERK method (2.4), we need the following result.

Lemma 2.1.

[24] If $g_{h}$ is Lipschitz-continuous with a constant $\ell_{g}>0$ and $\kappa\geq 2\ell_{g}$ , then

\displaystyle\big{\langle}u-v,\tfrac{1}{2}L_{\kappa}(u+v)-g_{\kappa}(v)\big{\rangle}\geq E[u]-E[v]+\tfrac{1}{2}(\kappa-2\ell_{g})\big{\|}u-v\big{\|}^{2},

where the energy $E$ is defined in (1.2).

Our theoretical framework contains three main steps:

(1)

Compute difference coefficients: we introduce a class of difference coefficients, for $i=2,\cdots,s$ ,

	$\displaystyle\underline{a}_{i,i}:=a_{i,i}\quad\text{and}\quad\underline{a}_{i,j}:=a_{i,j}-a_{i-1,j}\quad\text{for $1\leq j\leq i-1$;}$
	$\displaystyle\underline{\hat{a}}_{i,i-1}:=\hat{a}_{i,i-1}\quad\text{and}\quad\underline{\hat{a}}_{i,j-1}:=\hat{a}_{i,j-1}-\hat{a}_{i-1,j-1}\quad\text{for $2\leq j\leq i-1$.}$

It is not difficult to derive from (2.4) that

\displaystyle M_{h}^{-1}\delta_{\tau}U^{n,i+1}=

\displaystyle\,\tau\sum_{j=1}^{i}\underline{a}_{i+1,j+1}\sum_{\ell=1}^{j}L_{\kappa}\delta_{\tau}U^{n,\ell+1}+{\tau}\sum_{j=1}^{i}\underline{\hat{a}}_{i+1,j}\left[L_{\kappa}U^{n,1}-g_{\kappa}(U^{n,j})\right]

(2.5)

for $1\leq i\leq s_{\mathrm{I}}$ . The associated Butcher difference tableaux read

\begin{array}[]{c|ccccc}c_{1}=0&0&&&&\\ c_{2}&\underline{a}_{21}&\underline{a}_{22}&&&\\ c_{3}&\underline{a}_{31}&\underline{a}_{32}&\underline{a}_{33}&&\\ \vdots&\vdots&\vdots&\ddots&\ddots&\\[2.0pt] c_{s}=1&\underline{a}_{s,1}&\underline{a}_{s,2}&\cdots&\underline{a}_{s,s-1}&\underline{a}_{s,s}\\[2.0pt] \hline\cr&0&0&\cdots&0&0\end{array},\quad\begin{array}[]{c|ccccc}\hat{c}_{1}=0&0&&&&\\ \hat{c}_{2}&\underline{\hat{a}}_{21}&0&&&\\ \hat{c}_{3}&\underline{\hat{a}}_{31}&\underline{\hat{a}}_{32}&0&&\\ \vdots&\vdots&\vdots&\ddots&\ddots&\\[2.0pt] \hat{c}_{s}=1&\underline{\hat{a}}_{s,1}&\underline{\hat{a}}_{s,2}&\cdots&\underline{\hat{a}}_{s,s-1}&0\\[1.0pt] \hline\cr\\[-8.0pt] &0&0&\cdots&0&0\end{array}.

(2)

Determine DOC kernels and differential form: we introduce the so-called discrete orthogonal convolution (DOC) kernels $\underline{\theta}_{k,j}(z)$ with respect to the explicit coefficient $\underline{\hat{a}}_{i+1,j}$ , cf. [22, 23, 24, 25],

\displaystyle\underline{\theta}_{k,k}:=\frac{1}{\underline{\hat{a}}_{k+1,k}}\quad\text{and}\quad\underline{\theta}_{k,j}:=-\sum_{\ell=j+1}^{k}\underline{\theta}_{k,\ell}\frac{\underline{\hat{a}}_{\ell+1,j}}{\underline{\hat{a}}_{j+1,j}},\quad\text{$1\leq j\leq k-1$}

(2.6)

for $k=1,2,\cdots,s_{\mathrm{I}}$ . It is easy to check the following discrete orthogonal identity,

\displaystyle\sum_{\ell=j}^{m}\underline{\theta}_{m,\ell}\underline{\hat{a}}_{\ell+1,j}\equiv\delta_{m,j}\quad\text{for $1\leq j\leq m\leq s_{\mathrm{I}}$},

(2.7)

where $\delta_{m,j}$ is the Kronecker delta symbol with $\delta_{m,j}=0$ if $j\neq m$ . Multiplying the third term of (2.5) by the DOC kernels $\underline{\theta}_{k,i}$ , and summing $i$ from 1 to $k$ , one can exchange the summation order and apply the discrete orthogonal identity (2.7) to find that

	$\displaystyle\tau\sum_{i=1}^{k}\underline{\theta}_{k,i}\sum_{j=1}^{i}\underline{\hat{a}}_{i+1,j}$	$\displaystyle\,\left[L_{\kappa}U^{n,1}-g_{\kappa}(U^{n,j})\right]=\tau\sum\limits_{j=1}^{k}\left[L_{\kappa}U^{n,1}-g_{\kappa}(U^{n,j})\right]\sum_{i=j}^{k}\underline{\theta}_{k,i}\underline{\hat{a}}_{i+1,j}$
	$\displaystyle=$	$\displaystyle\,\tau L_{\kappa}U^{n,1}-\tau g_{\kappa}(U^{n,k})=\tau L_{\kappa}U^{n,k+1}-\tau g_{\kappa}(U^{n,k})-\tau\sum_{\ell=1}^{k}L_{\kappa}\delta_{\tau}U^{n,\ell+1}$
	$\displaystyle=$	$\displaystyle\,\tfrac{\tau}{2}L_{\kappa}(U^{n,k+1}+U^{n,k})-\tau g_{\kappa}(U^{n,k})-\tau\sum_{\ell=1}^{k}\big{(}1-\tfrac{1}{2}\delta_{k,\ell}\big{)}L_{\kappa}\delta_{\tau}U^{n,\ell+1}$

for $1\leq k\leq s_{\mathrm{I}}$ . Similarly, by multiplying the second term of (2.5) by the DOC kernels $\underline{\theta}_{k,i}$ , and summing $i$ from 1 to $k$ , one has

\displaystyle\tau\sum_{i=1}^{k}\underline{\theta}_{k,i}\sum_{j=1}^{i}\underline{a}_{i+1,j+1}\sum_{\ell=1}^{j}L_{\kappa}\delta_{\tau}U^{n,\ell+1}=\tau\sum\limits_{j=1}^{k}\sum_{\ell=1}^{j}L_{\kappa}\delta_{\tau}U^{n,\ell+1}\sum_{i=j}^{k}\underline{\theta}_{k,i}\underline{a}_{i+1,j+1}

for $1\leq k\leq s_{\mathrm{I}}$ . With the help of the above two equalities, we have an equivalent form of the IERK method (2.4),

	$\displaystyle\sum\limits_{\ell=1}^{k}\Big{(}\underline{\theta}_{k,\ell}-{\tau}M_{h}L_{\kappa}\sum_{j=\ell}^{k}\sum_{i=j}^{k}\underline{\theta}_{k,i}\underline{a}_{i+1,j+1}+{\tau}M_{h}L_{\kappa}-\tfrac{1}{2}{\tau}M_{h}L_{\kappa}\delta_{k,\ell}\Big{)}\delta_{\tau}U^{n,\ell+1}$
	$\displaystyle=\tau M_{h}\left[\tfrac{1}{2}L_{\kappa}(U^{n,k+1}+U^{n,k})-g_{\kappa}(U^{n,k})\right]\quad\text{for $1\leq k\leq s_{\mathrm{I}}$.}$

Thus we obtain the following differential form of the IERK method (2.4)

\displaystyle\sum_{\ell=1}^{k}d_{k,\ell}(\tau M_{h}L_{\kappa})\delta_{\tau}U^{n,\ell+1}=\tau M_{h}\left[\tfrac{1}{2}L_{\kappa}(U^{n,k+1}+U^{n,k})-g_{\kappa}(U^{n,k})\right]

(2.8)

for $1\leq k\leq s_{\mathrm{I}}$ , where the elements $d_{k,\ell}$ are defined by $d_{k,\ell}(z):=0$ for $\ell>k$ , and

\displaystyle d_{k,\ell}(z):=\underline{\theta}_{k,\ell}-z\sum_{j=\ell}^{k}\sum_{i=j}^{k}\underline{\theta}_{k,i}\underline{a}_{i+1,j+1}+z-\frac{z}{2}\delta_{k,\ell}\quad\text{for $1\leq\ell\leq k\leq s_{\mathrm{I}}$}.

(2.9)

The associated lower triangular matrix $D:=(d_{k,\ell})_{s_{\mathrm{I}}\times s_{\mathrm{I}}}$ is called the differentiation matrix.

Let $E_{s_{\mathrm{I}}}:=(1_{i\geq j})_{s_{\mathrm{I}}\times s_{\mathrm{I}}}$ be the lower triangular matrix full of element 1, and let $I_{s_{\mathrm{I}}}$ be the identity matrix of the same size as $A_{\mathrm{I}}$ . One has

(\underline{\hat{a}}_{i+1,j})_{s_{\mathrm{I}}\times s_{\mathrm{I}}}=E_{s_{\mathrm{I}}}^{-1}A_{\mathrm{E}},\quad(\underline{\theta}_{k,\ell})_{s_{\mathrm{I}}\times s_{\mathrm{I}}}=A_{\mathrm{E}}^{-1}E_{s_{\mathrm{I}}}\quad\text{and}\quad(\underline{a}_{i+1,j+1})_{s_{\mathrm{I}}\times s_{\mathrm{I}}}=E_{s_{\mathrm{I}}}^{-1}A_{\mathrm{I}}.

Thus the above differentiation matrix $D$ defined in (2.9) can be formulated as

\displaystyle D(z)=D_{\mathrm{E}}-D_{\mathrm{EI}}z\quad\text{with}\quad D_{\mathrm{E}}:=A_{\mathrm{E}}^{-1}E_{s_{\mathrm{I}}},\quad D_{\mathrm{EI}}:=A_{\mathrm{E}}^{-1}A_{\mathrm{I}}E_{s_{\mathrm{I}}}-E_{s_{\mathrm{I}}}+\tfrac{1}{2}I_{s_{\mathrm{I}}}.

(2.10)

Always, if the symmetric part $\mathcal{S}(D;z):=\tfrac{1}{2}[D(z)+D(z)^{T}]$ is positive (semi-)definite, we say that the matrix $D(z)$ is positive (semi-)definite.

(3)

Establish stage energy dissipation law: this process is standard and we have the following result, which simulates the original energy dissipation law (1.3) at all stages.

Theorem 2.1.

If the differentiation matrix $D(z)$ in (2.10) is positive (semi-)definite for $z\leq 0$ , the IERK methods (2.4) with the stabilized parameter $\kappa\geq 2\ell_{g}$ preserve the original energy dissipation law (1.3) at all stages,

\displaystyle E[U^{n,j+1}]-E[U^{n,1}]\leq

\displaystyle\,\frac{1}{\tau}\sum_{k=1}^{j}\Big{\langle}M_{h}^{-1}\delta_{\tau}U^{n,k+1},\sum_{\ell=1}^{k}d_{k,\ell}({\tau}M_{h}L_{\kappa})\delta_{\tau}U^{n,\ell+1}\Big{\rangle}

(2.11)

for $1\leq j\leq s_{\mathrm{I}}$ , and in particular, by taking $j:=s_{\mathrm{I}}$ ,

\displaystyle E[u_{h}^{n}]-E[u_{h}^{n-1}]\leq\frac{1}{\tau}\sum_{k=1}^{s_{\mathrm{I}}}\Big{\langle}M_{h}^{-1}\delta_{\tau}U^{n,k+1},\sum_{\ell=1}^{k}d_{k,\ell}({\tau}M_{h}L_{\kappa})\delta_{\tau}U^{n,\ell+1}\Big{\rangle}\quad\text{for $n\geq 1$.}

Proof.

Making the inner product of the equivalent form (2.8) with $\frac{1}{\tau}M_{h}^{-1}\delta_{\tau}U^{n,k+1}$ and summing $k$ from $k=1$ to $j$ , one can find that

	$\displaystyle\frac{1}{\tau}\sum_{k=1}^{j}$	$\displaystyle\,\Big{\langle}M_{h}^{-1}\delta_{\tau}U^{n,k+1},\sum_{\ell=1}^{k}d_{k,\ell}({\tau}M_{h}L_{\kappa})\delta_{\tau}U^{n,\ell+1}\Big{\rangle}$
		$\displaystyle\,=\sum_{k=1}^{j}\Big{\langle}\delta_{\tau}U^{n,k+1},\tfrac{1}{2}L_{\kappa}(U^{n,k+1}+U^{n,k})-g_{\kappa}(U^{n,k})\Big{\rangle}$

for $1\leq j\leq s_{\mathrm{I}}$ . Lemma 2.1 yields the following energy dissipation law at each stage

\displaystyle E[U^{n,j+1}]-E[U^{n,1}]-

\displaystyle\,\frac{1}{\tau}\sum_{k=1}^{j}\Big{\langle}M_{h}^{-1}\delta_{\tau}U^{n,k+1},\sum_{\ell=1}^{k}d_{k,\ell}({\tau}M_{h}L_{\kappa})\delta_{\tau}U^{n,\ell+1}\Big{\rangle}\leq 0

for $1\leq j\leq s_{\mathrm{I}}$ . It completes the proof. ∎

For $1\leq j\leq s_{\mathrm{I}}$ , let $D_{j}:=D[1:j,1:j]$ be the $j$ -th sequential sub-matrix of the differentiation matrix $D(z)$ . The above stage energy dissipation law (2.11) can be formulated as

\displaystyle E[U^{n,j+1}]-E[U^{n,1}]\leq\frac{1}{\tau}\Big{\langle}M_{h}^{-1}\delta_{\tau}\vec{U}_{n,j+1},D_{j}({\tau}M_{h}L_{\kappa})\delta_{\tau}\vec{U}_{n,j+1}\Big{\rangle}\quad\text{for $1\leq j\leq s_{\mathrm{I}}$.}

(2.12)

The involved vector $\delta_{\tau}\vec{U}_{n,j+1}:=(\delta_{\tau}U^{n,2},\delta_{\tau}U^{n,3},\cdots,\delta_{\tau}U^{n,j+1})^{T}$ and the matrix $D(\tau M_{h}L_{\kappa})$ can be formulated as

\displaystyle D(\tau M_{h}L_{\kappa})=D_{\mathrm{E}}\otimes I+D_{\mathrm{EI}}\otimes(-\tau M_{h}L_{\kappa}),

(2.13)

where $I$ is the identity matrix of the same size as $L_{\kappa}$ and $\otimes$ represents the Kronecker product. Since the matrix $-\tau M_{h}L_{\kappa}$ is symmetric and positive semi-definite, the properties of Kronecker product and the formula (2.13) arrive at the following corollary.

Corollary 2.1.

If the two $s_{\mathrm{I}}\times s_{\mathrm{I}}$ matrices $D_{\mathrm{E}}$ and $D_{\mathrm{EI}}$ are positive (semi-)definite, the results of Theorem 2.1 are valid.

2.2 Average dissipation rate

Theorem 2.1 shows that the IERK method (2.4) is unconditionally energy stable if the differentiation matrix $D(\tau M_{h}L_{\kappa})$ is positive semi-definite, that is, all eigenvalues of the symmetric part $\mathcal{S}(D;\tau M_{h}L_{\kappa})$ are nonnegative. A necessary condition is that the average eigenvalue is nonnegative,

\overline{\lambda}\big{[}\mathcal{S}(D;\tau M_{h}L_{\kappa})\big{]}=\frac{1}{s_{\mathrm{I}}m_{L_{\kappa}}}\mathrm{tr}\big{[}D(\tau M_{h}L_{\kappa})\big{]}\geq 0,

where $m_{L_{\kappa}}$ is the size of $L_{\kappa}$ and $\mathrm{tr}\left(D\right)$ is the trace of $D$ . By comparing the discrete energy dissipation law (2.12) with the continuous counterpart (1.3), the overall energy dissipation rate of the energy $E[u_{h}^{n}]$ could be roughly estimated by the average eigenvalue. Following the idea in [24], we will use the average dissipation rate, defined by

\displaystyle\mathcal{R}:=\overline{\lambda}\big{[}\mathcal{S}(D;\tau M_{h}L_{\kappa})\big{]}=\frac{1}{s_{\mathrm{I}}m_{L_{\kappa}}}\mathrm{tr}\big{[}D(\tau M_{h}L_{\kappa})\big{]},

(2.14)

to examine the energy dissipation behaviors of IERK methods. By using the definitions (2.9) and (2.6), one can compute the diagonal elements of the matrices $D_{\mathrm{E}}$ and $D_{\mathrm{EI}}$ ,

	$\displaystyle(D_{\mathrm{E}})_{k,k}=$	$\displaystyle\,\underline{\theta}_{k,k}=\tfrac{1}{\underline{\hat{a}}_{k+1,k}}=\tfrac{1}{\hat{a}_{k+1,k}},$
	$\displaystyle(D_{\mathrm{EI}})_{k,k}=$	$\displaystyle\,\underline{\theta}_{k,k}\underline{a}_{k+1,k+1}-\tfrac{1}{2}=\tfrac{\underline{a}_{k+1,k+1}}{\underline{\hat{a}}_{k+1,k}}-\tfrac{1}{2}=\tfrac{a_{k+1,k+1}}{\hat{a}_{k+1,k}}-\tfrac{1}{2}\quad\text{for $1\leq k\leq s_{\mathrm{I}}$.}$

Thus, by using the Kronecker product form (2.13) and the following property

\mathrm{tr}\big{[}D_{\mathrm{E}}\otimes I+D_{\mathrm{EI}}\otimes(-\tau M_{h}L_{\kappa})\big{]}=\mathrm{tr}(D_{\mathrm{E}})\mathrm{tr}(I)+\mathrm{tr}(D_{\mathrm{EI}})\mathrm{tr}(-\tau M_{h}L_{\kappa}),

it is easy to obtain the following result.

Lemma 2.2.

If the two matrices $D_{\mathrm{E}}$ and $D_{\mathrm{EI}}$ are positive (semi-)definite, then the average dissipation rate of the IERK method (2.4) is nonnegative, that is,

\displaystyle\mathcal{R}=\frac{1}{s_{\mathrm{I}}}\mathrm{tr}(D_{\mathrm{E}})+\frac{1}{s_{\mathrm{I}}}\mathrm{tr}(D_{\mathrm{EI}})\tau\overline{\lambda}_{\mathrm{ML}}=\frac{1}{s_{\mathrm{I}}}\sum_{k=1}^{s_{\mathrm{I}}}\frac{1}{\hat{a}_{k+1,k}}+\frac{1}{s_{\mathrm{I}}}\sum_{k=1}^{s_{\mathrm{I}}}\Big{(}\frac{a_{k+1,k+1}}{\hat{a}_{k+1,k}}-\frac{1}{2}\Big{)}\tau\overline{\lambda}_{\mathrm{ML}}\geq 0,

(2.15)

where the average eigenvalue $\overline{\lambda}_{\mathrm{ML}}=\overline{\lambda}_{\mathrm{ML}}(\kappa,h)>0$ of the symmetric, positive definite matrix $-M_{h}L_{\kappa}$ is determined by the stabilized parameter $\kappa$ , the method of spatial approximation and the grid spacing $h$ as well.

Since the average dissipation rate (2.14) is defined via the similarity between the energy dissipation law (1.3) and the discrete counterpart (2.12), we say that an IERK method is a “good” candidate to preserve the original energy dissipation law (1.3) if the average dissipation rate $\mathcal{R}\geq 0$ and is as close to 1 as possible within properly large range of $\tau\overline{\lambda}_{\mathrm{ML}}$ . The specific expression of $\overline{\lambda}_{\mathrm{ML}}$ will depend on the phase field model involved; but always, the larger the value of $\kappa$ , the greater the value of $\overline{\lambda}_{\mathrm{ML}}$ , while the smaller the spatial step $h$ , the greater the value of $\overline{\lambda}_{\mathrm{ML}}$ . At the same time, the values of $\frac{1}{s_{\mathrm{I}}}\mathrm{tr}(D_{\mathrm{E}})$ and $\frac{1}{s_{\mathrm{I}}}\mathrm{tr}(D_{\mathrm{EI}})$ are entirely determined by the IERK method itself, regardless of the time-step size $\tau$ and the spatial step $h$ adopted by potential users. Throughout this paper, we aim to design some “good” IERK methods by minimizing the values of $\frac{1}{s_{\mathrm{I}}}\mathrm{tr}(D_{\mathrm{EI}})$ and $\frac{1}{s_{\mathrm{I}}}\mathrm{tr}(D_{\mathrm{E}})$ as far as possible to make the associated average dissipation rate $\mathcal{R}$ as close to 1 as possible within a given number of implicit stages.

Obviously, the average dissipation rate $\mathcal{R}$ is only a rough indicator for describing discrete energy behavior of IERK method, that is, it would be a qualitative rather than quantitative indicator, since it ignores the computational accuracy of the algorithm. As remarked in [24], one can not expect that the long-time dynamics of the energy can be understood by such a scalar alone. Nonetheless, Lemma 2.2 provides us a very simple criterion to evaluate the overall energy dissipation rate of an IERK method so that one can choose proper parameters in some parameterized IERK methods or compare different IERK methods with the same time accuracy.

As an example, consider the first-order IERK (IERK1) method [14, p. 383],

\begin{array}[]{c|c}\mathbf{c}&A\\ \hline\cr\\[-10.0pt] &\mathbf{b}^{T}\end{array}=\begin{array}[]{c|cc}0&0&0\\ 1&1-\theta&\theta\\ \hline\cr&1-\theta&\theta\end{array}\;,\quad\begin{array}[]{c|c}\hat{\mathbf{c}}&\widehat{A}\\ \hline\cr\\[-10.0pt] &\hat{\mathbf{b}}^{T}\end{array}=\begin{array}[]{c|cc}0&0&0\\ 1&1&0\\ \hline\cr&1&0\end{array}\;.

(2.16)

Obviously, $A_{\mathrm{I}}=(\theta)$ and $A_{\mathrm{E}}=(1)$ such that the differentiation matrix $D^{(1)}(z)=1-z\left(\theta-\tfrac{1}{2}\right)$ is positive semi-definite for $z\leq 0$ provided the weighted parameter $\theta\geq\frac{1}{2}$ according to Corollary 2.1. Lemma 2.2 says that the average dissipation rate

\displaystyle\mathcal{R}^{(1)}(\theta)=1+\left(\theta-\tfrac{1}{2}\right)\tau\overline{\lambda}_{\mathrm{ML}}\quad\text{for $\theta\geq\tfrac{1}{2}$}.

(2.17)

To enlarge the choice of $\tau\overline{\lambda}_{\mathrm{ML}}$ , one can choose $\theta=1/2$ such that the average dissipation rate, $\mathcal{R}^{(1)}(\tfrac{1}{2})=1$ , is independent of $\tau\overline{\lambda}_{\mathrm{ML}}$ . We say that the stabilized Crank-Nicolson type scheme

\displaystyle\delta_{\tau}u_{h}^{n}=\tau M_{h}\left[\tfrac{1}{2}L_{\kappa}(u_{h}^{n}+u^{n-1})-g_{\kappa}(u^{n-1})\right]\quad\text{for $n\geq 1$,}

is a “good” candidate to preserve the original energy dissipation law (1.3) unconditionally.

Here and hereafter, the superscript $(p)$ is always used to indicate the formal order of the method, that is to say, $D^{(p)}$ and $\mathcal{R}^{(p)}$ denote the associated differential matrix and the average dissipation rate, respectively, of a $p$ -th order IERK method.

3 Discrete energy laws of IERK2 methods

3.1 Lobatto-type IERK2 methods

Second-order methods require two implicit stages, $s_{I}=2$ . Consider the 3-stage IERK methods that satisfy the canopy node condition and the two order conditions for the first-order accuracy,

\begin{array}[]{c|c}\mathbf{c}&A\\ \hline\cr\\[-10.0pt] &\mathbf{b}^{T}\end{array}=\begin{array}[]{c|ccc}0&0&&\\ c_{2}&c_{2}-a_{22}&a_{22}&\\ 1&1-a_{32}-a_{33}&a_{32}&a_{33}\\ \hline\cr\\[-10.0pt] &1-a_{32}-a_{33}&a_{32}&a_{33}\end{array}\;,\quad\begin{array}[]{c|c}\hat{\mathbf{c}}&\widehat{A}\\ \hline\cr\\[-10.0pt] &\hat{\mathbf{b}}^{T}\end{array}=\begin{array}[]{c|ccc}0&0&&\\ c_{2}&c_{2}&0&\\ 1&1-\hat{a}_{32}&\hat{a}_{32}&0\\ \hline\cr\\[-10.0pt] &1-\hat{a}_{32}&\hat{a}_{32}&0\end{array}\;.

We should determine five independent coefficients. By the two order conditions for second-order accuracy, see the second line in Table 1, one has three independent coefficients to be determined. In detail, from the stand-alone condition for explicit part, $\hat{\mathbf{b}}^{T}\mathbf{c}=\tfrac{1}{2}$ , one has $\hat{a}_{32}=\tfrac{1}{2c_{2}}$ . From the stand-alone condition for implicit part, $\mathbf{b}^{T}\mathbf{c}=\tfrac{1}{2}$ , one has $a_{32}c_{2}+a_{33}=\tfrac{1}{2}$ such that $a_{32}=\frac{1-2a_{33}}{2c_{2}}$ . According to Lemma 2.2 and the formula (2.15), we always choose

a_{22}=2c_{2}^{2}a_{33}>0

such that the trace $\mathrm{tr}(D_{\mathrm{EI}}^{(2,1)})$ can attain the minimum value, that is,

\displaystyle\mathrm{tr}\big{[}D_{\mathrm{EI}}^{(2,1)}(c_{2},a_{33})\big{]}=\frac{a_{22}}{c_{2}}+2c_{2}a_{33}-1=4c_{2}a_{33}-1\quad\text{for $a_{33},c_{2}>0$.}

(3.1)

We obtain the following two-parameter class of second-order IERK (IERK2-1) methods with the associated Butcher tableaux, also see [26],

\begin{array}[]{c|ccc}0&0&&\\ c_{2}&c_{2}-2c_{2}^{2}a_{33}&2c_{2}^{2}a_{33}&\\[2.0pt] 1&1-\frac{1}{2c_{2}}+\frac{a_{33}(1-c_{2})}{c_{2}}&\frac{1-2a_{33}}{2c_{2}}&a_{33}\\ \hline\cr\\[-10.0pt] &1-\frac{1}{2c_{2}}+\frac{a_{33}(1-c_{2})}{c_{2}}&\frac{1-2a_{33}}{2c_{2}}&a_{33}\end{array}\;,\quad\begin{array}[]{c|ccc}0&0&&\\ c_{2}&c_{2}&0&\\[2.0pt] 1&1-\tfrac{1}{2c_{2}}&\tfrac{1}{2c_{2}}&0\\ \hline\cr\\[-10.0pt] &1-\tfrac{1}{2c_{2}}&\tfrac{1}{2c_{2}}&0\end{array}\;.

(3.2)

Now we are to determine the two independent coefficients $c_{2}$ and $a_{33}$ such that the IERK2-1 methods (3.2) are “good” candidates to preserve the original energy dissipation law (1.3) unconditionally. By the definition (2.10), one has

\displaystyle D_{\mathrm{E}}^{(2,1)}(c_{2}):=A_{\mathrm{E}}^{-1}E_{2}=\begin{pmatrix}\frac{1}{c_{2}}&0\\[4.0pt] 2c_{2}+\frac{1}{c_{2}}-2&2c_{2}\end{pmatrix}.

The determinant of the symmetric part

\mathrm{Det}\big{[}\mathcal{S}(D_{\mathrm{E}}^{(2,1)};c_{2})\big{]}=\big{(}c_{2}+\sqrt{2}-1+\tfrac{1}{2c_{2}}\big{)}\big{(}-c_{2}+\sqrt{2}+1-\tfrac{1}{2c_{2}}\big{)}.

It is easy to check that

\displaystyle 0.228788\approx\tfrac{1+\sqrt{2}-\sqrt{1+2\sqrt{2}}}{2}\leq c_{2}\leq\tfrac{1+\sqrt{2}+\sqrt{1+2\sqrt{2}}}{2}\approx 2.18543

(3.3)

is sufficient to ensure the positive semi-definiteness of $D_{\mathrm{E}}^{(2,1)}(c_{2})$ . Now compute the lower triangular matrix $D_{\mathrm{EI}}^{(2,1)}$ defined by (2.10),

\displaystyle D_{\mathrm{EI}}^{(2,1)}(c_{2},a_{33}):=A_{\mathrm{E}}^{-1}A_{\mathrm{I}}E_{2}-E_{2}+\tfrac{1}{2}I_{2}=\begin{pmatrix}2c_{2}a_{33}-\frac{1}{2}&0\\[4.0pt] -2a_{33}(2c_{2}^{2}-2c_{2}+1)&2c_{2}a_{33}-\frac{1}{2}\end{pmatrix}.

The positive semi-definiteness of $D_{\mathrm{EI}}^{(2,1)}(c_{2},a_{33})$ requires $a_{33}\geq\frac{1}{4c_{2}}$ and

	$\displaystyle\mathrm{Det}\big{[}\mathcal{S}(D_{\mathrm{EI}}^{(2,1)};c_{2},a_{33})\big{]}=$	$\displaystyle-a_{33}^{2}\left(2c_{2}^{2}+1\right)(2c_{2}^{2}-4c_{2}+1)-2a_{33}c_{2}+\frac{1}{4}$
	$\displaystyle=$	$\displaystyle(2c_{2}^{2}+1)(4c_{2}-2c_{2}^{2}-1)\Big{(}a_{33}-\tfrac{1}{2(2c_{2}^{2}+1)}\Big{)}\Big{(}a_{33}-\tfrac{1}{2(4c_{2}-2c_{2}^{2}-1)}\Big{)}\geq 0,$

in which we set $4c_{2}-2c_{2}^{2}-1>0$ , or $1-\frac{\sqrt{2}}{2}<c_{2}<1+\frac{\sqrt{2}}{2}$ . Otherwise, the fact $a_{33}\geq\frac{1}{4c_{2}}>\frac{1}{4c_{2}^{2}+2}$ arrives at the negative determinant, $\mathrm{Det}\big{[}\mathcal{S}(D_{\mathrm{EI}}^{(2,1)};c_{2},a_{33})\big{]}<0$ . Thus, by using the restriction (3.3) and the fact $\frac{1}{4c_{2}^{2}+2}\leq\frac{1}{4c_{2}}\leq\frac{1}{2(4c_{2}-2c_{2}^{2}-1)}$ , one has the following conditions for the coefficient $a_{33}$ ,

\displaystyle a_{33}\geq\frac{1}{2(4c_{2}-2c_{2}^{2}-1)}=\frac{1}{4(c_{2}-1+\frac{\sqrt{2}}{2})(1+\frac{\sqrt{2}}{2}-c_{2})}\quad\text{for $1-\tfrac{\sqrt{2}}{2}<c_{2}<1+\tfrac{\sqrt{2}}{2}$.}

(3.4)

By using Corollary 2.1 and Lemma 2.2, we have the following theorem.

Theorem 3.1.

In simulating the gradient flow system (2.3) with $\kappa\geq 2\ell_{g}$ , the two-parameter IERK2-1 methods (3.2) with the parameter setting (3.4) preserve the original energy dissipation law (1.3) unconditionally at all stages in the sense that

\displaystyle E[U^{n,j+1}]-E[U^{n,1}]\leq\frac{1}{\tau}\Big{\langle}M_{h}^{-1}\delta_{\tau}\vec{U}_{n,j+1},D_{j}^{(2,1)}(c_{2},a_{33};{\tau}M_{h}L_{\kappa})\delta_{\tau}\vec{U}_{n,j+1}\Big{\rangle}\quad\text{for $1\leq j\leq 2$,}

where the differentiation matrix $D^{(2,1)}$ is defined by

D^{(2,1)}(c_{2},a_{33};{\tau}M_{h}L_{\kappa}):=D_{\mathrm{E}}^{(2,1)}(c_{2})\otimes I+D_{\mathrm{EI}}^{(2,1)}(c_{2},a_{33})\otimes(-\tau M_{h}L_{\kappa}).

The associated average dissipation rate

\displaystyle\mathcal{R}^{(2,1)}(c_{2},a_{33})=\tfrac{1}{2}\mathrm{tr}\big{[}D_{\mathrm{E}}^{(2,1)}(c_{2})\big{]}+\tfrac{1}{2}\mathrm{tr}\big{[}D_{\mathrm{EI}}^{(2,1)}(c_{2},a_{33})\big{]}\tau\overline{\lambda}_{\mathrm{ML}}=c_{2}+\tfrac{1}{2c_{2}}+(2c_{2}a_{33}-\tfrac{1}{2})\tau\overline{\lambda}_{\mathrm{ML}}.

(3.5)

The parameter setting (3.4) implies that one can minimize the trace $\mathrm{tr}\big{[}D_{\mathrm{EI}}^{(2,1)}(c_{2},a_{33})\big{]}$ by choosing the lower bound of $a_{33}$ , that is,

\displaystyle\mathrm{tr}\big{[}D_{\mathrm{EI}}^{(2,1)}(c_{2},a_{33})\big{]}\geq

\displaystyle\,\frac{c_{2}}{4c_{2}-2c_{2}^{2}-1}-\frac{1}{2}\quad\text{for $1-\tfrac{\sqrt{2}}{2}<c_{2}<1+\tfrac{\sqrt{2}}{2}$.}

The minimum value is attained by setting $c_{2}=\frac{\sqrt{2}}{2}$ due to the fact $\frac{\,\mathrm{d}}{\,\mathrm{d}c_{2}}\big{(}\frac{c_{2}}{4c_{2}-2c_{2}^{2}-1}\big{)}=\frac{2c_{2}^{2}-1}{(2c_{2}^{2}-4c_{2}+1)^{2}}$ . In this case, the two-parameter IERK2-1 methods (3.2) reduce into the one-parameter IERK2 (IERK2-2) methods with the parameter $a_{33}\geq\tfrac{1+\sqrt{2}}{4}$ ,

\begin{array}[]{c|ccc}0&0&&\\ \frac{\sqrt{2}}{2}&\frac{\sqrt{2}}{2}-a_{33}&a_{33}&\\[1.0pt] 1&\frac{\sqrt{2}-1+(2-\sqrt{2})a_{33}}{\sqrt{2}}&\frac{1-2a_{33}}{\sqrt{2}}&a_{33}\\[1.0pt] \hline\cr\\[-10.0pt] &\frac{\sqrt{2}-1+(2-\sqrt{2})a_{33}}{\sqrt{2}}&\frac{1-2a_{33}}{\sqrt{2}}&a_{33}\end{array}\;,\quad\begin{array}[]{c|ccc}0&0&&\\ \frac{\sqrt{2}}{2}&\frac{\sqrt{2}}{2}&0&\\[1.0pt] 1&\frac{2-\sqrt{2}}{2}&\frac{\sqrt{2}}{2}&0\\[1.0pt] \hline\cr\\[-10.0pt] &\frac{2-\sqrt{2}}{2}&\frac{\sqrt{2}}{2}&0\end{array}\;.

(3.6)

As expected, from the perspective of average dissipation rate, the optimal IERK2-2 methods maintain a fixed diagonal ratio between the implicit and explicit parts, that is, $\frac{a_{k+1,k+1}}{\hat{a}_{k+1,k}}=\sqrt{2}a_{33}$ for $1\leq k\leq 2$ . The associated average dissipation rate

\displaystyle\mathcal{R}^{(2,2)}(\tfrac{\sqrt{2}}{2},a_{33})=\sqrt{2}+\sqrt{2}(a_{33}-\tfrac{\sqrt{2}}{4})\tau\overline{\lambda}_{\mathrm{ML}}\quad\text{for $a_{33}\geq\tfrac{1+\sqrt{2}}{4}$}.

(3.7)

Note that, the 3-stage Lobatto-type IERK2 method in [6, Section 2] would be not suitable for gradient flow problems because the associated differentiation matrix is not positive definite.

3.2 Radau-type IERK2 methods

For completeness, we also briefly consider the Radau-type DIRK methods, also known as ARS-type IERK methods [2, 11, 31], for the implicit part by assuming $\mathbf{a}_{1}=\mathbf{0}$ . In this case, one has three independent coefficients, see the associated Butcher tableaux as follows,

\begin{array}[]{c|c}\mathbf{c}&A\\ \hline\cr\\[-10.0pt] &\mathbf{b}^{T}\end{array}=\begin{array}[]{c|ccc}0&0&&\\ c_{2}&0&c_{2}&\\ 1&0&1-a_{33}&a_{33}\\ \hline\cr\\[-10.0pt] &0&1-a_{33}&a_{33}\end{array}\;,\quad\begin{array}[]{c|c}\hat{\mathbf{c}}&\widehat{A}\\ \hline\cr\\[-10.0pt] &\hat{\mathbf{b}}^{T}\end{array}=\begin{array}[]{c|ccc}0&0&&\\ c_{2}&c_{2}&0&\\ 1&1-\hat{a}_{32}&\hat{a}_{32}&0\\ \hline\cr\\[-10.0pt] &1-\hat{a}_{32}&\hat{a}_{32}&0\end{array}\;.

From $\hat{\mathbf{b}}^{T}\mathbf{c}=\tfrac{1}{2}$ , one has $\hat{a}_{32}=\tfrac{1}{2c_{2}}$ . The condition $\mathbf{b}^{T}\mathbf{c}=\tfrac{1}{2}$ leads to $a_{33}=\frac{1-2c_{2}}{2(1-c_{2})}$ . We obtain the following $c_{2}$ -parameterized class of Radau-type IERK (IERK2-Radau) methods [2]

\begin{array}[]{c|ccc}0&0&&\\ c_{2}&0&c_{2}&\\ 1&0&\frac{1}{2(1-c_{2})}&\frac{1-2c_{2}}{2(1-c_{2})}\\ \hline\cr\\[-10.0pt] &0&\frac{1}{2(1-c_{2})}&\frac{1-2c_{2}}{2(1-c_{2})}\end{array},\quad\begin{array}[]{c|ccc}0&0&&\\ c_{2}&c_{2}&0&\\ 1&1-\tfrac{1}{2c_{2}}&\tfrac{1}{2c_{2}}&0\\ \hline\cr\\[-10.0pt] &1-\tfrac{1}{2c_{2}}&\tfrac{1}{2c_{2}}&0\end{array}.

(3.8)

Obviously, the restriction (3.3) is sufficient to ensure the positive semi-definiteness of $D_{\mathrm{rad},\mathrm{E}}^{(2)}$ . We also apply the definition (2.10) to find

\displaystyle D_{\mathrm{rad},\mathrm{EI}}^{(2)}(c_{2}):=A_{\mathrm{E}}^{-1}A_{\mathrm{I}}E_{2}-E_{2}+\tfrac{1}{2}I_{2}=\begin{pmatrix}\frac{1}{2}&0\\[4.0pt] 0&\frac{4c_{2}^{2}-3c_{2}+1}{2(c_{2}-1)}\end{pmatrix}.

The positive semi-definiteness of $D_{\mathrm{rad},\mathrm{EI}}^{(2)}(c_{2})$ requires $c_{2}>1$ . By using Corollary 2.1 and Lemma 2.2, we have the following theorem.

Theorem 3.2.

In simulating the gradient flow system (2.3) with $\kappa\geq 2\ell_{g}$ , the $c_{2}$ -parameterized IERK2-Radau methods (3.8) with the parameter $1<c_{2}\leq\tfrac{1+\sqrt{2}+\sqrt{1+2\sqrt{2}}}{2}$ preserve the original energy dissipation law (1.3) unconditionally at all stages in the sense that

\displaystyle E[U^{n,j+1}]-E[U^{n,1}]\leq\frac{1}{\tau}\Big{\langle}M_{h}^{-1}\delta_{\tau}\vec{U}_{n,j+1},D_{\mathrm{rad},j}^{(2)}(c_{2};{\tau}M_{h}L_{\kappa})\delta_{\tau}\vec{U}_{n,j+1}\Big{\rangle}\quad\text{for $1\leq j\leq 2$,}

where the associated differentiation matrix $D_{\mathrm{rad}}^{(2)}$ is defined by

D_{\mathrm{rad}}^{(2)}(c_{2};{\tau}M_{h}L_{\kappa}):=D_{\mathrm{rad},\mathrm{E}}^{(2)}(c_{2})\otimes I+D_{\mathrm{rad},\mathrm{EI}}^{(2)}(c_{2})\otimes(-\tau M_{h}L_{\kappa}).

The associated average dissipation rate

\displaystyle\mathcal{R}_{\mathrm{rad}}^{(2)}(c_{2})=\tfrac{1}{2c_{2}}+c_{2}+\big{(}2c_{2}+1+\tfrac{1}{c_{2}-1}\big{)}\tau\overline{\lambda}_{\mathrm{ML}}\quad\text{for $1<c_{2}\leq\tfrac{1+\sqrt{2}+\sqrt{1+2\sqrt{2}}}{2}$}.

(3.9)

To enlarge the possible choices of $\tau\overline{\lambda}_{\mathrm{ML}}$ , we choose $c_{2}=1+\frac{\sqrt{2}}{2}$ such that $\mathrm{tr}\big{[}D_{\mathrm{rad},\mathrm{EI}}^{(2)}(c_{2})\big{]}$ attains the minimum value and $\mathcal{R}_{\mathrm{rad}}^{(2)}(1+\tfrac{\sqrt{2}}{2})=2+(3+2\sqrt{2})\tau\overline{\lambda}_{\mathrm{ML}}$ . However, it is always larger than the average dissipation rate $\mathcal{R}^{(2,2)}(\tfrac{\sqrt{2}}{2},a_{33})$ of the IERK2-2 methods (3.6) for $\tfrac{1+\sqrt{2}}{4}\leq a_{33}\leq\tfrac{8+7\sqrt{2}}{4}$ . Under the same stabilization strategy, the IERK2-2 methods (3.6) would be more competitive than the IERK2-Radau methods (3.8), at least, from the perspective of average dissipation rate.

To improve the average dissipation rate of Radau-type IERK2 methods, one can consider four-stage procedure such that there are nine independent coefficients to be determined. An example is the four-stage Radau-type IERK2 procedure in [31, Section 4], which is constructed by imposing the simplifying setting $A_{\mathrm{I}}=A_{\mathrm{E}}$ and has a better average dissipation rate $R_{\mathrm{rad},S}^{(2)}=\tfrac{3}{2}+\tfrac{1}{2}\tau\overline{\lambda}_{\mathrm{ML}}$ . It is to mention that, the 3-stage Radau-type IERK method in [2] would not be suitable for the gradient flow problems since it does not satisfy the assumptions of Corollary 2.1.

3.3 Tests of IERK2 methods

Example 1.

Consider the Cahn-Hilliard model, $\partial_{t}u=\partial_{xx}(-\epsilon^{2}\partial_{xx}u-u+u^{3})+f$ , subject to the initial data $u_{0}=\sin x$ on $\Omega=(0,2\pi)$ with the interface parameter $\epsilon=0.2$ . The source term $f$ is set by choosing the exact solution $u=e^{-t}\sin x$ . Always, the spatial operators are approximated by the Fourier pseudo-spectral approximation with 256 grid points.

At first, we test the convergence of IERK2-1 (3.2), IERK2-2 (3.6) and IERK2-Radau (3.8) methods by choosing the final time $T=1$ and the stabilized parameter $\kappa=4$ . Figure 1 lists the $L^{\infty}$ norm error $e(\tau):=\max_{1\leq n\leq N}\left\|u_{h}^{n}-u(t_{n})\right\|_{\infty}$ for the three classes of IERK2 methods on halving time steps $\tau=2^{-k}/10$ for $0\leq k\leq 9$ . As expected, the IERK2-1 (3.2), IERK2-2 (3.6) and IERK2-Radau (3.8) methods are second-order accurate in time. It seems that the different parameters for the IERK2-1 (3.2) and IERK2-2 (3.6) methods would arrive at different numerical precision; while the IERK2-Radau methods (3.8) with two different parameters $c_{2}=1.5$ and $c_{2}=2$ generates almost the same solution. Also, it can be observed that the error of IERK2-2 method for $a_{33}=\tfrac{1+\sqrt{2}}{4}\approx 0.6036$ is always the smallest for the same time-step size $\tau$ .

Refer to caption — (a) IERK2-1 (3.2) with $c_{2}=1$

Example 2.

Consider the Cahn-Hilliard model, $\partial_{t}u=\partial_{xx}(-\epsilon^{2}\partial_{xx}u-u+u^{3})$ on $\Omega=(-\pi,\pi)$ with the interface parameter $\epsilon=0.1$ , subject to the initial data

\displaystyle u_{0}=

\displaystyle\,\tfrac{1}{3}\tanh(2\sin x)-\tfrac{1}{10}e^{-23.5(\left|x\right|-1)^{2}}+e^{-27(\left|x\right|-4.2)^{2}}+e^{-38(\left|x\right|-5.4)^{2}}.

The reference solution is generated with a small time-step size $\tau=10^{-3}$ by the IERK2-1 method (3.2) for the parameters $c_{2}=1$ and $a_{33}=1$ .

Figure 2 depicts the final solution $u_{h}^{N}$ at $T=150$ for $\tau=0.01,\kappa=2$ , and the discrete energy $E[u_{h}^{n}]$ generated by the IERK2-1 methods (3.2) for three different scenes: (i) $\tau=0.01,\kappa=2$ ; (ii) $\tau=0.05,\kappa=2$ and (iii) $\tau=0.05,\kappa=3$ . As predicted by Theorem 3.1, the original energy curves in all cases always monotonically decreasing. Also, the numerical dissipation rates of original energy seem quite different for three different choices of the parameter $a_{33}$ . For three different discrete scenes in Figure 2 (b)-(d), the energy curves for $a_{33}=0.5$ are always closest to the reference energy, while the energies for $a_{33}=2$ are the farthest away. More interestingly, these phenomena can be predicted by (3.5), or $\mathcal{R}^{(2,1)}(1,a_{33})=\tfrac{3}{2}+(2a_{33}-\tfrac{1}{2})\tau\overline{\lambda}_{\mathrm{ML}}$ , which suggests that the discrete energy gradually deviates from the continuous energy as the method parameter $a_{33}$ increases.

It seems that, the above differences between discrete energy curves and continuous energy shown in Figure 2 (b)-(d) can be explained by the different precision of numerical solutions. Actually, Figure 1 (a) shows that the solution for the case $a_{33}=0.5$ is a bit more accurate than that for $a_{33}=2$ although both of them are second-order accurate. Maybe, this would not be the whole story. In Figure 2 (b)-(d), we also include the discrete energy generated by the first-order IERK1 method (2.16) with $\theta=1/2$ , which has the average dissipation rate $\mathcal{R}^{(1)}(\tfrac{1}{2})=1$ . More surprisingly, the energy curve generated by the first-order IERK1 method with $\theta=\tfrac{1}{2}$ is even closer to the reference one than some IERK2-1 schemes especially when the time-step size is properly large. Although we are unsure of the complete mechanism behind it, they suggest that the selection of method parameters in IERK methods is at least as important as the selection of different IERK methods, if not more important. Obviously, the choice $a_{33}=\frac{1}{2}$ is a good choice for IERK2-1 methods (3.2), at least, for this example.

We also run the one-parameterized IERK2-2 (3.6) and IERK2-Radau (3.8) methods up to the time $T=150$ , and present the discrete energy curves in Figures 3 and 4, respectively, for three different scenes: (i) $\tau=0.01,\kappa=2$ ; (ii) $\tau=0.05,\kappa=2$ and (iii) $\tau=0.05,\kappa=3$ . As predicted by our theory, the discrete energies in all cases always monotonically decrease, and the discrete energies generated by the same method gradually deviate from the continuous energy as the associated average dissipation rate deviates from 1 (reminding that the average dissipation rate is an increasing function with respect to the step size $\tau$ and the parameter $\kappa$ ), cf. the formulas (3.7) and (3.9). For this example, one can find that the case $a_{33}=\tfrac{1+\sqrt{2}}{4}\approx 0.6036$ of IERK2-2 methods (3.6), and the case $c_{2}=\frac{3}{2}$ of IERK2-Radau methods (3.8) seem better than other parameters for the same time-step size and the stabilized parameter. Again, from the difference between discrete energy and continuous energy (for the same time-step size), the IERK1 scheme (2.16) with $\theta=\frac{1}{2}$ is also comparable to some second-order IERK2 schemes, although it only has first-order time accuracy.

Also, in Figure 5, we compare the energy behaviors of three “good” schemes, namely, the IERK2-1 method for $c_{2}=1$ and $a_{33}=\frac{1}{2}$ , the IERK2-2 method for $a_{33}=\tfrac{1+\sqrt{2}}{4}\approx 0.6036$ , and the IERK2-Radau method for $c_{2}=\frac{3}{2}$ , in previous tests shown in Figures 2-4. We see that, for this example, the discrete energy generated by the IERK2-2 scheme is the closest to the continuous energy; while the discrete energy generated by the IERK2-Radau method is the farthest from the continuous one. Experimentally, they confirm our theoretical predictions in previous subsection.

Table 2: Average dissipation rate of second-order IERK methods.

Type	Method	Average dissipation rate	Best choice
Lobatto-type	3-stage IERK2-1 (3.2)	$c_{2}+\tfrac{1}{2c_{2}}+(2c_{2}a_{33}-\tfrac{1}{2})\tau\overline{\lambda}_{\mathrm{ML}}$	$a_{33}=\frac{1}{2(4c_{2}-2c_{2}^{2}-1)}$
	3-stage IERK2-2 (3.6)	$\sqrt{2}+\sqrt{2}(a_{33}-\tfrac{\sqrt{2}}{4})\tau\overline{\lambda}_{\mathrm{ML}}$	$a_{33}=\tfrac{1+\sqrt{2}}{4}$
	3-stage IERK2 [6]	NPD^∗
Radau-type	3-stage IERK2-Radau (3.8)	$\tfrac{1}{2c_{2}}+c_{2}+\big{(}2c_{2}+1+\tfrac{1}{c_{2}-1}\big{)}\tau\overline{\lambda}_{\mathrm{ML}}$	$c_{2}=1+\frac{\sqrt{2}}{2}$
	4-stage IERK2 [31]	$\tfrac{3}{2}+\tfrac{1}{2}\tau\overline{\lambda}_{\mathrm{ML}}$
	3-stage IERK2 [2]	NPD

^{* NPD means there exists a $z_{0}<0$ such that the associated differential matrix $D(z_{0})$ is not positive (semi-)definite.}

To end this section, we present the average dissipation rates and practical parameter choices in Table 2, which summarizes the above numerical results, the theoretical results in Theorems 3.1-3.2, and some results of existing IERK algorithms in the literature.

4 Discrete energy laws of IERK3 methods

4.1 IERK3 methods with four stages

Third-order methods require three implicit stages, $s_{I}=3$ , at least. We should determine eleven independent coefficients in the following four-stage IERK methods that satisfy the canopy node condition and the two order conditions for first-order accuracy,

\begin{array}[]{c|c}\mathbf{c}&A\\ \hline\cr\\[-10.0pt] &\mathbf{b}^{T}\end{array}=\begin{array}[]{c|cccc}0&0&&&\\ c_{2}&c_{2}-a_{22}&a_{22}&\\ c_{3}&c_{3}-a_{32}-a_{33}&a_{32}&a_{33}\\ 1&1-a_{42}-a_{43}-a_{44}&a_{42}&a_{43}&a_{44}\\ \hline\cr\\[-10.0pt] &1-a_{42}-a_{43}-a_{44}&a_{42}&a_{43}&a_{44}\end{array},\begin{array}[]{c|c}\hat{\mathbf{c}}&\widehat{A}\\ \hline\cr\\[-10.0pt] &\hat{\mathbf{b}}^{T}\end{array}=\begin{array}[]{c|cccc}0&0&&&\\ c_{2}&c_{2}&0&&\\ c_{3}&c_{3}-\hat{a}_{32}&\hat{a}_{32}&0&\\ 1&1-\hat{a}_{42}-\hat{a}_{43}&\hat{a}_{42}&\hat{a}_{43}&0\\ \hline\cr\\[-10.0pt] &1-\hat{a}_{42}-\hat{a}_{43}&\hat{a}_{42}&\hat{a}_{43}&0\end{array}\;.

From Table 1, eight order conditions for second- and third-order accuracy are required such that there are three independent coefficients. Considering the stand-alone conditions for explicit part, $\hat{\mathbf{b}}^{T}\widehat{A}\mathbf{c}=\tfrac{1}{6}$ , $\hat{\mathbf{b}}^{T}\mathbf{c}^{.2}=\tfrac{1}{3}$ and $\hat{\mathbf{b}}^{T}\mathbf{c}=\tfrac{1}{2}$ , we find that $\hat{a}_{42}=\frac{3c_{3}-2}{6c_{2}(c_{3}-c_{2})}$ , $\hat{a}_{43}=\frac{2-3c_{2}}{6c_{3}(c_{3}-c_{2})}$ and $\hat{a}_{32}=\frac{c_{3}(c_{3}-c_{2})}{c_{2}(2-3c_{2})}$ with the independent variables $c_{2}$ and $c_{3}$ . Here we have excluded the case $c_{2}(c_{3}-c_{2})=0$ since the case $c_{2}=c_{3}$ can not arrive at our aim. From the stand-alone conditions for implicit part, ${\mathbf{b}}^{T}\mathbf{c}^{.2}=\tfrac{1}{3}$ and ${\mathbf{b}}^{T}\mathbf{c}=\tfrac{1}{2}$ , one has ${a}_{42}=\frac{6{a}_{44}c_{3}-6{a}_{44}-3c_{3}+2}{6c_{2}(c_{2}-c_{3})}$ and ${a}_{43}=\frac{6{a}_{44}c_{2}-6{a}_{44}-3c_{2}+2}{6c_{3}(c_{3}-c_{2})}.$ Thus the coupling condition $\mathbf{b}^{T}\widehat{A}\mathbf{c}=\tfrac{1}{6}$ arrives at $\tfrac{{a}_{44}c_{2}}{3c_{2}-2}=0$ and then

a_{44}=0,

that is, the procedure at the fourth stage is purely explicit. Then it follows that ${a}_{42}=\hat{a}_{42}$ and ${a}_{43}=\hat{a}_{43}$ . It means that the coupling condition $\hat{\mathbf{b}}^{T}A\mathbf{c}=\tfrac{1}{6}$ is equivalent to the stand-alone condition for implicit part $\mathbf{b}^{T}A\mathbf{c}=\tfrac{1}{6}$ . Therefore we have four independent coefficients.

By choosing $c_{2}=\frac{1}{3}$ , $c_{3}=\frac{2}{3}$ and $a_{33}=\frac{1}{3}$ , we obtain the following $a_{22}$ -parameterized IERK3 methods with the associated Butcher tableaux

\displaystyle\begin{array}[]{c|c}\mathbf{c}&A\\ \hline\cr\\[-10.0pt] &\mathbf{b}^{T}\end{array}=\begin{array}[]{c|cccc}0&0&&&\\ \tfrac{1}{3}&\tfrac{1}{3}-a_{22}&a_{22}\\[2.0pt] \tfrac{2}{3}&\tfrac{1}{3}&0&\tfrac{1}{3}\\[2.0pt] 1&\tfrac{1}{4}&0&\tfrac{3}{4}&0\\[1.0pt] \hline\cr\\[-10.0pt] &\tfrac{1}{4}&0&\tfrac{3}{4}&0\end{array}\,,\quad\begin{array}[]{c|c}\hat{\mathbf{c}}&\widehat{A}\\ \hline\cr\\[-10.0pt] &\hat{\mathbf{b}}^{T}\end{array}=\begin{array}[]{c|cccc}0&0\\ \tfrac{1}{3}&\tfrac{1}{3}&0\\[2.0pt] \tfrac{2}{3}&0&\tfrac{2}{3}&0\\[2.0pt] 1&\tfrac{1}{4}&0&\tfrac{3}{4}&0\\[1.0pt] \hline\cr\\[-10.0pt] &\tfrac{1}{4}&0&\tfrac{3}{4}&0\end{array}\;.

(4.19)

By the definition (2.9), one can get

\displaystyle D_{\mathrm{E}}^{(3)}:=\begin{pmatrix}3&0&0\\[3.0pt] \tfrac{3}{2}&\tfrac{3}{2}&0\\[3.0pt] \tfrac{1}{3}&\tfrac{4}{3}&\tfrac{4}{3}\end{pmatrix},\quad D_{\mathrm{EI}}^{(3)}:=\begin{pmatrix}3a_{22}-\tfrac{1}{2}&0&0\\[3.0pt] -\tfrac{1}{2}&0&0\\[3.0pt] -a_{22}&0&-\tfrac{1}{2}\end{pmatrix}.

Simple computing yields that the matrix $D_{\mathrm{E}}^{(3)}$ is positive definite. Also, it is easy to find that $\mathrm{Det}\big{[}\mathcal{S}(D_{\mathrm{EI},1}^{(3)};a_{22})\big{]}=3a_{22}-\frac{1}{2}$ , $\mathrm{Det}\big{[}\mathcal{S}(D_{\mathrm{EI},2}^{(3)};a_{22})\big{]}=-\tfrac{1}{16}$ and $\mathrm{Det}\big{[}\mathcal{S}(D_{\mathrm{EI}}^{(3)};a_{22})\big{]}=\tfrac{1}{32}$ . That is, the associated differentiation matrix $D^{(3)}(z)=D_{\mathrm{E}}^{(3)}-D_{\mathrm{EI}}^{(3)}z$ can not be positive semi-definite for all $z<0$ , and the sufficient condition of Theorem 2.1 is not fulfilled. Although the third-order methods (4.19) may maintain the discrete energy dissipation law (1.3) under certain time-step conditions, they may not be stabilized to preserve the energy dissipation law unconditionally no matter how large the stabilization parameter $\kappa$ we set in the current stabilization strategy (2.3).

As the numerical evidence, we run the IERK3 methods (4.19) with $a_{22}=1$ , 2 and 3 for Example 2 with a fixed stabilized factor $\kappa=4$ . The associated discrete energy curves, depicted in Figure 6, for three different time-steps $\tau=0.01$ , $0.005$ and $0.001$ can not maintain the decaying property. It is to mention that, the differentiation matrix $D(z)=D_{\mathrm{E}}-D_{\mathrm{EI}}z$ can not be positive semi-definite for all $z<0$ if there are one zero diagonal element in the coefficient matrix $A$ of the implicit part. For example, if $a_{k_{0}+1,k_{0}+1}=0$ , the formula (2.9) gives the diagonal element $d_{k_{0},k_{0}}=\tfrac{1}{\hat{a}_{k_{0}+1,k_{0}}}+\tfrac{z}{2}$ , which cannot be always non-negative for all $z<0$ . Then the condition of Theorem 2.1 can not be fulfilled because the diagonal elements of a positive semi-definite matrix must be non-negative. As mentioned before, four-stage IERK3 methods always have a zero diagonal element, $a_{44}\equiv 0$ . In this sense, we say that there are no four-stage IERK3 methods fulfill the condition of Theorem 2.1. Note that, similar situations can be found in [6, Section 2] and [26, Sections 8 and 9], where many IERK methods up to fourth-order accuracy have zero diagonal elements.

4.2 IERK3 methods with five stages

We consider the five-stage third-order IERK methods with the following Butcher tableaux

\begin{array}[]{c|ccccc}0&0&&&&\\ c_{2}&c_{2}-a_{55}&a_{55}&&\\ c_{3}&c_{3}-a_{32}-a_{55}&a_{32}&a_{55}&\\ c_{4}&c_{4}-a_{42}-a_{43}-a_{55}&a_{42}&a_{43}&a_{55}\\ 1&1-a_{52}-a_{53}-a_{54}-a_{55}&a_{52}&a_{53}&a_{54}&a_{55}\\ \hline\cr\\[-10.0pt] &1-a_{52}-a_{53}-a_{54}-a_{55}&a_{52}&a_{53}&a_{54}&a_{55}\end{array},\quad\begin{array}[]{c|ccccc}0&0&&&&\\ c_{2}&c_{2}&0&&&\\ c_{3}&c_{3}-c_{2}&c_{2}&0&&\\ c_{4}&c_{4}-\hat{a}_{42}-c_{2}&\hat{a}_{42}&c_{2}&0&\\ 1&1-\hat{a}_{52}-\hat{a}_{53}-c_{2}&\hat{a}_{52}&\hat{a}_{53}&c_{2}&0\\ \hline\cr\\[-10.0pt] &1-\hat{a}_{52}-\hat{a}_{53}-c_{2}&\hat{a}_{52}&\hat{a}_{53}&c_{2}&0\end{array},

where we assume $a_{k+1,k+1}=a_{55}$ and $\hat{a}_{k+1,k}=c_{2}$ for $1\leq k\leq 4$ to reduce the degree of freedom so that there are five independent variables. Such simplifying assumptions are set to maintain a fixed diagonal ratio $\frac{a_{k+1,k+1}}{\hat{a}_{k+1,k}}$ between the implicit and explicit parts so that one can minimize the value of $\frac{1}{s_{\mathrm{I}}}\mathrm{tr}(D_{\mathrm{EI}})$ as far as possible according to the formula (2.15) of average dissipation rate.

After lots of trial and error processes, we keep one free variable $a_{55}$ by choosing the abscissas $c_{2}=\tfrac{4}{5}$ , $c_{3}=\tfrac{7}{5}$ , $c_{4}=\tfrac{6}{5}$ and the coefficient $a_{43}=-\tfrac{3}{5}$ , and obtain the five-stage IERK3 (IERK3-1) methods with the following Butcher tableaux

\displaystyle\begin{array}[]{c|ccccc}0&0&&&\\[3.0pt] \tfrac{4}{5}&\tfrac{4}{5}-a_{55}&a_{55}\\[3.0pt] \tfrac{7}{5}&\tfrac{3}{5}-\tfrac{5a_{55}}{16}&\tfrac{4}{5}-\tfrac{11a_{55}}{16}&a_{55}\\[3.0pt] \tfrac{6}{5}&\tfrac{977a_{55}}{4032}-\tfrac{473}{10080}&\tfrac{18617}{10080}-\tfrac{5009a_{55}}{4032}&-\tfrac{3}{5}&a_{55}\\[3.0pt] 1&a_{51}&a_{52}&a_{53}&a_{54}&a_{55}\\[3.0pt] \hline\cr\\[-10.0pt] &a_{51}&a_{52}&a_{53}&a_{54}&a_{55}\end{array}\,,\quad\begin{array}[]{c|ccccc}0&0&&&\\[3.0pt] \tfrac{4}{5}&\tfrac{4}{5}&0\\[3.0pt] \tfrac{7}{5}&\tfrac{3}{5}&\tfrac{4}{5}&0\\[3.0pt] \tfrac{6}{5}&\tfrac{10111}{10080}&-\tfrac{6079}{10080}&\tfrac{4}{5}&0\\[3.0pt] 1&\tfrac{313}{840}&\tfrac{131}{360}&-\tfrac{169}{315}&\tfrac{4}{5}&0\\[3.0pt] \hline\cr\\[-10.0pt] &\tfrac{313}{840}&\tfrac{131}{360}&-\tfrac{169}{315}&\tfrac{4}{5}&0\end{array}\,,

(4.34)

where $a_{51}=\tfrac{313}{840}-\tfrac{191a_{55}}{9590}$ , $a_{52}=\tfrac{131}{360}-\tfrac{797a_{55}}{4110}$ , $a_{53}=\tfrac{7087a_{55}}{14385}-\tfrac{169}{315}$ and $a_{54}=\tfrac{4}{5}-\tfrac{876a_{55}}{685}$ .

By the definitions in (2.10), one can write out the matrices $D_{\mathrm{E}}^{(3,1)}$ and $D_{\mathrm{EI}}^{(3,1)}$ as follows,

\displaystyle D_{\mathrm{E}}^{(3,1)}:=\begin{pmatrix}\tfrac{5}{6}&0&0&0\\[3.0pt] \tfrac{49}{30}&\tfrac{48}{35}&0&0\\[3.0pt] \tfrac{23}{9}&\tfrac{62}{21}&\tfrac{5}{2}&0\\[3.0pt] \tfrac{583}{1134}&\tfrac{851}{1323}&\tfrac{125}{126}&\tfrac{5}{3}\end{pmatrix},D_{\mathrm{EI}}^{(3,1)}:=\begin{pmatrix}\tfrac{5a_{55}-2}{4}&0&0&0\\[3.0pt] -\tfrac{35a_{55}}{64}&\tfrac{5a_{55}-2}{4}&0&0\\[3.0pt] \tfrac{21}{16}-\tfrac{130885a_{55}}{57344}&\tfrac{70715a_{55}}{32256}-\tfrac{7}{4}&\tfrac{5a_{55}-2}{4}&0\\[3.0pt] \tfrac{169}{192}-\tfrac{1213756279a_{55}}{659914752}&\tfrac{1301752603a_{55}}{1113606144}-\tfrac{169}{144}&\tfrac{67633a_{55}}{138096}&\tfrac{5a_{55}-2}{4}\end{pmatrix}.

Simple calculations can verify the positive semi-definiteness of $D_{\mathrm{E}}^{(3,1)}$ . For the matrix $D_{\mathrm{EI}}^{(3,1)}(a_{55})$ , We have $\mathrm{Det}\big{[}\mathcal{S}(D_{\mathrm{EI},1}^{(3,1)};a_{55})\big{]}=\tfrac{5a_{55}-2}{4}$ , $\mathrm{Det}\big{[}\mathcal{S}(D_{\mathrm{EI},2}^{(3,1)};a_{55})\big{]}=\tfrac{24375a_{55}^{2}}{16384}-\tfrac{5a_{55}}{4}+\tfrac{1}{4}$ , and

	$\displaystyle\mathrm{Det}\big{[}\mathcal{S}(D_{\mathrm{EI},3}^{(3,1)};a_{55})\big{]}=$	$\displaystyle\,-\tfrac{156124878125a_{55}^{3}}{266355081216}+\tfrac{4850385555625a_{55}^{2}}{2130840649728}-\tfrac{9210215a_{55}}{4718592}+\tfrac{969}{2048},$
	$\displaystyle\mathrm{Det}\big{[}\mathcal{S}(D_{\mathrm{EI},4}^{(3,1)};a_{55})\big{]}=$	$\displaystyle\,-\tfrac{7024871482286543654375a_{55}^{4}}{5079525965637831622656}+\tfrac{11933410165242198623465a_{55}^{3}}{2539762982818915811328}$
		$\displaystyle\,-\tfrac{26459898695651552391913a_{55}^{2}}{5079525965637831622656}+\tfrac{192673313809999a_{55}}{82103953784832}-\tfrac{1969849}{5308416}.$

One can check that the following condition

\displaystyle 0.717374\leq a_{55}\leq 1.74727

(4.35)

is sufficient for the positive semi-definiteness of $D_{\mathrm{EI}}^{(3,1)}$ . By using Corollary 2.1 and Lemma 2.2, we have the following theorem.

Theorem 4.1.

In simulating the gradient flow system (2.3) with $\kappa\geq 2\ell_{g}$ , the one-parameter IERK3-1 methods (4.34) with the parameter setting (4.35) preserve the original energy dissipation law (1.3) unconditionally at all stages in the sense that

\displaystyle E[U^{n,j+1}]-E[U^{n,1}]\leq\frac{1}{\tau}\Big{\langle}M_{h}^{-1}\delta_{\tau}\vec{U}_{n,j+1},D_{j}^{(3,1)}(a_{55};{\tau}M_{h}L_{\kappa})\delta_{\tau}\vec{U}_{n,j+1}\Big{\rangle}\quad\text{for $1\leq j\leq 4$,}

where the differentiation matrix $D^{(3,1)}$ is defined by

D^{(3,1)}(a_{55};{\tau}M_{h}L_{\kappa}):=D_{\mathrm{E}}^{(3,1)}\otimes I+D_{\mathrm{EI}}^{(3,1)}(a_{55})\otimes(-\tau M_{h}L_{\kappa}).

The associated average dissipation rate reads

\displaystyle\mathcal{R}^{(3,1)}(a_{55}):=\tfrac{5}{4}+\tfrac{5a_{55}-2}{4}\tau\overline{\lambda}_{\mathrm{ML}}\quad\text{for $0.717374\leq a_{55}\leq 1.74727$}.

(4.36)

Also, one can take $a_{43}$ as the free variable by choosing the abscissas $c_{2}=\tfrac{4}{5}$ , $c_{3}=\tfrac{7}{5}$ , $c_{4}=\tfrac{6}{5}$ and the diagonal element $a_{55}=\tfrac{18}{25}$ , and obtain the following one-parameter IERK3 (IERK3-2) methods with the associated Butcher tableaux

\begin{array}[]{c|ccccc}0&0\\[3.0pt] \tfrac{4}{5}&\tfrac{2}{25}&\tfrac{18}{25}\\[3.0pt] \tfrac{7}{5}&\tfrac{3}{8}&\tfrac{61}{200}&\tfrac{18}{25}\\[3.0pt] \tfrac{6}{5}&\tfrac{3a_{43}}{4}+\tfrac{7277}{12600}&-\tfrac{7a_{43}}{4}-\tfrac{1229}{12600}&a_{43}&\tfrac{18}{25}\\[3.0pt] 1&\tfrac{1030769}{2877000}&\tfrac{276523}{1233000}&-\tfrac{196127}{1078875}&-\tfrac{2068}{17125}&\tfrac{18}{25}\\[3.0pt] \hline\cr\\[-10.0pt] &\tfrac{1030769}{2877000}&\tfrac{276523}{1233000}&-\tfrac{196127}{1078875}&-\tfrac{2068}{17125}&\tfrac{18}{25}\end{array},\;\;\;\begin{array}[]{c|ccccc}0&0&&&\\[3.0pt] \tfrac{4}{5}&\tfrac{4}{5}&0\\[3.0pt] \tfrac{7}{5}&\tfrac{3}{5}&\tfrac{4}{5}&0\\[3.0pt] \tfrac{6}{5}&\tfrac{10111}{10080}&-\tfrac{6079}{10080}&\tfrac{4}{5}&0\\[3.0pt] 1&\tfrac{313}{840}&\tfrac{131}{360}&-\tfrac{169}{315}&\tfrac{4}{5}&0\\[3.0pt] \hline\cr\\[-10.0pt] &\tfrac{313}{840}&\tfrac{131}{360}&-\tfrac{169}{315}&\tfrac{4}{5}&0\end{array}\,.

(4.37)

Note that, the IERK3-2 methods (4.37) is quite different from the IERK3-1 methods (4.34) in the sense that the average dissipation rate $\mathcal{R}^{(3,2)}$ is independent of the parameter $a_{43}$ , while the average dissipation rate $\mathcal{R}^{(3,1)}(a_{55})$ of the latter is always dependent on the choice of $a_{55}$ .

By the definitions in (2.10), one has $D_{\mathrm{E}}^{(3,2)}:=D_{\mathrm{E}}^{(3,1)}$ and

\displaystyle D_{\mathrm{EI}}^{(3,2)}(a_{43}):=\begin{pmatrix}\tfrac{2}{5}&0&0&0\\[3.0pt] -\tfrac{63}{160}&\tfrac{2}{5}&0&0\\[3.0pt] -\tfrac{15a_{43}}{16}-\tfrac{128073}{143360}&\tfrac{5a_{43}}{4}+\tfrac{5183}{8960}&\tfrac{2}{5}&0\\[3.0pt] -\tfrac{845a_{43}}{1344}-\tfrac{6774788911}{8248934400}&\tfrac{845a_{43}}{1008}+\tfrac{264498203}{1546675200}&\tfrac{67633}{191800}&\tfrac{2}{5}\end{pmatrix}\,.

Simple calculations gives $\mathrm{Det}\big{[}\mathcal{S}(D_{\mathrm{EI},1}^{(3,2)};a_{43})\big{]}=\tfrac{2}{5}$ , $\mathrm{Det}\big{[}\mathcal{S}(D_{\mathrm{EI},2}^{(3,2)};a_{43})\big{]}=\tfrac{2483}{20480}$ ,

	$\displaystyle\mathrm{Det}\big{[}\mathcal{S}(D_{\mathrm{EI},3}^{(3,2)};a_{43})\big{]}=$	$\displaystyle-\tfrac{1055a_{43}^{2}}{8192}-\tfrac{2184569a_{43}}{14680064}-\tfrac{91452759}{6576668672},$
	$\displaystyle\mathrm{Det}\big{[}\mathcal{S}(D_{\mathrm{EI}}^{(3,2)};a_{43})\big{]}=$	$\displaystyle-\tfrac{14878761752461399a_{43}^{2}}{499921851934310400}-\tfrac{523061384625360827a_{43}}{17497264817700864000}-\tfrac{342702875992479487397}{48992341489562419200000}.$

It is easy to check that the following condition

\displaystyle-0.633312\leq a_{43}\leq-0.371114

(4.38)

is sufficient for the positive semi-definiteness of $D_{\mathrm{EI}}^{(3,2)}$ . Thus we have the following result.

Theorem 4.2.

In simulating the gradient flow system (2.3) with $\kappa\geq 2\ell_{g}$ , the one-parameter IERK3-2 methods (4.37) with the parameter setting (4.38) preserve the original energy dissipation law (1.3) unconditionally at all stages in the sense that

\displaystyle E[U^{n,j+1}]-E[U^{n,1}]\leq\frac{1}{\tau}\Big{\langle}M_{h}^{-1}\delta_{\tau}\vec{U}_{n,j+1},D_{j}^{(3,2)}(a_{43};{\tau}M_{h}L_{\kappa})\delta_{\tau}\vec{U}_{n,j+1}\Big{\rangle}\quad\text{for $1\leq j\leq 4$,}

where the differentiation matrix $D^{(3,2)}$ is defined by

D^{(3,2)}(a_{43};{\tau}M_{h}L_{\kappa}):=D_{\mathrm{E}}^{(3,2)}\otimes I+D_{\mathrm{EI}}^{(3,2)}(a_{43})\otimes(-\tau M_{h}L_{\kappa}).

The associated average dissipation rate is parameter-independent, that is,

\displaystyle\mathcal{R}^{(3,2)}:=\tfrac{5}{4}+\tfrac{2}{5}\tau\overline{\lambda}_{\mathrm{ML}}\quad\text{for $-0.633312\leq a_{43}\leq-0.371114$}.

(4.39)

4.3 Radau-type IERK3 methods with five stages

We also briefly discuss the five-stage Radau-type IERK3 methods by assuming $\mathbf{a}_{1}=\mathbf{0}$ . That is, consider the following Butcher tableaux that satisfy the canopy node condition and the two order conditions for first-order accuracy,

\begin{array}[]{c|ccccc}0&0&&&&\\ c_{2}&0&c_{2}&&\\ c_{3}&0&c_{3}-c_{2}&c_{2}&\\ c_{4}&0&c_{4}-a_{43}-c_{2}&a_{43}&c_{2}\\ 1&0&1-a_{53}-a_{54}-c_{2}&a_{53}&a_{54}&c_{2}\\ \hline\cr\\[-10.0pt] &0&1-a_{53}-a_{54}-c_{2}&a_{53}&a_{54}&c_{2}\end{array},\quad\begin{array}[]{c|ccccc}0&0&&&&\\ c_{2}&c_{2}&0&&&\\ c_{3}&c_{3}-\hat{a}_{32}&\hat{a}_{32}&0&&\\ c_{4}&c_{4}-\hat{a}_{42}-\hat{a}_{43}&\hat{a}_{42}&\hat{a}_{43}&0&\\ 1&1-\hat{a}_{52}-\hat{a}_{53}-\hat{a}_{54}&\hat{a}_{52}&\hat{a}_{53}&\hat{a}_{54}&0\\ \hline\cr\\[-10.0pt] &1-\hat{a}_{52}-\hat{a}_{53}-\hat{a}_{54}&\hat{a}_{52}&\hat{a}_{53}&\hat{a}_{54}&0\end{array},

where we set $a_{k+1,k+1}=c_{2}$ ( $1\leq k\leq 4$ ) to reduce the degree of freedom. According to Table 1, eight order conditions for second- and third-order accuracy are required so that there are four independent coefficients. This simplifying setting is also useful for practical applications since it provides the same iteration matrix for the systems of linear equations at all stages.

After lots of trial and error processes, we choose $c_{2}=\tfrac{4}{5}$ , $c_{3}=\tfrac{93}{200}$ and $c_{4}=\frac{171}{200}$ and obtain the $\hat{a}_{43}$ -parameterized Radau-type IERK3 (IERK3-Radau) methods with the Butcher tableaux

\displaystyle\begin{array}[]{c|ccccc}0&0&\\[3.0pt] \tfrac{4}{5}&0&\tfrac{4}{5}\\[3.0pt] \tfrac{93}{200}&0&-\tfrac{67}{200}&\tfrac{4}{5}\\[3.0pt] \tfrac{171}{200}&0&-\tfrac{9361649}{5132200}&\tfrac{241098}{128305}&\tfrac{4}{5}\\[3.0pt] 1&0&-\tfrac{5309}{11055}&\tfrac{9998}{7839}&-\tfrac{766}{1287}&\tfrac{4}{5}\\[3.0pt] \hline\cr\\[-10.0pt] &0&\tfrac{25}{1488}&\tfrac{4024}{1395}&-\tfrac{1397}{144}&\tfrac{39}{5}\end{array},\quad\begin{array}[]{c|ccccc}0&0&\\[3.0pt] \tfrac{4}{5}&\tfrac{4}{5}&0\\[3.0pt] \tfrac{93}{200}&\tfrac{10391}{32000}&\tfrac{4489}{32000}&0\\[3.0pt] \tfrac{171}{200}&\hat{a}_{41}&\hat{a}_{42}&\hat{a}_{43}&0\\[3.0pt] 1&\tfrac{2053}{11066}&\tfrac{3785983}{24466926}&\tfrac{20893310}{43373187}&\tfrac{1267730}{7120971}&0\\[3.0pt] \hline\cr\\[-10.0pt] &\tfrac{2053}{11066}&\tfrac{3785983}{24466926}&\tfrac{20893310}{43373187}&\tfrac{1267730}{7120971}&0\end{array},

(4.54)

where $\hat{a}_{42}=\tfrac{9690263}{12256000}-\tfrac{93\hat{a}_{43}}{160}$ and $\hat{a}_{41}=\tfrac{171}{200}-\hat{a}_{42}-\hat{a}_{43}$ .

By the definition (2.10), one can get the associated matrices

\displaystyle[D_{\mathrm{rad},\mathrm{E}}^{(3)}(\hat{a}_{43})]^{T}:=\begin{pmatrix}\tfrac{5}{4}&\tfrac{1135}{268}&\tfrac{200}{67}-\tfrac{4986267}{2052880\hat{a}_{43}}&\tfrac{114597584391147}{17436733668080\hat{a}_{43}}-\tfrac{2528991857}{339751640}\\[3.0pt] 0&\tfrac{32000}{4489}&\tfrac{18600}{4489}-\tfrac{7970976}{1719287\hat{a}_{43}}&\tfrac{183194079827616}{14603264447017\hat{a}_{43}}-\tfrac{67097015181}{5690839970}\\[3.0pt] 0&0&\tfrac{1}{\hat{a}_{43}}&\tfrac{7120971}{1267730}-\tfrac{22982641}{8493791\hat{a}_{43}}\\[3.0pt] 0&0&0&\tfrac{7120971}{1267730}\end{pmatrix},

and

\displaystyle D_{\mathrm{rad},\mathrm{EI}}^{(3)}(\hat{a}_{43}):=\begin{pmatrix}\tfrac{1}{2}&0&0&0\\[3.0pt] 0&\tfrac{46711}{8978}&0&0\\[3.0pt] 0&\tfrac{10391}{4489}-\tfrac{15730338}{8596435\hat{a}_{43}}&\tfrac{4}{5\hat{a}_{43}}-\tfrac{1}{2}&0\\[3.0pt] 0&\tfrac{361524711062658}{73016322235085\hat{a}_{43}}-\tfrac{94060069247}{14227099925}&\tfrac{476922}{3169325}-\tfrac{91930564}{42468955\hat{a}_{43}}&\tfrac{25314559}{6338650}\end{pmatrix}.

As done in the above subsection, one can check that, if $0.598442\leq\hat{a}_{43}\leq 1.05134$ , $D_{\mathrm{rad},\mathrm{E}}^{(3)}(\hat{a}_{43})$ and $D_{\mathrm{rad},\mathrm{EI}}^{(3)}(\hat{a}_{43})$ are positive semi-definite. Corollary 2.1 and Lemma 2.2 yield the following theorem.

Theorem 4.3.

In simulating the gradient flow system (2.3) with $\kappa\geq 2\ell_{g}$ , the one-parameter IERK3-Radau methods (4.54) for $0.598442\leq\hat{a}_{43}\leq 1.05134$ preserve the original energy dissipation law (1.3) unconditionally at all stages in the sense that

\displaystyle E[U^{n,j+1}]-E[U^{n,1}]\leq\frac{1}{\tau}\Big{\langle}M_{h}^{-1}\delta_{\tau}\vec{U}_{n,j+1},D_{\mathrm{rad},j}^{(3)}(\hat{a}_{43};{\tau}M_{h}L_{\kappa})\delta_{\tau}\vec{U}_{n,j+1}\Big{\rangle}\quad\text{for $1\leq j\leq 3$,}

where the differentiation matrix $D_{\mathrm{rad}}^{(3)}$ is defined by

D_{\mathrm{rad}}^{(3)}(\hat{a}_{43};{\tau}M_{h}L_{\kappa}):=D_{{\mathrm{rad}},\mathrm{E}}^{(3)}(\hat{a}_{43})\otimes I+D_{{\mathrm{rad}},\mathrm{EI}}^{(3)}(\hat{a}_{43})\otimes(-\tau M_{h}L_{\kappa}).

The associated average dissipation rate is

\displaystyle\mathcal{R}_{\mathrm{rad}}^{(3)}(\hat{a}_{43}):=\tfrac{1}{4}\big{(}\tfrac{1}{\hat{a}_{43}}+\tfrac{159293897563}{11381679940}\big{)}+\big{(}\tfrac{1}{5\hat{a}_{43}}+\tfrac{130839697713}{56908399700}\big{)}\tau\overline{\lambda}_{\mathrm{ML}}\quad\text{for $0.598442\leq\hat{a}_{43}\leq 1.05134$.}

(4.55)

To enlarge the possible choices of $\tau\overline{\lambda}_{\mathrm{ML}}$ , one can choose the parameter $\hat{a}_{43}=1$ such that

\mathcal{R}_{\mathrm{rad}}^{(3)}(1)=\tfrac{170675577503}{45526719760}+\tfrac{142221377653}{56908399700}\tau\overline{\lambda}_{\mathrm{ML}}\approx 3.74891+2.49913\tau\overline{\lambda}_{\mathrm{ML}}.

It is a little larger than the average dissipation rate $\mathcal{R}^{(3,2)}=\tfrac{5}{4}+\tfrac{2}{5}\tau\overline{\lambda}_{\mathrm{ML}}$ of the IERK3-2 methods (4.37). From the perspective of average dissipation rate, the IERK3-Radau methods (4.54) are more competitive than the existing Radau-type IERK3 method in [11, Example 5] with the associated average dissipation rate $\mathcal{R}_{\mathrm{rad},F}^{(3)}\approx 3.24727+16.3779\tau\overline{\lambda}_{\mathrm{ML}}$ , although they would be inferior to the one-parameter IERK3-1 (4.34) and IERK3-2 (4.37) methods.

To search a “better” IERK3 method, that is to say, to obtain a IERK method with the average dissipation rate closer to 1, one can consider higher-stage procedures. An example is the 7-stage IERK3 method in [31], which achieves a better average dissipation rate, $\mathcal{R}_{\mathrm{rad},S}^{(3)}=2+\tfrac{1}{2}\tau\overline{\lambda}_{\mathrm{ML}}$ , at the cost of two additional implicit stages. We also mentioned that, the 5-stage Radau-type IERK method in [2] does not satisfy the sufficient conditions of Corollary 2.1 so that it is not a proper candidate for the gradient flow system (1.2).

4.4 Tests for IERK3 methods

We use Example 1 to test the convergence of IERK3-1 (4.34), IERK3-2 (4.37) and IERK3-Radau (4.54) methods by choosing the final time $T=1$ and the stabilized parameter $\kappa=4$ . Figure 7 lists the $L^{\infty}$ norm error $e(\tau):=\max_{1\leq n\leq N}\left\|u_{h}^{n}-u(t_{n})\right\|_{\infty}$ for the three classes of IERK3 methods on halving time steps $\tau=2^{-k}/10$ for $0\leq k\leq 9$ . As expected, the IERK3-1 (4.34), IERK3-2 (4.37) and IERK3-Radau (4.54) methods are third-order accurate in time. It is observed that the different parameters for the IERK3-1 (4.34) and IERK3-2 (4.37) methods would arrive at different numerical precision; while the IERK3-Radau methods (4.54) with three different parameters $\hat{a}_{43}=0.6,0.8$ and $1$ generates almost the same solution.

Next, we address the discrete energy behaviors of three IERK3 methods for Example 2. Always, we use a small time-step size $\tau=10^{-3}$ to compute the reference solution and energy by the IERK3-2 method (4.37) for $a_{43}=-0.4$ . Figures 8 and 9 depict the discrete energies generated by the IERK3-1 (4.34) and IERK3-Radau (4.54) methods, respectively, for three different scenes: (i) $\tau=0.01,\kappa=2$ ; (ii) $\tau=0.05,\kappa=2$ and (iii) $\tau=0.05,\kappa=3$ . As predicted by our theory, the discrete energies in all cases always monotonically decrease, and the discrete energies generated by the same method gradually deviate from the continuous energy as the associated average dissipation rates deviate from 1, cf. the formulas (4.36) and (4.55). It is seen that there are always significant differences between the reference curve and the three discrete energy curves (corresponds to three different parameters) generated by the IERK3-1 methods (4.34); but such differences are not obvious for those of the IERK3-Radau methods (4.54).

The above differences between discrete energy curves and continuous energy shown in Figures 8 and 9 can be partly explained by the different precision of numerical solutions. Actually, Figure 7 (a) shows that the solution for the case $a_{55}=0.8$ of the IERK3-1 methods (4.34) is a bit more accurate than that for $a_{55}=1.7$ although both of them are third-order accurate. Figure 7 (c) shows that the solutions for the three cases $\hat{a}_{43}=0.6,$ 0.8 and $1$ of the IERK3-Radau methods (4.54) are very close. At the same time, this would not be the whole story. In Figures 8 and 9, we also include the discrete energy generated by the IERK2-2 method (3.6) for $a_{33}=\tfrac{1+\sqrt{2}}{4}\approx 0.6036$ , which seems to be the “best” scheme of second-order accuracy. One can observe that the discrete energy curve generated by the IERK2-2 method is even closer to the reference energy curve than the IERK3-1 method (4.34) with $a_{55}=1.7$ . They suggest that the selection of method parameters in IERK3 methods is at least as important as the selection of different IERK3 methods. For this example, the parameters $a_{55}=0.8$ and $\hat{a}_{43}=1$ are “good” choices for the IERK3-1 (4.34) and IERK3-Radau (4.54) methods, respectively.

We run the IERK3-2 methods (4.37) with three parameters $a_{43}=-0.6,-0.5$ and $-0.4$ up to the time $T=150$ , and exhibit the discrete energy curves in Figure 10 for three scenes: (i) $\tau=0.01,\kappa=2$ ; (ii) $\tau=0.05,\kappa=2$ and (iii) $\tau=0.05,\kappa=3$ . As predicted by Theorem 4.37, the discrete energies in all cases always monotonically decreasing; however, the dissipation rates of discrete energies are quite different from those in Figures 8 and 9. We see that the discrete energies for three parameters $a_{43}=-0.6,-0.5$ and $-0.4$ approach the reference curve in a staggered manner (fast and slow) at different time periods. For example, the energy curves for $a_{43}=-0.6$ are always closest to the reference energy in the third fast varying period, while they are not the closest one to the reference energy in the second fast varying period. This phenomenon cannot be directly understood by the parameter-independent dissipation rate, $\mathcal{R}^{(3,2)}=\tfrac{5}{4}+\tfrac{2}{5}\tau\overline{\lambda}_{\mathrm{ML}}$ , and remains mysterious to us. As expected, the IERK3-2 methods (4.37) perform better than the “best” second-order method (3.6).

Furthermore, the IERK3-2 methods (4.37) would be the “best” candidates among the three classes of IERK3 methods. In Figure 11, we compare the energy behaviors of three “good” schemes, namely, the IERK3-1 method for $a_{55}=0.8$ , the IERK3-2 method for $a_{43}=-0.6$ , and the IERK3-Radau method for $\hat{a}_{43}=1$ , in previous tests shown in Figures 8-10. The discrete energy curve generated by the IERK3-1 method is always closer to the reference energy curve than that of the IERK3-Radau method. Thanks to the minimum dissipation rate $\mathcal{R}^{(3,2)}=\tfrac{5}{4}+\tfrac{2}{5}\tau\overline{\lambda}_{\mathrm{ML}}$ among the three methods, the energy curve of the IERK3-2 method for $a_{43}=-0.6$ is always the closest one to the reference curve. By contrast, the discrete energy generated by the IERK2-2 method stays farthest away from the reference energy. Experimentally, they support our theoretical results.

To end this section, we summarize the theoretical results in Theorems 4.1-4.3, the numerical results in this section and results of some existing IERK algorithms in the literature in Table 3, which presents the average dissipation rates and some practical choices.

Table 3: Average dissipation rate of third-order IERK methods.

Type	Method	Average dissipation rate	Practical choice
Lobatto-type	5-stage IERK3-1 (4.34)	$\tfrac{5}{4}+\tfrac{5a_{55}-2}{4}\tau\overline{\lambda}_{\mathrm{ML}}$	$a_{55}=0.8$
	5-stage IERK3-2 (4.37)	$\tfrac{5}{4}+\tfrac{2}{5}\tau\overline{\lambda}_{\mathrm{ML}}$	$a_{43}=-0.6$
	4-stage IERK3 [6]	NPD
	5-stage IERK3 [26]	NPD
Radau-type	5-stage IERK3-Radau (4.54)	$\tfrac{1}{4\hat{a}_{43}}+3.49891+\big{(}\tfrac{1}{5\hat{a}_{43}}+2.29913\big{)}\tau\overline{\lambda}_{\mathrm{ML}}$	$\hat{a}_{43}=1.0$
	5-stage IERK3 [11]	$3.24727+16.3779\tau\overline{\lambda}_{\mathrm{ML}}$
	7-stage IERK3 [31]	$2+\tfrac{1}{2}\tau\overline{\lambda}_{\mathrm{ML}}$
	5-stage IERK3 [2]	NPD

5 Discrete energy laws of IERK4-A methods

5.1 IERK4-A methods

For six-stage IERK methods that satisfy the canopy node condition and the two order conditions for first-order accuracy, there are 29 coefficients to be determined by 26 order conditions up to fourth-order accuracy in Table 1. That is to say, there are three independent coefficients; however, we are not able to find any fourth-order IERK methods within six stages to fulfill the assumptions of Corollary 2.1, that is, the positive semi-definiteness of the two matrices $D_{\mathrm{E}}$ and $D_{\mathrm{EI}}$ .

Unfortunately, we also fail to find energy stable fourth-order IERK methods with seven stages that fulfill all order conditions in Table 1. Alternatively, one can present some approximately fourth-order IERK (IERK4-A) methods, where the method coefficients are computed by the least squares approximation of 18 order conditions for fourth-order accuracy with a given tolerance $\varepsilon_{\text{tol}}$ , such as $\varepsilon_{\text{tol}}=10^{-6}$ . It is to note that, for the IERK4-A methods satisfying the canopy node condition and 10 order conditions up to third-order accuracy, the truncation errors are about of $O(\varepsilon_{\text{tol}}\tau^{3}+\tau^{4})$ if the solution is sufficiently regular. That is to say, a lower bound $\tau\geq O(\varepsilon_{\text{tol}})$ will be required to attain the fourth-order accuracy of IERK4-A methods although such time-step restriction is rather weak in practical applications.

The first one of fourth-order accuracy with the tolerance $\varepsilon_{\text{tol}}=10^{-6}$ , called IERK4-A1 method, has the following Butcher tableaux

		$\displaystyle\begin{array}[]{c\|c}\mathbf{c}&A\\ \hline\cr\\[-10.0pt] &\mathbf{b}^{T}\end{array}=\begin{array}[]{c\|ccccccc}0&0\\[3.0pt] \frac{95341769}{200000000}&a_{21}&\tfrac{2400249}{2000000}\\[3.0pt] \frac{292103}{800000}&a_{31}&-\tfrac{173504613}{50000000}&\tfrac{2486}{625}\\[3.0pt] \frac{59556813}{200000000}&a_{41}&\tfrac{13944041}{20000000}&a_{43}&\tfrac{1326491}{1000000}\\[3.0pt] \frac{2580667}{5000000}&a_{51}&-\tfrac{92214113}{200000000}&-\tfrac{942329}{1000000}&\tfrac{11063869}{10000000}&\tfrac{1185669331}{2000000000}\\[3.0pt] \frac{150085929}{200000000}&a_{61}&-\tfrac{2010707}{2500000}&-\tfrac{151161011}{200000000}&a_{64}&-\tfrac{23680671}{40000000}&\tfrac{15240463}{20000000}\\[3.0pt] 1&a_{71}&\tfrac{188918701}{250000000}&-\tfrac{205244563}{500000000}&-\tfrac{463474901}{250000000}&a_{75}&a_{76}&\tfrac{1339351}{2000000}\\[3.0pt] \hline\cr\\[-10.0pt] &a_{71}&\tfrac{188918701}{250000000}&-\tfrac{205244563}{500000000}&-\tfrac{463474901}{250000000}&a_{75}&a_{76}&\tfrac{1339351}{2000000}\end{array}$		(5.13)
		$\displaystyle\begin{array}[]{c\|c}\hat{\mathbf{c}}&\widehat{A}\\ \hline\cr\\[-10.0pt] &\hat{\mathbf{b}}^{T}\end{array}=\begin{array}[]{c\|ccccccc}0&0\\[3.0pt] \frac{95341769}{200000000}&\hat{a}_{21}&0\\[3.0pt] \frac{292103}{800000}&\hat{a}_{31}&\tfrac{558887}{2000000}&0\\[3.0pt] \frac{59556813}{200000000}&\hat{a}_{41}&-\tfrac{13454127}{250000000}&\tfrac{121}{1000}&0\\[3.0pt] \frac{2580667}{5000000}&\hat{a}_{51}&\hat{a}_{52}&\tfrac{32018089}{200000000}&\tfrac{593979}{2000000}&0\\[3.0pt] \frac{150085929}{200000000}&\hat{a}_{61}&\tfrac{18031311}{200000000}&\hat{a}_{63}&\tfrac{1177953}{10000000}&\tfrac{313441}{1000000}&0\\[3.0pt] 1&\hat{a}_{71}&-\tfrac{23838287}{200000000}&\tfrac{31003}{160000}&\hat{a}_{74}&\hat{a}_{75}&\tfrac{70189993}{100000000}&0\\[3.0pt] \hline\cr\\[-10.0pt] &\hat{a}_{71}&-\tfrac{23838287}{200000000}&\tfrac{31003}{160000}&\hat{a}_{74}&\hat{a}_{75}&\tfrac{70189993}{100000000}&0\end{array}$		(5.26)

where the coefficients $a_{k,1}=c_{k}-\sum_{j=2}^{k}a_{k,j},\hat{a}_{k,1}=c_{k}-\sum_{j=2}^{k-1}\hat{a}_{k,j}$ for $2\leq k\leq 7$ and

$\displaystyle a_{43}=$	$\displaystyle\,-\tfrac{1585409690050693626522959}{10^{24}},$	$\displaystyle a_{64}=\tfrac{2344693633530028154195338}{10^{24}},$
$\displaystyle a_{75}=$	$\displaystyle\,\tfrac{12541381954770160599120029}{3627843526172490000000000},$	$\displaystyle a_{76}=-\tfrac{5666098495504076766418243}{2637342719402745375000000},$
$\displaystyle\hat{a}_{52}=$	$\displaystyle\,\tfrac{33045877519515636367149}{10^{25}},$	$\displaystyle\hat{a}_{63}=\tfrac{983273174147729070884395}{10^{25}},$
$\displaystyle\hat{a}_{74}=$	$\displaystyle\,\tfrac{952796123534851512831817}{1560502861592322600000000},$	$\displaystyle\hat{a}_{75}=-\tfrac{145732093978331037774401}{338092153983867000000000}.$

The associated matrices $D_{\mathrm{E}}^{(4,1)}$ and $D_{\mathrm{EI}}^{(4,1)}$ defined by (2.10) are positive definite with the eigenvalue vectors, respectively,

		$\displaystyle\big{(}12.4798,4.72381,3.05289,1.11175,0.530827,0.0238641\big{)}^{T}$
	and	$\displaystyle\big{(}14.4357,10.3404,2.47970,1.62679,1.21228,0.000779\big{)}^{T}.$

The corresponding average dissipation rate is

\mathcal{R}^{(4,1)}\approx 3.65382+5.01594\tau\overline{\lambda}_{\mathrm{ML}}.

The second one of fourth-order accuracy with the tolerance $\varepsilon_{\text{tol}}=10^{-6}$ , called IERK4-A2 method, has the following Butcher tableaux

		$\displaystyle\begin{array}[]{c\|c}{\mathbf{c}}&{A}\\ \hline\cr\\[-10.0pt] &{\mathbf{b}}^{T}\end{array}=\begin{array}[]{c\|ccccccc}0&0\\[3.0pt] \tfrac{429533}{1000000}&a_{21}&\tfrac{63137}{200000}\\[3.0pt] \tfrac{4785663}{10000000}&a_{31}&-\tfrac{917757}{1000000}&\tfrac{100379}{100000}\\[3.0pt] \tfrac{1182276}{1000000}&a_{41}&-\tfrac{1929}{1250}&a_{43}&\tfrac{15281}{20000}\\[3.0pt] \tfrac{915703}{1000000}&a_{51}&\tfrac{98637}{1000000}&\tfrac{196933}{1000000}&a_{54}&\tfrac{531}{625}\\[3.0pt] \tfrac{7336053}{10000000}&a_{61}&\tfrac{302663}{1000000}&\tfrac{5967}{125000}&\tfrac{150781}{1000000}&-\tfrac{156231}{125000}&\tfrac{142387}{100000}\\[3.0pt] 1&a_{71}&\tfrac{1843487}{10000000}&a_{73}&-\tfrac{129513}{1000000}&-\tfrac{820173}{2000000}&a_{76}&\tfrac{88673}{125000}\\[3.0pt] \hline\cr\\[-10.0pt] &a_{71}&\tfrac{1843487}{10000000}&a_{73}&-\tfrac{129513}{1000000}&-\tfrac{820173}{2000000}&a_{76}&\tfrac{88673}{125000}\end{array}\,,$		(5.39)
		$\displaystyle\begin{array}[]{c\|c}\hat{\mathbf{c}}&\widehat{A}\\ \hline\cr\\[-10.0pt] &\hat{\mathbf{b}}^{T}\end{array}=\begin{array}[]{c\|ccccccc}0&0\\[3.0pt] \tfrac{429533}{1000000}&\hat{a}_{21}&0\\[3.0pt] \tfrac{4785663}{10000000}&\hat{a}_{31}&\hat{a}_{32}&0\\[3.0pt] \tfrac{1182276}{1000000}&\hat{a}_{41}&\tfrac{4507}{6250}&\tfrac{1025153}{2000000}&0\\[3.0pt] \tfrac{915703}{1000000}&\hat{a}_{51}&\tfrac{587731}{2000000}&\hat{a}_{53}&\tfrac{371251}{2000000}&0\\[3.0pt] \tfrac{7336053}{10000000}&\hat{a}_{61}&\tfrac{82583}{200000}&-\tfrac{1415767}{10000000}&-\tfrac{33393}{200000}&\tfrac{119457}{250000}&0\\[3.0pt] 1&\hat{a}_{71}&\tfrac{28277}{100000}&\hat{a}_{73}&\hat{a}_{74}&-\tfrac{7683}{100000}&\tfrac{99459}{250000}&0\\[3.0pt] \hline\cr\\[-10.0pt] &\hat{a}_{71}&\tfrac{28277}{100000}&\hat{a}_{73}&\hat{a}_{74}&-\tfrac{7683}{100000}&\tfrac{99459}{250000}&0\end{array}\,,$		(5.52)

where the first coefficients $a_{k,1}=c_{k}-\sum_{j=2}^{k}a_{k,j}$ and $\hat{a}_{k,1}=c_{k}-\sum_{j=2}^{k-1}\hat{a}_{k,j}$ for $2\leq k\leq 7$ , and

$\displaystyle a_{43}=$	$\displaystyle\,\tfrac{219830108841087453607347}{2\times 10^{23}},$	$\displaystyle a_{54}=-\tfrac{4498694297454501655541833}{10^{25}},$
$\displaystyle a_{73}=$	$\displaystyle\,\tfrac{2298242610563399947}{4576990146963750000},$	$\displaystyle a_{76}=-\tfrac{556251214988653}{1754043394750312500},$
$\displaystyle\hat{a}_{32}=$	$\displaystyle\,\tfrac{1017529648895428446045183}{25\times 10^{23}},$	$\displaystyle\hat{a}_{53}=\tfrac{3161854834370097094143699}{1\times 10^{25}},$
$\displaystyle\hat{a}_{73}=$	$\displaystyle\,\tfrac{4063870960730480933}{25257881055233250000},$	$\displaystyle\hat{a}_{74}=\tfrac{422441222477275261}{6239843169579000000}.$

The associated matrices $D_{\mathrm{E}}^{(4,2)}$ and $D_{\mathrm{EI}}^{(4,2)}$ defined by (2.10) are positive definite with the eigenvalue vectors, respectively,

		$\displaystyle\big{(}8.31573,4.30483,1.87790,1.55745,0.673643,0.0000139\big{)}^{T}$
	and	$\displaystyle\big{(}5.12769,3.33053,1.57402,0.952105,0.0473756,0.0000135\big{)}^{T}.$

The corresponding average dissipation rate is $\mathcal{R}^{(4,2)}\approx 2.78826+1.83862\tau\overline{\lambda}_{\mathrm{ML}}.$

By using Corollary 2.1 and Lemma 2.2, we have the following theorem.

Theorem 5.1.

In simulating the gradient flow system (2.3) with $\kappa\geq 2\ell_{g}$ , the seven-stage IERK4-1 (5.13) and IERK4-2 (5.39) methods preserve the original energy dissipation law (1.3),

\displaystyle E[U^{n,j+1}]-E[U^{n,1}]\leq\frac{1}{\tau}\Big{\langle}M_{h}^{-1}\delta_{\tau}\vec{U}_{n,j+1},D_{j}^{(4,\mu)}({\tau}M_{h}L_{\kappa})\delta_{\tau}\vec{U}_{n,j+1}\Big{\rangle}\quad\text{for $1\leq j\leq 6$,}

where $\mu=1,2$ and the differentiation matrices $D^{(4,\mu)}$ are defined by

D^{(4,\mu)}({\tau}M_{h}L_{\kappa}):=D_{\mathrm{E}}^{(4,\mu)}\otimes I+D_{\mathrm{EI}}^{(4,\mu)}\otimes(-\tau M_{h}L_{\kappa}).

The associated average dissipation rates are

\displaystyle\mathcal{R}^{(4,1)}\approx 3.65382+5.01594\tau\overline{\lambda}_{\mathrm{ML}}\quad\text{and}\quad\mathcal{R}^{(4,2)}\approx 2.78826+1.83862\tau\overline{\lambda}_{\mathrm{ML}}.

(5.53)

The seven-stage IERK4-1 (5.13) and IERK4-2 (5.39) methods are considered only to show the existence of the energy stable methods with fourth-order time accuracy. As mentioned early, Liu and Zou [26] proposed several Lobatto-type fourth-order IERK methods; but these methods have zero diagonal entries and do not satisfy our assumptions in Corollary 2.1. Actually, in the literature, we do not find any Lobatto-type or Radau-type fourth-order IERK methods that the associated differentiation matrices are positive (semi-)definite.

5.2 Tests for IERK4-A methods

Example 1 is also used to test the convergence of IERK4-A1 (5.13) and IERK4-A2 (5.39) methods by choosing the final time $T=1$ with different stabilized parameters $\kappa=2,3$ . The $L^{\infty}$ norm errors $e(\tau)$ of the two IERK4-A methods on halving time steps $\tau=2^{-k}/10$ for $0\leq k\leq 9$ are shown in Figure 12. We see that the two IERK4-A methods can achieve fourth-order accuracy and the error of IERK4-A2 method is always smaller than the error of IERK4-A1 method.

The discrete energy behaviors of the two IERK4-A methods are examined by Example 2, in which the reference solution and energy are computed by the IERK4-A1 method with a small time-step size $\tau=10^{-3}$ . Figure 13 depicts the discrete energies generated by the IERK4-A1 and IERK4-A2 methods, for three different scenes: (i) $\tau=0.01,\kappa=2$ ; (ii) $\tau=0.05,\kappa=2$ and (iii) $\tau=0.05,\kappa=3$ . As predicted by Theorem 5.1 , the discrete energies in all cases always monotonically decrease. Also, the discrete energy obtained from the IERK4-A2 method is always closer to the reference energy than that of the IERK4-A1 method. It is predictable by the numerical errors in Figure 12, which shows that the IERK4-A2 method is always more accurate than the IERK4-A1 method, and the associated average dissipation rates in (5.53), that is, $\left|\mathcal{R}^{(4,2)}-1\right|<\left|\mathcal{R}^{(4,1)}-1\right|$ .

In Figure 13 (a)-(c), we also include the discrete energy generated by the IERK3-2 method (4.37) with $a_{43}=-0.6$ , which has the average dissipation rate $\mathcal{R}^{(3,2)}=\tfrac{5}{4}+\tfrac{2}{5}\tau\overline{\lambda}_{\mathrm{ML}}$ and is regarded as the “best” IERK3 method in this paper. More surprisingly, the energy curve generated by the IERK3-2 method (4.37) with $a_{43}=-0.6$ is always closer to the reference one than the two IERK4-A schemes. Actually, this phenomenon can not be explained by the precision of numerical solutions. As seen in Figure 12, the solution errors of IERK3-2 method are larger than the errors of IERK4-A1 and IERK4-A2 methods for the time-step sizes $\tau=0.01$ and $0.05$ . Obviously, the well preservation of the original energy dissipation law does not entirely depend on the numerical precision of solution although we are unsure of the complete mechanism behind it. One possible reason is that the average dissipation rates in (5.53) of the two IERK4-A schemes are relatively large.

This raises an important issue: finding some IERK4-A method or formally fourth-order IERK method that preserves the original energy dissipation law (1.3) and perform better than the “best” third-order method (4.37) with $a_{43}=-0.6$ . We have planned to continue exploring this issue and will present the relevant results in a forthcoming report.

6 Concluding remarks

We construct some parameterized IERK methods by adopting the stiffly-accurate Lobatto-type DIRK methods in the implicit part, and present a unified theoretical framework to examine the energy dissipation properties at all stages of IERK methods up to fourth-order accurate for gradient flow problems. The novel framework contains two main parts: one is the differential form and the associated differentiation matrix of an IERK method by using the difference coefficients of method and the so-called discrete orthogonal convolution kernels; the other is the average energy dissipation rate defined via the average eigenvalue of the differentiation matrix. The rough but simple concept of average energy dissipation rate seems very useful to evaluate the overall energy dissipation behaviors of the IERK methods, including choosing proper parameters in some parameterized IERK methods or comparing different IERK methods with the same accuracy. Our main contributions include:

(i)

Among the second-order IERK2-1 (3.2), IERK2-2 (3.6) and IERK2-Radau (3.8) methods preserving the original energy dissipation law (1.3), the one-parameter IERK2-2 method (3.6) with $a_{33}=\tfrac{1+\sqrt{2}}{4}$ would be the best one from both the average energy dissipation rate $\mathcal{R}^{(2,2)}(\tfrac{\sqrt{2}}{2},\tfrac{1+\sqrt{2}}{4})=\sqrt{2}+\tfrac{\sqrt{2}}{4}\tau\overline{\lambda}_{\mathrm{ML}}$ , and the numerical experiments in Subsection 3.2.
(ii)

Among the third-order IERK3-1 (4.34), IERK3-2 (4.37) and IERK3-Radau (4.54) methods preserving the original energy dissipation law (1.3) unconditionally, the one-parameter IERK3-2 method (4.37) with the parameter $a_{43}=-0.6$ would be the best one from both the average energy dissipation rate $\mathcal{R}^{(3,2)}=\frac{5}{4}+\tfrac{2}{5}\tau\overline{\lambda}_{\mathrm{ML}}$ , and the numerical tests in Subsection 4.4.
(iii)

Two approximately fourth-order IERK methods with the order-condition tolerance $\varepsilon_{\text{tol}}=10^{-6}$ are constructed and shown to preserve the original energy dissipation law (1.3) unconditionally. From the average dissipation rates and the numerical experiments in Subsection 5.2, the IERK4-A2 method (5.39) is better than the IERK4-A1 method (5.13).
(iv)

From the perspective of maintaining the original energy dissipation law (1.3), numerical tests show that the IERK1 method (2.16) with $\theta=\tfrac{1}{2}$ performs better than some of IERK2 methods; the IERK2-2 method (3.6) with $a_{33}=\tfrac{1+\sqrt{2}}{4}$ performs better than some of IERK3 methods and the IERK3-2 method (4.37) with $a_{43}=-0.6$ performs better than the two IERK4-A methods. They suggest that the time accuracy of an IERK method cannot fully characterize the degree to which it maintains the original energy dissipation law. In practice, the selection of method parameters in IERK methods is at least as important as the selection of different IERK methods, if not more important.

At the same time, our theory is far away from complete. There are many interesting issues that we have not yet addressed. Some of them are listed as follows:

(a)

We address two approximately fourth-order IERK methods that preserve the energy dissipation law (1.3) unconditionally; however, up to now, we can not construct any formally fourth-order IERK methods that preserve the energy dissipation law (1.3).
(b)

In order to maintain the continuous dissipation rate of the original energy as much as possible, the average dissipation rate, $\mathcal{R}=\frac{1}{s_{\mathrm{I}}}\mathrm{tr}(D_{\mathrm{E}})+\frac{1}{s_{\mathrm{I}}}\mathrm{tr}(D_{\mathrm{EI}})\tau\overline{\lambda}_{\mathrm{ML}}$ , is hoped to be as close to 1 as possible. It is seen from (2.17) that the value of $\mathrm{tr}(D_{\mathrm{EI}})$ can reach zero, while we always have $\frac{1}{s_{\mathrm{I}}}\mathrm{tr}(D_{\mathrm{E}})\geq 1$ for the presented eight IERK methods, cf. Tables 2-3. We wonder what the minimum value of $\frac{1}{s_{\mathrm{I}}}\mathrm{tr}(D_{\mathrm{E}})$ may take, and how to construct the corresponding IERK method.
(c)

As we have repeatedly stated, the average dissipation rate can only be used to compare different IERK algorithms with the same accuracy. Is it possible to reasonably integrate the accuracy information of the algorithm to form a more effective indicator so that the discrete energy dissipation behaviors of IERK methods with different accuracies can be compared?
(d)

Last but not least, the global Lipschitz continuity assumption in Lemma 2.1 of the nonlinear function $g_{h}$ is limited. For specific gradient flow problems, weakening the above assumption is theoretically important, and as the closely related issue, determining the minimum stabilized parameter $\kappa$ in (2.2) is also practically useful.

References

[1] G. Akrivis, B. Li and D. Li, Energy-decaying extrapolated RK-SAV methods for the Allen-Cahn and Cahn-Hilliard equations, SIAM J. Sci. Comput., 41:6 (2019), pp. A3703-A3727.
[2] U. Ascher, S. Ruuth and R.J. Spiteri, Implicit-explicit Runge-Kutta methods for time dependent partial differential equations, Appl. Numer. Math., 25 (1997), pp. 151-167.
[3] S. Boscarino, L. Pareschi and G. Russo, A unified IMEX Runge-Kutta approach for hyperbolic systems with multiscale relaxation, SIAM J. Numer. Anal., 55:4 (2017), pp. 2085-2109.
[4] S. Boscarino and G. Russo, On a class of uniformly accurate IMEX Runge-Kutta schemes and applications to hyperbolic systems with relaxation, SIAM J. Sci. Comput., 31:3 (2009), pp. 1926-1945.
[5] A. Cardone, Z. Jackiewicz, A. Sandu and H. Zhang, Extrapolated implicit-explicit Runge-Kutta methods, Math. Model. Anal., 19:1 (2014), pp. 18-43.
[6] G. J. Cooper and A. Sayfy, Additive Runge-Kutta methods for stiff ordinary differential equations, Math. Comp., 40:161 (1983), pp. 207-218.
[7] G. Dimarco and L. Pareschi, Asymptotic preserving implicit-explicit Runge-Kutta methods for nonlinear kinetic equations, SIAM J. Numer. Anal., 51:2 (2013), pp. 1064-1087.
[8] Q. Du, L. Ju, X. Li and Z. Qiao, Maximum bound principles for a class of semilinear parabolic equations and exponential time-differencing schemes, SIAM Rev., 63 (2021), pp. 317-359.
[9] M. Fasi, S. Gaudreault, K. Lund and M. Schweitzer, Challenges in computing matrix functions, arXiv:2401.16132v1, 2024.
[10] Z. Fu and J. Yang, Energy-decreasing exponential time differencing Runge-Kutta methods for phase-field models, J. Comput. Phys., 454 (2022), 110943.
[11] Z. Fu, T. Tang and J. Yang, Energy diminishing implicit-explicit Runge-Kutta methods for gradient flows, Math. Comput., 2024, doi: 10.1090/mcom/3950.
[12] Y. Gong, J. Zhao and Q. Wang, Arbitrarily high-order unconditionally energy stable schemes for thermodynamically consistent gradient flow models, SIAM J. Sci. Comput., 42:1 (2020), pp. B135-B156.
[13] Y. Gui, C. Wang and J. Chen, IMEX-RK methods for Landau-Lifshitz equation with arbitrary damping, arXiv: 2312.15654, 2023.
[14] W. Hundsdorfer and J.G. Verwer, Numerical Solution of Time-Dependent Advection Diffusion-Reaction Equations, Springer Ser. Comput. Math. 33, Springer-Verlag, Berlin, 2003.
[15] M. Hochbruck and A. Ostermann, Exponential integrators, Acta Numerica., 19 (2010), pp. 209-286.
[16] G. Izzo and Z. Jackiewicz, Highly stable implicit-explicit Runge-Kutta methods, Appl. Numer. Math., 113 (2017), pp. 71-92.
[17] M. Jiang, Z. Zhang and J. Zhao, Improving the accuracy and consistency of the scalar auxiliary variable (SAV) method with relaxation, J. Comput. Phys., 456 (2022), 110954.
[18] S. Jin, Runge-Kutta methods for hyperbolic conservation laws with stiff relaxation terms, J. Comput. Phys., 122:1 (1995), pp. 51-67.
[19] A. Kanevsky, M.H. Carpenter, D. Gottlieb and J.S. Hesthaven, Application of implicit-explicit high order Runge-Kutta methods to discontinuous-Galerkin schemes, J. Comput. Phys., 225:2 (2007), pp. 1753-1781.
[20] C.A. Kennedy and M.H. Carpenter, Additive Runge-Kutta schemes for convection-diffusion-reaction equations, Appl. Numer. Math., 44 (2003), pp. 139-181.
[21] C.A. Kennedy and M.H. Carpenter, Diagonally implicit Runge-Kutta methods for ordinary differential equations: a review, Technical Memorandum: NASA/TM 2016-219173, 2016.
[22] Z. Li and H.-L. Liao, Stability of variable-step BDF2 and BDF3 methods, SIAM J. Numer. Anal., 60:4 (2022), pp. 2253-2272.
[23] H.-L. Liao, T. Tang and T. Zhou, Positive definiteness of real quadratic forms resulting from variable-step L1-type approximations of convolution operators, Sci. China. Math., 67:2 (2024), pp. 237-252.
[24] H.-L. Liao and X. Wang, Average energy dissipation rates of explicit exponential Runge-Kutta methods for gradient flow problems, Math. Comput., 2024, doi: 10.1090/mcom/4015.
[25] H.-L. Liao and Z. Zhang, Analysis of adaptive BDF2 scheme for diffusion equations, Math. Comput., 90 (2021), pp. 1207-1226.
[26] H. Liu and J. Zou, Some new additive Runge-Kutta methods and their applications, J. Comput. Appl. Math., 190 (2006), pp. 74-98.
[27] C. Moler and C. Van Loan, Nineteen dubious ways to compute the exponential of a matrix, twenty-five years later, SIAM Rev., 45:1 (2003), pp. 3-49.
[28] H. Olsson and G. Söderlind, The approximate Runge-Kutta computational process, BIT, 40:2 (2000), pp. 351-373.
[29] J. Shen, J. Xu and J. Yang, The scalar auxiliary variable (SAV) approach for gradient flows, J. Comput. Phys., 353 (2018), pp. 407-416.
[30] J. Shin, H.-G. Lee and J.-Y. Lee, Convex splitting Runge-Kutta methods for phase-field models, Comput. Math. Appl., 73:11 (2017), pp. 2388-2403.
[31] J. Shin, H.-G. Lee and J.-Y. Lee, Unconditionally stable methods for gradient flow using convex splitting Runge-Kutta scheme, J. Comput. Phys., 347 (2017), pp. 367-381.
[32] A. M. Stuart and A. R. Humphries, Dynamical systems and numerical analysis, Cambridge University Press, New York, 1998.
[33] X. Zhong, Additive Semi-implicit Runge-Kutta methods for computing high-speed nonequilibrium reactive flows, J. Comput. Phys., 128 (1996), pp. 19-31.