Some Notes on Summation by Parts Time Integration Methods

TL;DR

This paper investigates the stability and properties of summation by parts (SBP) time integration methods, revealing their connections to Runge-Kutta schemes and clarifying the conditions under which they possess desirable stability features.

Contribution

It provides new proofs and counterexamples to clarify the assumptions, relations, and stability properties of SBP schemes and their connection to Runge-Kutta methods.

Findings

A necessary technical assumption for SBP schemes is identified.

Not all Runge-Kutta schemes with SBP stability are derived from SBP schemes.

Classical collocation methods on Radau and Lobatto nodes are SBP schemes.

Abstract

Some properties of numerical time integration methods using summation by parts operators and simultaneous approximation terms are studied. These schemes can be interpreted as implicit Runge-Kutta methods with desirable stability properties such as $A$-, $B$-, $L$-, and algebraic stability. Here, insights into the necessity of certain assumptions, relations to known Runge-Kutta methods, and stability properties are provided by new proofs and counterexamples. In particular, it is proved that a) a technical assumption is necessary since it is not fulfilled by every SBP scheme, b) not every Runge-Kutta scheme having the stability properties of SBP schemes is given in this way, c) the classical collocation methods on Radau and Lobatto nodes are SBP schemes, and d) nearly no SBP scheme is strong stability preserving.

Equations

$$\forall t\in (0,T):u^{\prime }(t)=f(t,u(t)),u(0)=u_0,$$

$$MD+D^{T}M=t_R{t}_R^T-t_L{t}_L^T.$$

$$Du=f+\sigma M^{-1}t_{L}\bigl{(}u_{0}-t_{L}^{T}u\bigr{)}.$$

$$Du=f+M^{-1}t_{L}\bigl{(}u_{0}-t_{L}^{T}u\bigr{)}.$$

$$A=\frac{1}{T}(D+M^{-1}{t}_L{t}_L^T{)}^{-1}=\frac{1}{T}(MD+t_L{t}_L^T{)}^{-1}M,b=\frac{1}{T}M1,c=\frac{1}{T}({\tau }_1,\dots ,{\tau }_s{)}^T.$$

$$R(z)=1+z{b}^T(I-zA{)}^{-1}1=\frac{\det (I-zA+z1b^T)}{\det (I-zA)}$$

$$\forall t,u,v:⟨f(t,u)-f(t,v),u-v⟩\le \nu {\Vert u-v\Vert }^2,$$

$$D=\left(\begin{array}{cccc}-2& 1& 1& 0\\ -1& 0& 0& 1\\ -1& 0& 0& 1\\ 0& -1& -1& 2\end{array}\right),M=\frac{1}{4}\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right),t_L=\left(\begin{array}{c}1\\ 0\\ 0\\ 0\end{array}\right),t_R=\left(\begin{array}{c}0\\ 0\\ 0\\ 1\end{array}\right),$$

$$A=\frac{1}{48}\left(\begin{array}{cccc}27& -33-6\sqrt{6}& -3& 9+6\sqrt{6}\\ -7+2\sqrt{6}& 33& -9-2\sqrt{6}& -1\\ 7& 3+2\sqrt{6}& 33& -11-2\sqrt{6}\\ 21-6\sqrt{6}& 21& -21+6\sqrt{6}& 27\end{array}\right),b=\frac{1}{8}\left(\begin{array}{c}1\\ 3\\ 3\\ 1\end{array}\right),c=\frac{1}{3}\left(\begin{array}{c}0\\ 1\\ 2\\ 3\end{array}\right).$$

$$R(z)=\frac{\det (I-zA+z1b^T)}{\det (I-zA)}=\frac{12-18z+3z^2+3z^3}{12-30z+27z^2-11z^3+2z^4},$$

$$t_R=A^{-T}b=\frac{1}{16}{\left(-4+\sqrt{6},-\sqrt{6},12-\sqrt{6},8+\sqrt{6}\right)}^T.$$

$$A^{-1}1=(D+M^{-1}{t}_L{t}_L^T)1=M^{-1}{t}_L⟺t_L=M{A}^{-1}1.$$

$$t_L^{T}1=1,t_L^{T}c=0.$$

$$D=A^{-1}-M^{-1}{t}_L{t}_L^T=A^{-1}-A^{-1}1t_L^T.$$

