General linear methods - davidar/scholarpedia GitHub Wiki

Table of Contents Definition of general linear methods Classification of GLMs Diagonally implicit multistage integration methods Order and stage order conditions Absolute stability properties of GLMs Various approaches to the construction of DIMSIMs Examples of DIMSIMs GLMs with inherent Runge-Kutta stability Bibliography External links See also

__AUTOLINKER{attention,stability}

General linear methods (GLMs) for the numerical solution of ordinary differential equations (ODEs)

<math>\label{eq1}

\left\{ \begin{array}{ll} y'(x)=f(y(x)), & x\in [x_0,X], \\ y(x_0)=y_0, & \end{array} \right. </math>

are defined by

<math>\label{eq2}

\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. </math>

<math>n=0,1,\ldots,N\ ,</math> <math>Nh=X-x_0\ .</math> Here, <math>Y_i\ ,</math> <math>i=1,2,\ldots,s\ ,</math> are internal approximations to <math>y(x_{n-1}+hc_i)\ ,</math> <math>x_{n-1}=x_0+(n-1)h\ ,</math> <math>y(x)</math> is the solution to \eqref{eq1}, and <math>y_i^{[n]}\ ,</math> <math>i=1,2,\ldots,r\ ,</math> 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 \eqref{eq2} are determined by the coefficient matrix <math>A=[a_{ij}]</math> 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 <math>A</math> takes the following form

<math>

A= \left[ \begin{array}{ccccc} \lambda], </math> where <math>\lambda =0</math> for type <math>1</math> methods or <math>\lambda > 0</math> for type <math>2</math> methods. These methods are appropriate for non-stiff or stiff differential systems in a sequential computing environment. The type <math>3</math> and type <math>4</math> methods correspond to the matrix of the form

<math>

A=\textrm{diag}\big(\lambda,\lambda,\ldots,\lambda\big)=\lambda I, </math> where <math>\lambda =0</math> for type <math>3</math> methods or <math>\lambda > 0</math> for type <math>4</math> 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 <math>p=q=r=s\ ,</math> where <math>p</math> is the order, <math>q</math> is the stage order, <math>r</math> is the number of external stages and <math>s</math> is the number of internal stages. We will also assume that <math>U=I\ ,</math> where <math>I</math> is the identity matrix of dimension <math>r\ .</math> Moreover, we will assume that <math>V</math> is a rank one matrix of the form

<math>

V= \left[ \begin{array}{cccc} v_1], </math> with its only nonzero eigenvalue equal to one. This is equivalent to the condition <math>\sum_{i=1}^rv_i=1\ ,</math> and as the result this matrix is power bounded which ensures zero-stability of \eqref{eq2}. 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 \eqref{eq2} has order <math>p</math> and stage order <math>q=p</math> if and only if

<math>\label{eq3}

B=B_0-AB_1-VB_2+VA, </math>

where

<math>

\begin{array}{l} B_0=\left[\displaystyle\int_0^{1+c_i}\!\phi_j(x)dx/\phi_j(c_j)]_{i,j=1}^s, \end{array} </math>

<math>

\begin{array}{l} B_1=\left[\phi_j(1+c_i)/\phi_j(c_j)\right]_{i,j=1}^s, \end{array} </math>

<math>

\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} </math> and

<math>

\phi_j(x)=\displaystyle\prod_{k=1,k\neq j}^s (x-c_k). </math> The matrices <math>B_0\ ,</math> <math>B_1</math> and <math>B_2</math> can be precomputed for specific choices of the abscissa vector <math>c=[c_1,\ldots,c_s]^T</math> and the relation \eqref{eq3} plays an important role in the construction of DIMSIMs.

Applying the GLM \eqref{eq2} to the test equation <math>y'=\zeta y\ ,</math> where <math>\zeta</math> is a complex parameter, we obtain

<math>

y^{[n]}=M(z)y^{[n-1]}, </math> with <math>z=h\zeta</math> and the matrix <math>M(z)</math> defined by

<math>\label{eq4}

M(z)=V+zB(I-zA)^{-1}U. </math>

The matrix <math>M(z)</math> defined by \eqref{eq4} will be referred to as the stability matrix of \eqref{eq2} and the rational function defined by

<math>\label{eq5}

p(w,z)=\det (wI-M(z)) </math>

