General linear methods (GLMs) for the
numerical solution of ordinary differential equations (ODEs)
\[\tag{1}
\left\{
\begin{array}{ll}
y'(x)=f(y(x)), & x\in [x_0,X], \\
y(x_0)=y_0, &
\end{array}
\right.
\]
are defined by \[\tag{2} \left\{ \begin{array}{ll} Y_i=h\displaystyle\sum_{j=1}^sa_{ij}f(Y_j)+\displaystyle\sum_{j=1}^ru_{ij}y_j^{[n-1]}, & i=1,2,\ldots,s,\\ y_i^{[n]}=h\displaystyle\sum_{j=1}^sb_{ij}f(Y_j)+\displaystyle\sum_{j=1}^rv_{ij}y_j^{[n-1]}, & i=1,2,\ldots,r, \end{array} \right. \]
\(n=0,1,\ldots,N\ ,\) \(Nh=X-x_0\ .\) Here, \(Y_i\ ,\) \(i=1,2,\ldots,s\ ,\) are internal approximations to \(y(x_{n-1}+hc_i)\ ,\) \(x_{n-1}=x_0+(n-1)h\ ,\) \(y(x)\) is the solution to (1), and \(y_i^{[n]}\ ,\) \(i=1,2,\ldots,r\ ,\) are external stages. These methods were introduced by Burrage and Butcher (1980) (see also Burrage (1995), Butcher (1987, 2003), Hairer, Nørsett & Wanner (1993)). We also refer to a recent article by Butcher (2006) for an extensive review of many aspects of GLMs such as motivation for these formulas, order conditions, linear and non-linear stability, special families of methods, and order and stability barriers. GLMs include as special cases Runge-Kutta (RK) methods, linear multistep methods (LMMs), e.g. BDF methods, and predictor-corrector methods. As discussed in Butcher (2006) both RK methods LMMs have limitations and the class of GLMs offers new possibilities of constructing new formulas which attempt to combine the advantages of RK methods (large regions of stability) and LMMs (high stage order) at the same time avoiding the disadvantages of these methods (low stage order for RK formulas, small regions of stability for LMMs).
The implementation costs of (2) are determined by the coefficient matrix \(A=[a_{ij}]\) and depending on its structure GLMs can be are divided into four types which are appropriate for non-stiff or stiff differential systems in a sequential or parallel computing environment. In this review we will discuss the construction of GLMs for which the coefficient matrix \(A\) takes the following form \[ A= \left[ \begin{array}{ccccc} \lambda & 0 & 0 & \ldots & 0 \\ a_{21} & \lambda & 0 & \ldots & 0 \\ a_{31} & a_{32} & \lambda & \ldots & 0 \\ \vdots & \vdots & \ddots & \ddots & \vdots \\ a_{s1} & a_{s2} & \ldots & a_{s,s-1} & \lambda \end{array} \right], \] where \(\lambda =0\) for type \(1\) methods or \(\lambda > 0\) for type \(2\) methods. These methods are appropriate for non-stiff or stiff differential systems in a sequential computing environment. The type \(3\) and type \(4\) methods correspond to the matrix of the form \[ A=\textrm{diag}\big(\lambda,\lambda,\ldots,\lambda\big)=\lambda I, \] where \(\lambda =0\) for type \(3\) methods or \(\lambda > 0\) for type \(4\) methods. These formulas are appropriate for non-stiff or stiff differential systems in a parallel computing environment.
In the remainder of the paper we will restrict our attention mainly to the methods such that \(p=q=r=s\ ,\) where \(p\) is the order, \(q\) is the stage order, \(r\) is the number of external stages and \(s\) is the number of internal stages. We will also assume that \(U=I\ ,\) where \(I\) is the identity matrix of dimension \(r\ .\) Moreover, we will assume that \(V\) is a rank one matrix of the form \[ V= \left[ \begin{array}{cccc} v_1 & v_2 & \ldots & v_s \\ v_1 & v_2 & \ldots & v_s \\ \vdots & \vdots & \ddots & \vdots \\ v_1 & v_2 & \ldots & v_s \end{array} \right], \] with its only nonzero eigenvalue equal to one. This is equivalent to the condition \(\sum_{i=1}^rv_i=1\ ,\) and as the result this matrix is power bounded which ensures zero-stability of (2). The resulting methods are then called diagonally implicit multistage integration methods (DIMSIMs). This class of GLMs was introduced by Butcher (1993a) (see also Butcher (1993b,1994)) and further investigated by Butcher, Chartier & Jackiewicz (1997,1999) and Butcher & Jackiewicz (1993,1996,1997,1998,2001a,2001b,2002,2003, 2004a,2004b). In Butcher & Jackiewicz (1993,1996,1998) Butcher, Jackiewicz & Mittelmann (1997), and Jackiewicz & Mittelmann (1999) various approaches to the construction of such methods are described and in Butcher, Chartier & Jackiewicz (1997,1999), Butcher & Jackiewicz (1997,2001,2002,2003), Butcher & Podhaisky (2006) and Jackiewicz (2002,2005) implementation aspects are discussed such as local error estimation, changing stepsize using the Nordsieck technique, the construction of continuous interpolants, and the solution of the resulting systems of nonlinear equations for implicit methods by simplified Newton iterations. We believe that the algorithms based on these methods have great potential for practical use.
It was proved in Butcher (1993a) that the method (2) has order \(p\) and stage order \(q=p\) if and only if \[\tag{3} B=B_0-AB_1-VB_2+VA, \]
where \[ \begin{array}{l} B_0=\left[\displaystyle\int_0^{1+c_i}\!\phi_j(x)dx/\phi_j(c_j) \right]_{i,j=1}^s, \end{array} \]
\[ \begin{array}{l} B_1=\left[\phi_j(1+c_i)/\phi_j(c_j)\right]_{i,j=1}^s, \end{array} \]
\[ \begin{array}{l} B_2=\left[\displaystyle\int_0^{c_i}\!\phi_j(x)dx/\phi_j(c_j)\right]_{i,j=1}^s, \end{array} \] and \[ \phi_j(x)=\displaystyle\prod_{k=1,k\neq j}^s (x-c_k). \] The matrices \(B_0\ ,\) \(B_1\) and \(B_2\) can be precomputed for specific choices of the abscissa vector \(c=[c_1,\ldots,c_s]^T\) and the relation (3) plays an important role in the construction of DIMSIMs.
Applying the GLM (2) to the test equation \(y'=\zeta y\ ,\) where \(\zeta\) is a complex parameter, we obtain \[ y^{[n]}=M(z)y^{[n-1]}, \] with \(z=h\zeta\) and the matrix \(M(z)\) defined by \[\tag{4} M(z)=V+zB(I-zA)^{-1}U. \]
The matrix \(M(z)\) defined by (4)
will be referred to as the stability matrix of
(2) and the rational function defined by
\[\tag{5}
p(w,z)=\det (wI-M(z))
\]
as the stability function of the method (2). This function determines the absolute stability properties of the general linear method (2). The region of absolute stability is defined as the set where \(M(z)\) is power-bounded, i.e. \(\{z\in\mathbb{C}: p(w,z)=0\Rightarrow |w|\le 1\}\ .\) A general linear method is called A-stable if its region of absolute stability contains the left half of the complex plane.
The methods presented in Butcher & Jackiewicz (1993,1996,1998), Butcher, Jackiewicz & Mittelmann (1997), and Jackiewicz & Mittelmann (1999) were obtained by requiring that the GLM has some desirable stability properties (same region of absolute stability as the Runge-Kutta method of the same order for type \(1\) methods, A-stability for type \(2\) methods) and then solving the resulting systems of nonlinear equations for the remaining coefficients of the method. For low order methods \((1\le p\leq 3)\) we were able to generate and solve the resulting systems of nonlinear equations by symbolic manipulation packages such as MATHEMATICA or MAPLE. For order \(p=4\) MATHEMATICA was still able to generate the nonlinear systems in symbolic form but was no longer powerful enough to provide approximations to their solutions. We solved these systems instead with the aid of numerical algorithms based on a homotopy approach. We used continuation programs from PITCON in Rheinboldt & Burkardt (1983a,1983b), ALCON in Deuflhard, Fiedler & Kunkel (1987), and HOMEPACK in Watson, Billups & Morgan (1987), and examples of methods obtained in this way are listed in Butcher & Jackiewicz (1996). For higher orders \((p\geq 5)\) it was no longer possible to generate the corresponding systems of nonlinear equations in manageable form by symbolic manipulation packages and a different approach to the construction of such methods was needed. In Butcher & Jackiewicz (1998) we described an approach based a variant of the Fourier series method and on least squares minimization using the Levenberg-Marquardt algorithm (cf. Levenberg (1944)). Using this algorithm we obtained type \(1\) and type \(2\) methods of order \(p=5\) and \(p=6\) with desired stability properties. In Butcher, Jackiewicz & Mittelmann (1997) we used state-of-the-art optimization methods, particularly variable-model trust-region least-squares algorithms, to construct type \(1\) and type \(2\) GLMs of order \(p=7\) and \(p=8\ .\) This algorithm was further refined in Jackiewicz & Mittelmann (1999).
\[ A=\left[ \begin{array}{cc} 0 & 0 \\ 2 & 0 \end{array} \right], \quad U=\left[ \begin{array}{cc} 1 & 0 \\ 0 & 1 \end{array} \right], \quad B=\left[ \begin{array}{cr} \frac{5}{4} & \frac{1}{4} \\ \frac{3}{4} & -\frac{1}{4} \end{array} \right], \quad V=\left[ \begin{array}{cc} \frac{1}{2} & \frac{1}{2} \\ \frac{1}{2} & \frac{1}{2} \end{array} \right]. \ :\]
This method has the same region of absolute stability as an explicit two-stage Runge-Kutta method of order two.
\[ A=\left[ \begin{array}{cc} \frac{2-\sqrt{2}}{2} \\ \frac{6+2\sqrt{2}}{7} & \frac{2-\sqrt{2}}{2} \end{array} \right], \quad U=\left[ \begin{array}{cc} 1 & 0 \\ 0 & 1 \end{array} \right], \quad B=\left[ \begin{array}{cc} \frac{73-34\sqrt{2}}{28} & \frac{-5+4\sqrt{2}}{4} \\ \frac{87-48\sqrt{2}}{28} & \frac{-45+34\sqrt{2}}{28} \end{array} \right], \quad V=\left[ \begin{array}{cc} \frac{3-\sqrt{2}}{2} & \frac{-1+\sqrt{2}}{2} \\ \frac{3-\sqrt{2}}{2} & \frac{-1+\sqrt{2}}{2} \end{array} \right]. \ :\]
This method is L-stable. i.e. it is A-stable and \(\lim_{z\to\infty}\lambda_{\mbox{max}}(M(z))=0\ .\)
\[ A=\left[ \begin{array}{cc} 0 & 0 \\ 0 & 0 \end{array} \right], \quad U=\left[ \begin{array}{cc} 1 & 0 \\ 0 & 1 \end{array} \right], \quad B=\left[ \begin{array}{cr} -\frac{3}{8} & -\frac{3}{8} \\ -\frac{7}{8} & \frac{9}{8} \end{array} \right], \quad V=\left[ \begin{array}{cc} -\frac{3}{4} & \frac{7}{4} \\ -\frac{3}{4} & \frac{7}{4} \end{array} \right]. \ :\]
This method has the interval of absolute stability equal to \([-\frac{4}{3},0]\ .\)
\[ A=\left[ \begin{array}{cc} \frac{3-\sqrt{3}}{2} & 0 \\ 0 & \frac{3-\sqrt{3}}{2} \end{array} \right], \quad U=\left[ \begin{array}{cc} 1 & 0 \\ 0 & 1 \end{array} \right], \quad B=\left[ \begin{array}{cc} \frac{18-11\sqrt{3}}{4} & \frac{-12+7\sqrt{3}}{4} \\ \frac{22-13\sqrt{3}}{4} & \frac{-12+9\sqrt{3}}{4} \end{array} \right], \quad V=\left[ \begin{array}{cc} \frac{3-2\sqrt{3}}{2} & \frac{-1+2\sqrt{3}}{2} \\ \frac{3-2\sqrt{3}}{2} & \frac{-1+2\sqrt{3}}{2} \end{array} \right]. \ :\]
This method is L-stable.
A closely related class of GLMs with so called inherent Runge-Kutta stability (IRKS) was introduced recently by Butcher and Wright (2003) and Wright (2002a,2002b). These methods have the property that the stability function \(p(w,z)\) defined by (5) takes the form \[ p(w,z)=w^p\big(w-R(z)\big), \] where \(R(z)\) is the stability function of a Runge-Kutta method of order \(p\ .\) In contrast to DIMSIMs these methods with \(p=q\) and \(r=s=p+1\) can be constructed using only linear operations, Butcher and Wright (2003). The search for optimal formulas in this class and the design of algorithms based on these methods is the topic of current research.
Internal references