# Construction of Nordsieck Second Derivative General Linear Methods with Inherent Quadratic Stability.

1 Introduction

In the past 50 years, great efforts for obtaining the numerical solution of stiff ordinary differential equations (ODEs) have been done. For stiff ODEs good accuracy and some reasonably wide region of absolute stability, in the best situation A-stability, are required. Among the linear multistep methods, backward differentiation formulae (BDF)  are regularly used for the numerical solution of stiff initial value problems. Many modifications of BDF were introduced to improve the stability characteristics of BDF such as EBDF (extended BDF) , MEBDF (modified EBDF) , MF-MEBDF (Matrix free MEBDF) , A-BDF method  and A-EBDF method . On the other hand, implicit Runge-Kutta (IRK) methods have a special role in the numerical solution of stiff problems. Although there exist A-stable Runge-Kutta methods of arbitrarily high order, but the implementation cost of full implicit Runge-Kutta methods is high. So, many subclasses of IRK schemes, such as diagonally implicit Runge-Kutta (DIRK) and singly implicit Runge-Kutta (SDIRK) schemes have been developed to reduce implementation cost. The main restrictions of these methods are their relatively low order and suffering from the order reduction phenomenon when applied to stiff ODEs.

General linear methods (GLMs) as a comprehensive extension of these traditional methods were introduced by Butcher in 1966 . GLMs opened up the possibility of obtaining essentially new methods which are neither these two categories nor slight variations of these methods (see more [8,9,29]).

To construct methods with high orders and satisfactory stability properties, methods were introduced where use second derivatives of the solution [16,18, 21,27]. Although GLMs include linear multistep methods, Runge-Kutta and many other standard methods, but they don't cover second derivative methods. So, GLMs were extended to second derivative general linear methods (SGLMs) by Butcher and Hojjati in  and studied more by Abdi and Hojjati [2,3,4].

The general form of SGLM for the numerical solution of initial value problem of the form

y'(x)= f (y(x)), f : [R.sup.m] [right arrow] [R.sup.m], x [member of] [[x.sub.0], [bar.x]], y([x.sub.0]) = [y.sub.0],

where the function f is assumed to be sufficiently smooth and m is the dimensionality of the system, is

[mathematical expression not reproducible], (1.1)

where n =1, 2, ... , N, Nh = [bar.x] - [x.sub.0], h is the stepsize, [cross product] is the Kronecker product of two matrices and [I.sub.m] stands for the identity matrix of dimension m. Here, A, [bar.A] [member of] [R.sup.sxs], U [member of] [R.sup.sxr], B, [bar.B] [member of] [R.sup.rxs], and V [member of] [R.sup.rxr] are six coefficients matrices of the SGLM. Also, p and q are respectively order and stage order of the method, r is the number of input and output approximations, and s is the number of internal stages. The vector [Y.sup.[n]] = [[[Y.sup.[n].sub.i]].sup.s.sub.i=1] is an approximation of stage order q to the vector y([x.sup.n-1] + ch) = [[y([x.sup.n-1] + [c.sub.i]h)].sup.s.sub.i=1] with c = [[[c.sub.1] [c.sub.2] ... [c.sub.s]].sup.T] as the abscissa vector, and the vectors f([Y.sup.[n]]) = [[f([Y.sup.[n]i])].sup.s.sub.i=1] and g([Y.sup.[n]]) = [[g([Y.sup.[n].sub.i])].sup.s.sub.i=1] denote the first and second derivative stage values, where g(*) = f'(*)f(*). The vectors [y.sup.[n-1]] = [[[y.sup.[n-1].sub.i]].sup.r.sub.i=1] and [y.sup.[n]] = [[[y.sup.[n].sub.i]].sup.r.sub.i=1] stand as the input and output vectors at step number n, respectively. Construction and the main features of SGLMs including pre-consistency, consistency, stability and types of these methods have been discussed in [3,10]. To construct efficient SGLMs for stiff ODEs with lower implementation cost, we will always assume that the coefficient matrices A and [bar.A] have the form

[mathematical expression not reproducible]

Also, to ensure that the method (1.1) is zero-stable , here we will assume that the coefficient matrix V has the form

[mathematical expression not reproducible]. (1.2)

The stability behavior of these methods is considered using the standard linear test problem y' = [xi]y , where [xi] is a complex parameter with negative real part. Applying method (1.1) to this test problem, the stability matrix is achieved by

M(z) = V + (zB + [z.sup.2][bar.B])[(I - zA - [z.sup.2][bar.A]).sup.-1]U,

