Superconvergent interpolants for Gaussian collocation solutions of mixed order BVODE systems

1.   Introduction

The numerical solution of boundary value ODEs (BVODEs) arises as a key task in a wide variety of applications such as the modelling of optimal control of electric vehicles ^[31], tropical fruit heat treatment ^[16], black holes ^[24], enzyme kinetics ^[30], and metal uptake in biointerphases ^[9].

The general form for a BVODE system with separated boundary conditions, written in first order form is,

However, most BVODEs do not arise naturally in this form. It is common for higher derivatives to be present and in fact almost all of the problems considered in the applications cited above and in the well known BVODE collection in ^[4] include equations in which higher derivatives appear. These are referred to as mixed order systems. A standard form for a mixed order system consisting of $n$ differential equations of orders $m_1, m_2, \dots, m_n$, respectively, is

where

where $y_j^{(\ell)}(t)$ is the $\ell$th derivative of the $j$th solution component, $y_j(t)$.

The high quality Fortran COLSYS/COLNEW collocation software package, ^[2,3,6] has been widely used for the numerical solution of BVODEs for several decades. More recently, interfaces to the package have been developed for several computing platforms such as Scilab, R, and Python. The papers cited at the beginning of this section are examples of investigations where COLSYS/COLNEW is used directly or through one of these platforms.

COLSYS/COLNEW is able to directly handle BVODEs having the more general mixed order form (1.2). While it is straightforward to convert a BVODE system of the form (1.2) to the form (1.1), there are advantages to treating the mixed order form directly; these include user convenience, efficiency improvements, and higher continuity for some of the approximate solution components; see, e.g., ^[4] and ^[26], for further discussion. See, in particular, ^[4], Page 252, where the authors comment that, for the application of collocation directly to the BVODE mixed order system form rather than the corresponding first order system form, "the saving in computation and storage is substantial". See also, Table 7.3.1 of ^[12], where the computation times for COLSYS/COLNEW applied to a standard test problem are presented; the costs for the solution of the mixed order system form are approximately one third of those associated with the first order system form.

COLSYS/COLNEW employs an iterative, adaptive error control strategy in order to obtain a continuous solution approximation with an associated error estimate that satisfies a user-specified tolerance. The solver assumes that the solution components, $\{ y_j(t) \}_{j = 1}^n$, are approximated by piecewise polynomials of degree $k+m_j$, which are determined by imposing boundary and continuity conditions (the approximation to $y_j(t)$ will be required to have $C^{(m_j-1)}$-continuity), and by imposing collocation conditions that require the approximate solution to satisfy the BVODE at a set of $k$ collocation points on each subinterval of a mesh that partitions the problem domain. The derivatives of the $\{ y_j(t) \}_{j = 1}^n$ solution components are approximated by the corresponding derivatives of the piecewise polynomials that approximate $\{ y_j(t) \}_{j = 1}^n$.

On each subinterval, the collocation points are the images of the set of $k$ Gauss points mapped onto $[0, 1]$. It can be shown - see, e.g., ^[4] - that the $\underline z(t)$ values at the mesh points are substantially more accurate than elsewhere. Let $h = { \max_{i = 1}^{N} h_i}$, where $h_i = t_{i+1} - t_i$, and where $\{ t_i \}_{i = 0}^N$ are the mesh points that partition the problem interval. Then, under appropriate assumptions ^[4], the global error associated with $z_j(t_i)$ is $O(h^{2k})$, for all components of $\underline z(t_i)$. On the other hand, at a non-mesh point, $t$, the global error associated with the component of $\underline z(t)$ equal to $y_j^{(\ell)}(t), j = 1, \dots, n, \ell = 0, \dots, m_j-1$, is $O(h^{k+m_j-\ell})$, provided that $k+m_j-\ell \leq 2k$. Therefore, when $2k > k+m_j-\ell$, the mesh point values are more accurate than the approximations at non-mesh points, and the former are said to be superconvergent.

An essential point is that COLSYS/COLNEW controls an error estimate for the continuous collocation solution rather than an error estimate for the mesh point solution values. Due to the lower accuracy of the collocation solution at non-mesh point values, it takes substantially more computational time and memory for COLSYS/COLNEW to obtain a satisfactory solution than it would if the continuous approximate solution had the same high accuracy that the mesh point values have. More mesh iterations are required to obtain a satisfactory continuous collocation solution and the meshes upon which the solution is constructed have to be much finer, which implies more subintervals, which implies higher computational and storage costs.

1.1.   Superconvergent interpolants for collocation solutions

As explained above, a primary issue that impacts substantially on the efficiency of the computation performed by COLSYS/COLNEW is that the continuous collocation solution has substantially lower accuracy than that of the mesh point collocation solution values. A key idea in efforts to address this issue have focused on developing interpolants for all components of $\underline z(t)$ that have an error that is $O(h^{2k})$ (the mesh point error) for any $t$ in the problem domain. Such interpolants are referred to as superconvergent interpolants (SCIs).

There have been a number of papers that have investigated the development of SCIs for collocation solutions. To our knowledge, the first such effort was ^[33], which considered the development of SCIs based on superconvergent mesh point values from several adjacent subintervals; difficulties are reported for highly non-uniform meshes. Subsequently, ^[34] and ^[20] considered SCIs based on secondary collocation; this involves performing an additional collocation computation, similar in cost to the original collocation computation, and efficiency is therefore a concern.

The paper, ^[13], describes how to develop SCIs for first order BVODE systems, (1.1), based on the use of continuous Runge-Kutta (CRK) methods - see, e.g., ^[28,29,36]. The key idea is to embed information from the collocation solution within a CRK method which is constructed to have the same superconvergent accuracy as the mesh point collocation values have. The primary contribution of that paper is the derivation of CRK methods that provide the basis for SCIs for collocation solutions for $k$, the number of collocation points, equal to 1, 2, 3, and 4, leading to SCIs of orders 2, 4, 6, and 8, respectively.

