Enhancing Accuracy of Runge–Kutta-Type Collocation Methods for Solving ODEs

: In this paper, a new class of Runge–Kutta-type collocation methods for the numerical integration of ordinary differential equations (ODEs) is presented. Its derivation is based on the integral form of the differential equation. The approach enables enhancing the accuracy of the established collocation Runge–Kutta methods while retaining the same number of stages. We demonstrate that, with the proposed approach, the Gauss–Legendre and Lobatto IIIA methods can be derived and that their accuracy can be improved for the same number of method coefﬁcients. We expressed the methods in the form of tables similar to Butcher tableaus. The performance of the new methods is investigated on some well-known stiff, oscillatory, and from the literature.


Introduction
A special type of initial value problems, treated in this work, involves highly oscillatory systems and problems known as stiff systems. In the field of numerical analysis, scientific computing, and engineering, many such problems exhibit these properties [1][2][3]. Ordinary differential equations (ODEs) describing these systems are ill-conditioned in a computational sense and challenging to solve [1]. For such problems, implicit Runge-Kutta (RK) methods are a convenient and frequent choice. The methods gained prominence after Butcher showed that the s-stage fully implicit Gaussian RK method of order 2s is A-stable for all s [4].
A subset of implicit RK methods involves collocation methods. They yield a continuous approximation which makes them suitable for problems where globally continuously differentiable functions are required [5]. Despite being known for decades, the methods are still under continuous development. Publications examine the application of the methods to new problems, the derivation of new strategies, or the development of new procedures for solving nonlinear equations [3]. Three basic types of implicit RK methods, each based on different Gauss-Legendre-type quadrature formulae, are Gauss-Legendre quadrature formulae, Gauss-Radau quadrature formulae, and Gauss-Lobatto quadrature formulae [6].
Perez-Rodriguez et al. [7] applied an implicit RK method of Radau IIA type to semidiscretized, multidimensional partial differential equations (PDEs) of advection-diffusionreaction type. To mitigate the stiffness of these systems, a special iteration technique was applied to solve the resulting implicit relations encountered in the integration steps. Ramos and Garcia-Rubio [8] constructed a new family of implicit methods on the basis of the s-degree Chebyshev approximation using the Chebyshev-Gauss-Lobatto points. Similarly, Vigo-Aguiar and Ramos [9], in addition to constructing an A-stable iRK method, proposed an embedding technique for changing the step size using the Runge-Kutta method and a linear multistep one. Teh and Yaacob [10] recently developed several implicit RK methods, two of them being Gauss-Kronrod methods with five stages and order 8 accuracy. These are based on five-point Gauss-Kronrod quadrature formula. They also developed four implicit RK methods on the basis of Gauss-Kronrod-Radau quadrature formulae [11] and four seven-stage 10th-order implicit RK methods using seven-point Gauss-Kronrod-Lobatto quadrature formula [6]. Liu and Sun [12] studied the Lobatto III-type RK methods. They extended the existing Lobatto methods using the W-transformation and constructed general Lobatto III methods with three parameters of stages two to four. Implicit RK methods play a vital role in solving systems of differential algebraic equations (DAEs) [13]. For example, Jay [14] analyzed the application of RK methods to DAEs of index-2 in Hessenberg form and (see [15]) applied the s-stage Lobatto IIIA-Lobatto IIIB method as a partitioned RK method to integrate Hamiltonian systems with holonomic constraints. A particularly important type of implicit RK methods for stiff and DAE problems involves the Radau IIA methods, implemented in the well-known code for stiff and DAE systems, RADAU5 [16].
In this work, we present a new class of collocation methods for numerical integration of ODEs by basing the derivation on their integral form. The approach resembles the well-known Galerkin method, which is based on the so-called weak formulation of the governing differential equations; nevertheless, we could not find any references involving applying the presented idea to initial value problems and deriving collocation methods. With the approach, we present how the accuracy of RK-type collocation methods (more specifically, the Gauss-Legendre and Lobatto IIIA method) can be increased while retaining the number of method coefficients. In addition to proposing some method instances of the new class of methods, we also demonstrate how the Gauss-Legendre and Lobatto IIIA methods can also be derived with the proposed approach.
In this paper, some theoretical preliminaries are given first (Section 2), followed by a general derivation of the method (Section 3). Particular method families are then presented in Section 4 where their accuracy and linear stability are also investigated. The performance of methods is analyzed using several numerical examples in Section 5. The discussion and conclusion then follow in Sections 6 and 7.