where z = [xi]h. The characteristic polynomial of M(z), as the stability function

[p.sup.*](w, z) = det(wI - M(z))

is a polynomial of degree r with respect to w and its coefficients are rational functions with respect to z. To investigate stability properties of the methods corresponding to this function it is usually more convenient to work with the polynomial obtained by multiplying stability function by its denominator, i.e.,

p(w, z) = [(1 - [lambda]z - [mu][z.sup.2]).sup.s] det(wI - M(z)).

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.

If the stability function has only one nonzero root, then the method is said to possess Runge-Kutta stability (RKS) property. In  Abdi and Hojjati introduced a subclass of SGLMs as SDIMSIMs (second derivative diagonally implicit multistage integration methods) and constructed methods of this subclass with RKS property. They obtained order barriers for some types of SGLMs which have RKS property in [3, 4] and constructed SDIMSIMs of all types up to order 4 in . Also, the construction of Nordsieck SGLMs up to order four and their implementation in a variable stepsize environment have been investigated in . It is desirable that SGLMs have RKS property, but it is complicated task because it needs to solve large polynomial equations systems of high degree, especially for the methods with a large r and s parameters. In , a special case of SGLMs as A-[bar.A]-V methods was introduced which made easier to reach RKS property. Similar to what has been done in [8, 11, 12, 13, 33, 34] for GLMs with inherent Runge-Kutta stability (IRKS), SGLMs with inherent Runge-Kutta stability (SIRKS) were introduced in  in which the coefficients matrices satisfy some relationships to ensure the methods have RKS property.

In , SGLMs of high order with s = 2 were investigated which have quadratic stability (QS) property. In this paper, we are going to relax the concept of SIRKS to the concept of inherent quadratic stability (IQS). This property was first introduced in  for two-step Runge-Kutta (TSRK) methods and then presented for GLMs [6,14] which guarantees the stability function has only two nonzero roots. Using this approach, we solve fewer equations in comparison with methods based on SIRKS, which makes construction to be easier and gains some additional free parameters. We use these parameters in order to construct methods with 'nice' coefficient, L-stability property and reasonably small error constant.

The rest of the paper is organized as follows. In Section 2, the order and stage order conditions for the Nordsieck SGLMs are represented and some sufficient conditions are obtained which guarantee quadratic stability property. In Section 3, practical construction of the methods with IQS property is discussed. Finally, some numerical experiments of the constructed methods are given in Section 4 to verify the theoretical results.

2 SGLMs with IQS property

In this section, we first recall the order and stage order conditions for the Nordsieck SGLMs. Then, we present interrelations between the matrices which ensure that SGLMs have QS property.

2.1 Order and stage order conditions

To formulate order and stage order conditions SGLMs in Nordsieck form in the case p = q = s = r - 1, we assume that the components of the input vector, for the next step, satisfy

[y.sup.[n-1].sub.i] = [h.sup.i-1][y.sup.(i-1)] ([x.sub.n-1])+ O([h.sup.p+1]), i = 1, 2, ..., r.

The method (1.1) in the Nordsieck form has stage order q and order p if the components of the internal stages and the output vector satisfy

[mathematical expression not reproducible].

Order and stage order conditions are derived in the following theorem.

Theorem 1.  An SGLM in Nordsieck form has order p and stage order q = P if and only if

exp(cz) = zA exp(cz) + [z.sup.2][bar.A] exp(cz) + UZ + O([z.sup.p+1]), (2.1)

exp(z)Z = zB exp(cz) + [z.sup.2][bar.B] exp(cz) + VZ + O([z.sup.p+1]), (2.2)

where Z is the p + 1 dimensional vector with component number i equal to [z.sub.i-1] and exp(cz) is the component-by-component exponential of cz.

An equivalent condition for order p and stage order q = p is that U and V are related to A, A, B and B by

U = Cp - ACp Kp - [bar.A][C.sub.p] [K.sup.2.sub.p], (2.3)

V = [E.sub.p] - B[C.sub.p] [K.sub.p] - [bar.B][C.sub.p] [K.sup.2.sub.p], (2.4)

where [C.sub.p] is the scaled Vandermonde matrix [mathematical expression not reproducible] is the shifting matrix defined by [K.sub.p] = [0 [e.sub.1] ... [e.sub.p]] with [e.sub.j] as the jth unit vector and

[mathematical expression not reproducible]

