On the Rayleigh-Taylor instability for the two coupled fluids

1.   Introduction

Arguably, one of the most impactful numerical methods for solving continuous-time optimal control problems (CTOCPs) in the 20th century has been direct pseudospectral (PS) methods, which can accurately reduce CTOCPs into optimization problems of standard forms that can be easily treated by using typical optimization methods. The key success of these methods lies in their ability to converge to sufficiently smooth solutions with exponential rates by using relatively coarse mesh grids. PS methods are considered to be "one of the biggest technologies for solving PDEs" that were largely developed about a half century ago since the pioneering works of Orszag ^[1] and Patterson Jr and Orszag ^[2]. They have been continuously refined and extended in later decades to solve many problems in various scientific areas that were only tractable by these techniques. Perhaps, one of the brightest moments in the course of their development appeared on March 3, 2007, when the International Space Station completed a 180-degree maneuver without using any propellant by tracking an attitude trajectory developed with PS optimal control theory, saving NASA one million USD ^[3]. The development of PS methods continues to be very active, and the progress in this research area has been remarkable in recent years; see the works in ^[4,5,6], for a few references on the solution of optimal control problems using various schemes of PS methods.

PS methods are closely related to the popular class of spectral methods, but they expand the solutions in terms of their grid point values by means of interpolation in lieu of global and, usually, orthogonal basis polynomials. Such a nodal representation is extremely useful in the sense that the solution values are immediately available at the collocation points once the full discretization is implemented, as the the governing equations are satisfied pointwise in the physical space, whereas modal representations require a further step of computing the modal approximation after calculating the coefficients of the basis function expansions ^[7]. This places PS methods at the front of highly accurate methods that are particularly easy to apply to equations with variable coefficients and nonlinearities ^[8]. Clear expositions of spectral and PS methods exhibiting a wide range of outlooks on the subject include, to mention a few, the books ^[9,10,11].

Two of the most common alternatives to direct spectral and PS methods are indirect methods and parameterization methods, which include control parameterization and state and control parameterization. As their names suggest, only the control variable is parameterized in a control parameterization method, and the differential equations are solved via numerical integration, while both state and control variables are parameterized in state and control parameterization methods, and the differential equations are converted into algebraic constraints. Typical examples of control parameterization methods include shooting methods and multiple shooting methods. On the other hand, an indirect method requires first the derivation of the often complicated first-order necessary optimality conditions, which include the adjoint equations, the control equations and the transversality conditions, before carrying out the numerical discretization. In fact, direct PS methods often present more computational efficiency and robustness over these classical methods on solving optimal control problems due to many reasons; to mention a few, (i) multiple shooting methods do not handle problems with singular arcs appropriately without a priori information on the structure of the trajectories ^[12]; (ii) direct shooting methods are often intensive computationally and associated with sensitivity issues; in particular, the ability to successfully use a direct shooting method declines as the number of variables grows large ^[13]; (iii) the boundary value problem resulting from the necessary conditions of optimality are extremely sensitive to initial guesses in indirect methods ^[14]; (iv) in contrast to indirect methods and direct shooting methods, direct PS methods do not require a priori knowledge of the active and inactive arcs for problems with inequality path constraints ^[15]; (v) the user does not have to be concerned with the adjoint variables or the switching structures to determine the optimal control in direct PS methods ^[16]; (vi) direct PS methods show much bigger convergence radii than either the indirect methods or direct shooting methods, as they are much less sensitive to the initial guesses ^[17]; (vii) direct PS methods are often memory minimizing when performed using orthogonal collocation approximations, as they result in finite-dimensional nonlinear programming problems (NLPs) with considerably lower dimensions compared to other competitive methods in the literature ^[18]; (viii) direct PS methods perform numerical differentiation through constant operators that can be stored for certain sets of collocation points, and invoked later when implementing the computational algorithms; and (ix) there is no need for the Lyapunov function to construct the asymptotically stabilizing control ^[6]. The reader may consult ^[19] for a further review on spectral and PS methods and their advantages against other classical methods.

A robust variant of PS methods is the class of integral PS (IPS) methods (that is, PS integration methods), which is closely related to PS methods, but it requires an initial step of reformulating the dynamical system equations in their integral form first before the collocation phase starts; thus, it avoids the degradation of precision that is often caused by numerical differentiation processes ^[20,21]. The integral reformulation can be performed by either a direct integration of the dynamical system equations if they have constant coefficients, or by approximating the solution's highest-order derivative involved in the problem by using a nodal finite series in terms of its grid point values and then solving for those grid point values before successively integrating back in a stable manner to obtain the sought solution grid point values. It has been shown in recent decades through a large number of publications that IPS methods often exhibit faster convergence rates, produce higher accuracy and are much less sensitive to the types of collocation points used than the usual PS methods; see ^[23,24] to mention a few. Historically, the spectral approximation of the integral form of differential equations was put forward in the 1960s by Clenshaw and Curtis ^[25] in the spectral space, and by El-Gendi ^[26] in the physical space.

Among the many classes of CTOCPs, infinite-horizon optimal control problems (IHOCs) and optimal control problems defined on sufficiently large intervals have attracted a lot of research interest due to the extent of their applications in economics, engineering, computer science, business and management science, bio-medicine, aerospace, energy, etc.; see ^[27,28,29] to mention a few. Some classical results on the existence of solutions for IHOCs can be found in ^{[30,31,32,33]}. One of the most general and well-known results on the existence of solutions to IHOCs was proved by Balder ^[34] using the notion of uniform integrability. Sufficient conditions for the existence of a finitely optimal solution for a class of nonlinear IHOCs were derived by Carlson ^[35] under minimal convexity and seminormality conditions. An existence and uniqueness theorem for a class of IHOCs was proved by Wang ^[36] under certain conditions. Existence and uniqueness results for a class of linear-quadratic, convex IHOCs in weighted Sobolev spaces for the state and weighted Lebesgue spaces for the control were obtained by Pickenhain ^[37]. A recent extension to the existence results of Balder ^[34] to the case in which the integral functional is understood as an improper integral was proved by Besov ^[38] using the notion of uniform boundedness of pieces of the objective functional that was proposed earlier by Dmitruk and Kuz'kina ^[39]. Aseev ^[40] derived some sufficient conditions for the existence and boundedness of optimal controls for a class of generally nonlinear IHOCs without necessarily having a bounded set of control constraints. Basco and Frankowska ^[41] obtained some existence and uniqueness results for weak solutions of the nonautonomous Hamilton–Jacobi–Bellman equation associated with a class of IHOCs for the class of lower semicontinuous functions vanishing at infinity and under certain conditions of controllability. The most important and well-known necessary conditions of optimality were first derived by Halkin ^[42]; see also [33,Theorem 2.3].

While many direct PS methods appeared in the literature for solving finite-horizon optimal control problems (FHOCs), we could only find a few works on deterministic IHOCs governed by integer-order differential equations using this class of methods. In particular, we recognize the Legendre-Gauss (LG) and Legendre-Gauss-Radau (LGR) PS methods in ^[43,44,45], the flipped Legendre-Gauss-Radau PS method (FRPM) in ^[46] and the transformed Legendre spectral method in ^[47]. Although Legendre polynomials are commonly used in PS methods designed to solve IHOCs, we shall explore in our work the possibility of whether we could achieve better accuracy and convergence rates using Gegenbauer polynomials (that is, ultraspherical polynomials). There are a number of reasons that prompt us to consider this family of polynomials as a viable alternative to perform discretizations of IHOCs. (i) First, observe that Gegenbauer polynomials include both Chebyshev and Legendre polynomials as part of its bigger family, so all theoretical and experimental results on Gegenbauer polynomials directly apply to Chebyshev and Legendre polynomials by definition. (ii) Being a part of Gegenbauer polynomials allows us to apply any of the Chebyshev and Legendre polynomials with a single selection of the Gegenbauer parameter (index) $\alpha$, as one can simply set $\alpha = 0$ or $1/2$ in your code, thus giving us more flexibility. (iii) Gegenbauer polynomials are very useful in eliminating the Gibbs phenomenon and recovering the spectral accuracy up to the discontinuity points ^[48,49,50]. (iv) One measure for assessing the quality of spectral and PS methods in numerical discretizations is concerned with how large is the number of terms that are required in a spectral/PS expansion to achieve a certain level of accuracy. Of course, the smaller the number of terms, the more efficient the method is in terms of speed and computational complexity. In applications like CTOCPs, this property leads to optimization problems on the small scale, which can be solved very quickly with reduced computational work at a concrete level by using modern optimization software ^[18,21,23]. Now, with this being mentioned, it is important to realize that Chebyshev and Legendre polynomials are usually optimal for large spectral/PS expansions under the Chebyshev and Euclidean norms, respectively, but they are not necessarily optimal for a small/medium range; this is an observation that was proven numerically in a number of papers for certain polynomial and rational interpolations and collocations; see ^[8,20,51] to mention a few, which give us another reason to apply Gegenbauer polynomials as proper basis polynomials that may provide faster convergence rates. (v) A stability analysis conducted in ^[24] and grounded in the Lebesgue constant for polynomial interpolations in Lagrange-basis form based on flipped Gegenbauer-Gauss-Radau (FGGR) points showed that the Lebesgue constant is not minimal for Chebyshev polynomials, but, rather, was minimal for Gegenbauer polynomials associated with negative $\alpha$ values. This analysis proved, with no doubt, that some Gegenbauer polynomials with negative $\alpha$ values could be more plausible for employment in basis-form polynomial interpolation/collocation for short-medium-range mesh grid sizes. This observation is consistent with the earlier work of Light ^[52], who proved in the late 1970's that the Chebyshev and Legendre projection operators cannot be minimal as the norms of Gegenbauer projection operators increase monotonically with $\alpha$ for small expansions.

In light of the above arguments, we were motivated in this work to develop a novel direct IPS method for solving IHOCs using Gegenbauer polynomials, and to study its convergence. To this end, we derive some accurate and numerically stable Gegenbauer-Gauss-Radau (GGR)-based rational interpolation formulas and describe two useful computational algorithms for constructing their barycentric weights. We also show how to derive the associated quadrature formulas required for numerical integration in time. We shall then use these numerical instruments to approximate the optimal state and control variables after transforming the IHOC into an FHOC in integral form (FHOCI) by means of certain parametric maps and rational collocation. During the course of our paper presentation, we shall also try to investigate a number of interesting relevant questions to our work. For instance, which parametric map is more suitable for GGR-based rational collocations? How should we choose the Gegenbauer parameter values to carry out collocations in practice? A "poor" choice of $\alpha$ can largely ruin the accuracy of the numerical scheme, while a "good" choice can, in many cases, furnish superb approximations with higher accuracy than those enjoyed by Chebyshev and Legendre polynomials for sufficiently smooth functions by using relatively coarse mesh grids. Do PS and IPS methods based on Chebyshev, Legendre and Gegenbauer polynomials generally converge to the solutions of IHOCs for large mesh sizes? If they do not, then what are the causes? Through rigorous stability, error and convergence analyses, we shall prove that such methods often converge at an exponential rate to near exact solutions by using relatively small mesh grids, but they usually diverge for fine meshes under certain parametric maps.

