Efficient General Linear Methods of High Order with Inherent Quadratic Stability

We search for general linear methods with s internal stages and r = s+1 external stages of order p = s + 1 and stage order q = s. We require that stability function of these methods has only two non-zero roots. This is achieved by imposing the so-called inherent quadratic stability conditions. Examples of such general linear methods which are Aand L-stable up to the order p = 8 and stage order q = p− 1 are derived.


Introduction
Consider the initial value problem for an autonomous system of ordinary differential equations of the following form y (t) = f y(t) , y(t 0 ) = y 0 , t ∈ [t 0 , T ], where the function f : R m → R m is sufficiently smooth and y 0 ∈ R m is a given initial value.
The i , i = 1, 2, . . . , s, are approximations of stage order q to y(t n−1 + c i h), and y [n] i , i = 1, 2, . . . , r, are approximations of order p to the linear combinations of the derivatives of y at the grid point t n . Different representations of GLMs are discussed in the monograph [7].
We assume that the matrices A and V have the form a s−1,1 a s−1,2 · · · λ a s,1 a s,2 · · · a s,s−1 λ so that the resulting methods have low implementation costs and are zerostable. In [6] the authors investigated the class of Nordsieck methods of order p = r = s+1 and stage order q = p or q = p−1. Since the order of these methods is p = s+1 these formulas are more efficient than diagonally implicit multistage integration methods (DIMSIMs) with parameters p = q = r = s investigated in [8,9,10,11,13,22,24], or GLM with inherent Runge-Kutta stability (IRKS) property with parameters p = q = r −1 = s−1 investigated in [12,14,22,25,26]. However, direct approach to construction of these methods presented in [6] was limited to order p = 2, 3, and 4. In [2] (compare also [3,4,15,16,17]) the approach to construct Nordsieck methods via inherent quadratic stability (IQS) was presented. It means that some algebraic conditions on coefficient matrices of the methods are imposed to guarantee that the stability function has only two non-zero roots. Therefore we can limit the linear stability analysis to the quadratic part of the stability function.
In present paper we combine both results from [6] and [2] to overcome technical difficulties and construct the families of A-and L-stable methods with parameters p = q + 1 = r = s + 1 and IQS property. The order of constructed methods is p = 2, 3, . . . , 8.
The paper is organized as follow: in Section 2 we recall the order and stage order conditions derived in [6] together with representation formulas for methods of order p = s + 1 and stage order q = p − 1 and some aspects of error propagation analysis. In Section 3 we present the concept of inherent quadratic stability, and present results obtained in [2] and [3]. In Section 4 we combine the results presented in previous two sections to construct A-and L-stable methods of order p = s + 1 and stage order q = p − 1 with inherent quadratic property and give examples of such methods of order p = 2, 3, . . . , 8. Section 5 contains some concluding remarks.

Order conditions and representation formulas
We assume (compare, for example, [22]) that the components y [n−1] i of the input vector for the step from t n−1 to t n approximate the linear combinations of the scaled derivatives h k y (k) (t n−1 ) up to the order p = s + 1, i.e., for the next step from t n to t n+1 satisfy (2.1) with n − 1 replaced by n, i.e., 3) i = 1, 2, . . . , r, with the same constants q ik . Similarly as in [8,9,22] we can write the stage order and order conditions corresponding to (2.1), (2.2), and (2.3) in a compact form using functions of a complex variable z. We define the vectors q i by q i = [q 1,i · · · q r,i ] T , i = 0, 1, . . . , r, (2.4) and the vector functions w = w(z) and w = w(z) of a complex variable z by Now, the order conditions can be written in more convenient for use in symbolic manipulation software (Mathematica, Maple) form. The GLM (1.1) has stage order q = s or q = s + 1 and order p = s + 1 if and only if where e cz = [e c1z · · · e csz ] T .
We are interested in the class of GLMs (1.1) for which the vectors q 0 , q 1 , . . . , q s defined in (2.4) take the special form q 0 = e 1 , q 1 = e 2 , . . ., q s = e s+1 , where e 1 , e 2 , . . . , e s+1 is the canonical basis in the space R s+1 . Then the vector of external approximations take the form where z(t, h) is the Nordsieck vector defined by We will also call the class of resulting GLMs the Nordsieck methods.
Following presentation in [6] we introduce the matrices The Nordsieck method has stage order q = s if and only if while the Nordsieck method has order p = s + 1 if and only if where E = exp(K), and (compare [6]). The relation (2.6) is usually inverted in order to express the matrix B in terms of V and c. We partition the coefficient matrices B and V and the matrix E into the form , and E ∈ R s×s , and define the matrix Assuming that components of vector c are distinct we write the inverted relation (2.6) as and (2.10)

Error propagation analysis
Similarly as in [6] we assume that the components of the input vector y [n−1] at the grid point t n−1 have the form . . , r, with β 1 = δ 1 = 0, and that the components of the vector of internal approximations Y [n] have the form . . , s. Than the components of the vector y [n] at the next grid point t n = t n−1 + h satisfy the relations . . , r, for the same constants β i and δ i if and the vector η is defined by The matrices B, V and vectors b T and v T were defined in (2.8), and vectors β and δ are partitioned as Since the local discretization error le(t n ) of the method (1.1) is defined by we search for an approximations of terms h p+1 y (p+1) (t n ) and h p+1 ∂f It was proved in [6] that vectors ϕ, ψ, and ϕ, ψ satisfy the systems of equations (2.14) where we use the convention that (−1)! = ∞.