In ^[13], the SCI is computed in a post-processing step to augment the collocation solution computed by COLSYS/COLNEW. That is, COLSYS/COLNEW is first called with a given user tolerance and the code returns with a continuous collocation solution that has an error estimate that satisfies the tolerance. Then the routine that implements the SCI is called to obtain a continuous approximate solution that has an error that is of the same order of accuracy as that of the mesh point values. This means that the SCI is substantially more accurate than the user tolerance. What of course is more desirable is to have the code return an SCI with a corresponding error estimate that meets the requested tolerance. This would mean that COLSYS/COLNEW would compute a continuous collocation solution with an intermediate level of accuracy and then the SCI would be used to augment this solution to obtain a continuous approximate whose accuracy would meet the requested tolerance. This would lead to a more efficient computation because COLSYS/COLNEW would only be required to compute a collocation solution of intermediate accuracy.

This idea has been pursued in recent work, ^[1], which has involved the development of COLNEWSC, a modification of COLSYS/COLNEW that employs the SCI within the code, as the primary numerical solution that is returned to the user. On a given mesh, COLNEWSC first computes the collocation solution, and then from this collocation solution, an SCI is constructed. An error estimate for this SCI is also computed and used within an adaptive error control algorithm to refine the mesh if necessary. The goal of this solver is to obtain an SCI for which the associated error estimate satisfies the user tolerance. In ^[1], COLNEWSC has been shown to be substantially more efficient than COLSYS/COLNEW since it is able to obtain a sufficiently accurate numerical solution using much coarser meshes than those employed by COLSYS/COLNEW due to the additional accuracy delivered by the SCI on a given mesh.

However, the algorithm upon which COLNEWSC is based requires that the BVODE be expressed as a first order system. Thus, when COLSYS/COLNEW is applied directly to a mixed order BVODE system, the advantages associated with the introduction of an SCI are not available. And, as mentioned above, in a majority of cases, BVODEs arise in their natural form as mixed order systems, which means that they can be treated directly by COLSYS/COLNEW, but not by COLNEWSC.

1.2.   SCIs for mixed first/second order BVODE systems

In order to address this issue, the goal of this paper is to present a derivation of generalizations of CRK methods that can serve as the basis for SCIs for the most common subclass of (1.2) known as mixed first and second order BVODE systems. These systems have the form,

with boundary conditions

In this paper, we will derive schemes based on either Hermite interpolants or extensions of continuous Runge-Kutta-Nyström (CRKN) methods - see, e.g., ^[18,23]. These schemes can provide the basis for SCIs for continuous approximations of $\underline y_1(t), \underline y_2(t)$, and $\underline y_2'(t)$, that will have the same superconvergent order of accuracy as the mesh point values. We will also consider preliminary implementations of these SCIs that are used in a post-processing step, similar to what was done in ^[13].

The computational costs associated with the construction and evaluation of an SCI for either a first or mixed order BVODE system are quite small compared to the costs associated with obtaining the collocation solution upon which the SCI is based. The primary costs associated with the computation of the collocation solution involve the setup and solution of a sequence of almost block diagonal linear systems performed within the context of a Newton type iteration. On the other hand, the costs associated with constructing an SCI involve only a few function evaluations per subinterval. Later in this paper, we provide numerical results similar to those given in ^[13], where it was shown that the cost of constructing an SCI is only a small fraction of the cost of the computation associated with obtaining the collocation solution upon which the SCI is based.

The significance of the work described in this paper is that it will provide the basis for an extension of COLNEWSC to allow for the direct treatment of mixed first and second order BVODE systems. This will mean that this new version of COLNEWSC will return an SCI for which the corresponding error estimate satisfies the user tolerance. Because the SCI is substantially more accurate than the collocation solution, this will mean that COLNEWSC will be able to terminate more quickly using a coarser mesh than that which is used by COLSYS/COLNEW. This new version of COLNEWSC will be able to directly handle mixed first and second order systems just as COLSYS/COLNEW can, thereby taking advantage of the user convenience and efficiency that is achieved by COLSYS/COLNEW for such systems. As well, the approach described in this paper suggests how the treatment of (1.4) can be generalized to allow for the treatment of general mixed order systems (1.2).

The schemes described in this paper may also be relevant for Gaussian collocation software for boundary value differential-algebraic systems ^[5], where the challenge will be to construct SCIs for the algebraic solution components. The schemes from this paper may also be applicable to software where Gaussian collocation is employed for the spatial descritization of PDEs - see, e.g., ^[1,17,22,32] - that involve the second spatial derivative of the solution; it may be possible to adapt the schemes developed in this paper to obtain SCIs for the collocation solutions computed by these solvers.

Work that is related to the investigation discussed in this paper includes an alternative approach to developing SCIs based on boot-strapping ^[10]. The paper ^[14] describes the extension of this boot-strapping approach to general mixed order BVODE systems (1.2). Hermite-Birkhoff interpolants - see, e.g., ^[21] - are employed to provide SCIs for components of ${\underline z}(t)$ that approximate the solution components, $y_j(t), j = 1, \dots, n$. Continuous approximations for the derivatives of $y_j(t), j = 1, \dots, n$, are obtained by differentiating these SCIs; however, each differentiation leads to a decrease by one in the order of accuracy and therefore the approximation for a component of $\underline z(t)$ of the form $y_j^{(\ell)}(t)$ will have a global error that is only $O(h^{2k-\ell})$. Thus, in this approach, only the approximations to the $y_j(t)$ components will have the full superconvergent accuracy of the mesh point collocation solution values. As well, in ^[13] the boot-strapping approach is shown to require more evaluations of the right hand side of (1.1) than does the approach based on the use of CRK methods described in ^[13], which impacts negatively on the efficiency of these schemes. In more recent work, ^[35], the boot-strapping algorithm considered in ^[10] and ^[14] is extended to obtain SCIs for collocation solutions of Volterra integro-differential equations with delay.