as the stability function of the method \eqref{eq2}. This function determines the absolute stability properties of the general linear method \eqref{eq2}. The region of absolute stability is defined as the set where <math>M(z)</math> is power-bounded, i.e. <math>\{z\in\mathbb{C}: p(w,z)=0\Rightarrow |w|\le 1\}\ .</math> 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 <math>1</math> methods, A-stability for type <math>2</math> methods) and then solving the resulting systems of nonlinear equations for the remaining coefficients of the method. For low order methods <math>(1\le p\leq 3)</math> we were able to generate and solve the resulting systems of nonlinear equations by symbolic manipulation packages such as MATHEMATICA or MAPLE. For order <math>p=4</math> 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 <math>(p\geq 5)</math> 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 <math>1</math> and type <math>2</math> methods of order <math>p=5</math> and <math>p=6</math> 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 <math>1</math> and type <math>2</math> GLMs of order <math>p=7</math> and <math>p=8\ .</math> This algorithm was further refined in Jackiewicz & Mittelmann (1999).

• Type <math>1</math> method with <math>p=q=r=s=2</math> and <math>c=[0,1]^T\ ,</math> Butcher (1993a):

<math>

A=\left[ \begin{array}{cc} 0], \quad U=\left[ \begin{array}{cc} 1], \quad B=\left[ \begin{array}{cr} \frac{5}{4}], \quad V=\left[ \begin{array}{cc} \frac{1}{2}]. \ :</math>

This method has the same region of absolute stability as an explicit two-stage Runge-Kutta method of order two.

• Type <math>2</math> method with <math>p=q=r=s=2</math> and <math>c=[0,1]^T\ ,</math> Butcher (1993a):

<math>

A=\left[ \begin{array}{cc} \frac{2-\sqrt{2}}{2}], \quad U=\left[ \begin{array}{cc} 1], \quad B=\left[ \begin{array}{cc} \frac{73-34\sqrt{2}}{28}], \quad V=\left[ \begin{array}{cc} \frac{3-\sqrt{2}}{2}]. \ :</math>

This method is L-stable. i.e. it is A-stable and <math>\lim_{z\to\infty}\lambda_{\mbox{max}}(M(z))=0\ .</math>

• Type <math>3</math> method with <math>p=q=r=s=2</math> and <math>c=[0,1]^T\ ,</math> Butcher (1993a):

<math>

A=\left[ \begin{array}{cc} 0], \quad U=\left[ \begin{array}{cc} 1], \quad B=\left[ \begin{array}{cr} -\frac{3}{8}], \quad V=\left[ \begin{array}{cc} -\frac{3}{4}]. \ :</math>

This method has the interval of absolute stability equal to <math>[-\frac{4}{3},0]\ .</math>

• Type <math>4</math> method with <math>p=q=r=s=2</math> and <math>c=[0,1]^T\ ,</math> Butcher (1993a):

<math>

A=\left[ \begin{array}{cc} \frac{3-\sqrt{3}}{2}], \quad U=\left[ \begin{array}{cc} 1], \quad B=\left[ \begin{array}{cc} \frac{18-11\sqrt{3}}{4}], \quad V=\left[ \begin{array}{cc} \frac{3-2\sqrt{3}}{2}]. \ :</math>

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 <math>p(w,z)</math> defined by \eqref{eq5} takes the form

<math>

p(w,z)=w^p\big(w-R(z)\big), </math> where <math>R(z)</math> is the stability function of a Runge-Kutta method of order <math>p\ .</math> In contrast to DIMSIMs these methods with <math>p=q</math> and <math>r=s=p+1</math> 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.