The proposed method inherits all of the merits of the direct IPS method stated above; in addition, the current paper presents the following novel contributions to the literature of IHOCPs: (i) For a coarse collocation mesh grid size under certain parametric maps, we show that direct IPS methods often converge exponentially to near-optimal solutions of IHOCs when leveraged with a fast and accurate NLP solver. (ii) Despite the exceedingly accurate approximations achieved at coarse meshes, we prove that direct IPS methods based on Gegenbauer polynomials and certain parametric maps usually diverge as the number of collocation points grow large if the computations are carried out using floating-point arithmetic and the discretizations use a single mesh grid, regardless of whether they are of Gauss/Gauss-Radau (GR) type or equally spaced, which is a result that can be readily extended to include classical Jacobi polynomials in general. (iii) We derive some novel formulas for the construction of GGR-based rational interpolation weights, leading to a more numerically stable interpolation procedure, which we prefer to call "the switched rational (SR) interpolation, " due to the switching nature of the novel algorithm we developed to accurately calculate the barycentric weights. (iv) More novel formulas for the construction of GGR-based integration and differentiation matrices in barycentric-trigonometric forms are derived. (v) We show that Legendre polynomials are (near) optimal basis polynomials for GGR-based SR collocations over coarse meshes. (vi) We show that Gegenbauer polynomials associated with certain nonnegative index values are more apt for GGR-based SR interpolations over fine meshes. (vii) We prove that parametric logarithmic maps are more apt for the domain transformation of IHOCs than parametric algebraic maps for GGR-based SR collocations. (viii) To the best of our knowledge, this paper introduces the first direct IPS method for solving IHOCs using Gegenbauer polynomials and algebraic-logarithmic parametric maps.

The rest of the article is organized as follows. Sections 2 and 3 describe the IHOC under study and its transformation into an FHOCI via various parametric maps. Section 4 presents the discretization scheme of the FHOCI passing through the construction of the needed barycentric rational interpolants and their quadratures, and it closes with a setup of the IPS rational collocation at the GGR points and the optimality necessary conditions of the obtained NLP in Sections 4.1–4.3. Section 4.1 includes an analysis of the stability and sensitivity of the GGR-based rational interpolation/collocation developed in this paper. Rigorous error and convergence analyses are conducted in Section 5. Some divergence results of typical IPS collocation schemes of the FHOCI for fine meshes of the Gauss type using certain parametric maps are derived in Section 5. Simulation results are shown in Section 6, followed by some conclusions and future works in Section 7. The derivation of the barycentric rational formulas necessary for constructing the GGR-based differentiation matrix is shown in Appendix B. Two easy-to-implement pseudocodes of computational algorithms for computing the barycentric weights of our new SR interpolation method are described in Appendix C.

2.   Problem statement

Consider the following nonlinear, autonomous control system of ordinary differential equations:

which is subject to the following system of initial conditions:

where $\mathit{\boldsymbol{x}}(t) = \left[x_{1}(t), x_{2}(t), \ldots, x_{n_{x}}(t)\right]^{\top} \in \mathbb{{R}}^{n_{x}}, \mathit{\boldsymbol{u}}(t) = \left[u_{1}(t), u_{2}(t), \ldots, u_{n_{u}}(t)\right]^{\top} \in \mathbb{{R}}^{n_{u}}$, $\mathit{\boldsymbol{x}}_0 = [x_{1, 0}, x_{2, 0}, \ldots, x_{n_x, 0}]^{\top} \in \mathbb{{R}}^{n_{x}}$ is a constant specified vector, and $\mathit{\boldsymbol{f}}: \mathbb{{R}}^{n_{x}} \times \mathbb{{R}}^{n_{u}} \rightarrow \mathbb{{R}}^{n_{x}}: \mathit{\boldsymbol{f}} = [f_1, f_2, \ldots, f_{n_x}]^{\top}$. The problem is to find the optimal control $\mathit{\boldsymbol{u}}(t)$ and the corresponding state trajectory $\mathit{\boldsymbol{x}}(t)$ on the semi-infinite-domain $[0, \infty)$ that satisfy Eqs (2.1) and (2.2) while minimizing the functional

where $g: \mathbb{{R}}^{n_{x}} \times \mathbb{{R}}^{n_{u}} \rightarrow \mathbb{{R}}$. We assume that $\mathit{\boldsymbol{f}}$ and $g$ are generally nonlinear, continuously differentiable functions with respect to their arguments, and that the nonlinear IHOC (2.1)–(2.3) has a unique solution. In the rest of the article, for any row/column vector $\mathit{\boldsymbol{y}} = (y_i)_{1 \le i \le n}$ with $y_i \in \mathbb{{R}}\, \forall i$ and real-valued function $h: \mathit{\boldsymbol{\Omega}} \subseteq \mathbb{R} \to \mathbb{R}$, the notation $h(\mathit{\boldsymbol{y}})$ stands for a vector of the same size and structure of $\mathit{\boldsymbol{y}}$ such that $h(y_i)$ is the $i$th element of $h(\mathit{\boldsymbol{y}})$. Moreover, by $\mathit{\boldsymbol{h}}(\mathit{\boldsymbol{y}})$, we mean $[h_1(\mathit{\boldsymbol{y}}), \ldots, h_m(\mathit{\boldsymbol{y}})]^{\top}$ for any $m$-dimensional column vector function $\mathit{\boldsymbol{h}}$, with the realization that the definition of each array $h_i(\mathit{\boldsymbol{y}})$ follows the former notation rule for each $i$.

3.   Transformation of the IHOC into an FHOC

Given a differentiable, strictly monotonic mapping $T: [0, \infty) \to [-1, 1)$ defined by $T(\tau) = t$, one can transform the IHOC (2.1)–(2.3) into the following FHOC:

subject to

where $\mathit{\boldsymbol{f}}\left(\tilde{\mathit{\boldsymbol{x}}}(\tau), \tilde{\mathit{\boldsymbol{u}}}(\tau)\right) = \left[f_{1}\left(\tilde{\mathit{\boldsymbol{x}}}(\tau), \tilde{\mathit{\boldsymbol{u}}}(\tau)\right), \ldots, f_{n_{x}}\left(\tilde{\mathit{\boldsymbol{x}}}(\tau), \tilde{\mathit{\boldsymbol{u}}}(\tau)\right)\right]^{\top}$ and $\tilde{\boldsymbol{\eta}}(\tau) = \boldsymbol{\eta}\left(T(\tau)\right)\, \forall \boldsymbol{\eta} \in \lbrace \mathit{\boldsymbol{x}}, \mathit{\boldsymbol{u}} \rbrace$. To take advantage later of the well-conditioning of numerical integration operators during the collocation phase, we rewrite Eq (3.1b) in its integral formulation, as follows:

We refer to the FHOC described by Eqs (3.1a), (3.1c) and (3.1d) by the FHOCI. A wide variety of defining formulas exist for the mapping $T$. Five common defining formulas of such a mapping that occurred in the literature are as follows:

where $L \in \mathbb{{R}}^+$ is a scaling parameter that can stretch the image of the interval $[-1, 1]$ in the codomain $[0, \infty)$ as desired. An "optimal" choice of $L$ value can significantly improve the quality of the discrete approximations, as we shall demonstrate later in Section 6. We refer to $L$ as "the map scaling parameter." The parametric maps $T_{i, L}, i = 1, 2$ are often referred to as "algebraic" and "logarithmic" maps in the literature, respectively. Notice that the maps defined by Eqs (3.2c) and (3.2d) are special cases of the parametric maps $T_{i, L}, i = 1, 2$; in particular, $T_{1, 1} = \varphi_a, T_{2, 1} = \varphi_b$ and $T_{2, 2} = \varphi_c$. Since the value of either parametric map varies when $\alpha$ varies for arguments of type GGR points, it is more convenient in this work to denote them by $T_{i, L}^{(\alpha)}, i = 1, 2$ to emphasize this fact. In particular, we prefer to generalize Eqs (3.2a) and (3.2b) to implicitly allow for maps of a wide spectrum of $\alpha$ values, as follows:

Mesh-like surfaces of both parametric mappings are shown in Figures 1 and 2 for several values of $n, L$ and $\alpha$. Both figures show that the parametric mappings (i) increase monotonically for decreasing values of $\alpha$ while holding $n$ and $L$ fixed, (ii) increase monotonically for increasing values of $L$ while holding $n$ and $\alpha$ fixed and (iii) increase monotonically for increasing values of $n$ while holding $L$ and $\alpha$ fixed. Moreover, near $\tau = 1$, the rate of increase of ${T_{1, L}^{(\alpha)}}$ with respect to any of the arguments $n, L$ and $\alpha$ while holding the others fixed is much larger than that of ${T_{2, L}^{(\alpha)}}$, which grows very slowly. Loosely put, the stretching of the mesh grid near $\tau = 1$ is stronger for $T_{1, L}^{(\alpha)}$ than $T_{2, L}^{(\alpha)}$.

4.   Numerical discretization of the FHOCI

In this section, we provide a description of the proposed numerical discretization of the FHOCI using an IPS method based on Gegenbauer polynomials and GGR points.

4.1.   Barycentric rational interpolation at the GGR points

Let $\mathbb{{Z}}^{+}$ be the set of positive integers, $G_{n}^{(\alpha)}(\tau)$ be the $n$th-degree Gegenbauer polynomial with $\alpha > -1/2$ and $G_{n}^{(\alpha)}(1) = 1, \; \forall n \in \mathbb{{Z}}_{0}^{+} = \mathbb{{Z}}^{+}\cup \lbrace 0 \rbrace$, and $\mathbb{{S}}_{n} = \left\{\tau_{k}: \mathcal{{G}}_{n+1}^{(\alpha)}\left(\tau_k\right) = 0\text{ for }k = 0, \ldots, n, \text{ and } -1 = \tau_0 < \tau_1 < \ldots < \tau_n\right\}$ be the set of GGR nodes, where $\mathcal{{G}}_{n+1}^{(\alpha)}(\tau) = G_{n}^{(\alpha)}(\tau) + G_{n+1}^{(\alpha)}(\tau)$. The orthonormal Gegenbauer basis polynomials are defined by $\phi_{j}^{(\alpha)}(\tau) = G_{j}^{(\alpha)}(\tau)/\sqrt{\lambda_{j}}$, where