This paper is organized as follows. We show, in Section 2, that it is possible to reformulate collocation methods, applied to (1.4), in a form that will allow us to show that they fit within the framework of Runge-Kutta (RK) and Runge-Kutta-Nyström (RKN) methods. This is followed, in Section 3, by a review of RKN and CRKN methods, Hermite interpolants, and a generalization of the CRKN methods, known as continuous parameterized implicit Runge-Kutta-Nyström (CPIRKN) methods. In Section 4, the main section of the paper, we derive specific CPIRKN methods that can provide the basis for SCIs of orders 2, 4, 6, and 8, for collocation solutions of (1.4). In Section 5, numerical results are provided to demonstrate the improved accuracy, orders of convergence, and low computational costs of the SCIs. Our summary, conclusions, and suggestions for future work are provided in Section 6.

2.   Collocation methods in RK/RKN form

For the numerical solution of (1.4), and for the mesh, $\{ t_i \}_{i = 0}^N$, assume piecewise polynomials, i.e., collocation polynomials, $\hat {\underline u}_1(t)$ and $\hat {\underline u}_2(t)$, of degrees $k+1$ and $k+2$, respectively, which approximate $\underline y_1(t)$ and $\underline y_2(t)$, respectively. These collocation solutions must satisfy collocation conditions at $k$ Gauss points per subinterval; on the $i$th subinterval, $[t_i, t_{i+1}]$, these conditions have the form,

where $r = 1, \dots, k$, $\hat t_{r} = t_i + \rho_r h_i$, and $\rho_r$ is the image of the $r$th Gauss point mapped onto $[0, 1]$. It is also required that $\underline {\hat u}_1(t)$ have $C^0$ continuity and that $\underline {\hat u}_2(t)$ have $C^1$ continuity.

The mesh point superconvergence results ^[4] imply that,

and the non-mesh point convergence results ^[4] imply that,

for $\theta \in (0, 1)$. Therefore, for $k > 1$, the $\underline {\hat u}_1(t)$ and $\underline {\hat u}_2'(t)$ approximations at the mesh points will have a higher order of accuracy than elsewhere, and, for $k > 2$, all three of the $\underline {\hat u}_1(t)$, $\underline {\hat u}_2(t)$, and $\underline {\hat u}_2'(t)$ approximations at the mesh points will have a higher order of accuracy than elsewhere.

Let $\underline y_{1, i} = \underline {\hat u}_1(t_{i})$, $\underline y_{2, i} = \underline {\hat u}_2(t_{i})$, and $\underline y_{2, i}' = \underline {\hat u}_2'(\hat t_{i})$, and let $\underline {\hat u}_{1r} = \underline {\hat u}_1(\hat t_{r})$, $\underline {\hat u}_{2r} = \underline {\hat u}_2(\hat t_{r})$, and $\underline {\hat u}_{2r}' = \underline {\hat u}_2'(\hat t_{r})$. Then, following from the collocation and continuity conditions, it can be shown that (see ^[27]), on the $i$th subinterval,

and

where,

for $r = 1, \dots, k$, and,

for $r, j = 1, \dots, k$, where $\{ {\hat b}_j(\theta)\}_{j = 1}^k$ are the Lagrange interpolating polynomials for the abscissa set $\{ \rho_r \}_{r = 1}^k$.

Similarly, it can be shown that,

and

where

for $j = 1, \dots, k$. We note that the same collocation method is used for the approximation of $\underline y_1(t)$ and $\underline y'_2(t)$.

3.   CRKN methods, Hermite interpolants, and CPIRKN methods

3.1.   Discrete and continuous RK and RKN methods

Assuming a first order ODE system, (1.1), and given solution approximations, $\underline y_i$, at the point $t_i$, and $\underline y_{i+1}$, at the point $t_{i+1} = t_i + h_i$, an $s$-stage RK method relates these solution approximations through the equation,

where $\hat t_r = t_i + c_r h_i$ and where

for $r = 1, \dots, s$. The $\underline f(\hat t_1, \underline {\hat y}_1), \dots, \underline f(\hat t_s, \underline {\hat y}_s)$ values are called the stages of the RK method. An RK method of a desired order of accuracy is determined by requiring that the coefficients of the method, $\{c_r, b_r, a_{r, j}\}_{r, j = 1}^{s}$, satisfy a set of RK order conditions.

