\(^\ast \)Department of Mathematics, University of Benin, P.M.B 1154, Benin City, Nigeria,
e-mail: okunoghae01@yahoo.co.uk, mnoikhilo@yahoo.com.

1 Introduction

The representation of second derivative hybrid linear multistep methods (SD-HLMM) as SD-GLM and transformation to its Nordsieck form for easy implementation will be the main focus herein. Nordsieck method was introduced by Arnold Nordsieck in 1962 as an easy means of changing step size in LMM [ 1 ] . Variable step-size implementation is readily achieved for formulas of Nordsieck kinds and in fact, LMM can be represented by Nordsieck formulas, details can be found in [ 1 ] , [ 3 ] , [ 10 ] , [ 12 ] , [ 18 ] and [ 19 ] . Hybrid LMM which can be with one or more off-step points are aimed at bypassing the Dahlquist [ 7 ] order barrier for the conventional LMM in computing the numerical solution of the autonomous initial value problems (IVPs),

Examples of hybrid LMM are in [ 3 ] , [ 4 ] [ 13 ] and [ 15 ] . Indeed, it is somewhat difficult to implement the hybrid LMM in variable step-size because of the problems associated with attempt to put the input and output data in their standard representations at each step of approximations of the solution of 1. The problems associated with the direct implementation of the hybrid LMM make it less attractive compared to the Nordsieck GLM [ 3 ] , [ 4 ] . Indeed, over the years there is now a greater sophistication in the understanding of the GLM,

where \(n=1, 2,\ldots \); \(h=\left( x_{n+1}-x_{n}\right)\) is the stepsize. The GLM (2a) is readily represented also as,

where, \(k\) is the number of input quantities, \(s\) is the number of stages, \(Y=(Y_{1}, Y_{2},. . .,Y_{s})^{T}\) denote the stage values computed in step \(n\), then \(y^{[n-1]}=(y_{1}^{[n-1]}, y_{2}^{[n-1]}\cdots ,y_{k}^{[n-1]})^{T}\) represent the vector of input quantities at the start of step \(n\), while \(y^{[n]}=(y_{1}^{[n]}, y_{2}^{[n]}\cdots ,y_{k}^{[n]})^{T}\) is the output solution and denotes the vector quantities exported to the next step, \(n+1\). The \(F(Y)=(f(Y_{1}), f(Y_{2}),. . .,f(Y_{s}))^{T}\) is the stage derivatives computed in step \(n\) and the matrices defining the structure of the method are \(A=\{ a_{ij}\} \in \Re ^{(s\times s)}\), \(B=\{ b_{ij}\} \in \Re ^{(k\times s)}\), \(U=\{ u_{ij}\} \in \Re ^{(s\times k)}\), and \(V=\{ v_{ij}\} \in \Re ^{(k\times k)}\). The \(c=\left[c_{1},\ldots ,c_{s}\right]^{T}\) is a vector of the stage abscissas. The Nordsieck form of (2b) is

where,

with the Nordsieck vectors in the past and present being respectively,

Hybrid LMM can in general be easily represented in the GLM formalism as in (2b) and when it is possible to define Nordsieck vector in (3b), then it is readily implementable as a variable stepsize Nordsieck GLM in 2. The [ 4 ] has given a practical example of the transformation of hybrid LMM into the Nordsieck form 2 of a GLM in (2b).We are therefore motivated following [4] to consider the transformation of a certain class of hybrid second derivative LMM into its Nordsieck SD-GLM for the purpose of easy implementation in variable step-size. In this regard, the multi-derivative GLM is

in [11], [ 16 ] . This method engages directly the higher order derivatives of the solution of (1) in the implementation of (2b). In particular, the SD-GLM in the Butcher sense is

This class of methods introduced in [2] which incorporates second derivatives of the solution of 1, enhances (2) in order and stability wise for a comparable number of stages with respect to (2b). The terms in (5a) which extend (2) are the matrices \(A^{(1)}=\{ a_{ij}^{(1)}\} \in \Re ^{(s\times s)}\), \(A^{(2)}=\{ a_{ij}^{(2)}\} \in \Re ^{(s\times s)}\), \( B^{(1)}=\{ b_{ij}^{(1)}\} \in \Re ^{(k\times s)}\), \(B^{(2)}=\{ b_{ij}^{(2)}\} \in \Re ^{(k\times s)}\), \(U=\{ u_{ij}^{(j)}\} \in \Re ^{(s\times k)}\), \(V=\{ v_{ij}\} \in \Re ^{(k\times k)}\) and the second derivatives stages \(F'(Y)=(f'(Y_{1}), f'(Y_{2}),. . .,f'(Y_{s}))^{T}\) respectively. The picture of (5a) in Nordsieck form is now