They satisfy the discrete orthonormality relation

where $\varpi_j, j = 0, 1, \ldots, n$, are the corresponding Christoffel numbers of the GGR quadrature formula on the interval $[-1, 1]$ defined by

with

Given a set of $n + 1$ data points $\left\{(\tau_i, f_i)\right\}_{i = 0}^n$, the Gegenbauer polynomial interpolant $P_{n}f$ in Lagrange form is defined by

where $\mathcal{{L}}_{n, i}$ are the Lagrange polynomials given by

$P_{n}f$ can be evaluated fast and more stably by evaluating Lagrange polynomials through the "true" barycentric formula

which brings into play the barycentric weights $\xi_{i}, i = 0, \ldots, n$, given by

An interpolation in Lagrange form with Lagrange polynomials defined by Eq (4.6) is often referred to as "a barycentric rational interpolation." The barycentric weights associated with the GGR points can be expressed explicitly in terms of the corresponding Christoffel numbers through the following theorem.

Theorem 4.1. The barycentric weights for the GGR points are given by

Proof. Let $P_n^{(\alpha, \beta)}(\tau)$ be the Jacobi polynomial of degree $n$ and associated with the parameters $\alpha, \beta > -1$, as normalized by Szegö ^[56]. Through [56,Eq (4.5.4)] and [18,Eq (A.1)], we have

Therefore, $(1+\tau_i) P_{n+1}^{(\alpha-1/2, \, \alpha+1/2)}(\tau_i) = \mathcal{{G}}_{n+1}^{(\alpha)}(\tau_i)\, \forall i$, and Eqs (4.8a) and (4.8b) can be derived from [57,Theorem 3.6] by replacing both $\alpha$ and $\beta$ with $\alpha-1/2$. Recall that the GGR points cluster near $\pm 1$ as $n \rightarrow \infty$, so Eq (4.8b) may suffer from cancellation errors for values of $\tau_i$ sufficiently close to $1$. The following theorem provides two alternative trigonometric forms of Eq (4.8b).

Theorem 4.2. The barycentric weights corresponding to the interior GGR points are given by

Proof. Through the change of variables $\tau = \cos{\theta}$ and the double-angle rule for the cosine function, we find that

from which Eq (4.10a) is derived. Equation (4.10b) is established by realizing that

which completes the proof.

We refer to Eqs (4.10a) and (4.10b) by the trigonometric-barycentric weights. Equation (4.8b) is faster to compute and requires a smaller number of arithmetic operations compared with Eqs (4.10a) and (4.10b), but the latter two formulas may possibly produce smaller errors near $\tau = 1$, as we observed through numerical experiments. This suggests that it is better to perform Eq (4.8b) for all values of $\tau$, except when $\tau$ is sufficiently close to $1$, where we switch to the other trigonometric forms. To this end, we introduce "a switching parameter, " $0 < \varepsilon \ll 1$, at which the interchange of formulas is performed. The crossover value of $\varepsilon$, where it becomes more accurate to use the trigonometric form, will depend on the implementation; a prescription of this strategy is outlined in Algorithms 1 and 2. We refer to the explicit formulas used in Algorithms 1 and 2 to compute the barycentric weights by the "first and second switching formulas" of the barycentric weights for the GGR points, respectively. We also denote the errors in computing the barycentric weights using Eq (4.8b), Algorithm 1 and Algorithm 2, by $E_{i, n}^{(\alpha)}$, for $i = 1, 2$ and $3$, respectively. Figures 3 and 4 show comparisons between the three strategies, in terms of error, for certain value ranges of $\varepsilon, n$ and $\alpha$. While Algorithm 1 does not look promising against the usual Eq (4.8b) relative to the given input data, as clearly shown in Figure 3, Figure 4 manifests that Algorithm 2 is more numerically stable near $\varepsilon = 0.1$ and provides better approximations. We shall therefore use the latter algorithm with $\varepsilon = 0.1$ for the computation of the barycentric weights, and we refer to the rational interpolation with barycentric weights obtained through the second switching formulas as the "SR interpolation."

Remark 4.1. The near optimal value $\varepsilon = 0.1$ is based on extensive numerical experiments for candidate $\varepsilon$ values between $0$ and $1$ with no theoretical evidence.

Stability and sensitivity analyses of GGR-based SR interpolation/collocation

A valuable device for measuring the quality and numerical stability of polynomial interpolations is the Lebesgue constant, as it provides a measure of how close the interpolant of a function is to the best polynomial approximant of the function. The Lebesgue constant is also very useful in assessing the quality of approximate solutions obtained through collocation (that is, the collocated solutions), as their accuracy is related to the rate at which the Lebesgue constant increases.

Let $\left\| f \right\|_{ \mathbb{{S}}} = \sup{\lbrace |f(x)|:x \in \mathbb{{S}}\rbrace}$ be the uniform norm (or supremum norm) of a real-valued, bounded function $f$ defined on a set $\mathbb{{S}} \subseteq \mathbb{{R}}$. Suppose that $y_i, i = 0, \ldots, n$ and $y_{c, n}(\tau)$ denote the exact solution values at the GGR points $\tau_{i}, i = 0, \ldots, n,$ and the corresponding collocated solution, respectively. Let $p_n^*y$ and $p_{n}y$ be the best polynomial approximation to the exact solution $y$ on $[-1, 1)$ and the Lagrange interpolating polynomial of degree at most $n$ that interpolate the data set $\left\{ {\left({\tau_{i}, y_i} \right)} \right\}_{i = 0}^n$, respectively. Through the uniqueness of Lagrange interpolation, one can easily show that

where $\Lambda _n^{(\alpha)} = \mathop {\max }\nolimits_{-1 \le x < 1} \sum\limits_{i = 0}^n {\left| {\mathcal{{L}}_{n, i}(\tau)} \right|}$ denotes the Lebesgue constant associated with GGR-based SR interpolation. Therefore,

where $\delta y_{c, n}(\tau)$ is the difference between the Lagrange interpolating polynomial and the collocated solution. When ${\left\| {\delta {y_{c, n}}} \right\|_{[- 1, 1)}} \approx 0, \Lambda _n^{(\alpha)}$ roughly bounds the collocation error, that is, it nearly quantifies how much larger the collocation error ${\left\| {y - {y_{c, n}}} \right\|_{[- 1, 1)}}$ is compared to the smallest possible error, ${\left\| {y - p_n^*y} \right\|_{[- 1, 1)}}$, in the worst case. In this case, it is obvious from Eq (4.14) that, the smaller the Lebesgue constant, the better is the predicted collocated solution in the uniform norm. In other words, the collocation error is about at most a factor $1+\Lambda _n^{(\alpha)}$ worse than the best possible polynomial approximation. One can also clearly see that $\Lambda _n^{(\alpha)}$ depends on the location of the collocation points $\tau_i, \, i = 0, \ldots, n$, but not on the solution values $y_i, \, i = 0, \ldots, n$. Since the positions of the GGR points change as $\alpha$ and $n$ vary, we are interested in learning the apt choices of $\alpha$ that makes $\Lambda _n^{(\alpha)}$ as small as possible while holding $n$ fixed. This can provide some useful insight into how we should select the candidate range of the $\alpha$ values often used for collocations based on GGR points. In ^[24], our findings uncovered that $\Lambda _n^{(\alpha)}$ for FGGR-based polynomial interpolation in Lagrange basis form blows up as $\alpha \to -0.5$ and monotonically increases for increasing positive values of $\alpha$; see [24,Eqs (3.7) and (3.8)]. Moreover, it was noticed that $\Lambda _n^{(\alpha)}$ does not decrease monotonically for increasing negative values of $\alpha$, indicating that $\Lambda _n^{(\alpha)}$ is not minimal for Chebyshev polynomials, but, rather, attains its smallest value for the Gegenbauer polynomials associated with some negative values of $\alpha$; see [24,Figures 1 and 3]. It was roughly estimated in many earlier works through theoretical and numerical evidence that reasonably good $\alpha$ values for polynomial interpolation in basis form typically belong to the Gegenbauer collocation interval of choice, $I_{\varepsilon, r}^G$, defined by $I_{\varepsilon, r}^G = \left\{ {\alpha | - 1/2 + \delta \le \alpha \le r, 0 < \delta \ll 1, r \in[1,2]} \right\}$ for reasons pertaining to the stability and accuracy of numerical schemes employing Gegenbauer polynomials as basis polynomials; see ^{[7,20,21,24,58]} and the references therein.

On the other hand, it was discovered in a number of works that Lebesgue constants for rational interpolation at equally spaced nodes are much smaller than those associated with classical polynomial interpolation ^{[59,60,61,62,63]}. Moreover, the Lebesgue constant of Berrut's rational interpolation ^[64] at equidistant nodes is smaller than the Lebesgue constant for polynomial interpolation at Chebyshev nodes ^[65]. Figure 5 shows the surface of the Lebesgue constant for GGR-based rational interpolation characterized by Eqs (4.4) and (4.6) with barycentric weights obtained through the second switching formulas. The surface was constructed through least-squares approximation, and it is shown together with some of its cross sections with the vertical planes $\alpha = -0.49, -0.2, -0.1, 0, 0.5, 1, 1.5$. A number of remarks deserve to be made at this point. (i) First, notice the small size of the Lebesgue constant values for GGR-based rational interpolation compared with its values for the FGGR-based polynomial interpolation in basis form; see [24,Figures 1]. (ii) $\Lambda _n^{(\alpha)}$ does not blow up as $\alpha \to -0.5$, but, rather, remains bounded. (iii) The associated Lebesgue constant grows logarithmically in the number of collocation nodes. (iv) It is also interesting to see how the Lebesgue constant drops monotonically as $\alpha$ increases while holding $n$ fixed. This suggests that Legendre polynomials are generally more suited for GGR-based SR interpolations than Chebyshev polynomials. In fact, we can also observe from Figure 5 that Gegenbauer polynomials with increasing $\alpha$ values are associated with smaller Lebesgue constant values. This suggests that Gegenbauer polynomials with $\alpha > 1/2$ may also be more plausible to employ in SR interpolation for short/medium ranges of $n$ values. However, the work of Elgindy and Smith-Miles ^[20] manifests that Gegenbauer quadratures "may become sensitive to round-off errors for positive and large values of the parameter $\alpha$ due to the narrowing effect of the Gegenbauer weight function, " which drives the quadratures to become more extrapolatory with greater uncertainty in integral approximations; thus, the collocation is subject to a higher risk of producing meaningless results. In other words, ${\left\| {\delta {y_{c, n}}} \right\|_{[- 1, 1)}}$ may become large and the collocation error grows accordingly. It was also observed in ^[20] that the weight function ceases to exist near the boundaries $\tau = \pm 1$, and that its support is nonzero only on a subinterval centered at $\tau = 0$ for increasing values of $\alpha > 1$. If we refer to GGR-based collocations employing SR interpolations by the "SR collocations, " then this analysis suggests that, for a relatively large collocation mesh size, SR collocations are expected to produce higher-order approximations for nonnegative $\alpha$ values with apparently optimal $\alpha$ values within/near the "Gegenbauer SR collocation interval of choice (SRCIC)" $\Upsilon_{ \alpha_{c}^-, \alpha_{c}^+}^G = [\alpha_{c}^-, \alpha_{c}^+]: \alpha_{c}^- \approx 1/2, 1/2 \le \alpha_{c}^+ \le 1$; in addition, Gegenbauer polynomials with positive and large $\alpha$ values are generally not apt for SR interpolation/collocation. For small mesh sizes, however, there is no rule of thumb as to how should we select $\alpha$, since all Lebesgue constant curves converge to the same limit as $n \to 1$. The analysis in this section requires the assumption that the problem under study is well-conditioned. For sensitive problems, the interval of choice $\Upsilon_{ \alpha_{c}^-, \alpha_{c}^+}^G$ may change depending on the sources of sensitivity. In Section 5, we shall show that proper collocations of the FHOCI using any of the maps described by Eqs (3.3a) and (3.3b) entail shifting the right boundary, $\alpha_{c}^+$, of $\Upsilon_{ \alpha_{c}^-, \alpha_{c}^+}^G$ rightward as the mesh size grows large to reduce the divergence rate of the collocated solutions from the exact solutions.