$$C(\eta )$$

$$D(\zeta )$$

$$A^{T}M{c}^{q-1}=\frac{1}{q}M(1-c^q)⟺qM{c}^{q-1}=A^{-T}M(1-c^q),$$

$$qM{c}^{q-1}=(D^{T}M+t_L{t}_L^T)(1-c^q).$$

$$\forall v:q{v}^{T}M{c}^{q-1}=v^T(D^{T}M+t_L{t}_L^T)(1-c^q),$$

$$v^T{t}_L{t}_L^T(1-c^q)=v(0)(1-0^q)=v(0).$$

$$q{v}^{T}M{c}^{q-1}=q{\int }_0^{1}v(t)t^{q-1}dt$$

$$v^T{D}^{T}M(1-c^q)={\int }_0^1{v}^{\prime }(t)(1-t^q)dt=-v(0)+q{\int }_0^{1}v(t)t^{q-1}dt,$$

$$A{c}^{q-1}=\frac{1}{q}{c}^q⟺qM{c}^{q-1}=M{A}^{-1}{c}^q,$$

$$\forall v:q{v}^{T}M{c}^{q-1}=v^T(MD+t_L{t}_L^T)c^q,$$

$$\forall v:q{v}^{T}M{c}^{q-1}=v^T(-D^{T}M+t_R{t}_R^T)c^q.$$

$$v^T{t}_R{t}_R^T{c}^q=v(1)1^q=v(1).$$

$$q{v}^{T}M{c}^{q-1}=q{\int }_0^{1}v(t)t^{q-1}dt$$

$$-v^T{D}^{T}M{c}^q=-{\int }_0^1{v}^{\prime }(t)t^{q}dt=-v(1)+q{\int }_0^{1}v(t)t^{q-1}dt,$$

$$(D+M^{-1}{t}_L{t}_L^T)1=0+M^{-1}{t}_L=b_1^{-1}{t}_L.$$

$$\sum _{j=1}^s{a}_{i,j}(c_i-c_j)=\frac{1}{2}{c}_i^2,i\in \left\{1,\dots ,s\right\}.$$

$$\sum _{j=1}^s{a}_{i,j}=c_i,i\in \left\{1,\dots ,s\right\},\sum _{i=1}^s{b}_i{a}_{i,j}=b_j(1-c_j),j\in \left\{1,\dots ,s\right\}.$$

Full text

1 Known Results on SBP SAT Schemes

In order to solve an ordinary differential equation (ODE)

[TABLE]

a grid $0\leq\tau_{1}<\dots<\tau_{s}\leq T$ is introduced and the numerical solution is approximated pointwise as $u_{i}=u(\tau_{i})$ and $f_{i}=f(\tau_{i},u_{i})$. Summation by parts (SBP) operators can be defined as follows, cf. [10, 3, 2].

Definition 1.1**.**

An SBP operator of order $p\in\mathbb{N}$ on $[0,T]$ consists of

• •

a discrete operator $D$ approximating the derivative $Du\approx u^{\prime}$ with order of accuracy $p$,

• •

a symmetric and positive definite discrete mass/norm matrix $M$ approximating the $L^{2}$ scalar product $u^{T}Mv\approx\int_{0}^{T}u(t)v(t)\operatorname{d\!}t$,

• •

and interpolation vectors $t_{L},t_{R}$ approximating the interpolation to the boundary as $t_{L}^{T}u\approx u(0)$, $t_{R}^{T}u\approx u(T)$ with order of accuracy at least $p$, such that

[TABLE]

SBP operators mimic integration by parts discretely via the summation by parts property (2). An SBP time discretisation using a simultaneous approximation term (SAT) of (1) with parameter $\sigma\in\mathbb{R}$ is [8, 7, 1]

[TABLE]

Most stability results have been achieved for the choice $\sigma=1$, i.e.

[TABLE]

Hence, this discretisation will be considered in the following. The numerical solution at $t=T$ is given by $t_{R}^{T}u$, where $u$ solves (4). The interval $[0,T]$ can also be partitioned into multiple subintervals/blocks such that multiple steps of this procedure are used sequentially.

In order to guarantee that (4) can be solved for a dissipative linear scalar problem, the following assumption is introduced [8].

Assumption 1.2**.**

