Construction of Eﬃcient General Linear Methods for Non-Stiﬀ Diﬀerential Systems ∗

. This paper describes the construction of explicit general linear methods in Nordsieck form with inherent quadratic stability and large areas of the stability region. After satisfying order and inherent quadratic stability conditions, the remaining free parameters are used to ﬁnd the methods with large area of region of absolute stability. Examples of methods with p = q + 1 = s = r and p = q = s = r − 1 up to order 6 are given.


Introduction
This paper is concerned with the construction of a class of explicit general linear methods (GLM) for the numerical solution of a system of ordinary differential equations (ODEs): and the four integers: the order p, the stage order q, the number of external approximations r, and the number of stages or internal approximations s. In this paper it will be always assumed that s = r = p and q = s − 1, or p = q = s and r = s + 1. Moreover, the coefficient matrices A and V are supposed to have the form This representation of V implies that the GLM (1.2) is zero-stable, i.e., that the matrix V is power bounded. When the method (1.2) is applied to the basic test equation y = ξy, t ≥ 0, ξ ∈ C, the recurrence equation z [n] = M(z)z [n−1] , n ≥ 1 arises, with z = hξ, where the stability matrix is The stability function is defined as p(w, z) = det wI − M(z) , and the method is said to be absolutely stable in z if p(w, z) has roots of modulus less than 1. The set of the points z such that the method is stable, is defined as the region of absolute stability.
GLMs have a larger number of free parameters with respect to standard numerical methods, which can be suitably chosen in order to get efficient and stable methods. In recent years, classes of GLMs with optimal properties of efficiency and stability have been constructed and analyzed in the context of DIMSIMs (see for example [6] and also [5,15]), in the context of two-step Runge-Kutta methods [1,2,12,13,14,15] and in the context of GLMs in Nordsieck form [3,4,8,9,10,11,16,17]. A successful strategy for constructing high order and stable methods consists of imposing the order and the stage order conditions and then requiring that the stability function assumes a particular expression with desired stability properties. Such approach was first introduced in [9,17] and then developed in [8,16] where authors proposed GLM in Nordsieck form with inherent Runge-Kutta stability (IRKS), which guarantees that the stability function of the method has only one non-zero root. In particular, for explicit Nordsieck methods the stability function assumes the form with R(z, η) given by where η is a free parameter. In order to improve these results, the two-step Runge-Kutta methods with quadratic stability (QS), i.e. methods whose stability function has only two non-zero roots were introduced in [13,15] and successively developed in [3,10,11]. In [10,11] explicit Nordsieck methods of type (1.2) were proposed with p = q = s − 1, r = s with QS, i.e. methods with stability polynomial of type where p r−1 (z) = 1 + p r−1,1 z + p r−1,2 z 2 + · · · + p r−1,s z s , In particular, authors used the additional free parameters, compared to the IRKS methods, in order to maximize the area of the region of absolute stability. In such a way, the resulting methods have larger stability regions compared to IRKS methods. A-and L-stable Nordsieck methods with inherent quadratic stability of the type (1.5) were investigated in [3]. The aim of our research is now to complete the investigation started in [3,10,11], by constructing explicit GLM of Nordsieck type with r = s = p, q = s − 1 and with p = q = s, r = s + 1, with quadratic stability and maximum area of the region of absolute stability. The paper is organized as follows. In Section 2 we derive representation formulae for coefficient matrices using order conditions. In Section 3 we review the concept of inherent quadratic stability and investigate the structure of the matrix X involved in this definition. In Section 4 we describe the numerical search for methods with inherent quadratic stability and extensive region of absolute stability. In Sections 5 and 6 examples of such methods with p = q + 1 = s = r and p = q = s = r − 1 are given, respectively, and comparisons with the stability regions of Runge-Kutta methods are made. Last section contains some concluding remarks.

Representation Formulae
We recall the definitions of order and stage order for the GLM (1.2). Assume that (2.1) Put q k := [q 1,k , . . . , q r,k ] T . The method (1.2) has stage order q and order p if for the same parameters q i,k . Since for the GLM (1.2) z First we derive representation formulae for the coefficient matrices in the case p = r = s, q = p − 1. We have the following result [6] (see also [15]).
Let us define the matrices: Observe that the matrix F is equal to exp(K p+1 ) with the first column removed. We also define where I p is the identity matrix of dimension p, and partition W as The following result gives useful representation formulae for the matrices U and B depending on the abscissae vector c and on the coefficient matrices A and V .