For the second order ODE system, $\underline y''(t) = \underline f(t, \underline y(t), \underline y'(t))$, an $s$-stage RKN method relates solution and derivative approximations, $\underline y_{i+1}$ and $\underline y_{i+1}'$, at $t_{i+1}$, to solution and derivative approximations, $\underline y_i$ and $\underline y_i'$, at $t_i$, through the formulas,

and

where

and

for $r = 1, \dots, s$. The $\underline f(\hat t_1, \underline {\hat y}_1, \underline {\hat y}_1'), \dots, \underline f(\hat t_s, \underline {\hat y}_s, \underline {\hat y}_s')$ values are called the stages of the RKN method. An RKN method of a desired order of accuracy is determined by requiring that the coefficients of the method, $\{c_r, b_r, {b}_r', a_{r, j}, {a}_{rj}'\}_{r, j = 1}^{s}$, satisfy a set of RKN order conditions.

A comparison of the expressions for the collocation approximations given in the previous section with the general forms for the RK and RKN methods given above shows that the collocation methods are special cases of the RK and RKN methods. In particular, we note that the collocation formula, (2.8), (2.4), is a special case of the RK formula, (3.1), (3.2), and that the collocation formulas, (2.9), (2.5), (2.10), (2.6), are special cases of the RKN formulas, (3.3), (3.5), (3.4), (3.6).

A CRK method extends an RK method to provide a solution approximation at any point within the current subinterval. An $s$-stage CRK method approximates the solution, $\underline y(t)$, on the $i$th subinterval, $[t_i, t_{i+1}]$, by $\underline u(t)$ where, for $t \in [t_i, t_{i+1}]$, and $t = t_i+\theta h_i$,

where $\underline {\hat y}_r$ is defined as in (3.2) and $b_r(\theta)$ is a (weight) polynomial in $\theta$. A CRK method of a desired order of accuracy is determined by requiring that the coefficients of the method, $\{c_r, b_r(\theta), a_{r, j}\}_{r, j = 1}^{s}$, satisfy a set of continuous RK order conditions.

Similarly, CRKN methods, - see, e.g., ^[15,23,26], and references within - extend RKN methods to provide solution and derivative approximations at any point within the current subinterval. An $s$-stage CRKN method approximates the solution, $\underline y(t)$, on the $i$th step by $\underline u(t)$, where,

and approximates the derivative, $\underline y'(t)$, by $\underline {\bar u}(t)$, where,

with $\underline {\hat y}_r$ and $\underline {\hat y'}_r$ defined as in (3.6) and (3.5), and where $b_r(\theta)$ and $\bar b_r(\theta)$ are (weight) polynomials in $\theta$. A CRKN method of a desired order of accuracy is determined by requiring that the coefficients of the method, $\{c_r, b_r(\theta), {b}_r'(\theta), a_{r, j}, {a}_{rj}'\}_{r, j = 1}^{s}$, satisfy a set of continuous RKN order conditions.

In this paper, we will employ CRK and CRKN methods in order to obtain SCIs for the collocation solutions of mixed first and second order systems (1.4). We will use a CRK method to handle the first order part of the system and a CRKN method to handle the second order part of the system.

In order to simplify the derivation of these methods, we will require that the CRK and CRKN methods have the same abscissa, $\{c_r\}_{r = 1}^s$, that the $\{a_{rj}\}_{r, j = 1}^s$ coefficients of the CRK method be equal to the $\{{a}_{rj}'\}_{r, j = 1}^s$ coefficients of the CRKN method, and that the $\{b_r(\theta)\}_{r = 1}^s$ weight polynomials of the CRK methods be equal to the $\{\bar b_r(\theta)\}_{r = 1}^s$ weight polynomials of the CRKN method. We note that the collocation methods of the previous section are of this type.

Then, since the coefficients and weight polynomials used in the $\underline y_1(t)$ approximations will be the same as those used in the $\underline {y}_2'(t)$ approximations, there will be no need to provide the coefficients for the $\underline y_1(t)$ approximations separately; rather we can simply present the coefficients of the CRKN method with the understanding that the part of the CRKN formula that is used for the $\underline {y}_2'(t)$ approximation will also be used for the $\underline {y}_1(t)$ approximation.

Since COLSYS/COLNEW allows the user convenient access to the (already computed) collocation function evaluations, (2.7), the efficiency of the computation is improved when we embed these function evaluations within the CRKN methods we derive; i.e., the collocation function evaluations will be employed as stages of the CRKN methods. This is possible because, as we have seen, it is possible to represent a collocation method for a mixed first and second order system in terms of an RK method and an RKN method.

3.2.   Hermite interpolants

For $k = 1$ and 2, it is possible to employ Hermite interpolants as the basis for SCIs of the desired orders of accuracy.

Once a collocation solution for (1.4) is obtained, we will have, for the $i$th subinterval, $[t_i, t_{i+1}]$, the six superconvergent data values,

We can then compute the corresponding endpoint stages,

and then construct a Hermite interpolant, $\underline u_1(t) \approx \underline y_1(t)$, having the form

where $d_{10}(\theta)$, $d_{11}(\theta)$, $d_{12}(\theta)$, and $d_{13}(\theta)$ are the Hermite cubic polynomials. They are defined to give $\underline u_1(t_i) = \underline y_{1, i}$, $\underline u_1(t_{i+1}) = \underline y_{1, i+1}$, $\underline u_1'(t_i) = \underline f_{1, i}$, $\underline u_1'(t_{i+1}) = \underline f_{1, i+1}$. The polynomials, $d_{10}(\theta)$ and $d_{11}(\theta)$, satisfy the condition that, $d_{10}(\theta) = 1 - d_{11}(\theta)$. This Hermite interpolant has an interpolation error that is $O(h^4)$.

Similarly, we can compute the endpoint stages,

and form a Hermite interpolant, $\underline u_2(t) \approx \underline y_2(t)$, having the form,

where $d_{20}(\theta), \dots, d_{25}(\theta)$ are the Hermite quintic polynomials. They are defined to give $\underline u_2(t_i) = \underline y_{2, i}$, $\underline u_2(t_{i+1}) = \underline y_{2, i+1}$, $\underline u_2'(t_i) = \underline y_{2, i}'$, $\underline u_2'(t_{i+1}) = \underline y_{2, i+1}'$, $\underline u_2''(t_i) = \underline f_{2, i}$, $\underline u_2''(t_{i+1}) = \underline f_{2, i+1}$. The polynomials, $d_{10}(\theta)$, $d_{11}(\theta)$, $d_{22}(\theta)$, and $d_{23}(\theta)$, satisfy the conditions that, $d_{20}(\theta) = 1 - d_{21}(\theta)$ and that $d_{22}(\theta) = 1 - d_{23}(\theta)$.

We can also form a Hermite interpolant, $\underline {\bar u}_2(t) \approx \underline {y}_2'(t)$, of the form,

where $\bar d_{20}(\theta)$, $\bar d_{21}(\theta)$, $\bar d_{22}(\theta)$, and $\bar d_{23}(\theta)$ are the Hermite cubic polynomials. They are defined to give $\underline u_2(t_i) = \underline y_{2, i}$, $\underline u_2(t_{i+1}) = \underline y_{2, i+1}$, $\underline u_2'(t_i) = \underline f_{2, i}$, $\underline u_2'(t_{i+1}) = \underline f_{2, i+1}$. The polynomials, $\bar d_{20}(\theta)$ and $\bar d_{21}(\theta)$, satisfy the condition that, $\bar d_{20}(\theta) = 1 - \bar d_{21}(\theta)$. The Hermite interpolant, (3.14), has an interpolation error of $O(h^6)$ and the Hermite interpolant, (3.15), has an interpolation error that is $O(h^4)$.

We note that the interpolation conditions will imply that $\underline u_1(t)$ and $\underline {\bar u}_2(t)$ will have $C^1$-continuity and that $\underline u_2(t)$ will have $C^2$-continuity.

When $k = 1$ or 2, the superconvergent collocation mesh point approximations have errors that are $O(h^2)$ and $O(h^4)$, respectively, and therefore the Hermite interpolants, (3.12), (3.14), (3.15), can be used to provide SCIs. For $k = 3$, the superconvergent collocation mesh point approximations have errors that are $O(h^6)$, and thus (3.14) can be used to provide an SCI for the $\underline y_2(t)$ approximation but (3.12) is not accurate enough to provide an SCI for $\underline y_1(t)$ nor is (3.15) accurate enough to provide one for $\underline {y}_2'(t)$. For $k = 4$, none of the Hermite interpolants are accurate enough to provide SCIs. We will consider these points further in the next section.

3.3.   CPIRKN methods

For first order BVODE systems, (1.1), we have solution information associated with both endpoints of the $i$th subinterval, i.e., $\underline y_{i}$, $\underline y_{i+1}$. It is therefore useful to rewrite an RK method so that the stages have explicit dependence on the solution information from both endpoints. This leads to the methods known as parameterized implicit RK (PIRK) methods ^[11]; with the additional restriction that the each stage depend only on previously computed stages, we get the mono-implicit RK (MIRK) methods ^[7,8]. By introducing weight polynomials, one can obtain continuous PIRK (CPIRK) methods; a subclass of these, the continuous MIRK methods, have been considered in ^[25].

For second order BVODE systems, $\underline y''(t) = \underline f(t, \underline y(t), \underline y'(t))$, we have solution and derivative information associated with the endpoints of the $i$th subinterval, i.e., $\underline y_{i}$, $\underline {y'}_{i}$, $\underline y_{i+1}$, $\underline {y'}_{i+1}$. It is possible to rewrite the well-known RKN methods so that the stages have explicit dependence on solution and derivative information at both endpoints of each subinterval; this gives the parameterized implicit RKN (PIRKN) methods. With the further restriction that each stage must depend only on previously computed stages, we get the mono-implicit RKN (MIRKN) methods; see, e.g., ^[26] and references within. By introducing weight polynomials, one can obtain the continuous PIRKN (CPIRKN) methods; a subclass of these, the continuous MIRKN methods, have been considered in ^[26].

The general form of a CPIRKN method is the same as that of a CRKN method, (3.8), and the definition of $\hat t_r$ is the same as in (3.2), but the $\underline {\hat y}_1, \dots, \underline {\hat y}_s$ values and the $\underline {\hat y}_1', \dots, \underline {\hat y}_s'$ values are defined differently. Instead of the expressions for $\underline {\hat y}_r$ and $\underline {\hat y}_r'$ given in (3.6) and (3.5), we have

and

A CPRIKN method of a desired order of accuracy is determined by requiring that the coefficients of the method, $\{c_r, v_r, w_r, v_r', b_r(\theta), {b}_r'(\theta), x_{r, j}, {x}_{rj}'\}_{r, j = 1}^{s}$, satisfy a set of continuous order conditions. We represent the CPIRKN methods using a tableau of the form,

It is possible to rewrite the Hermite interpolants, (3.12), (3.14), and (3.15), in terms of CPIRK and CPIRKN schemes. To see this, first recall that,

satisfy (2.8), (2.9), and (2.10). We can therefore use these expressions to substitute for $\underline y_{1, i+1}$, $\underline y_{2, i+1}$, $\underline {y}_{2, i+1}'$ in (3.12), (3.14), and (3.15). These equations can then be rewritten, respectively, as

where $\hat d_{1j}(\theta) = \hat b_{1j}d_{11}(\theta)$,

where $\tilde d_{2j}(\theta) = \hat {b}_{2j}' \bar d_{21}(\theta)$, and

where $\hat d_{2j}(\theta) = \hat b_{2j}d_{21}(\theta) + \hat {b}_{2j}'d_{23}(\theta)$. These new forms for the Hermite interpolants, i.e., (3.17), (3.18), (3.19), now fit into the standard forms for CPIRK and CPIRKN schemes.

3.4.   CPIRKN methods for SCIs

In this subsection, we discuss the CPIRKN methods that we will use in order to develop SCIs for mixed first and second order systems when $k\geq3$. (For $k = 1$ or 2, recall that the Hermite interpolants can be employed as the basis for the SCIs.)

On each subinterval, these CPIRKN methods will employ three types of stages:

(ⅰ) the first two stages will be the endpoint stages, (3.11), (3.13), which are easily computed from the corresponding collocation solution and derivative values,

(ⅱ) the next $k$ stages will be the collocation stages, (2.7), which are already available from the collocation computation, and,

(ⅲ) the rest of the stages will be MIRKN stages; each such stage can depend on stages of types (ⅰ) and (ⅱ), as well as on previously computed MIRKN stages. These new MIRKN stages will have the form,

where $\hat t_r = t_i + c_r h_i$, and where

and

4.   Derivation of CPIRKN methods for SCIs

4.1.   Collocation with $k=1$

When $k = 1$, the collocation mesh point values have errors that are $O(h^2)$. For this value of $k$, the continuous collocation solutions also have errors that are $O(h^2)$. We can therefore improve upon the continuous collocation solutions only in terms of achieving higher continuity, which we can do by employing the Hermite interpolants (3.12), (3.14), (3.15). Representing these interpolants in the form of a CPIRKN method, we get a CPIRKN method with the tableau:

where

This scheme, in addition to having the same order as the collocation approximations at the mesh points, provides approximations for $\underline y_1(t)$ and $\underline {y}_2'(t)$ that are $C^1$-continuous, and an approximation for $\underline y_2(t)$ that is $C^2$-continuous (whereas the corresponding collocation polynomials have one order of continuity less in each case).

4.2.   Collocation with $k=2$

In this case, the mesh point collocation values have errors that are $O(h^4)$. The continuous collocation solution approximating $\underline y_2(t)$ will also have an error that is $O(h^4)$ but the continuous collocation solutions approximating $\underline y_1(t)$ and $\underline {y}_2'(t)$ will have errors that are only $O(h^3)$. The use of the Hermite interpolant, (3.14), allows us to improve the continuity of the $\underline y_2(t)$ approximation and we can use the Hermite interpolants, (3.12) and (3.15), to improve the order of accuracy and the continuity of the continuous approximations for $\underline y_1(t)$ and $\underline {y}_2'(t)$.

Expressing the Hermite interpolants in CPIRKN form gives a method with the tableau:

where

4.3.   Collocation with $k=3$

In this case, the mesh point collocation values have errors that are O($h^6$), while the continuous collocation solution approximating $\underline y_2(t)$ has an error that is only O$(h^5$) and the errors for the continuous collocation solutions approximating $\underline y_1(t)$ and $\underline {y}_2'(t)$ are only O($h^4$). While the use of the Hermite interpolant (3.14) can provide an SCI of the appropriate order for $\underline y_2(t)$, we cannot use the Hermite interpolants, (3.12) and (3.15), to obtain SCIs for $\underline y_1(t)$ and $\underline {y}_2'(t)$. We therefore consider the derivation of a CPIRKN scheme that will provide SCIs for $\underline y_1(t)$ and $\underline {y}_2'(t)$, with $C^1$-continuity, for which the errors are $O(h^6)$. Using the terminology associated with RK methods, this means that the CPIRKN must be of fifth order. (The definition of order for RK methods says that a $p$th order RK method has a local error that is $O(h^{p+1})$.)

As indicated earlier in this paper, the CPIRKN scheme will employ the two endpoint stages and the three collocation stages. For these five stages we observe that, for $q = 3$, the coefficients satisfy,

and

where $X'$ is the $s$ by $s$ matrix whose $r, j$th element is $x_{rj}'$, $X$ is the $s$ by $s$ matrix whose $r, j$th element is $x_{rj}$, $\underline v' = [v'_1, \dots, v'_s]^T$, ${\underline v} = [v_1, \dots, v_s]^T$, ${\underline w} = [w_1, \dots, w_s]^T$, $\underline c^0$ is an $s$ vector of 1's, and $\underline c^j = [c^j_1, \dots, c^j_s]^T$, for $j > 0$. The above conditions, (4.6), (4.7), are called the stage order conditions up to order $q$ - see ^[26].

A fifth order CPIRKN scheme satisfying these stage order conditions will have to satisfy the following order conditions involving the $\underline {\bar b}(\theta)$ polynomials:

and the single condition,

where $\underline {\bar b}(\theta) = [\bar b_1(\theta), \dots, \bar b_s(\theta)]^T$. Since there are six independent order conditions involving the $\underline {\bar b}(\theta)$ polynomials, this CPIRKN scheme will require a total of six stages. Since the scheme will use the two endpoint stages and the three collocation stages, one extra MIRKN stage will be required.

We will require that the new sixth stage have the maximum possible stage order, which turns out to be stage order six - that is, (4.6), (4.7), with $q = 6$. This forces $c_6 = \frac{1}{2} \pm \frac{\sqrt{10}}{10}$, with $v_6, w_6$, and $v_6'$ left free. We choose, arbitrarily, $c_6 = \frac{1}{2} - \frac{\sqrt{10}}{10}$, $v_6 = c_6$, $w_6 = -c_6$, and $v_6' = c_6$. We can then solve the above order conditions to obtain the ${\underline {\bar b}}(\theta)$ weight polynomials.

As mentioned earlier, the approximation for ${\underline y}_2(t)$ can be determined using Hermite interpolation, with the Hermite interpolant converted to CPIRKN form, using the approach described in the previous section.

The resultant scheme has the tableau:

where $v_6 = \frac{1}{2}-\frac{\sqrt{10}}{10}$, $w_6 = -\frac{1}{2}+\frac{\sqrt{10}}{10}$, and the values for $x_{6r}$ and ${x}_{6r}'$ for $r = 1, \dots, 5$, the expressions for $b_{r}(\theta)$, for $r = 1, \dots, 5$, and the expressions for $\bar b_{r}(\theta)$, for $r = 1, \dots, 6$, are given by,

4.4.   Collocation with $k=4$

When $k = 4$, the mesh point collocation values have errors that are O($h^8$). On the other hand, the continuous collocation solution approximating $\underline y_2(t)$ has an error that is only O$(h^6$) and the errors for the continuous collocation solutions approximating $\underline y_1(t)$ and $\underline {y}_2'(t)$ are only O($h^5$). Since none of the Hermite interpolants can provide SCIs of the desired order, we will derive a CPIRKN scheme that will provide an SCI approximating $\underline y_2(t)$, with $C^2$-continuity, for which the error is $O(h^8)$, and SCIs approximating $\underline y_1(t)$ and $\underline {y}_2'(t)$, with $C^1$-continuity, for which the errors are $O(h^8)$.

Since we want this CPIRKN scheme to have an error that is O($h^8$), the CPIRKN must be of seventh order. This CPIRKN scheme will employ the two endpoint stages and the four collocation stages. These six stages have coefficients that satisfy the stage order conditions, (4.6) and (4.7), for $q = 4$. The additional MIRKN stages will be required to satisfy the stage order conditions for $q = 6$.

A CPIRKN scheme satisfying the above stage order conditions will have to satisfy the following order conditions in order to achieve seventh order:

the three conditions involving $\underline {\bar b}(\theta)$,

and the single condition involving $\underline {b}(\theta)$,

where $\underline {b}(\theta) = [{ b}_1(\theta), \dots, { b}_s(\theta)]^T$.

Since there are a total of ten order conditions involving the $\underline {\bar b}(\theta)$ polynomials, it would appear then that this CPIRKN scheme will require a total of ten stages. However, it turns out that we can satisfy the order condition,

through a careful choice of the available free coefficients - see ^[27] for details. This allows us to reduce the number of order conditions from ten to nine which implies that the scheme will require only nine stages.

Since this scheme will include the two endpoint stages and four collocation stages, we will need three extra MIRKN stages. Applying the stage order conditions, (4.6), (4.7), with $q = 6$, to the MIRKN stages leads to the specification of most of the coefficients but there are several free coefficients that remain; we arbitrarily set $v_7 = c_7$, $v_8 = c_8$, $v_9 = c_9$, $w_7 = w_8 = w_9 = 0$, $x_{76} = x_{86} = x_{87} = x_{96} = x_{97} = x_{98} = 0$. The final step is to solve the nine remaining order conditions involving the $\underline{\bar b}(\theta)$ polynomials and the seven order conditions involving the $\underline b(\theta)$ polynomials in order to obtain the $\underline{\bar b}(\theta)$ and $\underline b(\theta)$ polynomials. Since there are only seven order conditions involving the $\underline b(\theta)$ polynomials, we can set $b_8(\theta) \equiv b_9(\theta) \equiv 0$ and use the seven order conditions involving the $\underline {b}(\theta)$ polynomials to determine $b_1(\theta), \dots, b_7(\theta)$.

The resultant scheme has the tableau:

where $c_3 = \frac{1}{2}-{\frac {1}{70}}\, \sqrt {525+70\, \sqrt {30}}$, $c_4 = \frac{1}{2}-{\frac {1}{70}}\, \sqrt {525-70\, \sqrt {30}},$ $c_5 = \frac{1}{2}+{\frac {1}{70}}\, \sqrt {525-70\, \sqrt {30}},$ $c_6 = \frac{1}{2}+{\frac {1}{70}}\, \sqrt {525+70\, \sqrt {30}},$ and $c_7 = v_7 = \frac{1}{2}+1/14\, \sqrt {7}.$

The expressions for the non-zero components of the $X$ and $X'$ matrices and the coefficients of the $\underline b(\theta)$ and $\underline{\bar b}(\theta)$ polynomials are too complicated to present exactly. Here we represent these values in terms of 20 digit decimal numbers.

The non-zeros of the matrix $X$ are

The non-zeros of the matrix $X'$ are

The polynomials, $b_1(\theta), \dots, b_7(\theta)$, are

and the polynomials, ${\bar b}_1(\theta), \dots, {\bar b}_9(\theta)$, are

5.   Numerical results

Here we present results for two test problems. Results for other test problems are given in ^[27].

Problem 1: The first problem consists of two nonlinear ODEs, one of third order, the other of second order; it has the form,

with boundary conditions

We choose $\gamma = 3.0$.

In order to apply the methods developed in this paper, we will rewrite this problem as one first order ODE and two second order ODEs; letting $z_1(x) = f(x), z_2(x) = f'(x), z_3(x) = g(x)$, we get,

with the boundary conditions,

This system has no known exact solution; in order to estimate the errors for the approximate solutions for this test problem, we computed a high-precision numerical solution to this problem using COLSYS/COLNEW.

Problem 2: The second problem consists of a coupled system of second order BVODEs obtained from the application of the transverse-method-of-lines with a fixed time step backward Euler method to discretize the time-dependent terms of the PDE:

where

The time step, $\Delta t$, is chosen to give a system of 20 second-order BVODEs; the $i$th BVODE has the form

where $z_i(x) \approx z(x, t_i)$ and $t_i = i\cdot \Delta t$. The corresponding boundary conditions are

Because it is a second order system, the methods discussed in this paper can be directly applied to this problem.

For each test problem and for each solution component, the initial solution approximations provided to COLSYS/COLNEW are a straight line through the boundary conditions, for components with associated boundary conditions, and zero otherwise.

5.1.   Convergence rates and maximum errors

In order to numerically verify the orders of convergence of the mesh point collocation values, the continuous collocation solutions, and the SCIs, we compute collocation solutions using COLSYS/COLNEW on a sequence of fixed uniform meshes. We consider meshes for which, $N$, the number of subintervals, is equal to 8, 16, 32, 64, and 128. (This required setting control parameters of COLSYS/COLNEW so that the mesh refinement capability was disabled.)

We provide numerical results for the cases $k = 3$ and 4, which are the ones where CPIRKN schemes are used to obtain the SCIs. We compute collocation solutions using COLSYS/COLNEW and then augment these to obtain the SCIs, based on the CPIRKN schemes derived earlier in this paper.

We consider numerical experiments in which we compute the maximum error over all solution components for a given evaluation point. The Mesh Point Solution error is the maximum error over all mesh point values. The errors for the Collocation Solution and for the Superconvergent Interpolant are the maximum errors of these continuous approximations over 10000 uniformly distributed sample points. Table 1 gives results for the $k = 3$ case, while Table 2 gives results for the $k = 4$ case.

The results demonstrate that the derived schemes achieve experimentally observable rates of convergence that are consistent with those predicted by the theory. We note that, for $k = 3$ and for the finest mesh, the error of the SCI is about three orders of magnitude smaller than that of the continuous collocation solutions. Similarly, for $k = 4$ and again for the finest mesh, the error of the SCI is about four orders of magnitude smaller than that of the continuous collocation solutions, for the finest mesh.

As well, we see that the error of the SCIs, for the finest meshes, is about 2 times larger, for $k = 3$, and about 7 times larger, for $k = 4$, than the error of the corresponding mesh point collocation solution values. There are several ways to improve the accuracy of the SCIs. We document these in the final section, as items for future work.

5.2.   Relative costs of the collocation and SCI computations

Let $n$ be the number of differential equations and $k$ be the number of collocation points. A typical computation requires that several intermediate collocation solutions be computed on a sequence of meshes of sizes $N_1, N_2, \dots, N_\ell$, where $\ell$ is the total number of meshes that are required. The total cost associated with obtaining the collocation solution is $O(N_1q_1(kn)^3 + N_2q_2(kn)^3 + \dots + N_\ell q_\ell(kn)^3)$, where $q_j$ is the number of Newton iterations required to obtain convergence for the computation of the $j$th intermediate collocation solution on the mesh of $N_j$ subintervals. On the other hand, the cost associated with constructing an SCI in a post-processing step is $O(N_\ell nk^2)$. From these costs, we can see that the cost of constructing the SCI in a post-processing step, compared to the cost of computing the collocation solution, is quite small.

Tables 3 and 4 give timing results in seconds for Problem 2, with $\omega = 100$, for $k = 3$ and 4, and for a range of tolerances, $tol = 10^{-2}, \dots, 10^{-10}$. We report the cost of the computation of the collocation solution with COLSYS/COLNEW and the cost of the call to the SCI-SETUP routine, which computes the extra stages that are needed for the SCI. We also report the cost of evaluating the collocation solution, using the COLSYS/COLNEW APPSLN routine, and the cost of evaluating the SCI, using the SCI-INTERP routine, at a 1000 uniformly spaced points across the problem domain. From these tables we see that the SCI setup cost is less than 0.5% of the cost of computing the collocation solution and that the cost of evaluating the SCI is about 2 to 3 times the cost of evaluating the collocation solution, with both of these costs being negligible compared to the cost of computing the collocation solution itself.

6.   Summary, conclusions, and future work

In this paper we have shown how to generalize the approach considered in ^[13] to develop SCIs for Gaussian collocation solutions of mixed first and second order BVODE systems. These new methods can substantially improve the accuracy of the returned approximate solution leading to significant gains in the efficiency.

The approach involves augmenting the superconvergent collocation mesh point values with interpolants that have the same superconvergent order of accuracy. For the low order cases ($k = 1, 2$), Hermite interpolants can be employed on each subinterval to obtain the SCIs. For the $k = 3$ and 4 cases, a more general approach based on the use of a generalization of the continuous Runge-Kutta-Nyström methods can be employed. These methods can be used to obtain SCIs of orders 2, 4, 6, and 8, corresponding to $k = 1, 2, 3,$ and 4. Thus the range of orders of SCIs provided in this paper is comparable to the range of orders provided by COLSYS/COLNEW.

Our numerical results show that

(ⅰ) it is possible to obtain SCIs that are substantially more accurate than the corresponding continuous collocation solutions,

(ⅱ) the SCIs have experimentally observable orders of accuracy that are consistent with the expected orders, based on the theoretical framework that was used to derive them, and,

(ⅲ) the computational cost associated with setting up an SCI is a small fraction of the cost associated with computing the collocation solution upon which the SCI is based.

As mentioned earlier, additional work could be done in order to improve the accuracy of the CPIRKN methods upon which the SCIs are based. One possibility could be to optimize of the choice of the free coefficients that arise during the derivations to attempt to minimize the magnitude of the leading order term in the error of each scheme. Another possibility could involve the use of CPIRKN methods that, for each $k$, are one order of accuracy higher than those derived in this paper. A third possibility could involve expressing the weight polynomials in a Barycentric Lagrange form rather than in terms of a monomial basis. The paper ^[19] shows that the use of a Barycentric Lagrange representation can lead to less interference from round-off error.

A major task for future work is to extend the work reported in ^[1] in order to incorporate the schemes described in this paper to allow COLNEWSC to directly treat mixed first and second order BVODE systems, leading to improved user convenience, numerical solutions with higher continuity, and a more efficient computation.

The collection of BVODEs given in ^[4] involves differential equations with solution derivatives of orders 1 to 4, and COLSYS/COLNEW is able to directly treat mixed order systems across this range of orders. This suggests another direction for future work which would involve extending the approach considered in this paper to derive methods leading to SCIs for mixed order systems involving derivatives of orders 1 to 4. The primary challenge will be to extend the CPIRKN methods to develop methods for differential equations in which derivatives of orders 3 and 4 appear. With this goal in mind, we have developed general forms for generalizations of CPIRKN methods for systems of differential equations involving derivatives of orders 3 and 4 - see the Appendix of this paper.

Finally, as mentioned earlier in this paper, because the boundary value differential-algebraic system solver COLDAE ^[5] is based on Gaussian collocation and the BACOL family of error control solvers for 1D PDEs (see ^[17,32] and references within) makes use of Gaussian collocation for the spatial discretization of the PDEs, it may be worthwhile to investigate if it is possible to extend the ideas from the current paper to develop SCIs for the collocation solutions computed by these solvers.

Acknowledgments

(ⅰ) Funding sources: NSERC, Saint Mary's University.

(ⅱ) The authors wish to thank the referees and the editor for many helpful comments.

Conflict of interest

All authors declare that they have no conflicts of interest in this paper.

Appendix

A.   Generalized CPIRK schemes for third order ODEs

In this subsection we give the general form for a generalized CPIRK method that can be applied directly to an ODE of the form,

to provide SCIs for this problem class. The general form for these methods is,

where

B.   Generalized CPIRK schemes for fourth order ODEs

The methods presented in this subsection are generalizations of the methods presented in the previous subsection to allow for the direct treatment of BVODEs of the form,

The general form for these methods is,

where