Theoretical Preliminaries
Consider the following initial value problem: with y(x n ) = y 0 , y ∈ R d , for which we seek an approximate solution in the interval [x n , x n + h], with h as the step size. Collocation methods utilize a polynomial, say u, of degree s to approximate the solution y. The polynomial needs to satisfy the initial condition u(x n ) = y 0 and the differential equation The numerical solution of the collocation method is defined by The general class of RK methods is defined as follows [17]: where a fully implicit RK method is obtained when the upper limit r in the sum defining the k i is equal to s. The collocation method is equivalent to the s-stage RK method, Equation (2), with the following coefficients: where l j (τ) is the Lagrange polynomial l j (τ) = ∏ l =j (τ − c l )/ c j − c l , j = 1, . . . , s [18]. The choice of coefficients a ij and b i , and points c 1 , . . . , c s defines a method instance. This work limits the arbitrary choice of c i points to either the roots of Legendre polynomials or a combination of Legendre polynomials of different stages. RK method subclasses based on these collocation points are Gauss-Legendre, Radau, and Lobatto [19]. The Gauss-Legendre quadrature employs a shifted Legendre polynomial of degree s on the interval [0, 1].
The method subclass accuracy is of order 2s and its stability function is the diagonal (s, s)-Padé approximation [20]. The "Radau I quadrature formula" defines c i as the zeros of a polynomial P * s (τ) + P * s−1 (τ), yielding c 1 = 0, while the "Radau II quadrature formula" yields c s = 1 on the basis of the polynomial P * s (τ) − P * s−1 (τ). The "Lobatto III quadrature formula" with both c 1 = 0 and c s = 1 originates from the polynomial P * s (τ) − P * s−2 (τ). The accuracy of the s-stage Radau IA and Radau IIA methods is of order 2s − 1, and their stability function is the (s − 1, s) subdiagonal Padé approximation. The accuracy of the Lobatto III class of method is of order 2s − 2, with the stability function being the diagonal (s − 1, s − 1) Padé approximation for Lobato IIIA and IIIB and (s − 2, s) for Lobatto IIIC [20].
In contrast to classical collocation methods, the novel method class is derived from the integral form of Equation (1).
where v(x) is the test function, arbitrary in general. The usage of the integral approach is, for instance, widely employed in finite element formulation. Its property is a functional approximation of the observed (primary) variable which needs to satisfy the governing (partial) differential equation in an integral fashion. The approximation function, denoted as u , is assumed to approximate the solution of the differential equation y . The system of equations determining the approximation polynomial u is constructed by substituting different test functions v(x) in Equation (5) and imposing the initial condition.
A new class of methods follows from Equation (5) by, firstly, employing different approximation polynomial functions for u . Functions are defined on the basis of roots of the Legendre polynomial P * s (τ) or its combinations (recall: P * s (τ) + P * s−1 (τ), P * s (τ) − P * s−1 (τ) or P * s (τ) − P * s−2 (τ)). The choice of u then determines the quadrature rule in integrating the left-hand side of Equation (5). Secondly, different quadrature formulae are employed for the integration of the right-hand side of Equation (5). This integration can be conducted by any quadrature rule (Gauss-Legendre, Radau I, Radau II, or Lobatto), independent of the choice for the left-hand side.