Theorem 2.
Assume that c i = c j , for any i = j and that the GLM (1.2) with r = s has order p = s and stage order q = s − 1. Then we have the following representation of the matrices U and B:

Proof. Expanding (2.2) into power series, it results
Observing that [q 0 , . . . , q p−1 ] = I p we obtain the first part of (2.5). Expanding (2.3) into power series, we get Vq 0 = q 0 , i.e. Ve 1 = e 1 , that is already satisfied by our choice of V, and also k l=0 which is equivalent to From our hypotheses C p is invertible, so the second representation formula (2.5) arises. Now we analyze the case p = q = s and r = s + 1. In order to establish representation formulae for the coefficient matrices, we partition B, V, and E p+1 as follows: where b T stands for the first row of B and 0 stands for vector or matrix of appropriate dimension, and obtain the following theorem: [see [3]] Assume that c i = c j , for any i = j and that the GLM (1.2) with r = s + 1 has order and stage order p = q = s. Then we have this representation of the matrices U and B:

Inherent Quadratic Stability
This concept, which was first introduced in [13,15] and successively developed in [3,10,11], means that there exists a matrix X ∈ R r×r such that Here, the relation P ≡ Q means that the matrices P and Q are identical with the exception of the first two rows. In [3,10] the following theorem was proved, asserting that IQS condition leads to quadratic stability. The following theorems analyze the structure of the matrix X appearing in the conditions (3.1), for the GLMs we consider in this paper.

Theorem 5.
For a GLM of type (1.2) with p = r = s and q = s − 1, the most general matrix X satisfying conditions (3.1) is of the form: Using the IQS conditions (3.1) we have Now, following the lines of the proof of Theorem 2 of [10] it results It results w(z) = φ(z) + q s z s , with φ(z) = [1 z . . . z s−1 ] T , and so: where we moved the O(z s+1 ) terms in the right-hand side of the equivalence. Since Equating to zero the coefficients of z 0 , z 1 , . . . , z s−2 , it follows that X − J must be zero except for the first two rows and the last column. Then, equating to zero the coefficient of z s−1 , it results x i,s = q i,s , i = 3, . . . , s. Thus the proof is completed.
In the case p = q = s and r = s + 1, the matrix X has the same structure as in the case r = s, p = q = s − 1, studied in [10]: Theorem 6. For a GLM of type (1.2) with p = q = s and r = s + 1, the most general matrix X satisfying conditions (3.1) is of the form: Proof. The proof follows along the lines of Theorem 2 of [10].

Search for Methods with IQS and Maximum Area of Absolute Stability Region
The search for GLMs (1.2) with IQS and the maximum area of the stability region consists of two main steps: • impose the order conditions and the IQS conditions; • find the method with the maximum area of the stability region, within the class of GLMs derived in the first step. In the case p = q = s, r = s + 1 we also fix abscissa vector c in advance. It follows from representation formulae (2.8) that matrices U and B can be expressed in terms of A, V, and c. As in the previous case we need to find in the matrix X only the elements [x 3,r , . . . , x r,r ] (see Theorem 6). We solve the system of equations arising from IQS conditions with respect to [x 3,r , . . . , x r,r ], all elements of matrix A and all elements of matrix V except the first row, i.e., v 1,2 · · · v 1,r . Therefore we obtain an s-parameters family of methods of order p = q = s with IQS depending on the vector m = [v 1,2 . . . v 1,r ].
The second step of the search is the same for both cases. The free parameters m of each family of GLMs are chosen by maximizing area(m), that is the area of the intersection of the absolute stability region with the negative half plane. This numerical search is carried our by solving the minimization problem −area(m) → min, by using the Matlab program fminsearch. The function area(m) was computed using integration in polar coordinates as explained in [12]. Examples of methods are given in Sections 5 and 6.
In the numerical search for GLMs with IQS and maximum area of the stability region, we need to compute the coefficients of the stability polynomial, and this operation is required a huge number of times. For the case p = r = s, q = s − 1, we adopted the Fourier series approach proposed in the context of DIMSIMs in [7] for the representation formula of polynomial coefficients. Moreover, in order to considerably reduce the time of computation, we implemented such formula by using the fast Fourier transform algorithm, as proposed in [10]. On the other hand, for the case p = q = s, r = s + 1 symbolic manipulation furnish simple and efficient representation formulae for the coefficients of the polynomial, thus we could avoid these two special techniques.