For $\sigma>\frac{1}{2}$, all eigenvalues of $D+\sigma M^{-1}t_{L}t_{L}^{T}$ have strictly positive real part.

The following characterisation of (4) as Runge-Kutta method has been developed in [1].

Theorem 1.3**.**

If assumption 1.2 is satisfied, (4) is equivalent to an implicit Runge-Kutta method with the following Butcher coefficients, where $1$ denotes also the vector $(1,\dots,1)^{T}\in\mathbb{R}^{s}$.

[TABLE]

The factor $\frac{1}{T}$ is needed since Runge-Kutta coefficients are normalised to the interval $[0,1]$.

In order to make this article sufficiently self-contained, some classical stability properties of Runge-Kutta methods will be recalled briefly, cf. [5, sections IV.3 and IV.12]. The absolute value of solutions of the scalar linear ODE $u^{\prime}(t)=\lambda u(t)$, $u(0)=u_{0}\in\mathbb{C}$, $\lambda\in\mathbb{C}$, cannot increase if $\operatorname{Re}\lambda\leq 0$. The numerical solution after one time step of a Runge-Kutta method with Butcher coefficients $A,b,c$ is $u_{+}=R(\lambda\,\Delta t)u_{0}$, where

[TABLE]

is the stability function of the Runge-Kutta method. The stability property is mimicked discretely as $\mathinner{\!\left\lvert u_{+}\right\rvert}\leq\mathinner{\!\left\lvert u_{0}\right\rvert}$ if $\mathinner{\!\left\lvert R(\lambda\,\Delta t)\right\rvert}\leq 1$.

Definition 1.4**.**

A Runge-Kutta method with stability function $R(z)$ is $A$-stable, if $\mathinner{\!\left\lvert R(z)\right\rvert}\leq 1$ for all $z\in\mathbb{C}$ with $\operatorname{Re}(z)\leq 0$. The method is $L$-stable, if it is $A$-stable and $\lim_{z\to\infty}R(z)=0$.

Hence, $A$-stable methods are stable for every time step $\Delta t>0$ and $L$-stable methods damp out stiff components corresponding to $\lambda=-x$ with large $x\in\mathbb{R}$ sufficiently fast.

Another classical stability property is connected with possibly nonlinear problems (1) in Hilbert spaces satisfying a one-sided Lipschitz condition

[TABLE]

where $\nu\in\mathbb{R}$ is the one-sided Lipschitz constant of $f$. This condition gives some bounds on the growth rate of the difference between two solutions. In particular, the distance between two solutions cannot increase if $\nu\leq 0$.

Definition 1.5**.**

A Runge-Kutta method is $B$-stable, if the contractivity condition (7) with $\nu\leq 0$ implies $\mathinner{\!\left\lVert u_{+}-v_{+}\right\rVert}\leq\mathinner{\!\left\lVert u_{0}-v_{0}\right\rVert}$ for all $\Delta t>0$.

The following stability properties have been obtained in [7, 1].

Theorem 1.6**.**

Suppose that assumption 1.2 holds. Then, the SBP SAT scheme (4) is $A$- and $L$-stable. If the mass matrix $M$ is diagonal, the scheme is also $B$-stable.

2 Assumptions and Algebraic Stability

In this section, the new results of this short note concerning the necessity of assumption 1.2 and the necessity of an SBP SAT form for stability properties guaranteed by Theorem 1.6 are presented.

2.1 Assumption on Eigenvalues of $D+\sigma M^{-1}t_{L}t_{L}^{T}$

Assumption 1.2 has been proved for classical second order SBP operators in [8] and for SBP operators on Gauss, Radau, and Lobatto quadrature nodes in [9]. It has been examined numerically for other classical finite difference SBP operators in [8]. Since assumption 1.2 holds for all known SBP SAT schemes investigated in [8, 7, 1, 9], it is interesting to know whether it follows from properties of SBP operators.

Theorem 2.1**.**

There are SBP operators that do not satisfy assumption 1.2.

Proof.

Consider the operators

[TABLE]