• Burrage, K. (1995) Parallel and Sequential Methods for Ordinary Differential Equations, Clarendon Press, Oxford.
• Burrage, K. and Butcher, J.C. (1980) Non-linear stability of a general class of differential equation methods, BIT 20, 185–203.
• Butcher, J.C. (1987) The Numerical Analysis of Ordinary Differential Equations. Runge-Kutta and General Linear Methods, John Wiley, New York.
• Butcher, J.C. (1993a) Diagonally-implicit multi-stage integration methods, Appl. Numer. Math. 11, 347–363.
• Butcher, J.C. (1993b) General linear methods for the parallel solution of ordinary differential equations, World Series in Applicable Analysis 2, 99–111.
• Butcher, J.C. (1994) A transformation for the analysis of DIMSIMs, BIT 34, 25–32.
• Butcher, J.C. (2003) Numerical Methods for Ordinary Differential Equations, John Wiley, Chichester.
• Butcher, J.C. (2006) General linear methods, Acta Numerica 15, 157-256.
• Butcher, J.C., Chartier, P. and Jackiewicz, Z. (1997) Nordsieck representation of DIMSIMs, Numerical Algorithms 16, 209–230.
• Butcher, J.C., Chartier, P. and Jackiewicz, Z. (1999) Experiments with a variable-order type<math>1</math> DIMSIM code, Numerical Algorithms 22, 237–261.
• Butcher, J.C. and Jackiewicz, Z. (1993) Diagonally implicit general linear methods for ordinary differential equations, BIT 33, 452–472.
• Butcher, J.C. and Jackiewicz, Z. (1996) Construction of diagonally implicit general linear methods of type 1 and 2 for ordinary differential equations, Appl. Numer. Math. 21, 385–415.
• Butcher, J.C. and Jackiewicz, Z. (1997) Implementation of diagonally implicit multistage integration methods for ordinary differential equations, SIAM J. Numer. Anal. 34, 2119–2141.
• Butcher, J.C. and Jackiewicz, Z. (1998) Construction of high order diagonally implicit multistage integration methods for ordinary differential equations, Appl. Numer. Math. 27, 1–12.
• Butcher, J.C. and Jackiewicz, Z. (2001) A reliable error estimation for DIMSIMs, BIT 41, 656–665.
• Butcher, J.C. and Jackiewicz, Z. (2002) Error estimation for Nordsieck methods, Numerical Algorithms 31, 75–85.
• Butcher, J.C. and Jackiewicz, Z. (2003) A new approach to error estimation for general linear methods, Numer. Math. 95, 487–502.
• Butcher, J.C. and Jackiewicz, Z. (2004a) Construction of general linear methods with Runge-Kutta stability properties, Numerical Algorithms 36, 53–72.
• Butcher, J.C. and Jackiewicz, Z. (2004b) Uncoditionally stable general linear methods for ordinary differential equations, BIT 44, 557–570.
• Butcher, J.C., Jackiewicz, Z. and Mittelmann, H.D. (1997) Nonlinear optimization approach to the construction of general linear methods of high order, J. Comput. Appl. Math. 81, 181–196.
• Butcher, J.C. and Podhaisky, H. (2006) On error estimation in general linear methods for stiff ODEs, Appl. Numer. Math. 56, 345–357.
• Butcher, J.C. and Wright, W.M. (2003) The construction of practical general linear methods, BIT 43, 695–721.
• Deuflhard, P., Fiedler, B. and Kunkel, P, (1987) Efficient numerical path following beyond critical points, SIAM J. Numer. Anal. 24, 912–927.
• E. Hairer, E., Nørsett, S.P. and Wanner, G. (1993) Solving Ordinary Differential Equations I. Nonstiff Problems, Springer-Verlag, Berlin, Heidelberg, New York.
• Jackiewicz, Z. (2002) Implementation of DIMSIMs for stiff differential systems, Appl. Numer. Math. 42, 251–267.
• Jackiewicz, Z. (2005) Construction and implementation of general linear methods for ordinary differential equations: a review, J. Sci. Comput. 25, 29–49.
• Jackiewicz, Z. and Mittelmann, H.D. (1999) Exploiting structure in the construction of DIMSIMs, J. Comput. Appl. Math. 107, 233–239.
• Levenberg, K. (1944) A method for the solution of certain problems in least squares, Quart. Appl. Math. 2, 164–168.
• Rheinboldt, W.C. and Burkardt, J.V. (1983a) A locally-parametrized continuation process, ACM Trans. Math. Software 9, 215–235.
• Rheinboldt, W.C. and Burkardt, J.V. (1983b) Algorithm <math>596\ :</math> a program for a locally-parametrized continuation process, ACM Trans. Math. Software 9, 236–241.
• Watson, L.T., Billups, S.C. and Morgan, A.P. (1987) HOMPACK: a suite of codes for globally convergent homotopy algorithms, ACM Trans. Math. Software 13, 281–310.
• Wright, W.M. (2002a) General linear methods with inherent Runge-Kutta stability, Ph.D. thesis, The University of Auckland, New Zealand.
• Wright, W.M. (2002b)Explicit general linear methods with inherent Runge-Kutta stability, Numer. Algorithms 31, 381–399.

Internal references

• Bill Gear (2007) Backward differentiation formulas. Scholarpedia, 2(8):3162.
• Lawrence F. Shampine and Skip Thompson (2007) Initial value problems. Scholarpedia, 2(3):2861.
• Kendall E. Atkinson (2007) Numerical analysis. Scholarpedia, 2(8):3163.
• John Butcher (2007) Runge-Kutta methods. Scholarpedia, 2(9):3147.

• Zdzislaw Jackiewicz's website

• Initial Value Problems
• Numerical Analysis
• Ordinary Differential Equations
• Runge-Kutta Methods

Category:Numerical Analysis Category:Dynamical Systems