Implementation of explicit Nordsieck methods with inherent quadratic stability

The present paper deals with the implementation in a variable-step algorithm of general linear methods in Nordsieck form with inherent quadratic stability and large stability regions constructed recently by Braś and Cardone. Various implementation issues such as rescale strategy, local error estimation, step-changing strategy and starting procedure are discussed. Some numerical experiments are reported, which show the performances of the methods and make comparisons with other existing methods.


Introduction
We consider initial value problem for the system of ordinary differential equations (ODEs) in the autonomous form: where f : R d → R d is a sufficiently smooth function. Consider the uniform grid t n = t 0 + nh, n = 0, 1, . . . , N, h = (T − t 0 )/N . The general linear method (GLM) in Nordsieck form with coefficient matrices A ∈ R s×s , U ∈ R s×r , B ∈ R r×s , V ∈ R s×s and abscissa vector c ∈ R s is defined by 2) n = 1, 2 . . . , N , where I is the identity matrix of dimension d, s is the number of internal stages and r is the number of input and output approximations. Y i = h i−1 y (i−1) (t n ) + O(h p+1 ), i = 1, 2, . . . , r. From the theory of order and stage order conditions for GLMs, the following theorem follows. We omit the proof, and the reader can refer, for example to [19]. where e cz = [ e c1z e c2z · · · e csz ] T and Z = [ 1 z · · · z r−1 ] T .
By suitable series expansion of order conditions (1.3) and (1.4), algebraic conditions on the coefficient matrices can been derived, as done in [4,9].
Here we consider GLMs in Nordsieck form with s = r−1, and the coefficient matrices A and V are assumed to be of the form a s,1 a s,2 · · · a s,s−1 0 , so the resulting method is explicit and zero-stable. Nordsieck representation of DIMSIMs was discussed in [6] and that of two-step Runge-Kutta methods was derived in [1]. In these paper we deal with the implementation issues related to the methods with p = q = s = r − 1, p = 1, 2, . . . , 5 constructed in [4], which satisfy the property of inherent quadratic stability (IQS), compare [3,4,9,10,11]. This property guarantees that the stability polynomial of the GLM takes the form thus its stability properties depend on the quadratic polynomial Such a property is useful for the practical derivation of highly stable methods (e.g. A-and L-stable) in the implicit case (compare [3,11,19]), and methods with large stability regions in the explicit one case (see [4,9,10]).
The paper is organized as follows. Section 2 is concerned with the variable stepsize formulation of Nordsieck methods and a reliable error estimation, Section 3 deals with some issues necessary to implement such methods in a variable stepsize environment, such as a technique to update the vector of external approximations, a step control strategy and starting procedures. Section 4 provides the numerical evidence originated by implementing Nordsieck GLMs with fixed and variable stepsize, while Section 5 contains some conclusions and further developments of this research.

Variable step-size formulation of the methods
For the variable step formulation of the methods, it is useful to set where y n is an approximation to y(t n ) and z [n] is an approximation to z(t n , h n ), h n = t n − t n−1 , with Then GLM (1.2), on the nonuniform grid t 0 < t 1 < · · · < t N , t N ≥ T , can be formulated as (compare [19]) , I stands for the identity matrix of dimension d.

Error propagation
The following result analyzes the error propagation from one step to another (compare [19], Sec. 8.3, [5,8]). We report the proof, which follows the lines of Theorem 2.1 of [5], for sake of completeness.
Theorem 2. Assume that the input quantities to the current step from t n−1 to t n = t n−1 + h n satisfy where y(t) is the solution to the differential system (1.1) and require that with the same vector β. Here, z(t n , h n ) is the Nordsieck vector corresponding to the solution y(t) of the initial value problem y(t n ) = y n .
Then it follows that (2.3) holds if Proof. From the hypothesis that q = p, it follows that for a certain vector σ. Therefore, Substituting (2.5), (2.6) and (2.7) into (2.3), and using the localizing assumptions given by (2.2), it follows By expanding in Taylor series y(t n−1 + ch n ) and y (t n−1 + ch n ) around t n−1 , equating terms of (2.8) of order up to p + 1 and considering that terms of order up to p cancel out by order conditions, from (2.8) and (2.9), we derive ξ as in (2.4). By expanding in Taylor series y(t n ) around t n−1 , and considering that and successively equating terms of order up to p + 1 in (2.9) (and considering again that terms of order up to p cancel out by order conditions), we get the coefficient E as in (2.4). Last, by expanding in Taylor series h k n y (k) (t n ) around t n−1 into z(t n , h n ), it follows with t p defined as in (2.4) and Substituting (2.12) in (2.10), equating terms of order up to p + 1, we obtain β as in (2.4).