The transformation matrix \(T\) is as in equation (3a). The stability matrix of (5a) is

and the stability polynomial is given by

It is readily obvious that (5a) and (5b) exhibit the same stability characteristics. The same for (2b) and 2. An example of the explicit SD-GLM (5a) is from [ 14 ] ,

The interval of absolute stability of this SD-GLM [ 14 ] is \(\left(-12, 0 \right)\). Again, see [ 2 ] and [ 14 ] for some further examples of the SD-GLM in (5a). Sequential SD-GLM have been derived in [17] and inherent Runge-Kutta stable methods are proposed in [18].
A direct generalisation of GLM of Butcher and O’Sullivan [ 4 ] incorporating multiple hybrid points is

with the last stage and also serving as an output method is given as

where \(f_{n+c_{i}-1}=f(y(x_{n-1}+c_{i}h))\) and \(x_{n-1}=x_{n}-h\). A new SD-GLM extension of this class of methods becomes

with last stage as the method,

This class of methods is of interest in the next section. This will be implemented as first same as last (FSAL), see [3]. To construct an explicit or diagonally implicit SD-GLM from (6a), (6b), and (7a), (7b) set the flag \(\ell \) accordingly as,

By setting \(\ell =0\), \(k=2\), \(s=3\) and the coefficients \(\lbrace a_{ji}^{(1)}\rbrace _{j=0}\), \(\lbrace u_{ji}^{(j+2)}\rbrace _{j=0}^{k}\), \(\lbrace a_{ji}^{(2)}\rbrace _{j=0}^{k}\), \(\lbrace v_{3i}(t)\rbrace _{j=0}^{k}\) and \(\lbrace a_{ji}^{(2)}(t)\rbrace _{j=0}^{k}\) in (6a) and (6b) respectively to zero, [4] using the Cauchy integral formula we obtain the hybrid LMM

By (3a) and (3b) the Nordsieck form of this becomes,

The abscissa vector for the GLM in 16 is \(c = \left(\frac{8}{15}, 1, 1\right)^{T}\). The \(T\) transformation in this representation is

The interval of absolute stability of 16 is \((-1.231, 0)\), see Fig. 1 for its stability plot. The region of absolute stability is the closed region in the left half plane (LHP) of the complex plane in Fig. 1.

Fig. 1. Region of absolute stability (closed region in the LHP) of the method in (9).

An implicit counterpart of 16 is

The \(c = \big(\frac{8}{15}, \frac{2}{5}, 1\big)^{T}\) is the abscissa vector for the implicit counterpart of the method in 16 with the angle of absolute stability, \(\alpha = 77^{\circ }\) and having 17 as its \(T\).

This paper is organized as follows. Section 2 is on the transformation of SD-HLMM(7) into its Nordsieck form which is illustrated for step numbers \(k=1, 2\). In section 3 we derive the SD-HLMM (7) of step number \(k\leq 3\), while section 4 presents the variable step-size implementation. The section 5 presents how the RK-starters for the methods can be obtained and with examples. Section 6, presents results of numerical experiments.

2 The SD-HLMM (7a,b)

Consider an extension of the methods of [5] provided in (7a) and (7b) to SD-GLM. The explicit stage predictors are,

and

for computing the off-step points at \(x_{n+c_{1}-1}\) and \(x_{n+c_{2}-1}\) respectively. The last stage that will also serve as the output method is,