General Derivation of the Method
Let us introduce in Equation (5) For the approximate solution of y , we consider a collocation polynomial u of stage s.
where k j represents discrete values of u (x n + τ h) at collocation points c j , k j = u x n + c j h , and l j (τ) represents the Lagrange polynomials. The collocation points c j , j = 1, . . . , s can generally be chosen arbitrarily. As discussed in the preliminaries, this work limits their choice to the roots of P * s (τ) or any of its combinations. We continue without setting their numerical values until the next section where certain options are investigated.
By considering Equation (6) in Equation (5), we obtain substituting the integrated approximation function.
Unknowns in Equation (8) are k j , j = 1, . . . , s, the discrete values of the approximation function u . Equations that determine the coefficients are generated by substituting different test functions in Equation (8). For the test function, we choose the Lagrange polynomials of degree s (as is common in finite element method).
The left-hand side of Equation (8) assumes the form of the standard quadrature formula, with l j (τ) and v i (τ) defined by Equations (7) and (10), respectively. The right-hand side of Equation (8) is approximated as where F j represents discrete values of f (x n + τ h, u(x n + τ h)) at collocation points c j , F j = f x n + c j h, u x n + c j h . Functions l j (τ) are the Lagrange polynomials.
where collocation points are c j , j = 1, . . . , s. As in the case of collocation points c j , the numerical values of c j are arbitrary in general and there is no requirement for c j = c j . In this work, we limit their values to the roots of P * s (τ) and some combinations, which we analyze in the next section.
The right-hand side of Equation (8) can thus be expressed as By employing Equations (6) and (10), the right-hand side of Equation (8) finally becomes ∑ s j=1 f x n + c j h, y 0 + h ∑ s m=1 k m a jm q ij , i = 1, . . . , s where with l j (τ) and v i (τ) defined by Equations (13) and (10), respectively.
The final form of Equation (8) is given as with the coefficients specifying a method instance summarized as follows: Equation (16) determines the discrete values of the approximation function k j , j = 1, . . . , s and can be linear or nonlinear. The end-point of the interval y(x n + h) is approximated as with coefficients b j defined with the last line in Equation (17). The numerical values of the coefficients summarized in Equation (17) are found after specifying the collocation points c j , j = 1, . . . , s, and c j , j = 1, . . . , s (Equations (7) and (13), respectively).

Particular Method Families: Their Accuracy and Stability
The behavior of a method instance is determined by the choice of the collocation points c j , j = 1, . . . , s, and c j , j = 1, . . . , s, arbitrary in general. In this work, we investigate the option of selecting the Gauss-Legendre and Lobatto points. We denote method instances according to the quadrature formula used, e.g., Gs|G s employs s-stage Gauss points for c j , j = 1, . . . , s, and s-stage Gauss points for c j , j = 1, . . . , s.
The coefficients of a method instance, Equation (17), are listed in the manner of Butcher [21] as follows: Method families that we introduce in this work are presented in Table 1. Methods are chosen due to their favorable accuracy (e.g., Gs|Gs+1 and Ls|Gs+1) or to demonstrate the option of enhanced integration accuracy of the right-hand side of the ODEs (for instance, by comparing methods Gs|Gs and Gs|Gs+1 or methods Ls|Ls and Ls|Ls+1). Method families of interest are also eLs|Gs and eLs|Gs+1. Their definition only differs in the computation of the first coefficient k 1 , i.e., explicitly by using the initial condition y 0 . Their labels include an "e" before the quadrature rule designation.
The coefficient arrays of the methods in Table 1, obtained for s = 2, 3, 4, are reported in the Supplementary Materials (in the form of analytical expressions), whereas Appendix A provides tableaus for methods Gs|Gs, Ls|Ls, and Ls|Ls+1. For convenience in practical usage, we provide the numerical values of coefficients in Appendix B. In Section 4.1, we prove that the Gs|Gs, and Ls|Ls method families are identical to the Gauss-Legendre and the Lobatto IIIA methods, respectively. The accuracy and linear stability of the methods are analyzed in Section 4.2.