Local error estimation
It follows from Theorem 2 that the local error is We look for estimation of the quantity h p+1 14) The following result holds (compare [5,19]). We report the proof for completeness.
Theorem 3. Consider the GLM (2.1) of order p and stage order q = p and assume that f is sufficiently smooth. The vectors ϕ ∈ R s and ψ ∈ R r−1 in (2.14) satisfy the linear system Proof. Similarly as in [19], we will use expansion: (2.16) By using (2.16) and (2.6) we get: (2.17) We substitute (2.11) on the left of (2.14), use localizing assumptions (2.2) and (2.17) on the right of (2.14), and we obtain Now, by equating the coefficients of terms of order up to p + 1, the thesis follows.
In our cases, p = q = s = r − 1, thus the linear system (2.15) consists of p + 1 equations, so there are p − 1 free parameters. For p ≥ 2, we require that the contribution of both terms in the last equation is the same, i.e., and the number of free parameters is reduced to p − 2. From order p = 3 we fix the remaining free parameters to zero, i.e., Different choices of ψ 3 , ..., ψ p are possible as well, see for example [2]. Now we compute vectors ϕ, ψ for the methods with p = q = s = r − 1, p = 1, 2, . . . , 5, constructed in [4].

Rescale strategy
After the step t n is performed by (2.1), we have computed z [n] ≈ z(t n , h n ). In order to perform next step, we need as a new input vector z [n] ≈ z(t n , h n+1 ), with h n+1 = t n+1 − t n . The most natural way to get z [n] is to rescale the vector with D(δ) = diag(δ, δ 2 , . . . , δ s ), and δ = h n+1 /h n . Other strategies are also available, see for example [19].

Implementation issues
This section is devoted to the description of the issues we have considered for the implementation of our methods in a variable stepsize environment.