From another perspective, it is interesting to mention that the Lebesgue constant is also a useful instrument in observing how the collocated solutions change as the input data are varied. By closely following the convention in ^[24], suppose that $\tilde y_{i}, i = 0, \ldots, n,$ and $\tilde y_{c, n}(\tau)$ are the perturbed solution values due to round-off or input data errors and the perturbed collocated solution, respectively. Moreover, assume that $\tilde p_{n}y(\tau)$ is the Lagrange interpolating polynomial of degree at most $n$ that interpolates the data set $\left\{ {\left({\tau_{i}, \tilde y_{i}} \right)} \right\}_{i = 0}^n$. Then, we have

where $\delta \tilde y_{c, n}(\tau)$ is the difference between the perturbed Lagrange interpolating polynomial and the perturbed collocated solution. When ${\left\| {\delta y_{c, n} - \delta \tilde y_{c, n}} \right\|_{[-1, 1)} }$ is relatively small, $\Lambda _n^{(\alpha)}$ nearly quantifies the larger size of the perturbation error of the collocated solution, ${\left\| {y_{c, n} - \tilde y_{c, n}} \right\|_{[-1, 1)} }$, compared to the maximum possible perturbation error of the solution at the collocation points, $\mathop {{\text{max}}}\limits_{0 \le i \le n} \left| {y_i - \tilde y_i} \right|$, or to the maximum solution perturbation error, ${\left\| {y - \tilde y} \right\|_{[- 1, 1)}}$, in the worst case.

4.2.   Barycentric GGR-based integration matrix (GRIM) and quadratures

Consider a real-valued function $f$ defined on the interval $[-1, 1]$ and its GGR-Based SR interpolation given by Eqs (4.4) and (4.6) and the second switching formulas of the barycentric weights. Following the work in ^[66], the formulas needed to construct the nonzero rows of the barycentric GRIM can be derived by integrating Eq (4.4) on the successive intervals $\left[-1, \tau_{j}\right], j = 1, \ldots, n$, to obtain

where $f_i = f(\tau_i)\, \forall i$. With the change of variable

we can rewrite Eq (4.16) as

where

Since the polynomials $\mathcal{{L}}_{n, i}\left(t; -1, \tau_{j}\right), i = 0, \ldots, n$, are of degree $n$, the integrals $\int_{-1}^{1} \mathcal{{L}}_{n, i}(t; -1, \tau_{j}) dt$ can be computed exactly using an $N = \left\lceil(n+1)/2\right\rceil$-point LG quadrature, where $\left\lceil. \right\rceil$ denotes the ceiling function. Let $\lbrace \bar{\tau}_{k}, \bar{\varpi}_{k}\rbrace_{k = 0}^{N}$ be the set of LG quadrature nodes and weights, respectively, where

and $L'_{N+1}$ denotes the derivative of the $(N+1)$st-degree Legendre polynomial $L_{N+1}$. Then,

Hence, Eqs (4.18) and (4.21) yield the GGR-SR quadrature rule

where $\mathit{\boldsymbol{f}} = \left[f_{0}, f_{1}, \ldots, f_{n}\right]^{\top}$ and $q_{j, i}, i, j = 0, \ldots, n$, are the elements of the first-order barycentric GRIM $\mathbf{{Q}}$ given by

We denote the $j$th row of $\mathbf{{Q}}$ by $\mathbf{{Q}}_{j}\, \forall j$. The derivation of the formulas required to construct the GGR-based differentiation matrix (GRDM) in barycentric form is described in Appendix B.

4.3.   IPS rational collocation of the FHOC at the GGR points

Let $\mathit{\boldsymbol{\tau}}_{n} = \left[\tau_{0}, \tau_{1}, \ldots, \tau_{n}\right]^{\top}, \tilde{\mathit{\boldsymbol{x}}}_i = \tilde{\mathit{\boldsymbol{x}}}\left(\tau_{i}\right), \tilde{\mathit{\boldsymbol{u}}}_i = \tilde{\mathit{\boldsymbol{u}}}\left(\tau_{i}\right), \mathit{\boldsymbol{f}}_i = \mathit{\boldsymbol{f}}(\tilde{\mathit{\boldsymbol{x}}}_i, \tilde{\mathit{\boldsymbol{u}}}_i), f_{i, j} = f_{i}(\tilde{\mathit{\boldsymbol{x}}}_j, \tilde{\mathit{\boldsymbol{u}}}_j)\, \forall i, j$, and

Then, collocating Eq (3.1d) at the GGR nodes yields

where

$\mathbf{{F}} = \left[\mathbf{{F}}_{1, n}, \ldots, \mathbf{{F}}_{n_{x}, n}\right], \mathbf{{F}}_{i, n} = \left[f_{i, 0}, \ldots, f_{i, n}\right]^{\top}\, \forall i, \mathit{\boldsymbol{\mathit{1}}}_{n}$ is the all-ones column vector of size $n$ and $[., .], \text{"vec"}, \circ$ and $\otimes$ denote the horizontal matrix concatenation, the vectorization of a matrix, the Hadamard product, and the Kronecker product, respectively. Let $\tau_{n + 1} = 1$ and define ${{\mathbf{{Q}}}_{n + 1}} = \left({{q_{n+1, i}}} \right)_{0 \le i \le n}: {q_{n+1, i}} = \sum_{k = 0}^{N} \bar{\varpi}_{k} \mathcal{{L}}_{n, i}\left(\bar{\tau}_{k}\right)\, \forall i$; then, the discrete cost functional $J$ can be approximated numerically using the LG quadrature, as follows:

where $\mathit{\boldsymbol{g}} = \left[g_{0}, \ldots, g_{n}\right]^{\top}: g_{i} = g\left(\tilde{\mathit{\boldsymbol{x}}}_i, \tilde{\mathit{\boldsymbol{u}}}_i\right)\, \forall i$. Hence, the FHOCI (3.1a), (3.1c) and (3.1d) can now be converted into a NLP in which the goal is to minimize the discrete cost functional given by Eq (4.27), as subject to the nonlinear system of equations given by Eq (4.25). If we define the image of the collocation points set $\mathbb{{S}}_n$ under the transformation $T$ by $\mathbb{{S}}_n^T = \{t_k: t_{k} = T\left(\tau_{k}\right), \; k = 0, \ldots, n\}$, and denote $\mathit{\boldsymbol{x}}(t_i)$ and $\mathit{\boldsymbol{u}}(t_i)$ by $\mathit{\boldsymbol{x}}_i$ and $\mathit{\boldsymbol{u}}_i\, \forall i$, respectively, then the NLP can be solved using well-developed optimization software for the unknowns $\tilde{\mathit{\boldsymbol{x}}}_i = \mathit{\boldsymbol{x}}_i, \, i = 1, \ldots, n$, and $\tilde{\mathit{\boldsymbol{u}}}_j = \mathit{\boldsymbol{u}}_j, \, j = 0, \ldots, n$. The approximate optimal state and control variables can then be calculated at any point $t \in [0, \infty)$ through the PS expansions

where $\mathit{\boldsymbol{t}}_{n} = \left[t_{0}, t_{1}, \ldots, t_{n}\right]^{\top}$ and $\boldsymbol{\mathcal{{L}}}_{n}(t) = \left[\mathcal{{L}}_{n, 0}(t), \mathcal{{L}}_{n, 1}(t), \ldots, \mathcal{{L}}_{n, n}(t)\right]^{\top}.$ In the special case, when $T = T_{1, L}^{(\alpha)}$, one can easily show that the NLP can be written as follows:

subject to

where $\oslash$ denotes the Hadamard division, $(\mathit{\boldsymbol{v}})_{r} = \underbrace { \mathit{\boldsymbol{v}} \circ \mathit{\boldsymbol{v}} \circ \ldots \circ \mathit{\boldsymbol{v}}}_{r\text{ times}}$, for any vector $\mathit{\boldsymbol{v}}$, and

Furthermore, when $T = T_{2, L}^{(\alpha)}$, the NLP can be formulated as follows:

subject to

where

We refer to the NLPs described by Eqs (4.29a), (4.29b), (4.31a) and (4.31b) by NLP1 and NLP2, respectively. We also refer to the present collocation method as the "GGR-IPS" method; the acronyms "GGR-IPS1" and "GGR-IPS2" stand for the GGR-IPS method performed using the parametric maps $T_{1, L}^{(\alpha)}$ and $T_{2, L}^{(\alpha)}$, respectively, while "GGR-IPS12" stands for the GGR-IPS method performed using either maps $T_{1, L}^{(\alpha)}$ and $T_{2, L}^{(\alpha)}$. Notice that the Lagrangian associated with the NLP described by Eqs (4.25) and (4.27) is defined by