The relation (2.4) can be inverted in order to express the matrix B in terms of V, [bar.B] and c. Let us partition the matrix [C.sub.p][K.sub.p] and [C.sub.p][K.sup.2.sub.p] as follows

[mathematical expression not reproducible].

Now, assuming that the components of vector c are distinct, we can write the relation (2.4) as

[mathematical expression not reproducible]. (2.5)

2.2 Criteria for IQS

After satisfying the order and stage order conditions, relations (2.3) and (2.5), the stability function of (1.1) takes the form

p(w, z) = [(1 - [lambda]z - [mu][z.sup.2]).sup.s][w.sup.r] - [p.sub.r-1](z)[w.sup.r-1] + ... + [(-1).sup.r][p.sub.0](z),

where [p.sup.i](z), i = 0,1, ..., r - 1 are the polynomials of degree 2s with coefficients that depend on elements of A, [bar.A], [bar.B], V and c. We will try to construct the SGLMs with the stability function to be

p(w,z) = [w.sub.r-2][((1 - [lambda]z - [mu][z.sup.2]).sup.s][w.sup.2] - [[??].sub.1](z)w + [[??].sub.0](z)) (2.6)

with a root w = 0 of multiplicity r - 2 and a quadratic polynomial. To do this, a direct idea is to compute the stability function and solve system of nonlinear equations [p.sub.j](z) = 0, j = 0, 1, ..., r -3. But, in general, it is a very complicated task, especially for methods with large number of stages s and r. However, it is possible to find interrelations between the coefficients matrices which will guarantee this property.

To investigate the Nordsieck methods with IQS, we first introduce equivalence relation between matrices of the same dimensions. We say that two matrices D and E are equivalent, denoted by "D [equivalent to] E", if they are equal except for the first two rows.

Definition 1. A Nordsieck SGLM with p = q = s = r - 1 and coefficients matrices A, [bar.A], B, [bar.B], U, and V defined by (1.2) has IQS property if there exists a matrix X [member of] [R.sup.(s+1)x(s+1)] such that

B[bar.A] [equivalent to] X[bar.B], BA + [bar.B] [equivalent to] XB, (2.7)

BU [equivalent to] XV - VX. (2.8)

We have the following result.

Theorem 2. Assume that the SGLM (1.1) with p = q = s = r - 1 has IQS property and that the matrices [I.sub.s] - zA - [z.sup.2][bar.A] and [I.sub.s+1] - zX are nonsingular. Then its stability function takes the form (2.6).

Proof. The relations (2.7) are equivalent to

B(I - zA - [z.sup.2][bar.A]) [equivalent to] (I - zX)(B + z[bar.B]),

which by non-singularity of the matrix I - z A - [z.sup.2][bar.A], it follows that

B [equivalent to] (I - zX)(B + z[bar.B])[(I - z A - [z.sup.2][bar.A]).sup.-1]. (2.9)

For investigating characteristic polynomial of the stability matrix M(z), we consider the matrix related to M(z) by similarity transformation. Using (2.8) and (2.9), and assuming that I - zX is nonsingular, it follows that

[mathematical expression not reproducible],

or

[mathematical expression not reproducible],

where [M.sub.1,1](z) is a 2 x 2 matrix and [M.sub.1,2](z) is a 2 x (s - 1) matrix. Thus, the matrix M(z) has only two non-zero eigenvalues. This completes the proof.

The next theorem investigates the structure of the matrix X appearing in the IQS conditions.

Theorem 3. For an SGLM with p = q = s = r - 1, the general form of matrix X satisfying IQS conditions is

[mathematical expression not reproducible].

Proof. Multiply both sides of the relation (2.1) by zB and (2.2) by I - zX, and by adding the resulting equations, we get

[mathematical expression not reproducible].

From the IQS conditions, it follows that

(exp(z)I - V)(I - zX)Z [equivalent to] O([z.sup.p+1]). (2.10)

Let us partition the matrix V and the vector D(z) := (I - zX)Z as follows

[mathematical expression not reproducible],

where

[mathematical expression not reproducible],

where 0 stands for the zero matrix of dimension (r - 2) x 2, [D.sub.1](z) contains the first two elements of D(z) and [D.sub.2](z) contains the rest of its elements. From (2.10) it follows that

[mathematical expression not reproducible].

Since the matrix exp(z)[I.sub.r-2] - [V.sub.2,2] is invertible, thus we can consider the matrix [bar.V] defined by

[mathematical expression not reproducible].

We have

[mathematical expression not reproducible]

and it follows that