on the uniform grid with four nodes $0,\frac{1}{3},\frac{2}{3},1$ in $[0,1]$. The SBP property (2) is satisfied, $t_{L}$ and $t_{R}$ are exact, and $D$ is a first order accurate SBP derivative operator. However, $(D+\sigma M^{-1}t_{L}t_{L}^{T})u=0$ for $u=(0,-1,1,0)^{T}$. Thus, zero is an eigenvalue of $D+\sigma M^{-1}t_{L}t_{L}^{T}$ for all $\sigma\in\mathbb{R}$. ∎

2.2 Algebraic Stability

Many stability properties such as $A$- and $B$-stability are satisfied if the following algebraic criterion is fulfilled by the coefficients of a Runge-Kutta method [5, Theorem 12.4].

Definition 2.2**.**

A Runge-Kutta method with Butcher coefficients $A,b,c$ is algebraically stable, if $\forall i\colon b_{i}\geq 0$ and the matrix $\operatorname{diag}(b)A+A^{T}\operatorname{diag}(b)-bb^{T}$ is positive semidefinite.

It has been noted in [1] that an SBP SAT scheme (4) with diagonal $M$ is algebraically stable, since the nodes $\tau_{i}$ are pairwise distinct, i.e. the corresponding Runge-Kutta method is nonconfluent. In that case, $B$- and algebraic stability are equivalent [5, Corollary 12.14]. This can also be proved directly, cf. [1, Theorem 5.8].

It is interesting to know whether all Runge-Kutta methods with stability properties guaranteed by Theorem 1.6 can be constructed as SBP SAT schemes. Since those schemes are $L$-stable, the classical Gauss collocation schemes (which are not $L$-stable) cannot be constructed in this way, cf. [1]. However, there is

Theorem 2.3**.**

Consider a Runge-Kutta method and the statements

1. i)

The Runge-Kutta method is $A$-, $L$-, $B$-, and algebraically stable with pairwise distinct nodes $c_{i}\in[0,1]$, only positive quadrature weights $b_{i}$, and invertible matrix $A$. 2. ii)

The Runge-Kutta method is given via Theorem 1.3 by SBP SAT schemes (4) with at least first order accurate operators satisfying assumption 1.2.

Theorem 1.6 and the preceeding discussion show that ii) and “$M$ is diagonal” imply i). However, i) does not imply ii).

Proof.

The following example has been constructed using the $W$-transformation [5, Sections IV.5, IV.13, and IV.14]. Consider the Runge-Kutta method with coefficients

[TABLE]

Then, the algebraic stability matrix $\operatorname{diag}(b)A+A^{T}\operatorname{diag}(b)-bb^{T}$ has the eigenvalues $\frac{5}{8}$, $\frac{3}{8}$, and zero (twofold). Hence, the Runge-Kutta method is algebraically stable (because $b_{i}>0$ is satisfied additionally) and therefore also $A$- and $B$-stable. Its stability functions

[TABLE]

fulfils $\lim_{z\to\infty}R(z)=0$. Thus, the scheme is also $L$-stable.

It suffices to consider $T=1$. If the scheme is given by an SBP SAT method (4) via Theorem 1.3, $b=M1$ and $A=(D+M^{-1}t_{L}t_{L}^{T})^{-1}$. The SBP property (2) yields $A^{-1}=-M^{-1}D^{T}M+M^{-1}t_{R}t_{R}^{T}$. Because of consistency, $D1=0$ and $t_{R}^{T}1=1$. Hence, $1^{T}MA^{-1}=t_{R}^{T}$. Inserting $M1=b$ results in

[TABLE]

Similarly, consistency of $D$ and $t_{L}$ implies

[TABLE]

$t_{R}$ defined by (11) is first order accurate, i.e. $t_{R}^{T}1=1$ and $t_{R}^{T}c=1$. The same accuracy of $t_{L}$ requires

[TABLE]

Because of (12), $D$ can be written as

[TABLE]

Since $M\in\mathbb{R}^{4\times 4}$ should be symmetric, it is determined by ten real parameters, e.g. $M_{11}$, $M_{12}$, $M_{13}$, $M_{14}$, $M_{22}$, $M_{23}$, $M_{24}$, $M_{33}$, $M_{34}$, $M_{44}$. $t_{R}$ is given explicitly by (11), $t_{L}$ depends linearly on $M$ via (12), and $D$ is given via an affine-linear function of $M$ in (14).