The \(c=\left(c_{1}, c_{2}\right)\in \left[0, 1\right]\), is the stage abscissas vector, the continuous coefficients \(v_{1,1}(t)\) and \(v_{1,2}(t)\) are assumed to be one and zero respectively. The first and the second stage derivatives are \(F(Y)=\left[f(y_{n+c_{1}-1})\quad f(y_{n+c_{2}-1})\right]^{T}\) and \(F'(Y)=\left[f'(y_{n+c_{1}-1})\quad f'(y_{n+c_{2}-1})\right]^{T}\) respectively. Here \(f_{n+c_{i}-1}=f(y(x_{n-1}+c_{i}h)), i=1(1)s\), \(f'_{n+c_{i}-1}=f'(y(x_{n-1}+c_{i}h)), i=1(1)s\) and \(x_{n-1}=x_{n}-h\). The order of the algorithms in 18, 19 and 20 are \(3k-1\), \(3k+1\) and \(3k+2\) respectively.

Similarly, the equivalent methods of the implicit case are

and

with the last stage as an output method,

The prospect of these classes of SD-HLMM is the fact that the explicit 20 and implicit 21 methods have a common \(T\) transformation (3a) for a fixed \(k\) in their Nordsieck GLM. Moreover, the explicit method 18-20 can serve as a starter for the implicit one in 19-21. The SD-HLMM can be transformed into Nordsieck methods. To get started, the SD-GLM representation of the SD-HLMM 18, 19, 20 when \(k=2\) is

while the implicit SD-GLM partitioned matrix of the SD-HLMM 19, 20,21 when \(k=2\) is

where \(B^{(1)}\), \(B^{(2)}\), \(U\), \(V\) are of the same structure as in (17a). The stages, stage derivatives, the inputs and the output approximations for the implicit method in (17a, b) are:

respectively. To achieve our aim we seek a transformation of (17a) into the form

using the Nordsieck vectors of (3b). Similarly, the implicit case (17b) is

The \(UT\), \(T^{-1}B^{(1)}\), \(T^{-1}B^{(2)}\) and \(T^{-1}VT\) are of similar form to that in (18a). To arrive at this the stage values and the stage derivatives at step number \(n\) are

Thus, the \((3+9)\times (3+3+9)\) partitioned matrix of the SD-GLM (17a) is

In the same manner,

for the implicit method in (18b), where \(y_{n-3}\), \(hf_{n-3}\) and \(h^{2}f'_{n-3}\) are the additional input data with coefficients zeros, which will be included to write the method (19a,b) in this \((3+9)\times (3+3+9)\) partitioned matrix, following the ideas of [ 4 ] . The incoming and the corresponding outgoing vectors are now,

The \(T\) obtained through Taylor’s expansion of \(y^{[n]}\) in 22 about \(x_{n}\) using (3a) is

Carrying out the operations \(UT\), \(T^{-1}B^{(1)}\), \(T^{-1}B^{(2)}\), and \(T^{-1}VT\) on (19a,b) yields the Nordsieck algorithms in (18a,b). The transformations above are illustrated by practical examples in the next section.

3 The transformation of SD-HLMM to SD-Nordsieck methods

First, consider the derivation of the SD-HLMM in 18 for \(k=1\) and its transformation to the explicit Nordsieck methods (18a). For the derivation let the polynomial solution of 1 be

to set up the following interpolatory system of linear equations,

where

for the algorithm in 20. For \(k=1\) in 18 the 25 reduces to

Solving this resulting system of linear equations for \(a_{j}, j=0(1)2\) and substituting the resulting values into 18 with \(c_{1}=(x-x_{n-1})/h\) yields the following coefficients:

A predictor exists for any choice of \(c_{1}\). Consider for example, fixing \(c_{1}=\frac{1}{2}\) into 27 which gives thus

Similarly, setting \(N=4\) in 24 yields 25. Solving 25 for \(a_{j}, j=0(1)4\) and substituting the result into 24 lead to

These coefficients in 29 show that an infinite number of output methods 20 exists for step number \(k=1\). Inserting \(c= \big[\frac{1}{2}, 1\big]\) into 29 yields

Following the approach in (17)–22 the Nordsieck form of 28,30 is

The transformation matrix \(T\) in 31 is

The eigenvalues of the stability matrix of (29) are

and the fact that only one of them is not identically zero means that the method will behave like a Runge-Kutta method with stability function \(R(z)= 1+z+\frac{z^2}{2}+\frac{z^3}{6}+\frac{z^4}{24}\), see [ 6 ] . The interval of absolute stability of the algorithm in 31 is thus \([-2.785, 0]\). For the stability plot, see Fig. 2.

Similarly, the abscissa vector \(c=[\frac{2}{3},1]^{T}\) in 27 and 29 gives another fourth order SD-HLMM for step number \(k=1\) as

The Nordsieck form of this is

The transformation matrix \(T\) in 34 is as in 32. The eigenvalues of stability matrix of 34 are \(\lbrace {1+z+\frac{z^2}{2}+\frac{z^3}{6}+\frac{z^4}{72}, 0, 0, 0, 0, 0, 0\rbrace }\), see [19]. There is only one non-zero eigenvalue and is thus RK is nearly stable because of the term \(\frac{z^4}{72}\). The interval of absolute stability of the method in (32) is \([-3.117, 0]\), see
Fig. 2.

Fig. 2. Region of absolute stability of the Nordsieck methods in (29) and (32).

The coefficients of the first predictor of the SD-HLMM18 for step number \(k=2\) are:

If \(c_{1}\neq -1\), a family of predictors exist and \(0\leq c_{1} \leq 1\). Inserting \(c_{1}=\frac{2}{3}\) into the coefficients and substituting the resulting values into 18 gives the first predictor scheme

The relevant equations for the second predictor 19 using 24 are thus obtained and following the above process leads to the coefficients in (A) of the appendix. Fixing \(c_{1}=\frac{2}{3}\) and \(c_{2}= 1\) into (A) gives 19 as

The coefficients in terms of \((c, t)\) are very large expressions hence we fix the abscissa value of the output method as \(t=1\) as defined in section 2 to give the coefficients of the output method in 20 in terms of the other abscissa values \(c= \left[c_{1}, c_{2}\right]\). The expressions are in (B) of the appendix. Now inserting \(c=[\frac{2}{3}, 1]\) into (B) and substituting the resulting values into 20 gives

Composing the methods in 36, 37 and 38 and reorganizing their data in SD-GLM format as in (19a) yields

where

Putting all these details together as in (18) with \(T\) in 23 we find the coefficient matrices of the Nordsieck method to be

with \(k=2\), where \(UT\), \(T^{-1}B^{(1)}, T^{-1}B^{(2)}, T^{-1}VT\) are in appendix C. By verification \(\operatorname {Det}(wI-\overline{V})=w^{k-1}(w-1)\) thus guaranteeing zero stability of the Nordsieck method in (19a). The stability polynomial of 40 is

The interval of absolute stability is \(\left(-0.2342, 0\right)\), see Fig. 4. The stability interval of explicit Nordsieck methods as with explicit LMM generally are usually very severely limited compared to their implicit counterpart; this is compensated by the flexibility in their automatic stepsize implementation allowed by such methods.

From 35, 37 and (B) with \(c=[\frac{1}{2}, 1, 1]^{T}\), the SD-HLMM for \(k=2\) is

Its Nordsieck GLM form is

where \(UT\), \(T^{-1}B^{(1)}, T^{-1}B^{(2)}, T^{-1}VT\) are in appendix D. To distinguish between the wide range of special subcases within 40 and 41, with a view to making a choice for practical use, we consider the stability region for the given abscissae vector \(c=[\frac{1}{2}, 1, 1]^T\). The stability polynomial \(\det (wI-M(z))\) of 41 is

The interval of absolute stability is \(\left(-0.4, 0\right)\), see Fig. 3.

Fig. 3. Region of absolute stability of the Nordsieck methods in (38)
and (40).

4 Variable stepsize implementation of the Nordsieck GLM in (5b)

Take note that Section 3 has shown how 18-20 can practically be represented as the Nordsieck method in (5b) and found in (18a,b). An effective implementation of (18a,b) requires step-size changing and rescaling of the vector of the external approximations as it becomes necessary, see [ 17 ] . The variable stepsize form of (5b) is

where the rescaling matrix is given by

and

where \(n= 1, 2,\ldots ,\quad TOL=10^{-m};\quad m\gg 1\). 44 is the step-size required to advance from step \(n-1\) to \(n\), if the last step was successful, while \(h_{n-1}\) represents the step-size in the last attempt either a successful or a failed step, \(p\) denotes the order of the corrector formula and \(TOL\) is the allowable error estimate specified by the user. The \(\Delta _{n}\) in 44 is the local error estimate with the condition \(\| \Delta _{n}\| {\lt} TOL\). If it passes this condition, the step is successful and the next step-size will be obtained using the step-size increment formula in 44. We compute the local error estimate using a linear combination of the stage derivatives of the current step only. That is,

where \(\Omega _{p+1}\) is the error constant of the output point of the Nordsieck method. Expanding 45 using the Taylor series, it follows that the values of \(\Theta _{i}\) and \(\psi _{i}\) satisfy the resulting system of equations

see [12], page 86-87. Consider for example, the method 31 of order \(p=4\), \(\Omega _{5}=\frac{1}{720}\) with abscissa vector \(c=[\frac{1}{2}, 1]^{T}\) in 44. Solving the resulting system of linear equations in 46 yields no value for the variables \(\Theta _{1}\), \(\Theta _{2}\), \(\psi _{1}\) and \(\psi _{2}\). No local error estimate exist for the methods in 34, 40 and 41 from 45, therefore, an alternative approach to estimate the error is adopted. To control the error in 31 and 34, \(N=3\) in 24 to obtain a method of order 3, one order less than the corrector. The resulting error estimators for 31\((j=1)\) and 34 \((j=2)\) will therefore be based on \(\| y_{n}^{e_{j}}-y_{1}^{[n]}\| \), \(j=1, 2\) with

Similarly, for the Nordsieck methods in 40 \((j=3)\) and 41 \((j=4)\) the methods adopted to control the error can be:

The \(\Delta _{n}\) in 44 is now \(\| y_{n}^{e_{j}}-y_{1}^{[n]}\| \), \(j=3, 4\) and is the error estimate for the proposed Nordsieck methods in 40, and 41. This is defined as the difference between the output methods of proposed Nordsieck methods in 40, and 41, and the alternative scheme used to control the error. As a result of the asymptotic behaviour of \(\| y_{n}^{e_{j}}-y_{1}^{[n]}\| , j= 3, 4\) we use \(r\) in 44 to increase or decrease the error predicted in the next step by multiplying it by \(h_{n-1}\), (see [ 3 ] , pp. 130 - 131, Algorithm 271\(\alpha \)). The availability of initial data \(y_{i}^{[n-1]}\) facilitates easy adjustment of stepsize. The Nordsieck vectors (3b) in 42 are the outgoing and incoming approximation at step number \(n\) respectively. The application of 42 to the scalar test problem \(y'(x)=0\), yields

The zero stability property of the method 42 is determined by the eigenvalues of the matrix \(D(r)T^{-1}VT\), see [ 12 ] and [ 17 ] . Consider for example the matrix \(T^{-1}VT\) in (29a) given by

The eigenvalues of \(D(r)T^{-1}VT\) from (29a) are found to be \(\{ 1,0,0,0,0,0\} \). We obtain similarly the eigenvalues of the methods in 40 and 41 as
\(\lbrace 1,\frac{1}{3}, 0,0 ,0,0,0,0,0 \rbrace \) and \(\lbrace 1,\frac{1}{2}, 0,0 ,0,0,0,0,0 \rbrace \) respectively. This means that the Nordsieck methods are zero-stable.

5 The RK and SD-RK starters

To start up the methods (5b), we find an approximation to the Nordsieck vector

by use of a RK starter

In the search of a high order starter for a given stage number, the SD-RK starter

comes in handy to advance the solution at \(x_{0}\) through one single step to \(x_{0}+h\). The vectors of \(\overline{U}\) and \(\overline{V}\) in (48a) and (48b) are \(\left(1,1,\ldots ,1\right)^{T},\) and \(\left(1,0,\ldots ,0\right)^{T}\) respectively. In the Nordsieck vector, the first and second components of 48 are known from the initial values defining the ODE 1, but we can determine the other components using a Runge-Kutta method with multiple outputs. To obtain (48a) we use the internal stages of the method given by

and for (48b) we use

where \(Y_{i}=y(x_{0}+\widehat{c}_{i}h)\). Expanding (49a) by Taylor series and equating the powers of \(h\) yields the first stage order condition as

for the second stage we have

while the third stage gives

and the fourth stage order conditions are

Also, the fifth stage order conditions are

Let \(\widehat{c}=[0, \frac{1}{4}, \frac{1}{2}, \frac{3}{4}, 1]^{T}\) in 51-55 and solving the resulting system of linear equations we obtain the value \(\overline{a}_{ij}'s\), of matrix \(\overline{A}\).

Similarly, the components of the output method in (48a) are obtained using

The \(t=0\), gives the output value at \(x=x_{0}\). If we set \(t=1\) the output values will be for the first point \(x_{1}=x_{0}+h\). Expanding 56 by Taylor series and equating both sides in powers of \(h\), the first component \(y_{1}^{[0]}=y(x_{0}+th)\), yields

The system of linear equations for the second component \(y_{2}^{[0]}=hy'(x_{0}+th)\) is

The third component \(y_{3}^{[0]}=h^{2}y''(x_{0}+th)\) gives

In the same manner, the fourth component \(y_{4}^{[0]}=h^{3}y'''(x_{0}+th)\) of the output method (48a) gives

Similarly, the fifth component \(y_{5}^{[0]}=h^{4}y^{(4)}(x_{0}+th)\) yields

while the sixth component \(y_{6}^{[0]}=h^{5}y^{(5)}(x_{0}+th)\) gives

The GLM representation of (48b) is defined by

respectively. Again, let \(\widehat{c}=[0, \frac{1}{4}, \frac{1}{2}, \frac{3}{4}, 1]^{T}\) and \(t=0\) in 57-62 and solving the resulting system of linear equations respectively gives \(\overline{b}_{ij}'s\) and their values are in matrix \(\overline{B}\) in (62), given here as,

At \(t=1\) we have

Following this, a starting method for the Nordsieck method in 40 and 41 is derived. The outputs of the scheme 63 are the starter for the explicit Nordsieck method in 31 and 34 respectively. As was noted in [ 12 ] , this idea of using a RK method with multiple outputs as starter for a Nordsieck method was introduced by Gear [ 8 ] . This can readily and analogously extended to SD-RK starters through (48b).

6 Numerical Experiments and Conclusion

To illustrate the application of the schemes, and for the purpose of comparison consider implementation of the methods

(i). order four Adams-Bashforth methods (ABF4) in [ 3 ] , page 116, Table 244(I),

(ii). order four Runge-Kutta method (RK4) in [ 3 ] , page 102, (235i),

(iii). order four second derivative Nordsieck method (SDNM4) in 34,

(iv). order four second derivative GLM with nearly Almost Runge-Kutta methods stability (SDNARK4) in [ 20 ] on the following DETEST problems:

Problem 1:

Problem 2:

Problem 3: A nonlinear chemical reaction.

In Tables 1–3, \(TOL\) is the tolerance, \(ns\) represent the number of total step, \(nfs\) is the number of failed step, \(nfe\) denotes the number of function evaluations, \(Error\) is the \(|y_{exact}-y_{n}|\), while \(\| \Delta _{n}\| _{\infty }\) represents the maximum error of the error estimate. The \(y_{exact}\) and \(y_{n}\) are the exact and the numerical solutions respectively. The results obtained by each of the fourth order methods for Problems 1–3 with \(Tol=10^{-j}, j=2,4,6\) are given in Tables 1–3.

Table 1 shows that the SDNM4 34 (i.e., the order 4 second derivative Nordsieck method) performs better than ABF4, RK4 and SDNARK4 in terms of accuracy and computational cost. The main advantage over other RK methods is that its stage order is two. In addition, the transformation of the SDHLMM 18–20 to Nordsieck method (18a) also enhanced the performance of the SDNM4 (32) in variable step-size implementation.

The numerical results given in Table 2 are the results of evidences for the very good performance of the proposed method compared with the methods considered for comparison. In terms of computational cost, the SDNM4 in 34 outperform the ABF4, RK4, and SDNARK4, for details, see Table 2.

The numerical results given in Table 3 confirm the improved performance of the SDNM4 over ABF4, RK4, and SDNARK4 in terms of accuracy and computational cost. In fact, in all the problems solved the SDNM4 has the smallest number of function evaluations when compared with other methods considered herein. However, the numerical results presented in Tables 1, 2 and 3 show the accuracy of the new scheme. Observe from Tables 1–3, as the tolerance becomes smaller, the better is the accuracy of the SDNM4 in 34; however, SDNARK4 at high tolerance gives poor response time.

Conclusively, we have described in this paper the construction of second derivative Nordsieck methods (5b) for nonstiff differential systems 1. This approach can be extended to a higher step numbers and off-step points in 18–20. We have shown thus that it is possible to extend the methods of [ 4 ] to multiple hybrid points with Nordsieck implementation.

Appendix A

Appendix B

Appendix C

7 Appendix D