D(z) = (I - zX)Z [equivalent to] O([z.sup.p+1]).

Since zJZ [equivalent to] Z, where J is transpose of matrix [K.sub.p], we have

z(J - X)Z [equivalent to] O([z.sup.p+1]),

which is equivalent to the relation (J - X)Z [equivalent to] O([z.sup.p]). Therefore X - J must be zero except for the first two rows and the last column. This concludes the proof.

3 Construction of IQS methods

In this section, we are going to construct L-stable SGLMs with IQS property. After applying the order conditions, relations (2.3) and (2.5), and applying IQS conditions, a number of parameters remain as the free parameters. We choose the remaining free parameters to construct L-stable methods (if the methods are A-stable) with nice coefficients and reasonably small error constant.

To control the error constant, we use the characteristic of the stability function p(w, z), i.e.,

[mathematical expression not reproducible],

where [C.sub.p+1] is the error constant of the method of order p. Assuming Astability, we search for the related subset of L-stable methods by solving the nonlinear system

[mathematical expression not reproducible].

The boundary locus method  enables us to make a search among the remaining free parameters corresponding to A-stable (and so L-stable) methods.

3.1 Methods with p=s=1

In this subsection, we construct Nordsieck SGLMs with p = q = s = r - 1 = 1 and IQS property. We set abscissa c = 1, [C.sub.2] = -[10.sup.-5] as the error constant of the method and use relations (2.3) and (2.5) for satisfying order condition. After applying order and stage order, error constant and L-stability conditions, the stability function takes the form

p(w, z) = (1 - [lambda]z - [mu][z.sup.2])[w.sup.2] - [[??].sub.1](z)w + [[??].sub.0](z),

where

[[??].sub.1](z) 150001/100000 z - 2z[lambda] - z[mu] +1, [[??].sub.0](z) = 50001/100000 z - z[mu] - z[lambda].

Now, we search among the free parameters that guarantee the L-stability property of the method. Using the boundary locus method, we find pairs of ([lambda], [mu]) values in domain [0, 2] x [-2,0] giving L-stability of method in Figure 1.

We select a single example characterized by [lambda] = 3/4 and [mu] = -1/5. The coefficients of the corresponding method are

[mathematical expression not reproducible].

3.2 Methods with p=s=2

In this subsection, we construct Nordsieck SGLMs with p = q = s = r - 1 = 2 and IQS property. We set the abscissa vector c = [[1/2 1].sup.T] and the error constant [C.sub.3] = -[10.sup.-5]. Applying the order and stage order, IQS, error constant and L-stability conditions, we obtain a 4-parameter family of the methods depending on [v.sub.12], [v.sub.23], [lambda] and [mu].