where $\mathit{\boldsymbol{T}}'_n = T'\left(\mathit{\boldsymbol{\tau}}_{n}\right)$ and $\mathit{\boldsymbol{r}} = \left[r_{10}, \ldots, r_{1n}, \ldots, r_{n_{x}0}, \ldots, r_{n_{x}n}\right]^{\top}$ is the vector of Lagrange multipliers. Therefore, the Karush–Kuhn–Tucker necessary conditions of optimality are given by

where the operators ${ \mathop {\nabla}\limits_{\tilde { \mathit{\boldsymbol{x}}}} = \left[\frac{\partial }{\partial \tilde{x}_{11}} \ldots \frac{\partial }{\partial \tilde{x}_{1n}} \ldots \frac{\partial }{\partial \tilde{x}_{n_{x}1}} \ldots \frac{\partial }{\partial \tilde{x}_{n_{x}n}} \right], \mathop {\nabla}\limits_{\tilde { \mathit{\boldsymbol{u}}}} = \left[\frac{\partial }{\partial \tilde{u}_{10}} \ldots \frac{\partial }{\partial \tilde{u}_{1n}} \ldots \frac{\partial }{\partial \tilde{u}_{n_{u}0}} \ldots \frac{\partial }{\partial \tilde{u}_{n_{u}n}} \right]}$, $\mathbf{{I}}_{n}$ is the identity matrix of size $n, \hat{\mathbf{{{F}}}} = \text{vec}\left(\mathbf{{F}}\right), \mathbf{{E}} = [\mathit{\boldsymbol{\mathit{0}}}_{n}^{\top}; \mathbf{{I}}_{n}]$ and $[.; .]$ denotes the vertical matrix concatenation.

5.   Error and convergence analyses

In this section, we derive the truncation error bounds for Eqs (4.25) and (4.27) and their convergence rates.

Theorem 5.1. Let $f \in C^{n+1}[-1, 1)$ be approximated by a Gegenbauer interpolant $P_{n}f$ based upon the GGR points set $\mathbb{{S}}_n$, as defined by Eq (4.4). Then, there exist some numbers $\xi_{i} \in (-1, 1), i = 0, \ldots, n,$ such that the truncation error of the GGR-SR quadrature rule of Eq (4.22) is given by

where $K_{n}^{(\alpha)} = 2^{n-1} {\frac{\Gamma(n+\alpha) \Gamma(2 \alpha+1)}{\Gamma(n+2 \alpha)\Gamma(\alpha+1)}}, \, n = 0, 1, \ldots$.

Proof. By definition, we can write

for some $\xi \in (-1, 1)$, where ${}_fE_n$ is the interpolation truncation error at the GGR points given by

The proof is established by realizing that $\mathcal{{G}}_{n+1}^{(\alpha)}(\tau) = K_{n+1}^{(\alpha)}\prod_{k = 0}^{n}\left(\tau-\tau_{k}\right)$ and integrating Eq (5.2) on $\left[-1, \tau_{i} \right)\, \forall i$.

The following result is a direct corollary of Theorem 5.1 by letting $\eta: [-1, 1) \to \mathbb{{R}}$ and $\psi_{k}: [-1, 1) \to \mathbb{{R}}$ such that $\eta(\tau) = T'(\tau) g(\tilde{\mathit{\boldsymbol{x}}}(\tau), \tilde{\mathit{\boldsymbol{u}}}(\tau))$ and $\psi_{k}(\tau) = T'(\tau) f_{k}\left(\tilde{\mathit{\boldsymbol{x}}}(\tau), \tilde{\mathit{\boldsymbol{u}}}(\tau)\right)$ for each $k = 1, \ldots, n_x$.

Corollary 5.1. The truncation errors of Eq (4.27) and each approximate equation of System (4.25),

of the Integral Constraints System (3.1d) at each point $\tau_j \in \mathbb{{S}}_n$ are given by

and

respectively, where $\zeta, \xi_{j} \in(-1, 1)\, \forall j$.

The following upper bounds on the truncation errors of Eqs (4.25) and (4.27) can be deduced from [24,Theorem 5.1].

Theorem 5.2. Let $\psi_{k} \in C^{n+1}[-1, 1)$ and $\left\| \psi_{k}^{(n + 1)} \right\|_{[-1, 1)} = A_{\psi_k, n} \in \mathbb{{R}}^{+}\, \forall k$ for some constant $A_{\psi_k, n}$ dependent on $n$ and $k$. Then, there exist some positive constants ${B_{1}^{(\alpha)}}$ and ${C_{1}^{(\alpha)}}$ dependent on $\alpha$ and independent of n such that the truncation errors of System (4.25) at each point $\tau_j \in \mathbb{{S}}_n$ are bounded by the following inequalities:

where $\xi_j \in (-1, 1)\, \forall j$ and

Moreover,

and

as $n \rightarrow \infty$, where $\mathop < \limits_ \sim$ means "less than or asymptotically equal to."

Theorem 5.3. Let $\eta \in C^{n+1}[-1, 1)$ and $\left\| \eta^{(n + 1)} \right\|_{[-1, 1)} = A_{\eta, n} \in \mathbb{{R}}^{+}$ for some constant $A_{\eta, n}$ dependent on $n$. Then, there exist some positive constants ${B_{2}^{(\alpha)}}$ and ${C_{2}^{(\alpha)}}$ dependent on $\alpha$ and independent of n such that the truncation error of Eq (4.27) is bounded by the following inequality:

where $\zeta \in (-1, 1)$. Moreover,

and

as $n \rightarrow \infty$.

Divergence of typical IPS collocation schemes of the FHOCI at any large mesh grid of Gauss type when $T = T_{1, L}^{(\alpha)}$ or $T = T_{2, L}^{(\alpha)}$

In this section, we derive some striking results regarding the convergence of typical collocation schemes of the FHOCI described by Eqs (3.1a), (3.1c) and (3.1d) when $T \in \{T_{1, L}^{(\alpha)}, T_{2, L}^{(\alpha)}\}$ and the mesh grid is large and of Gauss type. While the proof pertains to the FHOCI and employs the GGR points as the collocation points, it can be generalized to the usual form of the FHOC described by Eqs (3.1a)–(3.1c) and any large collocation points of Gauss type, causing it to become of considerably greater interest. We derive these interesting divergence results in the following two corollaries.

Corollary 5.2. Let $T \in \lbrace T_{1, L}^{(\alpha)}, T_{2, L}^{(\alpha)} \rbrace$, and suppose that $\exists\, \hat{k} \in \{1, \ldots, n_x\}: \psi_{{ \hat{k}}} \in C^{n}[-1, 1)$, $0 < {\left\| f_{{ \hat{k}}} \right\|_{[-1, 1)}} < \infty$ and $0 \le \left\| {\frac{d^{j}}{d \tau^{j}}}f_{{ \hat{k}}} \right\|_{[-1, 1)} < \infty \quad \forall j = 1, \ldots, n+1$; then, the upper truncation error bounds of the Approximated System (4.25) diverge at each collocation point as $n \rightarrow \infty$ for any map scaling parameter value $L$.

Proof. By the general Leibniz rule, the $(n+1)$st derivative of $\psi_{{ \hat{k}}}$ is given by

with

Let $T = T_{1, L}^{(\alpha)}$, and notice that $\left(T_{1, L}^{(\alpha)}\right)^{(m)}(\tau) = {\frac{2L(m)!}{(1-\tau)^{m+1}}}\, \forall m \in \mathbb{{Z}}^{+}$, which is a monotonically increasing function for increasing values of $\tau$ as can be clearly seen in Figure 6. Therefore, $A_{\psi_{ \hat{k}}, n} = O\left(\left\|\left(T_{1, L}^{(\alpha)}\right)^{(n+2)}\right\|_{[-1, 1)} \right)$. From Theorem 5.2, there exist some positive constants $\hat B_{1}^{(\alpha)}$ and $\hat C_{1}^{(\alpha)}$ dependent on $\alpha$ and independent on $n$ such that

and

from which we realize that the upper bound of $\left|{}_{\psi_{ \hat{k}}} {E_n}\right|$ at each collocation point $\tau_{j}$ diverges as $n\rightarrow \infty$. Consider now the case when $T = T_{2, L}^{(\alpha)}$. By a similar argument, notice first that $\left(T_{2, L}^{(\alpha)}\right)^{(m)}(\tau) = {\frac{L(m-1)!}{(1-\tau)^{m}}}\; \forall m \in \mathbb{{Z}}^{+}$ is also a monotonically increasing function for increasing values of $\tau$, as shown in Figure 6. Therefore, $A_{\psi_{ \hat{k}}, n} = O\left(\left\|\left(T_{2, L}^{(\alpha)}\right)^{(n+2)}\right\|_{[-1, 1)} \right)$. From Theorem 5.2, there exist some positive constant $\hat B_{2}^{(\alpha)}$ and $\hat C_{2}^{(\alpha)}$ dependent on $\alpha$ and independent on $n$ such that

and

from which we observe that the upper bound of $\left|{}_{\psi_{ \hat{k}}} {E_n}\right|$ at each collocation point $\tau_j$ diverges as $n \rightarrow \infty$.

Theorem 5.2 and Corollary 5.2 reveal an interesting fact. Under their assumptions, the proposed method is expected to converge with an exponential rate to near-optimal solutions for increasing $n$ values within a relatively small $n$ value range, as indicated by inequality (5.7), but, as $n$ grows large, the constant $A_{\psi_{ \hat{k}}, n}$ grows exponentially fast and ultimately dominates the error bounds when $n \to \infty$, as implied by the inequalities (5.16a)–(5.17b), regardless of how well we choose the map scaling parameter value $L$. In fact, the asymptotic results of Corollary 5.2 manifest that, for increasing large $n$ values, reducing the $L$ value abates the divergence of the approximations at the outset, but, as $n$ grows larger, this approach fails to cope with the soaring values of $n$ powers of the factors $1/(1-\tau_j)$ at the mesh points $\tau_j$ for sufficiently close mesh points values to $\tau = 1$; thus, divergence is inevitable.

While the above forward error analysis may be too pessimistic and may reject solutions that are sufficiently accurate, another concern arises when we analyze the sensitivity of NLP1 and NLP2 associated with the maps $T_{1, L}^{(\alpha)}$ and $T_{2, L}^{(\alpha)}$ to input data errors. Observe that both problems require the computations of the maps $T'_{1, L}$ and $T'_{2, L}$, which are ill-conditioned for arguments near $1$. In particular, suppose that $\tau \approx 1$ with a small perturbation $h$ to $\tau$. Then, the absolute errors in computing $T'_{1, L}(\tau)$ and $T'_{2, L}(\tau)$ are given by