Inherent Quadratic Stability
The stability properties of GLM (1.1) with respect to linear test equation [19] y = ξy, ξ ∈ C, are defined by stability function where z = ξh and stability matrix M(z) is given by Stability function p(w, z) is a rational function. To simplify analysis we multiply stability function p(w, z) by (1 − λw) s and we work with stability function denoted by the same symbol p which is a polynomial with respect to w whose coefficients are polynomials with respect to z. We recall (compare, for example [22]) that the method is said to be absolutely stable in z if p(w, z) has roots of modulus less than 1. The set of all z ∈ C such that methods is absolutely stable is called the region of absolute stability. The method is A-stable if its region of absolute stability contains a negative half plane.
The method is said to be L-stable, if it is A-stable and where ρ(M(z)) is the spectral radius of stability matrix M(z). For the discussion of algebraic stability property of general linear methods we refer to [5,18,20,21,22].
To investigate Nordsieck methods with so called quadratic stability function, i.e., stability function with only two nonzero roots, we introduce equivalence relation between matrices of the same dimensions. We say that the two matrices D and E are equivalent, denoted by D ≡ E, if they are equal except for the first two rows (compare [17,22]).
It was demonstrated in [3] that for GLM with p = q + 1 = s + 1 = r the most general matrix X satisfying conditions (3.1) is of the form 4 Search for Methods of Order p = s + 1 = r and Stage Order q = p − 1 In this section we will describe the search for Nordsieck methods of order p = s + 1 = r and stage order q = p − 1. It will be always assumed that vector c is given in advance. We begin with applying order representation formulas (2.5), (2.9) and (2.10) corresponding to order p = q = s = r − 1. Taking into account the form of matrix X given by (3.3), we solve IQS conditions (3.1) with respect to x ir , i = 3, 4, . . . , r, a ij , i = 2, 3 . . . , s, j = 1, 2, . . . , i, and v i,j , i = 2, 3, . . . , r − 1, j = i + 1, i + 2, . . . , r. We solve the remaining order condition (2.7) with respect to v 1r , q ir , i = 2, 3, . . . , r. We fix q 1r = 0, so y [n] 1 = y(t n ) + O(h s+2 ). Finally, we solve the system (compare (3.2)) p 1s = 0, p 0s = 0 (4.1) to ensure that A-stable methods are also L-stable. The remaining free parameters are used to find families of A-stable methods.
In the following subsections we will use the described above procedure to search for methods of order p = 3, 4, . . . , 8. Case p = 2 is an exception, since the stability polynomial has quadratic form without any additional assumptions. We also present sets of constants E and F , vectors β, δ, ξ and η for derived examples of methods. It was observed in [6], that systems (2.13) and (2.14) consist of s + 4 equations in 2s + 1 unknowns and have solutions for s ≥ 3. For s ≥ 4 we require that both terms in one before last equation of system (2.13) and last equation of system (2.14) have the same contributions, i.e., For s ≥ 5 we require that (4.2) is satisfied and, in addition, that the sums of squares of elements of vectors ϕ and ψ or ϕ and ψ, i.e., ϕ 2 1 + · · · + ϕ 2 s + ψ 2 1 + · · · + ψ 2 r or ϕ 2 1 + · · · + ϕ 2 s + ψ 2 1 + · · · + ψ 2 r are minimal (compare [1]).

Construction of methods with
The polynomial of the method with parameters p = r = 2 and s = q = 1 assumes quadratic form without any additional conditions. Similarly as in [6] we assume that q 2 = [0, 0] T , so y [n] = z(t n , h) + O(h 3 ). We derive coefficient matrices U and B from representation formulas, and solve the remaining order condition with respect to v 12 and c 1 , and obtain the one-parameter family of the Nordsieck methods depending on λ. The coefficients of these methods are c = 1 and We apply Schur theorem [23] (compare also [22] for illustrative examples) and obtain that these methods are A-stable if λ ≥ 1/2. There are no L-stable methods in this family. For λ = 1 2 we have For this method we have This method was derived in [6] and is presented here for the sake of completeness.

Construction of methods with
We fix in advance the abscissa vector c = [0, 1] T . After applying order and IQS conditions as described at the beginning of the section, we solve system (4.1) with respect to λ and v 1,2 . This system has two solutions λ = 1 6 (4 − √ 6), v 1,2 = 1 18 (−11 + 4 √ 6) and λ = 1 6 (4 + . We verify using Schur criterion that the method corresponding only to the second solution leads to A-stable, hence also L-stable, method. The quadratic part of the stability function is given by We note, that the method with the same value of λ parameter was constructed in [6] via quadratic stability, but these methods differ due to another approach and another choice of free parameters. The coefficient of this method are We fix in advance the abscissa vector c = [0, 1 3 , 2 3 , 1] T . After applying order and IQS conditions we solve system (4.1) with respect to v 1,3 and v 1,4 and obtain the two-parameter family of methods depending on λ and v 1,2 . The range of parameters λ and v 1,2 for which the methods are A-and Lstable is plotted in Figure 1. Again, this result was obtained by applying the We fix in advance the abscissa vector c = [0, 1 4 , 1 2 , 3 4 , 1] T . After applying order and IQS conditions we solve system (4.1) with respect to v 1,4 and v 1,5 and obtain the three-parameter family of methods depending on λ, v 1,2 and v 1,3 .
We use the Schur criterion to find A-, and hence, L-stable methods among this family. In Figure 2 we plotted the range of parameters λ and v 1,2 for some chosen values of parameter v 1,3 .
We use the Schur criterion to find A-, and hence, L-stable methods among this family. In Figure 3 we plotted the range of parameters λ and v 1,2 for some chosen values of parameters v 1,3 and v 1,4 .