Gauss-Legendre and Lobatto IIIA Method Families as Subclasses of the Novel Method Class
Let us show that (a) the Gs|Gs method family is identical to the Gauss-Legendre method and (b) the Ls|Ls method family is identical to the Lobatto IIIA method.
Proof. For c j = c j , j = 1, . . . , s, (s = s), the coefficients p ij = q ij in Equation (16) can be factored out, yielding exploiting q ij = p ij , i = j = 1, . . . , s. The system of Equation (19) is equivalent to the classical RK system in Equation (2) if and only if the coefficients a jm and c j in Equation (19) are identical to a ij and c j of a RK method, and the matrix P = p ij is invertible. Statement (a) is verified by reviewing the definition of the a jm coefficients (see Equation (17)). These can be defined in the exact same manner in the Gauss-Legendre and Lobatto IIIA methods (see Equation (2)) [18]. To verify statement (b), review the definition of p ij in Equation (17). In this case, the matrix P is a square matrix of size s × s. The nonzero determinant is a consequence of linearly independent equations (i.e., matrix rows), generated from test functions v i . These are, specifically, the Lagrange polynomials of degree s, Equation (10), with points ζ i as the roots of the Legendre polynomial. This is a sufficient condition for the matrix P to be invertible.

Accuracy and Linear Stability
The accuracy of the methods is investigated by integrating the following two test equations: where y e (x) is the exact solution. Equation A is a standard test equation for investigating the accuracy and linear stability of methods. Test equation B is proposed to demonstrate the impact of using different quadrature rules for the integration of the right-hand side of ODEs. The accuracies are determined by applying the methods to both test equations. The approximate and the exact solutions are expanded into a Taylor series about x n and subtracted at x n + h, yielding an error of the approximate solution. The stability analysis was performed using test equation A according to the established procedure. The same procedure was conducted for the Gauss-Legendre, Lobatto IIIA, and Lobatto IIIC method [13].
The results of accuracy and stability analysis are given in Table 2. The methods were grouped as follows by matching behavior:

Method
The results suggest that the behavior of a method is primarily determined by the integration rule of the right-hand side of the ODE, and that taking Gauss points of s stages is equivalent to using Lobatto points with s + 1 stages.
The effect of using a higher integration rule for the right-hand side of ODEs is evident upon comparing the accuracy of Ls|Ls vs. Ls|Ls+1, Gs|Ls vs. Gs|Ls+1, and Gs|Gs vs. Gs|Gs+1 methods ( Table 2). Their accuracy was improved by merely enhancing the integration of the right-hand side of the ODE while retaining the number of coefficients. This improvement in accuracy cannot be general since it was not observed in Ls|Gs vs. Ls|Gs+1, and Gs|Gs vs. Gs|Gs+1 for test equation A (Equation (20)).
The stability of the methods denoted with "e" (where the first coefficient is computed explicitly, eLs|Gs and eLs|Gs+1 in Table 2), is reduced to a (s, s − 1)-Padé approximation. This is the cost of having a smaller system of equations for computing the method's coefficients, having the size of s − 1 × s − 1. Unlike the other proposed method families, "e"-methods are not A-stable. However, an interesting feature of these methods is their accuracy. Their accuracy order was 2s − 1 in test equation A but increased to 2s for eLs|Gs and even 2s + 2 for eLs|Gs+1 in test equation B.