hence, the relative errors are ${\frac{2 |h|}{1-\tau}}$ and ${\frac{|h|}{1-\tau}}$, respectively, which blow up as $\tau \to 1$. Recall that GGR points cluster near $\pm 1$ as $n \to \infty$, so the sensitivity of the problem of calculating the maps' derivative functions $T'_{1, L}$ and $T'_{2, L}$ at arguments near $1$ increases for increasing values of $n$. For example, let $\tau = 0.9999999999999$ be an exact argument value and consider its approximation $\hat \tau = 0.9999999999998$ with a small perturbation of about $9.99 \times 10^{-14}$ to $\tau$. Then, the relative error in the input value is about $10^{-11} \%$. However, the relative errors in computing $T'_{1, L}(\tau)$ and $T'_{2, L}(\tau)$ are nearly $75 \%$ and $50 \%$. Hence, the relative changes in evaluating $T'_{1, L}(\tau)$ and $T'_{2, L}(\tau)$ are about $7.5$ and $5$ trillion times larger than the relative change in the input value in respective order. This example shows that increasing the mesh size shifts the positive collocation points closer and closer toward $\tau = 1$, and that wild ill-conditioning ultimately occurs, as the sensitivity of NLP2 progressively stiffens for arguments near $1$. Therefore, one should keep in mind that reducing $L$ may still improve the approximations for a certain range of $n$ values; nonetheless, this strategy is not prone to producing accurate approximations for relatively large values of $n$, in general, since both NLP1 and NLP2 are ill-conditioned near $\tau = 1$. It is noteworthy to mention here that this sensitivity of NLP1 and NLP2 near $\tau = 1$ is foreseen to relax or disappear if $g$ and $f_k\, \forall k$ decay exponentially fast such that $\mathop {\lim }\nolimits_{\tau \to 1} {T'(\tau) g\left(\tilde{\mathit{\boldsymbol{x}}}(\tau), \tilde{\mathit{\boldsymbol{u}}}(\tau)\right) } = 0$ and $\mathop {\lim }\nolimits_{\tau \to 1} {T'(\tau) \mathit{\boldsymbol{f}}(\tilde{\mathit{\boldsymbol{x}}}(\tau), \tilde{\mathit{\boldsymbol{u}}}(\tau))} = \mathit{\boldsymbol{\mathit{0}}}$. Under a similar proof to that of Corollary 5.2, one can derive the following second divergence result.

Corollary 5.3. Let $T \in \lbrace T_{1, L}^{(\alpha)}, T_{2, L}^{(\alpha)} \rbrace$, $\eta \in C^{n}[-1, 1)$, $0 < {\left\| g \right\|_{[-1, 1)}} < \infty$ and $0 \le \left\| {\frac{d^{j}}{d \tau^{j}}}g \right\|_{[-1, 1)} < \infty\, \forall j = 1, \ldots, n+1$; then, the upper truncation error bound of Eq (4.27) diverges as $n \rightarrow \infty$ for any map scaling parameter value $L$.

The present analysis begs another interesting question: Which map should we use if we desire to implement the proposed method? For small/medium ranges of $n$ values, the answer is a bit elusive; however, for large values, it seems that we have a crystal clear answer, as shown by the following corollary.

Corollary 5.4. A Gegenbauer-Gauss collocation of the FHOCI using the map $T_{1, L}^{(\alpha)}$ generally diverges faster than that achieved by applying the method lumped with the map $T_{2, L}^{(\alpha)}$ when $n \to \infty$.

Proof. The faster divergence exhibited using the map $T_{1, L}^{(\alpha)}$ compared with $T_{2, L}^{(\alpha)}$ as $n \to \infty$ can be easily justified by using the inequalities (5.16a)–(5.17b). Moreover, since $T'_{1, L}(\tau)$ grows faster than $T'_{2, L}(\tau)$ by a factor of $2/(1-\tau)$, which blows up as $\tau \to 1$, the ill-conditioning of $T'_{1, L}$ is clearly more severe than that of $T'_{2, L}$ for values of $\tau \approx 1$. Corollary 5.4 manifests that the map $T_{2, L}^{(\alpha)}$ is more likely a better choice than $T_{1, L}^{(\alpha)}$ for large $n$ values. We end this section by drawing the attention of the reader to the fact that integral reformulations of various mathematical models have received considerable attention in the literature because they often produce well-conditioned linear systems. While numerical quadratures and integration matrices are generally more stable than numerical differentiation operators and matrices, there is no strong reason to expect that standard PS collocations of the FHOC in its strong differential form, as obtained by using a single mesh grid of Gauss type and maps like $T_{1, L}^{(\alpha)}$ and $T_{2, L}^{(\alpha)}$, would exhibit any merits over the current method, and they would ultimately diverge for a large mesh grid size. These considerations lead naturally to the following interesting conjecture.

Conjecture 5.1. Classical Jacobi polynomial collocations of the FHOC in differential/integral form obtained through maps like $T_{1, L}^{(\alpha)}$ and $T_{2, L}^{(\alpha)}$ will likely diverge as the mesh size grows large if the computations are carried out using floating-point arithmetic and the discretizations use a single mesh grid, regardless of whether they are of Gauss/GR type or equally spaced. The former divergence case is a direct result of the present divergence analysis, while the latter case is due to Runge's phenomenon and the ill-conditioning of polynomial interpolation at equally spaced nodes as the degree of the polynomial grows.

6.   Numerical experiments

This section presents the results of some numerical experiments on three test examples, which demonstrate the accuracy and efficiency of the proposed GGR-IPS12 methods for small-medium-range mesh grid sizes and verifies the inevitable divergence as the mesh size grows large. All numerical experiments were carried out using MATLAB R2022a software installed on a personal laptop equipped with a 2.9 GHz AMD Ryzen 7 4800H CPU and 16 GB memory running on a 64-bit Windows 11 operating system. The NLPs obtained through the GGR-IPS12 methods were solved using either (i) MATLAB fmincon solver with the interior-point algorithm (fmincon-int) and sqp algorithm (fmincon-sqp), or (ii) the augmented Lagrange multiplier method ^[67,68] integrated with a modified Broyden–Fletcher–Goldfarb–Shanno method and a Chebyshev PS line search method ^[69], which is henceforth referred to as the "EALMM." It should be clearly understood by the reader when we coin the name of the current collocation method with any NLP solver that we are implementing them both to solve the FHOCI. For example, the acronym GGR-IPS12-EALMM stands for the GGR-IPS12 methods combined with the EALMM. In all numerical tests, the exact optimal state and control variables were calculated using MATLAB with $30$ digits of precision maintained in internal computations. The fmincon solver was carried out using the stopping criteria TolFun $=$ TolX $= 10^{-12}, 10^{-15}, 10^{-15}$ for Examples 1–3, respectively. The augmented Lagrange multiplier method was terminated when the lower bound on the change in the augmented Lagrangian function value during a step does not exceed $10^{-12}$. All experiments were conducted using the parameters values $L \in \{0.25(0.25)10\}$ and $\alpha \in \{-0.4(0.1)2\}$. Most of the numerical simulations were performed using two sets of initial guesses: $\mathbf{{\Omega}}_1 = \{(\tilde{\mathit{\boldsymbol{x}}}_{0}, \tilde{\mathit{\boldsymbol{u}}}_{0}): \tilde{\mathit{\boldsymbol{x}}}_{0} = \mathit{\boldsymbol{\mathit{1}}}_{n_x}, \tilde{\mathit{\boldsymbol{u}}}_{0} = \mathit{\boldsymbol{\mathit{1}}}_{n_u}\}$ and $\mathbf{{\Omega}}_2 = \{(\tilde{\mathit{\boldsymbol{x}}}_{0}, \tilde{\mathit{\boldsymbol{u}}}_{0}): \tilde{\mathit{\boldsymbol{x}}}_{0} = 0.5 \mathit{\boldsymbol{\mathit{1}}}_{n_x}, \tilde{\mathit{\boldsymbol{u}}}_{0} = 0.5 \mathit{\boldsymbol{\mathit{1}}}_{n_u}\}$; henceforth, $\mathbf{{\Omega}}_1 \cup \mathbf{{\Omega}}_2$ is denoted by $\mathbf{{\Omega}}$. Furthermore, by $AE_{J}$ and $MAE_{x, u}$, we mean the absolute error in the objective function value and the maximum absolute error of the state and control variables in respective order.

Example 1. Consider the IHOC (2.1)–(2.3) with $g\left(\mathit{\boldsymbol{x}}(t), \mathit{\boldsymbol{u}}(t)\right) = \left(\ln^{2} x(t) + u^{2}(t)\right)/2, f\left(\mathit{\boldsymbol{x}}(t), \mathit{\boldsymbol{u}}(t)\right) = x(t) \ln{x(t)}+x(t) u(t)$, and $\mathit{\boldsymbol{x}}_{0} = 2$. The exact state and control variables are

where

see ^[44]. The exact cost function ${J^*} = {(\ln 2)^2}{\mkern 1mu} \left({\sqrt 2 + 1} \right)/2 \approx 0.5799580911421756$ was rounded to $16$ significant digits. Through the change of variables

the IHOC described by Eqs (2.1)–(2.3) can be rewritten as an equivalent infinite-horizon linear-quadratic optimal control problem with $g\left(\mathit{\boldsymbol{z}}(t), \mathit{\boldsymbol{u}}(t)\right) = \left(z^2(t) + u^{2}(t)\right)/2, f\left(\mathit{\boldsymbol{z}}(t), \mathit{\boldsymbol{u}}(t)\right) = z(t) +u(t)$ and $\mathit{\boldsymbol{z}}_{0} = \ln 2$. We refer to the former and latter forms of the IHOC as Forms A and B, respectively. Form A of the example was previously solved in ^[44] by using LG- and LGR-PS methods and the three maps given by (3.2c) and (3.2d); the obtained NLPs were solved using the Sparse Nonlinear Optimizer (SNOPT) method ^[70,71]. The problem was solved later in ^[46] by using the FRPM ^[46] and the mapping $\zeta_c$:

which maps the scaled left half-open interval $(-1, 1]$ to the descending time interval $(\infty, 0]$. The obtained NLP was also solved using SNOPT. Table 1 shows a comparison between the LGR- and LG-PS methods, the FRPM and the GGR-IPS12-EALMM. The LGR- and LG-PS methods and the GGR-IPS12-EALMM were performed using the same initial guesses $\tilde{x}(\tau) = 2$ and $\tilde{u}(\tau) = \tau\, \forall \tau \in [-1, 1]$. Notice how the GGR-IPS2-EALMM generally enjoy superior stability properties and achieve higher-order approximations in this example for $n = 5(5)30$ compared with the other approaches, except for the LG-PS method, where they both achieve the same order of accuracy at $n = 30$. It is interesting here to recognize how the GGR-IPS2-EALMM defeats the LG-PS method for $n = 5(5)25$, although the latter employs a Gauss quadrature that is more accurate than the GR quadrature used by the former. One may connect the success of the GGR-IPS2-EALMM here to many reasons, namely, (i) the clever change of variables given by Eq (6.2) that converts the NLP into a linear-quadratic optimal control problem which can be collocated more accurately; (ii) the integral form of the system dynamics allows for gaining more digits of accuracy via numerical quadratures, which are well known for their numerical stability; (iii) the application of the highly accurate built-in Algorithm 2 to the current methods, as it applies the latest technology of SR interpolation; (iv) the parametric logarithmic map $T_{2, L}^{(\alpha)}$ that is favored over $T_{1, L}^{(\alpha)}$ for its slower growth and less sensitivity near $\tau = 1$; and (v) the map scaling parameter $L$, which permits faster convergence rates when "optimally" chosen. On the other hand, we observe that the errors of GGR-IPS12-EALMM generally decline gradually as the mesh grid size initially grow up to a certain limit, yet they bounce back beyond that limit as the mesh grid size continues to grow large, in agreement with the theoretical results of Section 5. It is interesting to see similar phenomena as in the control error profiles in ^[44] in the sense that (i) the control error plot of the LG-PS method does not appear as a (near) straight line in the shown log-scaled chart, but, rather, a convex-shaped curve, as it curves outward ([44,Figure 1(b)]), and (ii) the control error plot of the LGR-PS method suddenly increases at $n = 30$ much earlier before reaching the round-off plateau ([44,Figure 2(b)]). A closely related behavior to these observations can be found in ^[46], where we noticed that the control error history is also not linear in the shown log-scaled chart, but, rather, oscillate up and down at different levels again much earlier before reaching the round-off plateau ([46,Figure 5]). Another interesting remark lies in the smallest errors of the current methods; they were all recorded at/near $\alpha = 0.5$ with $\alpha \in \Upsilon_{0.5, 0.6}^G$, while several optimal values of $L$ were detected. Table 1 also shows the corresponding elapsed time (ET) to perform the GGR-IPS12-EALMM. The calculated execution times were measured multiple times and the shown data are the medians of the time measurements in seconds (s). Figures 7 and 8 show the plots of the exact state and control variables, in addition to their collocated solutions and absolute errors obtained by the GGR-IPS12 methods integrated with three NLP solvers using the same initial guesses set and several values of $n, L$ and $\alpha$. We can observe from the shown graphical data that the GGR-IPS2-EALMM generally achieves better accuracy and stability properties compared with the other methods.

Figures 9 and 10 show the $MAE_{x, u}$ of the GGR-IPS1-EALMM obtained by using $n = 10(10)50, \alpha = -0.4(0.1)2, L = 1(1)4$ and $(\tilde{\mathit{\boldsymbol{x}}}_{0}, \tilde{\mathit{\boldsymbol{u}}}_{0}) \in \mathbf{{\Omega}}$. It is interesting to observe here by visual inspection that, when holding $L$ fixed, the global minima of the error mesh surface plots occur near $\alpha = 0.5$ and the mesh surfaces rise up gradually as we move away, except when $\alpha \in \{-0.4, -0.3\}$, where sharp peaks may emerge suddenly for growing values of $n$. This suggests that Legendre polynomials seem to be an optimal choice among Gegenbauer basis polynomials when holding $L$ fixed, while Gegenbauer polynomials associated with $\alpha$ values near $-0.5$ may cause numerical instability as $n$ grows large. However, a different story emerges when the GGR-IPS2-EALMM is performed instead, as can be seen in Figures 11 and 12. Notice now that the errors appear to be monotonically decreasing for decreasing values of $\alpha$ when holding $L$ fixed at $1$ and $2$, and that the error surface shoots up as $n$ and $\alpha$ increases. Therefore, Gegenbauer polynomials with some negative $\alpha$ values seem optimal for relatively small values of $L$. Notice also that the sudden peaks observed before with the GGR-IPS1-EALMM in Figures 9 and 10 for $\alpha \in \{-0.4, -0.3\}$ and large $n$ values disappear. A further array of error mesh surface plots of the GGR-IPS1-EALMM are shown in Figures 13 and 14 for $n = 10(10)50, \alpha = -0.2, 0, 0.5, 1, L = 0.5(0.5)6$ and $(\tilde{\mathit{\boldsymbol{x}}}_{0}, \tilde{\mathit{\boldsymbol{u}}}_{0}) \in \mathbf{{\Omega}}$. While holding $\alpha$ fixed, there seems to be no general rule of thumb that can be drawn from the shown data. On the other hand, Figures 15 and 16 show the corresponding plots associated with the GGR-IPS2-EALMM, where the errors are very similar and can be clearly seen to surge as $L \to 0$, especially when $\alpha = 1$, but remain relatively small for $L = 2(0.5)6$.

Figures 17 and 18 show comparisons of the number of iterations required by the GGR-IPS12 methods when combined with three distinct NLP solvers using $(\tilde{\mathit{\boldsymbol{x}}}_{0}, \tilde{\mathit{\boldsymbol{u}}}_{0}) \in \mathbf{{\Omega}}$ and several $L$ and $\alpha$ values. Clearly, the integration of the GGR-IPS12 methods with the EALMM leads to a drastic reduction in the number of iterations in all cases. In fact, while the GGR-IPS12-EALMM often converged in only four/five iterations, other methods usually require many more iterations to converge; for example, the GGR-IPS12 methods combined with fmincon-int and fmincon-sqp took more than $200$ iterations to converge to the solutions of the problem when $(n, L, \alpha) = (48, 1, -0.2)$ and starting with any initial guess $(\tilde{\mathit{\boldsymbol{x}}}_{0}, \tilde{\mathit{\boldsymbol{u}}}_{0}) \in \mathbf{{\Omega}}$.

Table 2 shows the approximate cost function values obtained by the GGR-IPS2-EALMM for several parameter values. The fastest convergence was recorded at $\alpha = 0.5$ in all cases with $J \approx J_{16} = 0.579958091127$, which is in agreement with $J$ to $10$ significant digits. It is interesting to see through the tabulated data how the Gegenbauer polynomials with $\alpha \in \{-0.4, 0.25\}$ exhibit faster convergence rates than the Chebyshev polynomials (when $\alpha = 0$), capturing four correct significant digits as early as $n = 12$, whereas the Chebyshev polynomials are still lagging behind by one digit even when $n$ increases by two units. Gegenbauer polynomials with $\alpha = -0.2$ also performed better than the Chebyshev polynomials for $n \in \{14, 16\}$, while the poorest stability was that of Gegenbauer polynomials with $\alpha = 2$, scoring only one correct significant digit in all cases. Table 3 shows the smallest $MAE_{x, u}$ obtained by the GGR-IPS12-EALMM among the recorded errors for the parameter values $n = 10(10)80, \alpha = -0.4(0.1)2, L = 0.25(0.25)10$ and $(\tilde{\mathit{\boldsymbol{x}}}_{0}, \tilde{\mathit{\boldsymbol{u}}}_{0}) \in \mathbf{{\Omega}}$. The table shows the capacity of the GGR-IPS12-EALMM to achieve improved near-optimal solutions for increasing values of $n$ within a small/medium range of mesh grid size; however, the accuracy deteriorates beyond a certain limit as the mesh grid size grows larger, in agreement with the theoretical results proven in Section 5. The GGR-IPS2-EALMM is clearly superior to the GGR-IPS1-EALMM in terms of accuracy in all cases, and it is interesting here to see how the GGR-IPS1-EALMM diverges faster than the GGR-IPS2-EALMM for growing mesh sizes, as indicated earlier by the divergence analysis presented in Section 5. The best approximations obtained experimentally by the GGR-IPS2-EALMM were recorded at/near $\alpha = 0.5$ with $\alpha \in \Upsilon_{0.4, 0.6}^G$. The smallest errors of the GGR-IPS1-EALMM were also recorded at/near $\alpha = 0.5$ for $n = 5(5)55$ with $\alpha \in \Upsilon_{0.4, 0.6}^G$; however, the algorithm tends to favor larger positive $\alpha$ values beyond $n = 55$, where we noticed the travel of the right boundary, $\alpha_{c}^+$, of $\Upsilon_{0.4, \alpha_{c}^+}^G$ rightward from $\alpha_{c}^+ = 0.6$ into $\alpha_{c}^+ = 1.8$ as $n$ reaches $80$. This is no surprise. In fact, recall that $T_{1, L}^{(\alpha)}$ increases monotonically for decreasing values of $\alpha$ while holding $n$ and $L$ fixed, and we can observe from Figure 1 that the mapping escalates wildly as we continue decreasing the $\alpha$ values. The byproduct of this behavior is that increasing the $\alpha$ values moves the collocation points associated with large values of $t$ leftward, relocating them closer to regions where the solution changes rapidly. This graphical interpretation is consistent with the fact that the interior GGR points move monotonically toward the center of the interval $(-1, 1)$ as the parameter $\alpha$ increases ^[51,56]. From another perspective, the leftward movement of the collocation points near the right boundary $\tau = 1$ mitigates the effect of the ill-conditioning of $T'_{1, L}$ for arguments near $1$, as $T'_{1, L}$ is evaluated at mesh points that are gradually departing the vicinity of $\tau = 1$. This argument adds more tenability to the cause of using Gegenbauer polynomials as a basis polynomials for numerical collocations of FHOCIs obtained from IHOCs via $T_{1, L}^{(\alpha)}$ or $T_{2, L}^{(\alpha)}$ in the sense that, while Chebyshev and Legendre polynomials cease to downgrade the errors as the mesh grid size grows large, Gegenbauer polynomials have the additional advantage of being able to alleviate the growth rates of both $T_{1, L}^{(\alpha)}$ and $T_{2, L}^{(\alpha)}$ by increasing the $\alpha$ value whenever we wish, while sustaining the luxury of being able to apply either Chebyshev or Legendre polynomials. If we now turn our attention to the recorded $L$ values in the table, we can quickly spot that the smallest computed errors initially span a wide range of numerically optimal $L$ values; however, the solvers ultimately have a bias toward smaller values of $L$ in an attempt to damp the error in agreement with the forward error analysis presented in Section 5. Notice that the numerically optimal $L$ value for the GGR-IPS2-EALMM stays at $0.75$ for $n = 55(5)80$, while the corresponding values for the GGR-IPS1-EALMM occur at the smallest feasible $L$ value among the input range of experimental data. One may attribute this peculiar behavior of the solvers to the fact that the logarithmic map $T_{2, L}^{(\alpha)}$ increases at a much slower rate than that of the algebraic map $T_{1, L}^{(\alpha)}$; see Figures 1 and 2.

Example 2. Consider the IHOC given by Eqs (2.1)–(2.3) with $g\left(\mathit{\boldsymbol{x}}(t), \mathit{\boldsymbol{u}}(t)\right) = x_{1}^{2}(t)+x_{2}^{2}(t)/2+u^{2}(t)/4, \mathit{\boldsymbol{f}}\left(\mathit{\boldsymbol{x}}(t), \mathit{\boldsymbol{u}}(t)\right) = \left[x_{2}(t), 2x_{1}(t)-x_{2}(t)+u(t)\right]^{\top}$ and $\mathit{\boldsymbol{x}}_{0} = [-4, 4]^{\top}$. The exact state and control variables are

where

see ^[45,55,72]. This example is a linear-quadratic regulator problem with an optimal cost functional value $J^{*} = 19.85335656362790$, rounded to $16$ significant digits, as obtained in MATLAB by using the Symbolic Math Toolbox. Figure 19 shows the plots of the exact optimal state and control variables and their approximations obtained through GGR-IPS2-EALMM using some parameter values. Table 4 shows the $MAE_{x, u}$ of the LGR-PS method in ^[45] and the smallest corresponding $MAE_{x, u}$ and $AE_{J}$ pairs of the GGR-IPS2-EALMM at the collocation points, as obtained by using $(\tilde{\mathit{\boldsymbol{x}}}_{0}, \tilde{\mathit{\boldsymbol{u}}}_{0}) \in \mathbf{{\Omega}}$ and several values of $n$. The table also shows the corresponding ET taken to perform the GGR-IPS2-EALMM. The GGR-IPS2-EALMM proves again to be superior in terms of accuracy for a small range of mesh sizes, as it converges rapidly to near-optimal solutions at a much higher rate than that in ^[45]. However, the superb accuracy of the method starts to decline when $n$ grows larger, as anticipated earlier. Notice again here that the best accuracy in all cases was recorded at $\alpha = 0.5$ with SRCIC $\Upsilon_{0.5, 0.5}^G = \{0.5\}$.

Another comparison between GGR-IPS2-fmincon-int, GGR-IPS2-fmincon-sqp and the transformed LGR method in ^[47] is shown in Table 5. We can clearly see that the former two methods generally yield smaller $AE_{J}$ values. The rise and fall of accuracy as the mesh size grows is again peculiar in the observed approximations, which is in agreement with the presented divergence analysis in Section 5. Remarkably, a match with the exact $J^*$ to full machine precision was recorded as early as $n = 20$, indicating an exceedingly accurate numerical scheme with exponential convergence for coarse meshes. All of the smallest errors reported by the current methods occurred at $\alpha = 0.5$, except for $n \in \{90, 100\}$, where collocations at $\alpha = 0$ and $0.8$ furnished higher accuracy. On the other hand, the rounded errors in ^[47] decay to $1.35 \times 10^{-09}$ as soon as $n = 30$ when using the algebraic map, but they cease to vary any further for $n = 40(10)100$.

Table 6 shows a third comparison between GGR-IPS2-fmincon-int, GGR-IPS2-fmincon-sqp and the recent method of Mamehrashi and Nemati ^[73], which uses Laguerre functions and the Ritz spectral method. All errors produced by the current methods were smaller than those obtained in ^[73] under the condition of using the same number of approximation terms; notably, a large difference of nine and ten orders of magnitude was observed in favor of the current methods using as small as $n = 10$.

Example 3. Consider the IHOC given by (2.1)–(2.3) with

and $\mathit{\boldsymbol{x}}_{0} = [0, 0.8]^{\top}$, where $\mathbf{{Q}} = \mathbf{{R}} = \mathbf{{I}}_2$. This problem was solved earlier in ^[74] by using a piecewise Adomian decomposition method, and again recently in ^[73] by using Laguerre functions and the Ritz spectral method. Table 7 shows comparisons between these two methods and the GGR-IPS2-EALMM with $\alpha = 0.5, L = 1$ and $n = 1(1)10$. Because of the agreement of the approximations obtained by the GGR-IPS2-EALMM for $n = 9$ and $10$, we could reasonably expect the value of $J_{10}$ to be accurate to the places listed. It is remarkable here to notice that $J_1$ agrees with $J_{10}$ to three significant digits under the condition of using as little as two collocation nodes, which reflects the swift convergence rate achieved by the current method for coarse meshes. Notice also that the obtained $J_5$ value is smaller than the corresponding approximate optimal cost function values obtained using both rival methods. The state and control numerical solutions obtained by the GGR-IPS2-EALMM are shown in Figure 20.

7.   Conclusions and future work

Direct IPS methods for solving IHOCS using the logarithmic mapping $T_{2, L}^{(\alpha)}$ and the developed SR interpolation and quadrature formulas can produce excellent approximations to the optimal state and control variables for relatively small/medium mesh grids. However, this class of methods often suffer from numerical instability for fine meshes when endowed with any of the parametric maps $T_{i, L}^{(\alpha)}, i = 1, 2$; therefore, as the mesh size grows, they are not as useful as one might hope for computing the optimal state and control trajectories to within high precision. In fact, it has been shown in the current paper that two sources of difficulty arise in handling the horizon in IHOCs using a domain transformation that maps the infinite horizon to the finite horizon $[-1, 1)$ through the algebraic and logarithmic maps $T_{i, L}^{(\alpha)}, i = 1, 2$: (i) the exponential growth of the mappings surface slopes near the right boundary $\tau = 1$, which increases the truncation errors produced in the FHOCI discretization without bounds as $\tau \to 1$; and (ii) despite the fact that both mappings $T_{i, L}^{(\alpha)}, i = 1, 2$ have a singularity at $\tau = 1$ and we actually never evaluate them at the singularity, since the GGR collocation points are strictly less than 1, their derivatives are sensitive to input data errors for arguments near $\tau = 1$; thus, both NLP1 and NLP2 are ill-conditioned for $\tau \approx 1$. These theoretical facts, as well as the observed empirical data, are considerable reasons to say that typical direct spectral/PS and IPS methods based on classical Jacobi polynomials and the parametric maps $T_{i, L}^{(\alpha)}, i = 1, 2$ are foreseen to diverge as the mesh size grows large if the computations are carried out using floating-point arithmetic and the discretizations use a single mesh grid, regardless of whether they are of Gauss type or equally spaced, which is a result that contradicts the convergence claim of such methods using Legendre polynomials made in ^[44].

While Gegenbauer polynomials associated with certain nonpositive $\alpha$ values are well suited for FGGR-based polynomial interpolations in Lagrange basis form over fine meshes, as shown in ^[24]; this paper asserts that Gegenbauer polynomials associated with certain nonnegative $\alpha$ values are more apt for GGR-based SR interpolations over fine meshes. Moreover, for coarse mesh grids, Legendre polynomials are particularly (near) optimal basis polynomials for GGR-based SR collocations of FHOCIs, as argued in Section 4.1 and sustained through numerical simulations. On the other hand, Gegenbauer polynomials associated with certain positive values of $\alpha \in (1/2, 2]$ are optimal for IHOCI collocations over fine mesh grids, as they can largely slow down the exponential growth of both parametric maps $T_{i, L}^{(\alpha)}, i = 1, 2$, and their associated GGR collocation points are less dense near $\tau = 1$; thus, the sensitivity of computing $T'_{i, L}, i = 1, 2$ at arguments near $\tau = 1$ is significantly attenuated. The paper also shows that the parametric map $T_{1, L}^{(\alpha)}$ is more severely sensitive for $\tau \approx 1$ than $T_{2, L}^{(\alpha)}$, and that the family $\left\{\left(T_{1, L}^{(\alpha)}\right)^m\right\}_{m = 0}^\infty$ grows faster than $\left\{\left(T_{2, L}^{(\alpha)}\right)^m\right\}_{m = 0}^\infty$ as $\tau \to 1$; therefore, $T_{2, L}^{(\alpha)}$ is more apt for the domain transformation of IHOCs than $T_{1, L}^{(\alpha)}$ for collocation points of Gauss/GR type.

It is worthy to mention that direct IPS methods based on the proposed Gegenbauer SR collocation can exhibit faster convergence rates for coarse meshes by regulating the map scaling parameter $L$ and the Gegenbauer parameter $\alpha$. In light of the stability analysis conducted in Section 4.1, GGR-based SR collocations of well-conditioned problems are generally endorsed for $\alpha$ values within/near the SRCIC $\Upsilon_{1/2, 1}^G$; the current study also supports this rule of thumb for IHOCs that are converted into FHOCIs through the use of parametric maps $T_{i, L}^{(\alpha)}, i = 1, 2$ and then collocated at relatively coarse mesh grids. On the other hand, a rule of thumb for choosing the optimal map scaling parameter value, $L^*$, based on the derived error bounds in Section 5 and extensive numerical simulations performed in Section 6 suggests to choose $L$ within the range $(0, 1)$ for large mesh grids. However, the question of how can we find $L^*$ remains open for coarse meshes. An interesting direction for future work may involve a study of new mappings with slower growth rates and derivatives with less sensitivity to input data errors. Another interesting direction is to extend the current work to handle IHOCs subject to boundary value problems.

Conflict of interest

The authors declare no conflicts of interest regarding this paper.

Appendix

A.   Barycentric GRDM

To construct the barycentric GRDM, we follow the derivation presented in ^[23] and multiply both sides of Eq (4.6) by $x-\tau_{j}\, \forall j$ to render them differentiable at $x = \tau_{j}$ such that

Letting ${S(x) = \sum\limits_{k = 0}^{n}\frac{\xi_{k}\left(x-\tau_{j}\right)}{x-\tau_{k}}}$ and differentiating Eq (A.1) with respect to $x$ yields

Since ${S(\tau_{j}) = \xi_{j}, S'(\tau_{j}) = \sum\limits_{j \neq k}\frac{\xi_{k}}{ \tau_{j}-\tau_{k}}}$ and $\mathcal{{L}}_{n, i}(\tau_{j}) = 0\, \forall i \neq j$, the off-diagonal elements of the differentiation matrix $\mathbf{{D}} = (d_{j, i})_{0 \le j, i \le n}$ can be calculated by using the following formula:

For $i = j$, we have ${\sum\limits_{i = 0}^{n} \mathcal{{L}}_{n, i}\left(x\right) = 1}$, so ${\sum\limits_{i = 0}^{n} \mathcal{{L}}'_{n, i}\left(x\right) = 0}$ and

Hence, the derivative of a real-valued function $f \in C^{1}[-1, 1]$ can be approximated at the GGR points by using the following formula:

B.   Computational algorithms