Examples of Methods with p = s = r and q = s − 1
In this section we illustrate some examples of methods constructed by using the algorithm of numerical search described in the previous section. We will compare the stability regions of our methods with the regions of stability of the Runge-Kutta methods of the same order. The coefficients matrices will be given in Matlab rational format.

p
It is evident that in this case all methods possess the quadratic stability property. The GLM with the largest stability area has area equal to 9.0129, and error constant equal to 0.0833. The stability polynomial is In Fig. 1 we have plotted the stability region of this method and, for comparison, the stability region of explicit Runge-Kutta methods of order 2. In Fig. 2 we have plotted the stability region of this method and, for comparison, the stability region of explicit Runge-Kutta methods of order 3.

r
GLM method with IQS and the largest stability area has area equal to 18.3609, and error constant equal to −0.0071. The stability polynomial is  In Fig. 3 we have plotted the stability region of this method and, for comparison, the stability region of explicit Runge-Kutta methods of order 4. It is interesting to notice that the best quadratic polynomial, with respect to the In Fig. 4 we have plotted the stability region of this method and the stability region of an approximation of order 5 of the function exp(z).
6 Examples of Methods with s = p = q and r = s + 1 In this section we present results of search for GLMs of type (1.2) with p = q = s and r = s + 1, with IQS and maximum area of the stability region, up to order p = 6. To be more concise, we give in Matlab rational format the vector c and the coefficient matrices A and V, since matrices U and B can be computed with the representation formulae (2.8). We observe that these results cannot be improved by relaxing the IQS conditions, i.e. requiring simply that the method possesses the QS property, (as done for another class of GLMs in [10]). As a matter of fact, the IQS conditions give raise to a family of GLMs depending on s free parameters, used to maximize the area of the stability region. On the other hand a quadratic polynomial which satisfies the same conditions as the stability polynomial of a GLM of type (1.2) depends on s free parameters, as well. Therefore it is not surprising that the QS approach do not produce GLMs with larger regions of absolute stability.
For sake of clarity, we illustrate in detail how we compute matrices A and V (except for v 1,2 , . . . , v 1,r ). First we fix the vector c as composed by s equally spaced points in [0, 1]. Then we impose the order conditions (2.8) and the IQS conditions (3.1) with X of the form given in Theorem 6. In this way, It turns out that the matrix A has a very special structure with all non-zero elements equal to 1 p−1 , where p is the order of the method. The precise explanation of this phenomenon will be the subject of future work. The method coefficients are c = 1, In Fig. 5 we have plotted the stability region of this method and, for comparison, the stability region of explicit Runge-Kutta methods of order 1. In Fig. 6 we have plotted the stability region of this method and, for comparison, the stability region of explicit Runge-Kutta methods of order 2.
In Fig. 7 we have plotted the stability region of this method and, for comparison, the stability region of explicit Runge-Kutta methods of order 3.
In Fig. 9 we have plotted the stability region of this method and of an approximation of order 5 of the function exp(z). In Fig. 10 we have plotted the stability region of this method and of an approximation of order 6 of the function exp(z).

Concluding Remarks
We have presented the construction of explicit GLMs in Nordsieck form with p = r = s, q = p − 1 and with p = q = s, r = s + 1, with inherent quadratic stability. We derived the representation formulae form the coefficients matrices from the order conditions. We analyzed the structure of the matrix X related to the definition of IQS. Finally we carried out the construction of GLMs with IQS and maximum area of the stability region by solving suitable minimization problems. Examples of methods have been provided up to order 6. Comparisons with Runge-Kutta methods prove that the proposed methods have stability regions considerably larger than Runge-Kutta methods of the same order.
Future work will include the implementation of these methods and comparison with classical methods such as DIMSIMs and GLMs with IRKS. Moreover, we intend to further analyze the structure of the coefficient matrix A of GLMs in Nordsieck form with IQS, in the case p = q = s, r = s + 1.