Setting [v.sub.12] = [v.sub.23] = -0.01, we search for the free parameters [lambda] and [mu] to find L-stable methods. We have plotted pairs of (A, values in domain [0, 2] x [-2, 0] giving L-stability in Figure 2. An example of such a method corresponding to [lambda] = 3/5 and [mu] = -1/5 is

[mathematical expression not reproducible].

3.3 Methods with p=s=3

In this subsection, we construct Nordsieck L-stable SGLMs with p = q = s = r - 1 = 3, the abscissa vector c = [[[1/2] [3/4[ 1].sup.T], the error constant [C.sub.4] = -[10.sup.-5] and IQS property. We obtain a 9-parameter family of methods which by setting [a.sub.21] = [b.sub.12] = [b.sub.13] = [V.sub.13] = 0, [[bar.a].sub.21] = - 0.001, [[bar.a].sub.32] = - 0.01 and [v.sub.12] = 1/2, we search for the parameters [lambda] and [mu] to find L-stable methods. For such methods, the pairs of ([lambda], [mu]) values in domain [0, 2] x [ 2, 0] giving L-stability are shown in Figure 3.

The coefficients matrices of a single example of these methods are as follows

[mathematical expression not reproducible].

3.4 Methods with p=s=4

In this subsection, we construct Nordsieck L-stable SGLMs with p = q = s = r - 1 = 4, the abscissa vector c = [[[1/4] [1/2] [3/4] 1].sup.T], the error constant [C.sub.5] = -[10.sup.-5] and IQS property. We obtain a 13-parameter family of methods which by setting [a.sub.41] = [a.sub.42] = [[bar.a].sub.31] = ([[bar.a].sub.32] = ([[bar.a].sub.41] = [[bar.a].sub.43] = 0, [[bar.b].sub.12] = -0.003, [[bar.b].sub.13] = 0.002, [[bar.b].sub.14] = -0.01, [v.sub.12] = -0.001 and [v.sub.13] = -0.0031, we search for the parameters [lambda] and [mu] to find L-stable methods. For such methods, the pairs of ([lambda], [mu]) values in domain [0, 2] x [-2, 0] giving L-stable methods are shown in Figure 4.

We select a single example characterized by ([lambda], [mu]) = (3/5, - 9/50) which the coefficients of this method are

[mathematical expression not reproducible].

4 Numerical experiments

In this section, we present some results of numerical experiments to show efficiency of the constructed methods in Section 3.

At first, in order to validate the order of proposed methods, we consider the following problem:

Problem 1. The non-linear stiff system of ODEs

[mathematical expression not reproducible],

whose exact solution is [[[y.sub.1](x), [y.sub.2](x)].sup.T] = [[exp(-4x), exp(- x)].sup.T] and x [member of] [0,1].

We apply the methods of orders p =1, 2, 3, 4 with the fixed stepsizes h = [1/2.sup.k] for different integer values of k on Problem 1. In Table 1 and Table 2, we have listed the Euclidean norms of errors [parallel][e.sub.h]([bar.x])[parallel] at the endpoint of the integration and the observed order of convergence p computed by

p = log([parallel][e.sub.h]([bar.x])[parallel]/[parallel][e.sub.h/2]([bar.x])[parallel])/log(2),

where [parallel][e.sub.h]([bar.x])[parallel] and [parallel][e.sub.h/2]([bar.x])[parallel] are the errors corresponding to stepsizes h and h/2.

Also in Figure 5 we have plotted log([parallel][e.sub.h]([bar.x])[parallel]) versus log(h) and linear functions in terms of log(h) with a slope of p for methods of orders p = 1, 2, 3, 4. These numerical results confirm the theoretical order of convergence. Now in a similar way to the introduced strategy in , we are going to verify the capability and efficiency of the constructed method of order 4 in a variable stepsize environment.

The Nordsieck methods (1.1) in the variable stepsize mode takes the form

[mathematical expression not reproducible], (4.1)

for n = 1, 2, ..., N, where [h.sub.n] = [x.sub.n] - [x.sub.n-1]. Here [Y.sup.[n]] is an approximation of stage order q = p to the vector y([x.sub.n-1] + [ch.sub.n]) = [[y([x.sub.n-1] + [c.sub.i][h.sup.n)]].sup.s.sub.i=1], [y.sup.[n]] is an approximation of order p to the Nordsieck vector [[[h.sup.i-1.sub.n][y.sup.(i-1)]([x.sub.n])].sup.r.sub.i=1], and D([[delta].sub.n]) is the diagonal rescaling matrix defined by

D([[delta].sub.n]) := (1, [[delta].sub.n], [[delta].sup.2.sub.n], ..., [[delta].sup.p.sub.n]), [[delta].sub.n] = [h.sub.n]/ [h.sub.n-1].

The zero-stability properties of the methods determine by the eigenvalues of the matrix VD([[delta].sub.n]). Since the matrix VD([[delta].sub.n]) for all [[delta].sub.n], has a simple eigenvalue equal to 1 and eigenvalue zero of multiplicity p, it follows that the methods (4.1) are zero-stable for any variable stepsize pattern.

To determine the accuracy of the calculations and to achieve a suitable choice of stepsize for the next step, we must approximate the local truncation error. For this purpose, we approximate the principal local truncation error. So, est([x.sub.n]) [approximately equal to] [C.sub.p+1][h.sub.p+1][y.sub.(p+1])([x.sub.n]), can be calculated as an approximation to the local truncation error.

For the method of order 4, we use the linear combination of the form

est([x.sub.n]) = [C.sub.p+1]([s.summation over (i=1)] [[gamma].sub.pi][h.sup.2]g([Y.sup.[n].sub.i])),

where [[gamma].sub.41] = - 64, [[gamma].sub.42] = 192, [[gamma].sub.43] = - 192, [[gamma].sub.44] = 64.

To control the stepsize in the advancing from the step n to the step n +1, we use the following control

est([x.sub.n]) [less than or equal to] Rtol * max{[parallel][y.sub.n][parallel], [parallel][y.sub.n+1][parallel]} + Atol, (4.2)

where Atol and Rtol are given the absolute and relative tolerances. If the control (4.2) is not satisfied, the current step is repeated with the halved stepsize. Otherwise, the current step is accepted and the new stepsize is chosen according to the standard step control strategy as the following

[mathematical expression not reproducible].

In our numerical experiments, we have used Atol = Rtol = tol, [rho] = 0.95 and [DELTA] = 2.

We consider the following test problems:

Problem 2. (Chemical Akzo Nobel problem) The non-linear stiff system of ODEs 

[mathematical expression not reproducible],

with y(0) = [[0.437, 0.00123,0,0,0, 0.367].sup.T] and x [member of] [0, 180], where the [r.sub.i] and [F.sub.in] are auxiliary variables, given by

[mathematical expression not reproducible].

The values of the parameters [k.sub.1], [k.sub.2], [k.sub.3], [k.sub.4], K, klA, p(C[O.sub.2]) and H are

[mathematical expression not reproducible].

Problem 3. (Hires problem) The chemical reaction involving eight reactants was proposed 

[mathematical expression not reproducible],

with the initial values y(0) = [[1, 0, 0, 0, 0, 0, 0, 0.0057].sup.T], x [member of] [0, 321.8122].

Using the mentioned strategies, we have implemented the SGLM with IQS property of order 4, constructed in Subsection 3.4, in variable stepsize environment for solving Problems 2 and 3. To compare, we also present the results of numerical experiments of the L-stable Nordsieck GLM of order p = 4 given in  with the abscissa vector c = [[[1/4] [1/2] [3/4] 1].sup.T]. In Table 3 and Table 4, we have reported ns as the number of steps, nrs as the number of rejected steps, nJe as the number of Jacobian evaluations and ge as the global error at the end of the interval of integration for different given tolerances, to/.

The numerical results confirm the capability of the proposed methods to solve stiff problems.

5 Conclusions

Construction of SGLMs with quadratic stability property is desirable. To achieve this aim, we need to solve a large number of nonlinear equations. To remedy this situation, SGLMs with inherent quadratic stability (IQS) property is interesting. In this paper, we introduced SGLMs with IQS property, by finding some sufficient conditions on the coefficients matrices of the methods. Some L-stable IQS methods up to order 4 were constructed and confirmed by some numerical results. There are some lines of development that can be followed. One of them can be the construction of high-order implicit SGLMs with IQS property. It is our plan for future work.

https://doi.org/10.3846/13926292.2017.1269024

References

 A. Abdi. Construction of high-order quadratically stable second-derivative general linear methods for the numerical integration of stiff ODEs. J. Comput. Appl. Math., 303:218-228, 2016. https://doi.org/10.1016/jxam.2016.02.054.

 A. Abdi, M. Bras and G. Hojjati. On the construction of second derivative diagonally implicit multistage integration methods for ODEs. Appl. Numer. Math., 76:1-18, 2014. https://doi.org/10.1016/j.apnum.2013.08.006.

 A. Abdi and G. Hojjati. An extension of general linear methods. Numer. Algor., 57(2):149-167, 2011. https://doi.org/10.1007/sll075-010-9420-y.

 A. Abdi and G. Hojjati. Maximal order for second derivative general linear methods with Runge-Kutta stability. Appl. Numer. Math., 61(10):1046-1058, 2011. https://doi.org/10.1016Zj.apnum.2011.06.004.

 A. Abdi and G. Hojjati. Implementation of Nordsieck second derivative methods for stiff ODEs. Appl. Numer. Math., 94:241-253, 2015. https://doi.org/10.1016/j.apnum.2015.04.002.

 M. Bras. Nordsieck methods with inherent quadratic stability. Math. Model. Anal., 16(1):82-96, 2011. https://doi.org/10.3846/13926292.2011.560617.

 J.C. Butcher. On the convergence of numerical solutions to ordinary differential equations. Math. Comp., 20(93):1-10, 1966. https://doi.org/10.1090/S0025-5718-1966-0189251-X.

 J.C. Butcher. General linear methods for stiff differential equations. BIT, 41(2):240-264, 2001. https://doi.org/10.1023/A:1021986222073.

 J.C. Butcher. Numerical methods for ordinary differential equations. John Wiley & Sons, 2008.

 J.C. Butcher and G. Hojjati. Second derivative methods with RK stability. Numer. Algor., 40(4):415-429, 2005. https://doi.org/10.1007/s11075-005-0413-1.

 J.C. Butcher and Z. Jackiewicz. Construction of general linear methods with Runge-Kutta stability properties. Numer. Algor., 36(1):53-72, 2004. https://doi.org/10.1023/B:NUMA.0000027738.54515.50.

 J.C. Butcher and W.M. Wright. The construction of practical general linear methods. BIT, 43(4):695-721, 2003. https://doi.org/10.1023/B:BITN.0000009952.71388.23.

 J.C. Butcher and W.M. Wright. A transformation relating explicit and diagonally implicit general linear methods. Appl. Numer. Math., 44(3):313-327, 2003. https://doi.org/10.1016/S0168-9274(02)00149-6.

 A. Cardone and Z. Jackiewicz. Explicit Nordsieck methods with quadratic stability. Numer. Algor., 60(1):1-25, 2012. https://doi.org/10.1007/s11075-011-9509-y.

 J.R. Cash. On the integration of stiff systems of ODEs using extended backward differentiation formulae. Numer. Math., 34(3):235-246, 1980. https://doi.org/10.1007/BF01396701.

 J.R. Cash. Second derivative extended backward differentiation formulas for the numerical integration of stiff systems. SIAM J. Numer. Anal., 18(1):21-36, 1981. https://doi.org/10.1137/0718003.

 J.R. Cash. The integration of stiff initial value problems in ODEs using modified extended backward differentiation formulae. Comut. Math. Appl., 9(5):645-657, 1983. https://doi.org/10.1016/0898-1221(83)90122-0.

 R.P.K. Chan and A.Y.J. Tsai. On explicit two-derivative Runge-Kutta methods. Numer. Algor., 53(2):171-194, 2010. https://doi.org/10.1007/s11075-009-9349-1.

 D. Conte, R. D'Ambrosio and Z. Jackiewicz. Two-step Runge-Kutta methods with quadratic stability functions. J. Sci. Comput., 44(2):191-218, 2010. https://doi.org/10.1007/s10091-010-937-x.

 G. Dahlquist. A special stability problem for linear multistep methods. BIT, 3(1):27-43, 1963. https://doi.org/10.1007/BF01963532.

 W.H. Enright. Second derivative multistep methods for stiff ordinary differential equations. SIAM J. Numer. Anal., 11(2):321-331, 1974. https://doi.org/10.1137/0711029.

 A.K. Ezzeddine, G. Hojjati and A. Abdi. Sequential second derivative general linear methods for stiff systems. Bull. Iranian Math. Soc., 40(1):83-100, 2014.

 C. Fredebeul. A-BDF: a generalization of the backward differentiation formulae. SIAM J. Numer. Anal., 35(5):1917-1938, 1998. https://doi.org/10.1137/S0036142996306217.

 C. Gear. Simultaneous numerical solution of differential-algebraic equations. IEEE transactions on circuit theory, 18(1):89-95, 1971. https://doi.org/10.1109/TCT.1971.1083221.

 E. Hairer and G. Wanner. Solving Ordinary Differential Equations II: Stiff and differential-algebraic problems. Springer, Berlin, 2010.

 G. Hojjati, M.Y. Rahimi Ardabili and S.M. Hosseini. A-EBDF: an adaptive method for numerical solution of stiff systems of ODEs. Math. Comput. Simul., 66(1):33-41, 2004. https://doi.org/10.1016/j.matcom.2004.02.019.

 G. Hojjati, M.Y. Rahimi Ardabili and S.M. Hosseini. New second derivative multistep methods for stiff systems. Appl. Math. Model., 30(5):466-476, 2006. https://doi.org/10.1016/j.apm.2005.06.007.

 S.M. Hosseini and G. Hojjati. Matrix free MEBDF method for the solution of stiff systems of ODEs. Math. Comput. Model., 29(4):67-77, 1999. https://doi.org/10.1016/S0895-7177(99)00040-0.

 Z. Jackiewicz. General Linear Methods for Ordinary Differential Equations. John Wiley & Sons, 2009.

 J.D. Lambert. Computational methods in ordinary differential equations. Wiley, 1973.

 W.M. Lioen and J.J.B. de Swart. Test Set for Initial Value Problem Solvers. Centrum voor Wiskunde en Informatica, 1998.

 A. Movahedinejad, G. Hojjati and A. Abdi. Second derivative general linear methods with inherent Runge-Kutta stability. Numer. Algor., 73(2):371-389, 2016. https://doi.org/10.1007/s11075-016-0099-6.

 W.M. Wright. Explicit general linear methods with inherent Runge-Kutta stability. Numer. Algor., 31(1):381-399, 2002. https://doi.org/10.1023/A:1021195804379.

 W.M. Wright. General linear methods with inherent Runge-Kutta stability. PhD thesis, ResearchSpace@ Auckland, 2002.

Akram Movahedinejad (a), Gholamreza Hojjati (a) and Ali Abdi (a)

(a) Faculty of Mathematical Sciences, University of Tabriz 29 Bahman Boulevard 51666, Tabriz, Iran

E-mail(corresp.): ghojjati@tabrizu.ac.ir

E-mail: a_abdi@tabrizu.ac.ir

Received May 26, 2016; revised November 30, 2016; published online January 5, 2017

Caption: Figure 1. L-stable choices of ([lambda], [mu]) for the methods with p = s = 1.

Caption: Figure 2. L-stable choices of ([lambda], [mu]) for the methods with p = s = 2.

Caption: Figure 3. L-stable choices of ([lambda], [mu]) for the methods with p = s = 3.

Caption: Figure 4. L-stable choices of ([lambda], [mu]) for the methods with p = s = 4.

Caption: Figure 5. Numerical results for the SGLMs of order p = 1, 2, 3, 4.
```Table 1. Numerical results for IQS methods of orders p = 1 and
p = 2 applied to Problem
1.

k                  method of order p = 1

[parallel][e.sub.h]([bar.x])[parallel]    p

4     2.24 x [10.sup.-6]
5     1.19 x [10.sup.-6]                        0.91
6     6.10 x [10.sup.-7]                        0.96
7     3.10 x [10.sup.-7]                        0.98
8     1.56 x [10.sup.-7]                        0.99

k                 method of order p = 2

[parallel][e.sub.h]([bar.x])[parallel]    p

4     3.87 x [10.sup.-7]
5     9.76 x [10.sup.-8]                        1.99
6     2.45 x [10.sup.-8]                        1.99
7     6.16 x [10.sup.-9]                        1.99
8     1.55 x [10.sup.-9]                        1.99

Table 2. Numerical results for IQS methods of orders p = 3
and p = 4 applied to Problem 1.

k                 method of order p = 3

[parallel][e.sub.h]([bar.x])[parallel]   p

4    1.25 x [10.sup.-7]
5    1.62 x [10.sup.-8]                       2.95
6    2.08 x [10.sup.-9]                       2.97
7    2.67 x [10.sup.-10]                      2.96
8    3.45 x [10.sup.-11]                      2.95

k                method of order p = 4

[parallel][e.sub.h]([bar.x])[parallel]    p

4    6.44 x [10.sup.-8]
5    4.00 x [10.sup.-9]                        4.01
6    2.49 x [10.sup.-10]                       4.01
7    1.54 x [10.sup.-11]                       4.01
8    9.34 x [10.sup.-13]                       4.04

Table 3. Numerical results for Problem 2 solved by the method of
order 4 with ho = 10 3.

tol            Method   ns     nrs   nfe     nJe    ge

[10.sup.-4]    SGLM     47     12    670     438    6.17 x [10.sup.-5]
GLM      68     2     753     427    2.39 x [10.sup.-3]
[10.sup.-6]    SGLM     24     1     286     190    1.34 x [10.sup.-6]
GLM      158    2     1272    636    6.89 x [10.sup.-5]
[10.sup.-8]    SGLM     34     3     325     181    2.14 x [10.sup.-6]
GLM      455    2     3648    1824   3.69 x [10.sup.-6]
[10.sup.-10]   SGLM     64     4     536     268    1.42 x [10.sup.-9]
GLM      1402   1     11216   5608   6.59 x [10.sup.-8]

Table 4. Numerical results for Problem 3 solved by the method of order
4 with ho = 10 3.

tol            Method   ns     nrs   nfe     nJe     ge

[10.sup.-4]    SGLM     24     3     472     368     2.88 x [10.sup.-5]
GLM      141    12    1844    1236    2.54 x [10.sup.-3]
[10.sup.-6]    SGLM     35     5     723     567     2.90 x [10.sup.-6]
GLM      372    2     3295    1803    8.29 x [10.sup.-5]
[10.sup.-8]    SGLM     68     16    1050    718     6.09 x [10.sup.-8]
GLM      1139   1     9147    4591    6.06 x [10.sup.-6]
[10.sup.-10]   SGLM     142    17    1492    860     2.43 x [10.sup.-9]
GLM      3576   2     28652   14344   8.63 x [10.sup.-8]

```

Please Note: Illustration(s) are not available due to copyright restrictions.
COPYRIGHT 2017 Vilnius Gediminas Technical University
No portion of this article can be reproduced without the express written permission from the copyright holder.