The accuracy conditions (13) are linear in $t_{L}$ and hence linear in $M$. They can be used to eliminate two parameters, e.g. $M_{11}$ and $M_{12}$. Then, the SBP property (2) is a system of 16 equations that are quadratic in the parameters $M_{ij}$. This system can be solved uniquely, which has been verified using the function Reduce of Mathematica [11]. For this unique solution, one eigenvalue of $M$ is zero. Thus, $M$ is not positive definite, in contradiction to the assumptions. ∎

3 Classical Collocation Methods

In [1], it has been shown that the SBP SAT scheme with Lobatto quadrature on four nodes corresponds to the classical Lobatto IIIC method with $s=4$. It has been mentioned that this is similar for the Radau IA and Radau IIA schemes. However, to the authors knowledge, no general proof of this result has been given up to now. To prove it, the classical conditions

[TABLE]

will be used.

Theorem 3.1**.**

The SBP SAT scheme (4) using left Radau, right Radau, or Lobatto quadrature correspond to the classical Radau IA, Radau IIA, or Lobatto IIIC Runge-Kutta methods for all orders of accuracy.

Proof.

It suffices to consider the case $T=1$, i.e. the time interval $[0,1]$.

The weights and nodes of the left Radau quadrature (left endpoint [math] included) are the weights $b_{i}$ and nodes $c_{i}$ of the Radau IA method. The matrix $A$ of the Radau IA method is determined uniquely by the condition $D(s)$, i.e. $D(\zeta)$ with $\zeta=s$ in (16) [5, section IV.5]. Hence, it suffices to prove that the SBP SAT method satisfies $D(s)$, which can be written using $M=\operatorname{diag}(b)$ as

[TABLE]

where the exponentiation $c^{q}$ is performed pointwise. Inserting $A$ from (5) yields

[TABLE]

This is equivalent to

[TABLE]

where $v$ is any polynomial of degree $\leq s-1$, evaluated at the nodes $c_{i}$. Since the left endpoint [math] is included,

[TABLE]

The Radau quadrature is exact for polynomials of degree $\leq 2s-2$. Hence, for every $q\in\left\{1,\dots,s\right\}$,

[TABLE]

and (using integration by parts)

[TABLE]

proving $D(s)$.

The weights and nodes of the right Radau quadrature (right endpoint $1$ included) are the weights $b_{i}$ and nodes $c_{i}$ of the Radau IIA method. The matrix $A$ of the Radau IIA method is determined uniquely by the condition $C(s)$, i.e. $C(\eta)$ with $\eta=s$ in (15) [5, section IV.5]. Hence, it suffices to prove that the SBP SAT method satisfies $C(s)$, which can be written using $M=\operatorname{diag}(b)$ as

[TABLE]

where the exponentiation $c^{q}$ is again performed pointwise. Inserting $A$ from (5), this is equivalent to

[TABLE]

where $v$ is any polynomial of degree $\leq s-1$, evaluated at the nodes $c_{i}$. Using the SBP property (2), this can be rewritten as

[TABLE]

Since the right endpoint $1$ is included,

[TABLE]

Using the exactness of the Radau quadrature for polynomials of degree $\leq 2s-2$, for every $q\in\left\{1,\dots,s\right\}$,

[TABLE]

and (using integration by parts)

[TABLE]

proving $C(s)$.

Finally, the weights and nodes of the Lobatto quadrature (left and right endpoints $0,1$ included) are the weights $b_{i}$ and nodes $c_{i}$ of the Lobatto IIIC method. The matrix $A$ of the Lobatto IIIC method is determined uniquely by the condition $C(s-1)$ and $a_{i,1}=b_{1},i\in\left\{1,\dots,s\right\}$ [5, section IV.5]. Since the order of accuracy of the SBP operator is $s-1$, $C(s-1)$ is satisfied [1, Lemma 5.3]. This can also be proved using similar manipulations as above. Hence, it remains to show $a_{i,1}=b_{1},i\in\left\{1,\dots,s\right\}$. Since $D$ is exact for constants, $t_{L}=(1,0,\dots,0)^{T}$, and $M=\operatorname{diag}(b_{1},\dots,b_{s})$,

[TABLE]

Therefore, $(a_{i,1})_{i=1}^{s}=At_{L}=(D+M^{-1}t_{L}t_{L}^{T})^{-1}t_{L}=b_{1}1,$ proving $a_{i,1}=b_{1},i\in\left\{1,\dots,s\right\}$. ∎