Numerical Examples
We applied the novel and the established methods to several test problems and compared the results, finding that the novel methods offer either comparable or significantly improved performance. In this work, four examples are presented. Example 1 is proposed by the authors, while the remaining examples are well-known stiff and/or highly nonlinear test problems from the literature. In Example 1 (dynamic problem of a mass-spring system), we investigate the application of the methods to problems involving high-frequency oscillations. The behavior of the method with respect to the high initial stiffness of the system was analyzed in Example 2. In Example 3, the methods were examined for problems where the Lobatto IIIC methods might pose stability issues, and, in Example 4, the methods were investigated for a highly nonlinear system. The examples are integrated by Gs|Gs, Gs|Gs+1, Ls|Ls+1, Ls|Gs, and Ls|Gs+1 for stages s = 3, 4, which have orders of accuracy 4 and 6, respectively (Table 2). For comparison, we also analyzed the examples with the Lobatto IIIC method of the same orders of accuracy, 4 and 6.
Analyses were performed in the fixed step size environment for different step sizes h. The approximations were compared to the exact solutions by computing the error as the discrete L 2 norm, defined for the variable y as where y and y e are the approximate and the exact solutions, respectively. N is the number of intervals and x i represents the discrete values of the independent variable. All computations were conducted in Wolfram Mathematica 11. Nonlinear systems of equations were solved using the Newton-Raphson method with relative error tolerance 10 −15 .
The differential equations describe the dynamic behavior of a mass-spring system. The mass is subject to a harmonic excitation and is oscillating about its static equilibrium. The particular excitation used adds high-frequency oscillations to the natural response of the mass, requiring a small step size and making the problem computationally demanding. The approximation errors of different methods for the y 1 variable are reported in Figure 1 with squares and circles for methods of order 4 and 6, respectively. In the same manner, the results for methods eLs|Gs and eLs|Gs+1 are denoted by triangles.
In Figure 1, Ls|Gs and Gs|Gs yielded the same error and, likewise, the Gs|Gs+1 and Ls|Gs+1 pair, as observed for all analyzed cases. The results demonstrate that the behavior of the new method seems to be determined by the integration rule on the right-hand side of the ODE. Gs|Gs+1 and Ls|Gs+1 also performed significantly better than the other methods. Ls|Ls+1 performed slightly better than Lobatto IIIC. A significant finding is that the usage of a different quadrature rule for the right-hand side of the ODE improved the accuracy of a numerical method.
Methods eLs|Gs and eLs|Gs+1 produced equivalent results in the analyzed example, as evident in Figure 1. Although they are of higher order (5 and 7) than the rest, they employ the same number of coefficients per variable (except Lobatto IIIC, which has a coefficient more than Gs|Gs at a given order of accuracy). For example, Gs|Gs of order 4 (G2|G2) employs two coefficients per variable, matching eLs|Gs at order 5; G3|G3 is of order 6 and employs three coefficients per variable and so does eL4|G4 at order 7. Although these methods share the same number of coefficients per step, a significant advantage of the "e" methods is observed.
Methods eLs|Gs and eLs|Gs+1 were applied only to Example 1. They require a very small step size for solving the remaining problems because they are not A-stable.

Example 2
The problem proposed by Frank and van der Houwen [23], with exact solution y 1 (x) = e −2x , y 2 (x) = e −x is a common choice for benchmarking [2,24]. Its intentionally designed stiffness is confirmed by computing the Jacobian of the right-hand side at the initial point and determining the eigenvalues λ 1 = −1004, λ 2 = −1.00199 [2].
The computed errors of different methods for the y 1 variable are reported in Figure 2 with squares and circles for methods of order 4 and 6, respectively. In Figure 2, Ls|Ls+1 outperformed the rest. The problem also seems to be suitable for the Lobatto IIIC, which performed better than Gs|Gs, Ls|Gs, Gs|Gs+1, and Ls|Gs+1. The latter two outperformed the Gs|Gs. However, while Ls|Gs, Gs|Gs+1, Ls|Gs+1, and Gs|Gs have the same number of unknown coefficients per step, Lobatto IIIC has a coefficient more per variable.

Example 3
A problem presented by Vigo-Aguiar et al. [2], and Frank and van der Houwen [23] with the exact solution 1 ( ) = − , 2 ( ) = 0, and 3 ( ) = 1 − − . The same highly stiff Jacobian matrix is computed as in the famous Robertson problem (describing the kinetics of three species), only altered by the nonhomogeneous terms [23]. The computed errors for the 1 variable are reported in Figure 3 with squares and circles for methods of order 4 and 6, respectively.  Table 2.
athematics 2021, 9, x FOR PEER REVIEW Figure 3. The errors of different methods (as the 2 norm) vs. the step size in Example 3 (log-log scale). The method accuracy is reported in Table 2.
In this problem (see Figure 3), Lobatto IIIC was unstable at larger step sizes. Among the methods of order 4, L2|G3 and G2|G3 performed best, while L3|L4 seemed to be the most appropriate order 6 method.