Starting procedure
The implementation of Nordsieck methods needs the assessment of a suitable starting procedure for the derivation of the input vector associated to the initial stepsize. For simple problems, the derivatives contained in the input vector can be computed exactly from the given problem, since where the subscripts denote partial derivations with respect to y. Hence, the elements of the initial input vector z [0] can be exactly computed in terms of the elementary differentials of the problem. In more general cases, the employ of an additional starting method for the computation of the initial input vector is required. Following [19,21], we compute the vector z [0] by means of the formulas . . , p, being c = [c 1 , c 2 , . . . , c p ] T ∈ R p is a given abscissa vector for the starting method. The algebraic constraints on the coefficients of the starting procedure which guarantee the above accuracy requirements are (see [19,21] i, k = 1, 2, . . . , p. In [19], Chapter 8, A, B are computed for equally spaced abscissa c, with p = 1, . . . , 4. For p = 5, and equally spaced abscissa c, the tableau of the coefficient matrices of the starting method assume the form This starting method, similarly as those presented in [19,21], is implicit. However, Huang [18] proposed different starting formulae based on singly diagonally implicit GLMs.
Regarding the choice of the initial stepsize of integration h 0 , we assume that for given tolerance tol, being p the order of the method and · is the Euclidean norm. Here, N is an arbitrary number introduced in order to prevent the initial step form being too large; in our tests we set N = 100.

The necessity of a reliable error estimate
The derivation of a reliable error estimate is the first building block for the construction of an efficient strategy of control of the stepsize. In fact, after performing the current step with stepsize h n , we have compared the requested accuracy tol with the achieved error estimate, by checking if the following inequality holds based on the error estimate above described, in order to compute the next value of the stepsize accordingly. If the control (3.2) is not satisfied, the stepsize h n+1 used in the next step is halved. Otherwise, the new stepsize h n+1 is suitably chosen, according to one of the following step changing strategy.

Step control strategies
A standard step control strategy (see [16]) , which only depends on the estimate computed in the previous step, can often determine useless stepsize rejections. Thus, Gustafsson, Lundh and Söderlind [15] introduced a different stepsize control, the so-called PI stepsize control, based on control theory arguments [5,15,17]. The PI control involves the estimation of the local errors related to the two most recent subintervals of the discretization, i.e., where σ 1 and σ 2 are constant values which must be suitably chosen. As in [14,17] we have experimentally found some values for σ 1 and σ 2 which make the PI stepsize control competitive with the standard strategy: such values are where p the is order of the method.

Constant stepsize implementation
We apply methods of order p = 1, 2, . . . 5 to linear test equation (4.1), with λ = 40, T = 1 and fixed stepsize h = T /N , where N = 2 i+1 ·10, i = 1, 2, . . . , 10. We calculated the norm ||e h (T )|| of error at endpoint (here, error is the difference between exact solution and the one obtained by the method) and the effective order of convergence p eff given by The results are collected in Tab. 1, and they confirm the theoretical order of convergence. Table 1. Numerical results of method of order p = 1, 2, 3, 4, 5 for problem (4.1), T = 1, with constant stepsize. In Tab. 2 we compare the error of Nordsieck method with IQS of order p = 4 and method with IRKS of the same order, corresponding to η = 3/5 with error constant E = 1/300, whose coefficients are listed in [7]. We considered the test problem (4.1) on interval [0, 10] for some chosen values of λ and also test problem (4.2). The results show that our methods converge for a larger value of stepsize with respect to IRKS methods. We notice that for the same order our methods are also cheaper, since IRKS method has s = p + 1 internal stages (compare [19]), while our methods have s = p.

Variable stepsize implementation
In this section we show numerical results obtained by optimal methods found in [4] with p = q = s = r − 1, in a variable stepsize implementation. We adopted the rescale strategy described in Sec. 2.3 and we estimate the local error as described in Sec. 2.2.
As regards the estimation of the local error, we employ the error estimate described in Sec. 2.2 and compare it with the error computed as difference between the numerical solution and the exact one (when available), or between the numerical solution and the reference solution computed by applying the Matlab routine ode15s with the maximum precision (when the exact solution is not available).
Concerning the selection of the stepsize, we apply both strategies described in Sec. 3.3 and compare them in terms of their efficiency. Such comparison takes into account the number of rejected steps with respect to those accepted and the smoothness of the patter of the stepsize on the overall integration interval, helpful to test the step changing strategy also taking into account the stability properties of the method. In fact, the more the pattern is smooth, the more the values of the stepsize lie within the stability region and, consequently, the more the number of rejections of the stepsize is small. Fig. 1 and Fig. 2 report the results originated by applying the methods of order respectively 3 and 4 to (4.2), by employing both the standard and PI stepsize changing strategy. We observe that, in both cases, the error estimate is coherent with the true local error. Moreover, by employing the PI stepsize control, we observe a lower number of rejected steps which provides a lower computational cost: this advantage is commonly acknowledged in the literature (compare [17] and [14], where the variable stepsize implementation of two-step Runge-Kutta methods [12] based on modified collocation techniques [13] is presented).
Tab. 3-6 show the results achieved from the application of methods order from 1 up to 5 in variable stepsize environment, with several values of the tolerances (tol), for the numerical solution on (4.2) and (4.3). We compare global error (ge), total number of steps (ns), number of rejected steps (nrs), and total number of function evaluations (nf e). We observe that for more demanding tolerances high order methods have to be preferred. Moreover, in any reported case we observe the superiority of the PI stepsize controller in gaining the solution with a lower computational effort. Finally, Tab. 7 and Tab. 8 provide a comparison between our method of order 4 (named NORD4 method) with the embedded pairs of continuous explicit  Runge-Kutta methods of order 3 and 4 (named CERK43 method, see [2]). It results that NORD4 method is able to achieve the accuracy of CERK43 with a lower computational cost. The requested computational effort is even lower when the variable stepsize strategy is based on the PI controller.

Concluding remarks
We have treated the implementation of explicit GLMs in Nordsieck form with quadratic stability and large stability regions, constructed in [4]. We addressed the issues related to the implementation in a variable step algorithm, in particular: local error estimation, computation of the input vector for the new stepsize, stepsize changing, starting procedure. Numerical experiments carried on significative test examples proved that: • In a fixed stepsize algorithm, our methods converge with a larger stepsize with respect to IRKS methods.  • In a variable stepsize algorithm, our methods are more efficient than other existing methods.
• PI controller strategy considerably reduces the number of rejected steps, and therefore improves the efficiency of our methods.
The next step of our research will be the development of a variable step variable order algorithm. The design of a reliable order changing strategy is based on a fundamental building block which is provided by an accurate estimation of the higher order terms in the local error expansion. We aim to derive such an estimate by suitably extending the results reported in [8].