4 Strong Stability Preservation

Another desirable stability property of time integration methods is that they are strong stability preserving (SSP), i.e. that they preserve convex stability properties of the explicit Euler method [4].

Definition 4.1**.**

A numerical time integration method is called strongly stable for a given convex functional $\eta$ if $\eta(u_{+})\leq\eta(u_{0})$, possibly using some time step restriction of the form $0<\Delta t\leq\Delta t_{\mathrm{max}}$.

A numerical time integration method is called strong stability preserving with SSP coefficient $c>0$, if $\eta(u_{+})\leq\eta(u_{0})$ for all time steps $0<\Delta t\leq c\,\Delta t_{E}$ whenever the explicit Euler method is strongly stable for the convex functional $\eta$ and time steps $0<\Delta t\leq\Delta t_{E}$.

Typical convex functionals $\eta$ considered for SSP methods are the norm in a Hilbert space for dissipative operators or the total variation seminorm for semidiscretisations of scalar conservation laws.

Theorem 4.2**.**

No SBP SAT scheme (4) whose SBP operator has a diagonal norm matrix, satisfies assumption 1.2, and

1. a)

is either at least second order accurate 2. b)

or is at least first order accurate and contains at least one of the end points $0,1$ in the nodes $c_{i}$

can be strong stability preserving.

Proof.

An SSP scheme must satisfy $\forall i,j\in\left\{1,\dots,s\right\}\colon a_{i,j}\geq 0$ [4, Observation 5.2].

If the SBP operator is at least second order accurate, the corresponding Runge-Kutta method satisfies $C(2)$ [1, Lemma 5.3], i.e. $\sum_{j=1}^{s}a_{i,j}=c_{i}$ and $\sum_{j=1}^{s}a_{i,j}c_{j}=\frac{1}{2}c_{i}^{2}$ for $i\in\left\{1,\dots,s\right\}$. Subtracting the second equation from the first one multiplied by $c_{i}$ yields

[TABLE]

If $a_{i,j}$ were non-negative, the left hand side would be non-positive for $i=1$ (since $c_{j}\geq c_{1}$) and thus zero. Hence, the first row of $A$ would be zero, which is impossible, because $A$ is invertible.

If the SBP operator is at least first order accurate, the corresponding Runge-Kutta method satisfies $C(1)$ and $D(1)$ [1, Lemma 5.3 and Lemma 5.4], i.e.

[TABLE]

If the left endpoint $0=c_{1}$ is contained in the nodes, non-negativity of all $a_{i,j}$ and $C(1)$ imply $\forall j\in\left\{1,\dots,s\right\}\colon a_{1,j}=0$. Similarly, if the right endpoint $c_{s}=1$ is contained in the nodes, non-negativity of all $a_{i,j}$ and $D(1)$ imply $\forall i\in\left\{1,\dots,s\right\}\colon a_{i,s}=0$. But $A$ cannot have a zero row or column because it is invertible. ∎

Remark 4.3**.**

Classical finite difference SBP operators and those based on Radau or Lobatto quadrature include at least one endpoint and can thus not result in SSP schemes. The SBP SAT scheme (4) on two Gauss nodes does not contain an endpoint and has a first order accurate derivative operator. Nevertheless, the scheme is not SSP, since the corresponding matrix $A$ has a negative entry.

Example 4.4**.**

There is a first order accurate SBP operator with diagonal norm matrix not including any boundary node such that the resulting Runge-Kutta method given by Theorem 1.3 is SSP. Indeed, choose $T=1$ and

[TABLE]

The operators $D,t_{L},t_{R}$ are exact for polynomials of degree one, assumption 1.2 has been verified numerically for $\sigma\in(1/2,2)$, $A$ and $b$ have only non-negative entries, and the scheme is strong stability preserving with SSP coefficient $\approx 1.35$, computed using NodePy [6].

Acknowledgements

This work was supported by the German Research Foundation (DFG, Deutsche Forschungsgemeinschaft) under Grant SO 363/14-1. We acknowledge support by the German Research Foundation and the Open Access Publication Funds of the Technische Universität Braunschweig. The author would like to thank the anonymous reviewers for their helpful comments and valuable suggestions to improve this article.