Example 4
The problem presented by Vigo-Aguiar et al. [2], and Frank and van der Houwen [23]: with exact solution 1 ( ) = cos( ), 2 ( ) = 3 ( ) = sin( ). This is a more computationally intensive test problem due to strong nonlinearity. The approximation errors of different methods for the 1 variable are reported in Figure 4 with squares and circles for methods of order 4 and 6, respectively. The best was Ls|Ls+1, followed by Ls|Gs+1 and Gs|Gs+1, and then Gs|Gs.  Table 2.  Table 2.  Table 2.

Discussion
The results show that changing the quadrature rule for the left-hand side of differential Equation (4) had a minuscule effect on the results. However, increasing the order on the right-hand side, i.e., methods #s|#s+1, improved the accuracy at the same number of coefficients. This was demonstrated on test equations analytically, as well as on some numerical examples.
In contrast to the classical derivation of collocation RK-type methods, the new method class is derived from the integral form of the ODE. This approach allows integrating the left-hand and the right-hand side separately, enabling us to employ different quadrature rules for each, resulting in a multitude of method families. Proof is provided that using Gauss quadrature on both ODE sides results in the Gauss-Legendre method and, likewise, Lobatto quadrature returns the Lobatto IIIA method.
By examining the computational cost of the analyzed methods for the presented examples (data not shown), the Gauss-Legendre method somewhat outperformed the proposed methods whereas the Lobatto method was the slowest. However, despite producing computationally somewhat inefficient results, we still believe that the presented approach contributes to our understanding of RK-type collocation methods; it shows how their accuracy can be enhanced (for the same number of method's coefficients) and presents a new way of constructing them.
The proposed method class was found to be especially suitable for (high-frequency) oscillatory and stiff problems. We compared certain novel subclasses to implicit RKtype methods, i.e., the established collocation methods with the same framework and purpose, and we demonstrated the advantage of the new methods. Certain extensions and alternatives could outperform the proposed approach: the embedded formula collocation method [20,25], multivalue methods [26], parallel time integration methods [27], variable step size methods [28], block methods [29][30][31][32], etc. We focused on the popular RK-type methods to avoid an extensive comparison to all the possible alternatives. Tables similar to Butcher tableaus are provided for convenience and comparison to the RK methods. The behavior of methods was in this work demonstrated on some of the well-known ODEs from literature. These include stiff problems, as well as a nonlinear and highly oscillatory system. Systems where the accuracy of methods can behave differently, such as chaotic-type problems [33] or DAE systems [14], are beyond the scope of this work.
Although the implicit RK methods are known to be computationally expensive, much effort has been devoted to developing effective numerical algorithms and efficient computer software packages [3]. Well-known solvers based on the implicit RK approach are RADAU, STRIDE, DIRK, and SDIRK, popular for integrating stiff systems. Implicit RK methods can employ a much larger step size than the explicit ones. The relatively small interval of absolute stability makes the explicit RK methods unsuitable for stiff initial value problems [6]. The presented approach allows improving the implicit RK by increasing the integration accuracy of the right-hand side of the ODE at no additional computational cost.

Conclusions
A new class of methods for numerical integration of ODEs was proposed. The methods allow the usage of differing quadrature rules for the integration of the right-hand side of ODEs while using the same approximation function for the solution of the ODE. We demonstrated how the accuracy of Gauss-Legendre and Lobatto IIIA methods, which can also be derived from the proposed approach, can be increased at the same number of method stages.
In engineering and natural sciences, implicit RK-type methods are still widely used, especially for solving stiff systems [34,35]. As researchers are familiar with Butcher tableaus, the existing program codes can be enhanced in accuracy according to the proposed approach by forming an altered set of equations which is a rather small modification that retains the number of unknowns.