On Runge Kutta Methods Based on Indefinite Integral Approaches

1. Introduction

Differential equations are well-known among scientists as useful tools for scientific modeling. In this introductory discussion, we do not want to rehash what has already been said in numerous books and articles; see e.g. [1] and references therein, but rather to present widely recognized ideas in the terminology that will be used throughout the rest of this article.

Let us consider a system of ordinary differential equation (ODE) which can be stated in a compact form in the following form:

$$y^{\prime }={\displaystyle \frac{\mathrm{dy}\left(t\right)}{\mathrm{dt}}}=f(t,y,\lambda ),$$

(1)

i.e. $y:\mathbb{R}\to {\mathbb{R}}^n$ where n is the dimension of the solution space. The dash in (1) represents differentiation with respect to the independent variable, say t, which commonly represents time in dynamical systems. The dependent variables y are known as state variables, whereas $\lambda$ is a collection of parameters represented as a vector. When addressing the components of a vector differential equation, we use the term “system of differential equations” to refer to equation (1).

Using a differential equation to represent the development of a state vector y in a physical process poses the challenge of predicting future values of y based on the starting or initial value of y. The mathematical model is represented by a set of equations. Using a system of differential equations to represent the evolution of a state vector in a physical process presents the difficulty of forecasting future values based on the initial value $y\left(t_0\right)=y_0$. This system of equations in connection with (1) serves as the representation for the mathematical model in the initial state.

Every system of ordinary differential equations (1) is equivalent to a single differential equation including a higher derivative. The inverse of this assertion is true for higher-order differential equations, which can be solved to the highest order. As a result, we shall confine our discussion to systems of ordinary differential equations in the form (1).

For the system of first-order ordinary differential equations (1), we will provide numerical solution approaches based on implicit Runge-Kutta (RK) procedures. These approaches are primarily dependent on the presence of discrete representations of indefinite integrals ${\mathcal{J}}^{+}$. In other words, we shall demonstrate a significant numerical application of indefinite integrals to implicit Runge-Kutta techniques. We shall demonstrate that the discrete formulation of indefinite integrals leads directly to the numerical integration of initial value problems.

In this regard, there is nothing new; because indefinite integrals may be utilized to solve initial value problems via direct integration. What is new is that the Runge-Kutta techniques can now include the discretization of indefinite integrals in different forms. This generalizes to a large variety of indefinite integral representations using Sinc techniques or orthogonal polynomial (OP) roots [2,3]. We’ll utilize both Sinc and root based approaches.

The generalization stems from the finding that indefinite integration is the ancestor of the Runge-Kutta technique. We employ Sinc, Sinc-Exp, Lagrange and generalized Lagrangian interpolation for indefinite integration, which in this particular instance uses either Sinc functions or orthogonal polynomials as basis, respectively [3,4].

Runge and Kutta [5,6] developed a numerical solution process that allows for both explicit and implicit representations of calculation formulae. The explicit Runge-Kutta technique is simple to use, but the implicit Runge-Kutta approach needs a simultaneous system of equations to be solved at each step. Except for f being linear in y, they are nonlinear equations. For the first approach to discuss the relations, we use a scalar ODE to clearly demonstrate the links. Here, we will look at scalar initial value problems of the type

$$y^{\prime }\left(t\right)=f(t,y\left(t\right),\lambda ),\mathrm{with}y\left(t_0\right)=y_0,$$

(2)

where $y^{\prime }\left(t\right)=\mathrm{dy}/\mathrm{dt}$. Furthermore, we assume that they are non-stiff or stiff ODEs. In the context of stiff ODEs, Hairer and Wanner‘s work [8] has to be noted. Nonetheless, it will become clear that our method, which is based on the definition of indefinite integrals, can handle both stiff and non-stiff ODEs. The method is not new in and of itself; it dates back to the 1960s of the previous century. The Runge-Kutta technique can be reduced to indefinite integrals in a collocation process, as Guillou and Soulé discussed in a 1969 paper [9]. Wright’s publication [10] contained additional verification of this. This approach provides an estimate for both explicit and implicit Runge-Kutta methods presentations.

Some kinds of ordinary differential equations can be easily solved numerically using the Runge-Kutta technique. However, deriving high-order Runge-Kutta techniques can be a challenging undertaking. This is due to a variety of factors. Identifying the so-called order criteria is the first problem [1]. To make the error in the process of order for some digits, these nonlinear equations need to be met. Solving these equations is the second obstacle. In addition to being nonlinear, several heuristics and simplifying assumptions are typically used, and there is typically no unique solution. The problem of combinatorial explosion is another. To give an idea on the related problems, a twelfth-order approach has 7813 order criteria exploding the efficiency [1]. To avoid all these mentioned obstacles we take here the route via indefinite integration.

The article is structured as follows: Section 2 discusses the key features of indefinite integrals and Runge-Kutta techniques. Section 3 describes three methods for approximating indefinite integrals. Section 4 applies the methods described in Section 3 to implicit Runge-Kutta procedures. Section 5 discussed the stability problem of the implicit Runge-Kutta method. Section 6 uses examples to show the solution of stiff and non-stiff initial value problems. Section 7 summarizes the findings and addresses unresolved questions.

2. Methods

This section covers the principles of the implicit Runge-Kutta technique. The principles of indefinite integrals provide an important foundation for explaining the numerical technique, which is based on the collocation of a Volterra integral equation.

2.1. Indefinite Integrals

This section covers the numerical computation of indefinite integrals $\left({\mathcal{J}}^{\pm }f\right)$ for a function $f:(a,b)\to \mathbb{R}$. We use the notation ${\mathcal{J}}^{\pm }$ for the generic operators defined by the two following relations:

$$\left({\mathcal{J}}^{+}f\right)\left(x\right)={\int }_a^{x}f\left(t\right)dt,$$

(3)

and

$$\left({\mathcal{J}}^{-}f\right)\left(x\right)={\int }_x^{b}f\left(t\right)dt.$$

(4)

The two operators are linear i.e.

$${\mathcal{J}}^{\pm }\left(c_{1}f+c_{2}g\right)=c_1{\mathcal{J}}^{\pm }\left(f\right)+c_2{\mathcal{J}}^{\pm }\left(g\right).$$

(5)

Inspection provides the proof. Another characteristic of $\left({\mathcal{J}}^{\pm }f\right)$ is that the two operators are bounded [4]. In the following, we will investigate how the indefinite integral operators ${\mathcal{J}}^{\pm }$ may be represented by Sinc, Sinc-Gaussian, or Poly-Sinc methods, respectively. However, let us first discus the connection between indefinite integrals and Runge-Kutta methods.

2.2. Runge-Kutta Method

We start with the initial value problem (2) and represent it as an integral equation

$$y\left(t\right)=y\left(t_0\right)+{\int }_{t_0}^{t}f(\tau ,y\left(\tau \right),\lambda )d\tau .$$

(6)

To obtain a numerically approximate solution to this problem, we will choose a sequence of values $t_0=0,t_1,\dots ,t_n,\dots$ and we will define, step by step, approximations (with vector coefficients) $z_0\left(t\right)$,..., $z_n\left(t\right)$,... which will respectively approximate $y\left(t\right)$ on the successive intervals $\left(t_0,t_1\right)$, …, $\left(t_n,t_{n+1}\right)$, …. Here the $z_n\left(t\right)$ satisfies the following conditions:

1.

$z_n\left(t_n\right)=z_{n-1}\left(t_n\right)$, $z_0\left(0\right)=y_0$,

2.

${\mathrm{dz}}_n\left(t_i\right)/\mathrm{dt}=f\left(t_i,z_n\left(t_i\right),\lambda \right)$ with $t_i=t_n+\left(t_{n+1}-t_n\right){\tau }_i$ and $i=1,\dots \mathrm{s}$.

Next we replace the right hand side of (2) by the derivative in (6) and write

$$y\left(t\right)=y\left(t_n\right)+{\int }_{t_n}^t{y}^{\prime }\left(\tau \right)d\tau .$$

(7)

Using a specific basis $B_i\left(t\right)$ the derivative of y can be represented as an approximation in an Hilbert space as

$$y^{\prime }=\sum _{i=1}^s{k}_i{B}_i\left(t\right),$$

(8)

using an s-order approximation. Contrary to the standard notation, we call the approximation an s-order approximation to indicate the order of the approximating matrix $A^{+}$ which becomes clear in the following. This approximation results to the representation

$$z_n\left(t\right)=z_{n-1}\left(t_n\right)+{\int }_{t_n}^t\sum _{i=1}^s{k}_i{B}_i\left(\tau \right)d\tau =z_{n-1}\left(t_n\right)+\sum _{i=1}^s{k}_i{\int }_{t_n}^t{B}_i\left(\tau \right)d\tau ,$$

(9)

and subsequently after discretization of t we get

$$z_n\left(t_j\right)=z_{n-1}\left(t_n\right)+\sum _{i=1}^s{k}_i{\int }_{t_n}^{t_j}{B}_i\left(\tau \right)d\tau =z_{n-1}\left(t_n\right)+\sum _{i=1}^s{k}_i{A}_{i,j}^{+}.$$

(10)

Here $A^{+}$ denotes an $s\times s$ matrix of order s. This leads us to the following method: determining the $k_i$ using the system of implicit equations

$$k_i=f\left(t_0+c_{i}h,{\tilde{y}}_n+h\sum _{j=1}^s{A}_{i,j}^{+}{k}_j\right),i=1,2,\dots ,s,$$

(11)

and the discrete approximation results to

$${\tilde{y}}_{n+1}={\tilde{y}}_n+h\sum _{i=1}^s{b}_i{k}_i,$$

(12)

where ${\tilde{y}}_n=z_{n-1}\left(t_n\right)=z_n\left(t_n\right)$, ${\tilde{y}}_{n+1}=z_n\left(t_{n+1}\right)$, and $h=t_{n+1}-t_n$ with $n=0,1,2,\dots$.

The integration matrix $A^{+}$and the Runge-Kutta weights $b_i$ can be calculated by using integration as follows if the n discretization points $t_1,t_2,\dots ,t_n$ are known, i.e., either the Sinc points or the roots of classical orthogonal polynomials:

$$A_{i,j}^{+}={\int }_0^{t_i}{B}_j\left(\tau \right)d\tau ,$$

(13)

and

$$b_i={\int }_0^1{B}_i\left(\tau \right)d\tau ,\mathrm{with}i,j=1,2,\dots ,s,$$

(14)

with

$$c_i=\sum _{j=1}^s{A}_{i,j}^{+},$$

(15)

where $B_i\left(t\right)$ are either Sinc, Sinc-Exp functions, elementary Lagrange or elementary generalized Lagrange interpolation polynomials. Knowing this relations and the connection between discretization and the used base functions we will denote them as Sinc-RK, Sinc-Exp-RK, L-RK, and gL-RK methods, respectively. Integration in (13) and (14) is effectively based on Sinc quadrature [3,34]. Knowing the discretization vector $\mathit{V}\left(t\right)=\left(t_1,t_2,\dots ,t_n\right)$, the $s\times s$ integration matrix $A^{+}$, and the Runge-Kutta weights $\mathit{V}\left(b\right)=\left(b_1,b_2,\dots ,b_s\right)$ enables the use of the Butcher tableaux [1], which summarizes the determining elements of a Runge-Kutta technique:

(16)

The tableaux allow us to classify the Runge-Kutta technique as explicit or implicit. If matrix $A^{+}$ is completely filled with nonzero elements, it is an implicit method. If matrix $A^{+}$ can be reduced to a triangular matrix, an explicit method is provided. In the above mentioned cases, we will always have to deal with implicit methods in our different approaches since matrix $A^{+}$ is completely filled. Because matrix $A^{+}$ is a $s\times s$ matrix, which precisely defines the Runge-Kutta method, we can establish the method’s order using the dimension s. Contrary to commonly used nomenclature, we refer to s as the order of $A^{+}$ rather than the s-stage technique [1,8,11]. If we create $A^{+}$ with s discrete points $\mathit{V}\left(t\right)$, we have an s-order Runge-Kutta method. So we’ll utilize an implicit Runge-Kutta method of order s, which is derived from either a Sinc, Sinc-Exp, Lagrange, or a generalized Lagrange collocation.

3. Variants of Methods

In this section, we will look at the RK techniques’ discrete representations. The representation here refers to the discretization points and the precise basis elements used. We start with the Sinc approach, which takes translated Sinc functions as a base and discretizes over Sinc points. In the second technique, we utilize orthogonal polynomials as basis functions, with their roots used for discretization. In the final technique, we employ generalized Lagrange polynomials as a foundation, which are based on orthogonal polynomials, and discretize using Sinc points or orthogonal polynomial roots.

3.1. Sinc Methods

Let us start by providing two definitions that will be incredibly useful throughout the remainder of this article.

Definition 1 

(Basic Definitions.). Let $\mathcal{D}$ be a simply-connected domain having boundary $\partial \mathcal{D}$. Let a and b denote two distinct points of $\partial \mathcal{D}$, and let ϕ denote a conformal map of $\mathcal{D}$ onto ${\mathcal{D}}_d$, where ${\mathcal{D}}_d=\{w\in \mathbb{C}:|\mathfrak{I}\mathfrak{m}\left(w\right)|\le d\}$, such that $\varphi \left(a\right)=-\infty$, and $\varphi \left(b\right)=\infty$. Let $\psi ={\varphi }^{-1}$ denote the inverse map, and let Γ be defined by

$$\Gamma =\{z\in \mathbb{C}:z=\psi (u),u\in \mathbb{R}\}.$$

(17)

Given ϕ, ψ and a positive number h, let us set

$$z_k=z_k\left(h\right)=\psi \left(kh\right),k=0,\pm 1,\pm 2,\dots .$$

(18)

Let us also define ρ by

$$\rho \left(z\right)=e^{\varphi \left(z\right)}.$$

(19)

Let $1\le p\le \infty$, and let ${\mathcal{H}}^p\left(\mathcal{D}\right)$ denote the family of all functions F that are analytic in $\mathcal{D}$, such that

$${\mathcal{N}}_p(F,\mathcal{D})=\left\{\begin{array}{ccc}{\left({\int }_{\partial \mathcal{D}}{\left|F\left(z\right)\right|}^p\left|dz\right|\right)}^{1/p}<\infty & \mathit{if}& 1\le p<\infty \\ {\displaystyle \underset{z\in \mathcal{D}}{sup}}\left|F\left(z\right)\right|<\infty & \mathit{if}& p=\infty \end{array}.\right.$$

(20)

Corresponding to numbers α and β, let ${\mathit{L}}_{\alpha ,\beta }\left(\mathcal{D}\right)$ denote the family of all functions $F\in \mathit{Hol}\left(\mathcal{D}\right)$ for which there exists a constant C such that

$$\left|F\left(z\right)\right|\le {\displaystyle \frac{{C\left|\rho \left(z\right)\right|}^{\alpha }}{{(1+|\rho \left(z\right)\left|\right)}^{\alpha +\beta }}}$$

(21)

for all z in $\mathcal{D}$. Finally, let us define ${\mathit{L}}_{\alpha }\left(\mathcal{D}\right)$ by ${\mathit{L}}_{\alpha }\left(\mathcal{D}\right)={\mathit{L}}_{\alpha ,\alpha }\left(\mathcal{D}\right)$.

It is convenient to define yet another important family of functions, ${\mathit{M}}_{\alpha ,\beta }\left(\mathcal{D}\right)$ with $0<\alpha \le 1$, $0<\beta \le 1$, and $0<d<\pi$. The family ${\mathit{M}}_{\alpha ,\beta }\left(\mathcal{D}\right)$ consists of all those functions $F\in \mathit{Hol}\left(\mathcal{D}\right)\cap C\left(\overline{\mathcal{D}}\right)$ such that $G\in {\mathit{L}}_{\alpha ,\beta }\left(\mathcal{D}\right)$ where

$$G\left(z\right)=F\left(z\right)-{\displaystyle \frac{F\left(a\right)+\rho \left(z\right)F\left(b\right)}{1+\rho \left(z\right)}}.$$

(22)

Just as for ${\mathit{L}}_{\alpha ,\beta }\left(\mathcal{D}\right)$ and ${\mathit{L}}_{\alpha }\left(\mathcal{D}\right)$, we define ${\mathit{M}}_{\alpha }\left(\mathcal{D}\right)$ by ${\mathit{M}}_{\alpha }\left(\mathcal{D}\right)={\mathit{M}}_{\alpha ,\alpha }\left(\mathcal{D}\right)$.

Sinc methods rely on the ${\mathit{L}}_{\alpha }\left(\mathcal{D}\right)$ and ${\mathit{M}}_{\alpha }\left(\mathcal{D}\right)$ classes, which are easily recognized in applications. They also assure the quick convergence of the Sinc approximations presented in the next theorems, as well as the approximations utilized in subsequent sections and parts of this work (see also [2,3].

The second concept we’ll often employ is a set of basis functions for approximations. The next step in Sinc approximations is to define a Sinc basis, $\left\{w_i\right.$}, that uniformly approximates the functions F on the interval $(a,b)=\Gamma$. For this reason, we provide the following definition.

Definition 2 

(Sinc Basis.). Let the conformal maps ϕ be given. The density ρ is defined as in (19) for single exponential maps. Let the Sinc function $S(k,h)$ be $S(k,h)\circ (·)=sin\left(\pi \right(·-kh)/h)/(\pi ·)$ and the Sinc-Gaussian or Sinc-Exp $\mathfrak{S}(k,h,c)\circ (·)=S(k,h)\circ (·)exp\left(-c{((·-kh)/h)}^2\right)$ be given as orthogonal basis functions, then by using

$$B(k,h,c)\circ \left(x\right)=\left\{\begin{array}{cc}S(k,h)\circ \left(z\right)& \mathit{if}c=0\\ \mathfrak{S}(k,h,c)\circ \left(z\right)& \mathit{if}c\ne 0\end{array},\right.$$

(23)

we are able to define the set of basis functions as follows

$${\gamma }_j(x,c)=B(j,h,c)\circ \varphi \left(x\right),j=-M,\dots ,N$$

(24)

$$w_j(x,c)=\left\{\begin{array}{ccc}{\displaystyle \frac{1}{1+\rho }}+\sum _{l=-(M-1)}^N\psi (lh,-1,0){\gamma }_l(x,c)& \mathit{if}& j=-M\\ {\gamma }_j(x,c)& \mathit{if}& -M+1\le j\le N-1\\ {\displaystyle \frac{\rho }{1+\rho }}-\sum _{l=-M}^{N-1}\psi (lh,0,1){\gamma }_l(x,c).& \mathit{if}& j=N\end{array}\right.,$$

(25)

where $\psi (x,a,b)={\varphi }^{-1}$ is the inverse conformal map using the interval limits as additional parameters.

The functions $w_j$ defined in (25) meet the relation $w_j\left(z_k\right)=1$ if $j=k$, and $w_j\left(z_k\right)=0$ if $j\ne k$ otherwise. This is nothing more than the interpolation condition. Note the piecewise definition changes ${\gamma }_j$ only at the endpoints by interpolation to guaranty a nearly homogenous local error.

This notation allows us to define a row vector ${\mathit{V}}_m\left(\mathfrak{B}\right)$ of basis functions

$${\mathit{V}}_m\left(\mathfrak{B}\right)=\left(w_{-M},\dots ,w_N\right),$$

(26)

with $w_j$ defined as in (25). For a given vector ${\mathit{V}}_m\left(u\right)={\left(u_{-M},\dots ,u_N\right)}^T$ we now introduce the dot product as an approximation of the function $u\left(z\right)$ by

$$u\left(z\right)\approx V_m\left(\mathfrak{B}\right).V_m\left(u\right)=\sum _{k=-M}^N{u}_k{w}_k,$$

(27)

which is simply a Sinc or SG interpolation of a known function, $u\left(z\right)$ in a single exponential representation.

Based on this notation, we will introduce in the next few subsections the different integrals we need [3].

3.1.1. Sinc and Sinc-Exp Approximation of Functions

In our approximations we will use two types of basis function a Sinc based and a Sinc-Gaussian based approximation. In general a basis function $B(k,h,c)$ is used to represent the approximation of the function f as follows

$$C_{h,c,M,N}\left[f\right]\left(z\right)=\sum _{k=-M}^{N}f\left(z_k\right)B(k,h,c)\circ \varphi \left(z\right),$$

(28)

where the basis functions are classified using the c parameter as

$$B(k,h,c)\left(x\right)=\left\{\begin{array}{cc}S(k,h)\circ \left(z\right)& \mathrm{if}c=0\\ \mathfrak{S}(k,h,c)\circ \left(z\right)& \mathrm{if}c\ne 0\end{array}\right.,$$

(29)

with $S(k,h)\left(z\right)=sin\left(\pi \right(z/h-k\left)\right)/\left(\pi \right(z/h-k\left)\right)=\mathrm{Sin}\mathrm{c}(z/h-k)$ represents the Sinc function and $\mathfrak{S}(k,h,c)\left(z\right)=\mathrm{Sin}\mathrm{c}(z/h-k)exp\left(-c{(z/h-k)}^2\right)$ is the Sinc-Gaussian (S/G). The first type of approximation results to the representation

$$\begin{array}{ccc}\hfill C_{h,M,N}\left[f\right]\left(z\right)& =& \sum _{k=-M}^{N}f\left(z_k\right)\mathrm{Sin}\mathrm{c}\left({\displaystyle \frac{\varphi \left(z\right)}{h}}-k\right)\hfill \\ & =& \sum _{k=-M}^{N}f\left(z_k\right)S(k,h)\circ \varphi \left(z\right),\hfill \end{array}$$

(30)

using the basis set of orthogonal functions

$$B_S={\{S(k,h)\circ \varphi \left(z\right)\}}_{k=-M}^N,$$

(31)

where $\varphi \left(z\right)$ is a conformal map. As discussed by Schmeisser and Stenger [12] a Sinc approximation of a function f can be given in connection with a Gaussian multiplier in the following representation

$$C_{h,c,M,N}\left[f\right]\left(z\right)=\sum _{k=-M}^{N}f\left(z_k\right)\mathrm{Sin}\mathrm{c}\left({\displaystyle \frac{\varphi \left(z\right)}{h}}-k\right)exp\left(-c{\left({\displaystyle \frac{\varphi \left(z\right)}{h}}-k\right)}^2\right),$$

(32)

with $c\in {\mathbb{R}}^{+}$ a constant and $\varphi$ denoting a conformal map. This type of approximation allows to represent a function $f\left(z\right)$ on an arc $\Gamma$ with an exponential decaying accuracy [12]. As demonstrated in [12] the approximation works effective for analytic functions. The definition of the Sinc-Gaussian basis follows next

$$\begin{array}{ccc}\hfill B_{S/G}& =& {\left\{S(k,h)\circ \varphi \left(z\right)exp\left(-c{\left({\displaystyle \frac{\varphi \left(z\right)-kh}{h}}\right)}^2\right)\right\}}_{k=-M}^N\hfill \\ & =& {\{\mathfrak{S}(k,h,c)\circ \varphi \left(z\right)\}}_{k=-M}^N\hfill \end{array}$$

(33)

allowing us to write the approximation in a compact form as

$$C_{h,c,M,N}\left[f\right]\left(z\right)=\sum _{k=-M}^{N}f\left(z_k\right)\mathfrak{S}(k,h,c)\circ \varphi \left(z\right).$$

(34)

Here, $z_k=\psi \left(kh\right)={\varphi }^{-1}\left(kh\right)$ are the discrete points based on Sinc points $kh$. Note that the formulae also apply when single exponential conformal maps (SE) are replaced with double exponential maps (DE).

The following theorem states that both Sinc and Sinc-Gaussian approximations provide a bounded error.

Theorem 1 

(Sinc Approximation [13]. ). Let $f\in {\mathit{L}}_{\alpha ,\beta }\left(\mathcal{D}\right)$ for $\alpha >0$ and $\beta >0$, take M=[β N/α], where [x] denotes the greatest integer in x, and then set $m=M+N+1$. If $f\in {\mathit{M}}_{\alpha ,\beta }\left(\mathcal{D}\right)$, and if $h={(\pi d/\left(\beta N\right))}^{1/2}$ then there exists a positive constant $K_1$ and $k_1$ independent of N, such that

$${∥f\left(x\right)-C_{h,c,M,N}\left[f\right]\left(x\right)∥}_2=K_1{N}^{1/2}exp\left(-k_1{N}^{1/2}\right),$$

(35)

with $w_k$ the base function (see (25)). $K_1$ and $k_1$ are constants independent of N.

The proof of this theorem for the pure Sinc case is presented in [2], whereas [12] examines the $S/G$ case. The choice $h={(\pi d/\left(\beta N\right))}^{1/2}$ is near to optimal for an approximation in the space ${\mathit{M}}_{\alpha ,\beta }\left(\mathcal{D}\right)$, in the sense that the error bound in Theorem 1 cannot be significantly improved regardless of the basis [13]. It is also optimum in the sense that the Lebesgue measure achieves an optimal value lower than the Chebyshev approximations [14].

3.1.2. Sinc Indefinite Integral

This subsection describes how indefinite integrals can be numerically defined [13] and how these definitions are related to the definition of definite integrals. For collocating an indefinite integral and for obtaining explicit approximations of the functions $\left({\mathcal{J}}^{+}f\right)\left(x\right)$ and $\left({\mathcal{J}}^{-}f\right)\left(x\right)$ defined by

$$\left({\mathcal{J}}^{+}f\right)\left(x\right)={\int }_a^{x}f\left(t\right)dt,\mathrm{with}x\in (a,b),$$

(36)

and

$$\left({\mathcal{J}}^{-}f\right)\left(x\right)={\int }_x^{b}f\left(t\right)dt,\mathrm{with}x\in (a,b).$$

(37)

We use the following basic relations [13]. Let $\mathbb{Z}$ denote the set of all integers, and let $\mathbb{C}$ denote the complex plane. Let $\mathrm{Sin}\mathrm{c}\left(x\right)$ be given by $sin\left(\pi x\right)/\left(\pi x\right)$ and $e_k$ be defined by using

$${\sigma }_k={\int }_0^k\mathrm{Sin}\mathrm{c}\left(x\right)dx={\displaystyle \frac{1}{\pi }}\mathrm{Si}\left(\pi k\right),$$

(38)

with $\mathrm{Si}\left(x\right)$ the sine integral. This put us into position to write

$$e_k={\displaystyle \frac{1}{2}}+{\sigma }_k,k\in \mathbb{Z}.$$

(39)

If now $\varphi$ denotes a one–to–one transformation of the interval $(a,b)$ onto the real line $\mathbb{R}$, let h denote a fixed positive number, and let the Sinc points be defined on $(a,b)$ by $z_k={\varphi }^{-1}\left(kh\right)$, $k\in \mathbb{Z}$, where ${\varphi }^{-1}$ denotes the inverse function of the conformal map $\varphi$. Let M and N be positive integers, set $m=M+N+1$, and for a given function f defined on $(a,b)$, define a diagonal matrix $D\left(f\right)$ by $D\left(f\right)=\mathrm{diag}\left[f\left(z_{-M}\right),\dots ,f\left(z_N\right)\right]$. Let $I^{(-1)}$ be a square Töplitz matrix of order m having $e_{i-j}$, as its ${(i,j)}^{\mathrm{th}}$ entry, so that

$${\left(I^{(-1)}\right)}_{i,j}=e_{i-j},\mathrm{with}i,j=-M,\dots ,N.$$

(40)

Define square matrices $A_m^{+}$ and $A_m^{-}$ by

$$A_m^{+}=h{I}^{(-1)}D(1/{\varphi }^{\prime }),$$

(41)

$$A_m^{-}=h{\left\{I^{(-1)}\right\}}^{T}D(1/{\varphi }^{\prime }),$$

(42)

where the superscript “T” denotes the transpose. The collocated representation of the indefinite integrals (36) and (37) are thus given by

$$\begin{array}{ccc}\hfill \left({\mathcal{J}}^{+}u\right)\left(x\right)& \approx & {\mathcal{J}}_m^{+}u={\mathit{V}}_m\left(\mathfrak{B}\right)A_m^{+}{\mathit{V}}_m\left(u\right)\hfill \\ & =& h{\mathit{V}}_m\left(\mathfrak{B}\right)I^{(-1)}D(1/{\varphi }^{\prime }){\mathit{V}}_m\left(u\right),\hfill \end{array}$$

(43)

and

$$\begin{array}{ccc}\hfill \left({\mathcal{J}}^{-}u\right)\left(x\right)& \approx & {\mathcal{J}}_m^{-}u={\mathit{V}}_m\left(\mathfrak{B}\right)A_m^{-}{\mathit{V}}_m\left(u\right)\hfill \\ & =& h{\mathit{V}}_m\left(\mathfrak{B}\right){\left\{I^{(-1)}\right\}}^{T}D(1/{\varphi }^{\prime }){\mathit{V}}_m\left(u\right).\hfill \end{array}$$

(44)

These are collocated representations of the indefinite integrals defined in (36) and (37), respectively [3].

The following theorem then enables us to collocate (linear or non-linear, non-stiff or stiff) initial value problems over an interval or an arc.

Theorem 2 

(Indefinite Integral Approximation [3]. ). Let $f\in {\mathit{L}}_{\alpha ,\beta }\left(\mathcal{D}\right)$ for $\alpha >0$ and $\beta >0$, take M=[β N/α], where [x] denotes the greatest integer in x, and then set $m=M+N+1$. If $f\in {\mathit{M}}_{\alpha ,\beta }\left(\mathcal{D}\right)$, and if $h={(\pi d/\left(\beta N\right))}^{1/2}$ then there exists a positive constant $K_1$ and $k_1$ independent of N,such that

$${∥{\mathcal{J}}^{+}f-{\mathcal{J}}_m^{+}f∥}_2=K_1{N}^{1/2}exp\left(-k_1{N}^{1/2}\right),$$

(45)

and

$${∥{\mathcal{J}}^{-}f-{\mathcal{J}}_m^{-}f∥}_2=K_1{N}^{1/2}exp\left(-k_1{N}^{1/2}\right).$$

(46)

$K_1$ and $k_1$ are constants independent of N.

Theorem 2 states that in both cases of the operator $\left({\mathcal{J}}^{\pm }\right)$ the error decay is the same. In addition these are collocated representations of the indefinite integrals defined in (36) and (37), respectively.

3.1.3. Sinc-Exp Indefinite Integration

In this section we pose the question how to approximate indefinite integrals on a sub-domain of $\mathbb{R}$. The approximation will use our basis system introduced in Section 3.1.1. We will introduce in this section the second version of the basis systems; i.e. a Sinc-Gaussian basis. The Sinc Gaussian basis includes the Sinc basis by setting $c=0$, it turns out that for both basis systems, we are able to get an approximation converging exponentially. Specifically we are interested in indefinite integrals of the two types

$${\mathcal{J}}^{+}\left(f\right)={\int }_{-\infty }^{x}f\left(t\right)dt,$$

(47)

and

$${\mathcal{J}}^{-}\left(f\right)={\int }_x^{\infty }f\left(t\right)dt.$$

(48)

If the function f is approximated by one of the approximations given in Section 3.1.1, we write for ${\mathcal{J}}^{+}\left(f\right)$,

$$\begin{array}{ccc}\hfill {\mathcal{J}}^{+}\left(f\right)& =& {\int }_{-\infty }^{x}f\left(t\right)dt\approx {\mathcal{J}}_{h,c,M,N}\left[f\right]\left(x\right)\hfill \\ & =& {\int }_{-\infty }^x\sum _{k=-M}^{N}f\left(t_k\right)B(k,h,c)\circ \varphi \left(t\right)dt\hfill \\ & =& \sum _{k=-M}^{N}f\left(t_k\right){\int }_{-\infty }^{x}B(k,h,c)\circ \varphi \left(t\right)dt\hfill \\ & =& \sum _{k=-M}^{N}f\left(t_k\right)\left\{\begin{array}{cc}{\int }_{-\infty }^{x}S(k,h)\circ \varphi \left(t\right)dt& \mathrm{if}c=0\\ {\int }_{-\infty }^{x}S(k,h)\circ \varphi \left(t\right)exp\left(-c{\left({\displaystyle \frac{\varphi \left(t\right)-kh}{h}}\right)}^2\right)dt& \mathrm{if}c\ne 0\end{array}.\right.\hfill \end{array}$$

(49)

Scaling the variable $\xi =t/h$ and collocating the expression with respect to x, we end up with the representation

$${\mathcal{J}}^{+}\left(f\right)\approx h\sum _{k=-M}^{N}f\left(t_k\right)\left\{\begin{array}{cc}{\int }_{-\infty }^{x/h}\mathrm{Sin}\mathrm{c}(\varphi \left(\xi \right)-k)d\xi & \mathrm{if}c=0\\ {\int }_{-\infty }^{x/h}\mathrm{Sin}\mathrm{c}(\varphi \left(\xi \right)-k)exp\left(-c{(\varphi \left(\xi \right)-k)}^2\right)d\xi & \mathrm{if}c\ne 0\end{array}\right.$$

(50)

The collocation of the variable x delivers an integral $I_{j,k}^c$ which will be our target in the next steps.

$$I_{j,k}^c={\int }_{-\infty }^j{\displaystyle \frac{sin\left(\pi \right(x-k\left)\right)}{\pi (x-k)}}{e}^{-c{(x-k)}^2}dx.$$

(51)

Note with $c=0$ the integral simplifies to the expression $I_{j,k}^0={\int }_{-\infty }^j\mathrm{Sin}\mathrm{c}(x-k)dx$. The discrete approximation of the integral ${\mathcal{J}}^{+}\left(f\right)$ thus becomes

$${\mathcal{J}}^{+}\left(f\right)\approx h\sum _{k=-M}^{N}f\left(t_k\right){\displaystyle \frac{1}{{\varphi }^{\prime }\left(x_k\right)}}{I}_{j,k}^c={\mathcal{J}}_j^{+},$$

(52)

with $x_k=\psi \left(kh\right)$, and $-M\le j\le N$, and delivers the approximation of ${\mathcal{J}}^{+}\left(f\right)$ via

$${\mathcal{J}}^{+}\left(f\right)\approx {\mathcal{J}}_{h,c,M,N}\left[f\right]\left(x\right)=\sum _{j=-M}^N{\mathcal{J}}_j^{+}B(j,h,c)\circ \varphi \left(x\right).$$

(53)

For the discrete approximation, we need to know the value of $I_{j,k}^c$ to be able to find the approximation of the indefinite representation

$${\mathcal{J}}_j^{+}\left(f\right)=h\sum _{k=-M}^{N}f\left(t_k\right){\displaystyle \frac{I_{j,k}^c}{{\varphi }^{\prime }\left(x_k\right)}}=h\left[I_{j,k}^c\right]D(1/{\varphi }^{\prime })V\left(f\right).$$

(54)

The matrix $I_{j,k}^c$ can be written as

$$\begin{array}{ccc}\hfill I_{j,k}^c& =& {\int }_{-\infty }^j{\displaystyle \frac{sin\left(\pi \right(x-k\left)\right)}{\pi (x-k)}}{e}^{-c{(x-k)}^2}dx\hfill \\ & =& {\int }_{-\infty }^0{\displaystyle \frac{sin\left(\pi \right(x-k\left)\right)}{\pi (x-k)}}{e}^{-c{(x-k)}^2}dx+{\int }_0^j{\displaystyle \frac{sin\left(\pi \right(x-k\left)\right)}{\pi (x-k)}}{e}^{-c{(x-k)}^2}dx\hfill \\ & =& {\int }_{-\infty }^{j-k}{\displaystyle \frac{sin\left(\pi y\right)}{\pi y}}{e}^{-c{y}^2}dy=I_{j-k,0}^c,\hfill \end{array}$$

(55)

reducing the problem to the structure of a Töplitz matrix if the indefinite integral has some finite values. To get the values of the integrals we divide the integration domain into two parts

$$\begin{array}{ccc}\hfill I_{j-k,0}^c& =& {\int }_{-\infty }^{j-k}{\displaystyle \frac{sin\left(\pi y\right)}{\pi y}}{e}^{-c{y}^2}dy\hfill \\ & =& {\int }_{-\infty }^0{\displaystyle \frac{sin\left(\pi y\right)}{\pi y}}{e}^{-c{y}^2}dy+{\int }_0^{j-k}{\displaystyle \frac{sin\left(\pi y\right)}{\pi y}}{e}^{-c{y}^2}dy\hfill \\ & =& I_0^c+I_{j-k}^c.\hfill \end{array}$$

(56)

The two integrals deliver in straight forward integration the following

$$I_0^c=\left\{\begin{array}{cc}{\displaystyle \frac{1}{2}}& \mathrm{if}c=0\\ \mathrm{erf}\left({\displaystyle \frac{\pi }{2\sqrt{c}}}\right)& \mathrm{if}c\ne 0\end{array}\right.,$$

(57)

and

$$I_{j-k}^c=\left\{\begin{array}{cc}\mathrm{Si}(j-k)& \mathrm{if}c=0\\ \mathrm{sign}(j-k){\sigma }_{|j-k|}^c& \mathrm{if}c\ne 0\end{array}\right.,$$

(58)

with ${\sigma }_{k+1}^c={\sigma }_k^c+f_k^c$, and $f_k^c={\int }_0^1\mathrm{Sin}\mathrm{c}(y+k)exp\left(-c{(y+k)}^2\right)\mathrm{dy}$, thus ${\sigma }_n^c={\sum }_{k=0}^{n-1}{f}_k^c$. A similar procedure can be applied to the integral ${\mathcal{J}}^{-}\left(f\right)$. The difference is a transposition of the matrix $I_{j-k}^c$ in the representation of approximation. The following Theorem summarizes the results.

Theorem 3 

(Indefinite Integral Approximation.). If ϕ denotes a one–to–one transformation of the interval $(a,b)$ onto the real line $\mathbb{R}$, let h denote a fixed positive number, and let the Sinc points be defined on $(a,b)$ by $z_k={\varphi }^{-1}\left(kh\right)$, $k\in \mathbb{Z}$, where ${\varphi }^{-1}$ denotes the inverse function of the conformal map ϕ. Let M and N be positive integers, set $m=M+N+1$, and for a given function f defined on $(a,b)$, define the vector $\mathit{V}\left(f\right)=\left(f\left(z_{-M}\right),\dots ,f\left(z_N\right)\right)$, and a diagonal matrix $D\left(f\right)$ by $D\left(f\right)=\mathit{diag}\left[f\left(z_{-M}\right),\dots ,f\left(z_N\right)\right]$. Let $\mathit{V}\left(B\right)=\left(B_{-M},\dots ,B_N\right)$ the vector of basis functions, and let $I^{(-1)}$ be a square Töplitz matrix of order m having $I_{i-j}^c$, as its ${(i,j)}^{\mathrm{th}}$ entry, $i,j=-M,\dots ,N$,

$${\left(I^{(-1)}\right)}_{i,j}=I_{i-j,0}^c,\mathrm{with}i,j=-M,\dots ,N.$$

(59)

Define square matrices $A_m^{+}$ and $A_m^{-}$ by

$$A_m^{+}=h{I}^{(-1)}D(1/{\varphi }^{\prime }),$$

(60)

$$A_m^{-}=h{\left\{I^{(-1)}\right\}}^{T}D(1/{\varphi }^{\prime }),$$

(61)

where the superscript “T” denotes the transpose. Then the indefinite integrals (47) and (48) are approximated by

$${\mathcal{J}}^{+}\left(f\right)\approx {\mathcal{J}}_m^{+}f=\mathit{V}\left(\mathfrak{B}\right)A_m^{+}\mathit{V}\left(f\right),$$

(62)

and

$${\mathcal{J}}^{-}\left(f\right)\approx {\mathcal{J}}_m^{-}f=\mathit{V}\left(\mathfrak{B}\right)A_m^{-}{\mathit{V}}_m\left(f\right).$$

(63)

The error of this approximation was estimated for the pure Sinc case $c=0$ in [2] as

$$E_N={∥\mathcal{J}\left(f\right)-{\mathcal{J}}_{h,c,M,N}\left[f\right]\left(x\right)∥}_2\sim K_1{N}^{1/2}exp\left(-k_1{N}^{1/2}\right),$$

(64)

where $K_1$ and $k_1$ are constants independent of N.

Note the matrices $A_m^{+}$ and $A_m^{-}$ have eigenvalues with $\Re \left({\lambda }_i\right)>0$ which guaranties the convergence of the solution. A proof for the pure Sinc case for the matrix $I^{(-1)}$ was recently given by Han and Xu [15].

If we use the properties of the matrices $A_m^{+}$ and $A_m^{-}$ defined above we know that

$$A_m^{+}=h{I}^{(-1)}D(1/{\varphi }^{\prime })=h\left(I_0^c\widehat{I}+\left[I_{j-k}^c\right]\right)D(1/{\varphi }^{\prime }),$$

(65)

where $\widehat{I}$ is a $m\times m$ matrix filled with 1’s and $\left[I_{j-k}^c\right]$ is antisymmetric. The elements of the diagonal matrix $D=D(1/{\varphi }^{\prime })$ are all positive according to their definition. If we assume that $A_m^{\pm }$ is diagonalizable then for each eigenvalue ${\lambda }_i\in \mathbb{C}$ there exists a complex valued eigenvector $w_i\in {\mathbb{C}}^m$. Since the eigenvalues of $A_m^{\pm }$ are the same as $D^{1/2}{A}_m^{\pm }{D}^{-1/2}$ we can write for $A_m^{+}$ for example

$$\begin{array}{ccc}\hfill w^{*}{D}^{1/2}{A}_m^{+}{D}^{-1/2}w& =& h{w}^{*}{D}^{1/2}{I}^{(-1)}{D}^{1/2}w\hfill \\ & =& h{w}^{*}{D}^{1/2}\left(I_0^c\widehat{I}+\left[I_{j-k}^c\right]\right)D^{1/2}w\hfill \\ & =& {\lambda \parallel w\parallel }^2,\hfill \end{array}$$

(66)

which is equivalent to

$$\left.I_0^c{w}^{*}{D}^{1/2}\widehat{I}{D}^{1/2}w+w^{*}{D}^{1/2}\left[I_{j-k}^c\right]\right)D^{1/2}w={\displaystyle \frac{\lambda }{h}}{\parallel w\parallel }^2,$$

(67)

or

$$I_0^c\underset{W^{*}}{\underbrace{\left(\sum _{k=-M}^N{w}_k^{*}{d}_k^{1/2}\right)}}\left(\sum _{k=-M}^N{w}_k{d}_k^{1/2}\right)+w^{*}{D}^{1/2}\left[I_{j-k}^c\right]D^{1/2}w={\displaystyle \frac{\lambda }{h}}{\parallel w\parallel }^2,$$

(68)

which is

$$I_0^c{W}^{*}W+w^{*}{D}^{1/2}\left[I_{j-k}^c\right]D^{1/2}w={\displaystyle \frac{\lambda }{h}}{\parallel w\parallel }^2,$$

(69)

in other terms

$$I_0^c{\parallel W\parallel }^2+w^{*}{D}^{1/2}\left[I_{j-k}^c\right]D^{1/2}w={\displaystyle \frac{\lambda }{h}}{\parallel w\parallel }^2.$$

(70)

The real part of this expression delivers

$$=\Re \left({\displaystyle \frac{\lambda }{h}}\right)=\Re \left(I_0^c{\displaystyle \frac{{\parallel W\parallel }^2}{{\parallel w\parallel }^2}}\right)+\underset{=0}{\underbrace{\mathrm{Re}\left(w^{*}{D}^{1/2}\left[I_{j-k}^c\right]D^{1/2}w\right)}},$$

(71)

the second term vanishes because the matrix $\left[I_{j-k}^c\right]$ is antisymmetric and thus since $I_0^c>0$ and $h>0$ we have

$$\Re \left(\lambda \right)=h{I}_0^c{\displaystyle \frac{{\parallel W\parallel }^2}{{\parallel w\parallel }^2}}\ge 0.$$

(72)

If we assume that the norm is defined on the vector space, we can use the Rayleigh quotient of a matrix H given by $R_H\left(w\right)=w^{*}Hw/\left(w^{*}w\right)$ to bound the eigenvalues. According to [16,17] there exist a minimal and maximal eigenvalue defined by

$${\lambda }_m=\underset{w}{min}{R}_H\left(w\right),$$

(73)

and

$${\lambda }_M=\underset{w}{max}{R}_H\left(w\right).$$

(74)

Based on Grenander and Szegö theory [18], we may conclude that the eigenvalues of $A_m$ are strictly positive for finite matrix sizes.

3.2. Methods Using Orthogonal Polynomials

In this part, we will largely focus on approximation via interpolation at the roots $x_j$ of orthogonal polynomials ${\Pi }_n$. These polynomials include Legendre polynomials, Chebyshev polynomials, Gegenbauer polynomials, and Jacobi polynomials, which are initially defined on $(-1,1)$. We focus upon studying on the classical set of orthogonal polynomials. For information on other classes and their computations, refer to Gautschi [19]. The following presentation includes polynomials such as Laguerre- and Hermite polynomials defined on $(0,\infty )$ and $(-\infty ,\infty )$.

For example we may categorize the orthogonal polynomials as:

1.

For Legendre polynomials, $x_j$ are the n zeros of the Legendre polynomial $P_n\left(x\right)$, which are orthogonal on the interval $(a,b)=(-1,1)$ with regard to the weight function w, which is identically 1 on $(-1,1)$.

2.

For the case of Hermite polynomial interpolation, using the Hermite polynomials $H_n\left(x\right)$ that are orthogonal over $\mathbb{R}$ with respect to the weight function w , with $w\left(x\right)=exp\left(-x^2\right)$, and with $H_n\left(x_j\right)=0$ for $j=1,...,n$; and

3.

Other polynomials that are orthogonal with respect to a weight function, such as Jacobi polynomials, Gegenbauer polynomials, etc.

We should remark that the approach presented here, indefinite integrals via root-based Lagrangian interpolation, is novel.

Let ${\Pi }_n$ be a set of polynomials with degrees less than or equal to n, and $f\left(x\right)$ be a smooth entire or analytic function defined on the interval $(a,b)$. For the function $f\left(x\right)$, the interpolation polynomial $p_n\left(x\right)$ fits the following requirements.

$$p_{n-1}\in {\Pi }_{n,}$$

(75)

for n different root positions $x_j$ in $(a,b)$. The discrete points $x_1$, $\dots ,x_n\in (a,b)$ are taken as roots of orthogonal polynomials, as indicated before. The resulting polynomial of order $n-1$ meets the interpolation condition: $p_{n-1}\left(x_i\right)=f\left(x_i\right)$ for $i=1,\dots ,n$. The scenario is analogous to Lagrange polynomials.

A Lagrange polynomial approximation of order $n-1$ is:

$$f\left(x\right)\approx p_{n-1}\left(x\right)=\sum _{i=1}^{n}f\left(x_i\right){\ell }_i\left(x\right),$$

(76)

where

$${\ell }_i\left(x\right)=\prod _{\underset{j\ne i}{j=1}}^n{\displaystyle \frac{x-x_j}{x_i-x_j}}={\displaystyle \frac{{\omega }_n\left(x\right)}{{\omega }_n^{\prime }\left(x_i\right)\left(x-x_i\right)}}\mathrm{for}i=1,\dots ,n,$$

(77)

are the discrete values based Lagrange polynomials employing $x_i$ with

$${\omega }_n\left(x\right)=\left(x-x_1\right)\dots \left(x-x_{n-1}\right)\left(x-x_n\right)=\prod _{i=1}^n\left(x-x_i\right).$$

(78)

The Lagrange polynomials satisfy the elementary interpolation requirement, as shown by

$${\ell }_i\left(x_j\right)=\left\{\begin{array}{cc}1\hfill & j=i,\hfill \\ 0\hfill & j\ne i.\hfill \end{array}\right.$$

(79)

We can now estimate the global approximation error for an entire function f. The infimum of the absolute difference between f and the approximating $n^{\mathrm{th}}$ order polynomial $p_n$ is roughly given as

$$E_n=\underset{p\in {\Pi }_n}{inf}\left|f\left(x\right)-p_n\left(x\right)\right|\simeq {\displaystyle \frac{M\left(\rho \right)}{{\rho }^r}}.$$

(80)

If f is entire, then $\rho \simeq 2r$; e.g. if $f\left(x\right)=e^x$ then $M\left(\rho \right)=e^{2r}$ which delivers

$$\underset{p\in {\Pi }_n}{inf}\underset{-1<x<1}{sup}\left|exp\left(x\right)-p_n\left(x\right)\right|\simeq {\displaystyle \frac{e^{2r}}{{\left(2r\right)}^n}}.$$

(81)

The maximum occurs when we differentiate with respect to r and put the result equal to zero.

$${\partial }_r\left(\alpha {\displaystyle \frac{e^{2r}}{{\left(2r\right)}^n}}\right)=-2^{-n}{e}^{2r}n{r}^{-1-n}\alpha +2^{1-n}{e}^{2r}{r}^{-n}\alpha =0$$

(82)

Solving the equation with respect to n yields a relationship between r and n which is $n=2r$. The solution for r then delivers $r=n/2.$ Finally, this demonstrates that there is an error that decays as $E_n=\alpha n^{-n}{e}^n$. This relationship confirms that the error relation becomes

$$\underset{p\in {\Pi }_n}{inf}\left|f\left(x\right)-p_n\left(x\right)\right|\simeq {\displaystyle \frac{M\left(\rho \right)}{{\rho }^r}}=\alpha n^{-n}{e}^n=E_n.$$

(83)

The error thus comprises two important aspects. $E_n=t_1{t}_2$ with $t_1=e^n$, which is increasing; $t_2=n^{-n}$, which is decreasing. The maximum of the expression $E_n$ is assumed when $n=1$. Because $n^{-n}$ decays faster than the exponential function, the overall error $E_n=\alpha n^{-n}{e}^n$ as a function of n decays faster than the exponential function $e^{-n}$.

On the other hand, it is well known that if a function f is analytic in a simply connected domain D in the complex plane and if a closed interval $[a,b]$ is in the interior of D, then we can approximate f on $[a,b]$ via a polynomial of degree n for which the error approaches zero at a rate of $\mathcal{O}(exp(-\beta n\left)\right)$, where $\beta$ is a positive constant. If the function f on D is entire, we are able to predict an error rate of $\mathcal{O}\left(n^{-n}exp\left(\beta n\right)\right)$ as discussed above. If none of the two are correct, the usual error rate for polynomial approximation is $\mathcal{O}\left(n^{-\beta }\right)$. These three error rates will serve in the following as basis for the indefinite integral approximation based on Lagrange polynomials.

Based on those results and the Lagrange approximation provided above, let us look at the following two indefinite integrals.

$${\mathcal{J}}^{+}\left(wf\right)\left(x\right)={\int }_a^{x}w\left(x\right)f\left(x\right)dx,$$

(84)

and

$${\mathcal{J}}^{-}\left(wf\right)\left(x\right)={\int }_x^{b}w\left(x\right)f\left(x\right)dx,$$

(85)

It will be approximated with the aforementioned relationships. As a result, we offer the following Definition.

Definition 3 

( Root Based Indefinite Integrals.). Consider generalized Lagrange interpolation at distinct points $x_j$ on $(a,b)$, with $a<x_1<x_2<\dots <x_m<b$, which takes the form

$$f\left(x\right)\approx \sum _{k=1}^{n}f\left(x_k\right){\sigma }_k\left(x\right),$$

(86)

using the Lagrange interpolation for orthogonal polynomials ${\pi }_m\left(x\right)\in {\Pi }_m$

$${\sigma }_k\left(x\right)={\displaystyle \frac{{\pi }_m\left(x\right)}{\left(x-x_k\right){\pi }_m^{\prime }\left(x_k\right)}},\mathrm{for}k=1,\dots ,m.$$

(87)

Referring to integral (84) and (85), we are able to define the entries $A_{j,k}^{\pm }$ of the matrices $A^{\pm }$ by the matrix elements as

$$A_{j,k}^{+}={\int }_a^{x_j}w\left(\xi \right){\sigma }_k\left(\xi \right)d\xi ,$$

(88)

and

$$A_{j,k}^{-}={\int }_{x_j}^{b}w\left(\xi \right){\sigma }_k\left(\xi \right)d\xi ,$$

(89)

where w is the weight function related to ${\pi }_m\left(x\right)\in \left\{P_m\left(x\right),T_m\left(x\right),U_m\left(x\right),C_m^{\left(\mu \right)}\left(x\right),P_m^{(\alpha ,\beta )}\left(x\right),\right.$$\left.L_m^{\alpha }\left(x\right),H_m\left(x\right),\dots \right\}={\Pi }_m$ which is the set of classical orthogonal polynomials (see Olver et al., pp. 438 [20]), that is positive a.e. on $(a,b)$, and is such that the moments ${\int }_a^{b}w\left(x\right)x^{j}dx$ exists for every non-negative integer j. The function w, where $w\left(x\right)>0$ a.e. on $(a,b)$ defines a sequence of orthogonal polynomials ${\left\{{\pi }_n\right\}}_{n=2}^{\infty }$, for which all n zeros of ${\pi }_n$ are located on $(a,b)$.

Now, setting

$$\mathit{V}\left(f\right)={\left(f\left(x_1\right),\dots ,f\left(x_n\right)\right)}^T,$$

(90)

and

$$\mathit{L}\left(x\right)={\left({\ell }_1\left(x\right),\dots ,{\ell }_n\left(x\right)\right)}^T,$$

(91)

where ${\ell }_i\left(x\right)$ are the basic Lagrange polynomials (77). We are defining ${\mathcal{J}}_n^{\pm }$ by

$${\mathcal{J}}^{\pm }\left(wf\right)\left(x\right)\approx \left({\mathcal{J}}_n^{\pm }wf\right)\left(x\right)=\mathit{L}\left(x\right)A^{\pm }\mathit{V}\left(f\right),$$

(92)

we obtain a practicable approximation of an indefinite integral using standard Lagrange interpolation. Assuming our function f is either entire, analytic, or none of the two, we obtain the error decay of the approximation provided by either the relation (83), exponential, or $\mathcal{O}\left(n^{-\beta }\right)$, respectively.

3.3. Generalized Lagrange Methods

Aitken and Neville suggested a generalization of the Lagrange interpolation in 1932 [21,22], followed by Mühlbach [23]. The numerical computations were performed using an iterative technique. Occorsio and Russo [24] recently employed an iterative technique, while Moalemi presented a generalization in analytic form in her PhD thesis [25]. Cheney and Light discuss a more general form of Lagrange interpolation [26]. With prescribed nodes $x_0$,…, $x_n$ in X, where X is an arbitrary set, they define

$${\nu }_j\left(x\right)=\prod _{\underset{k\ne j}{k=0}}^n{\displaystyle \frac{\phi \left(x,x_k\right)}{\phi \left(x_j,x_k\right)}},\mathrm{for}j=0,1,2,\dots ,n,$$

(93)

and $\phi :X\times X\to \mathbb{R}$ is a function such that $\phi (x,y)=0$ if and only if $x=y$. Therefore,

$${\nu }_j\left(x_k\right)={\delta }_{j,k}=\left\{\begin{array}{cc}1& k=j,\\ 0& k\ne j,\end{array}\right.$$

(94)

and the interpolating polynomial is defined by $p_n\left(x\right)={\sum }_{j=0}^{n}f\left(x_j\right){\nu }_j\left(x\right)=\mathit{V}\left(f\right).\mathit{V}\left(N\right)$ using $\mathit{V}\left(f\right)={\left(f\left(x_0\right),\dots ,f\left(x_n\right)\right)}^T$ and $\mathit{V}\left(N\right)=\left({\nu }_0\left(x\right),\dots ,{\nu }_n\left(x\right)\right.$ interpolates the function f.

We will briefly discuss the method used in [25] and then apply it to the set of classical orthogonal polynomials in this section.

Consider $B_{n-1}={\left\{b_i\left(x\right)\right\}}_{i=1}^{n-1}$ be a set of continuous functions defined on $(a,b)$ so that $b_i\left(x_l\right)\ne b_i\left(x_j\right)$ for $l\ne j$ and $X={\left\{x_k\right\}}_{k=1}^n$ represent n unique discretization locations in $(a,b)$. The generalized Lagrange function $B_k\left(x\right)$ for $k=1,2,\dots \mathrm{n}$ is defined as follows:

$$B_k\left(x\right)=\prod _{i,j\in {\mathcal{I}}_{n,k}}{\displaystyle \frac{b_i\left(x\right)-b_i\left(x_j\right)}{b_i\left(x_k\right)-b_i\left(x_j\right)}},\mathrm{for}k=1,2,\dots ,n,$$

(95)

using the index set given by:

$$\begin{array}{ccc}\hfill {\mathcal{I}}_{n,k}& =& \left.\left\{\right(\rho ,\rho )\right|\rho =1,2,\dots ,k-1\wedge k>1\}\cup \hfill \\ & & \left.\left.\left\{\right(\rho ,\rho +1)\right|\rho =k,k+1,\dots ,n-1\wedge k<n\right\}.\hfill \end{array}$$

(96)

The interpolation condition is directly verified in (95) by

$$B_k\left(x_i\right)={\delta }_{k,i},k,i=1,2,\dots \mathrm{n}.$$

(97)

This implies that the collection of basic functions $B_{n-1}={\left\{B_k\left(x\right)\right\}}_{k=1}^n$ represents a collection of linearly independent functions on $(a,b)$. In other terms:

$$\mathrm{span}\left\{B_1\left(x\right),\dots ,B_n\left(x\right)\right\}={\Pi }_n,$$

(98)

creates an interpolation space that enables for representation:

$$f\left(x\right)=\sum _{k=1}^{n}f\left(x_k\right)B_k\left(x\right),$$

(99)

satisfying

$$f\left(x_i\right)=\sum _{k=1}^{n}f\left(x_k\right)B_k\left(x_i\right)=\sum _{k=1}^{n}f\left(x_k\right){\delta }_{k,i},$$

(100)

ensuring a unique interpolation function, therefore establishing a Chebyshev system [27].

To construct the basis $B_{n-1}={\left\{B_k\left(x\right)\right\}}_{k=1}^n$, it is first necessary to have a sequence of discrete numbers $X={\left\{x_k\right\}}_{k=1}^n$, which can be any numbers satisfying $x_k\ne x_l$ for $k,l=1,2,\dots ,n$. This set is used in conjunction with another set of functions $B_{n-1}={\left\{b_i\left(x\right)\right\}}_{i=1}^{n-1}$ to form the generalized Lagrange basis $B_k\left(x\right)$ using (95). It is important that the interplay between the two sets X and $B_{n-1}$ is independent of the discrete function values $f\left(x_k\right)$.

For the specific interpolation we will develop, we are employing the classical normal orthogonal polynomials. For $P_m\left(x\right)$, we put $B_{n-1}={\left\{P_k\left(x\right)\right\}}_{k=1}^{n-1}$. Where $P_k\left(x\right)\in \left\{T_k\left(x\right),U_k\left(x\right),P_k\left(x\right),C_k^{\left(m\right)}\left(x\right),P_k^{(\alpha .\beta )}\left(x\right),L_k^{\left(\alpha \right)}\left(x\right),H_k\left(x\right),...\right\}$. For the discretization points, we utilize the roots of the polynomial $P_n\left(x\right)=0$. Note that the distribution of roots on $(a,b)$ is not equidistant.

Thus equation (95) takes the form:

$$B_k\left(x\right)=\prod _{i,j\in {\mathcal{I}}_{n,k}}{\displaystyle \frac{P_i\left(x\right)-P_i\left(x_j\right)}{P_i\left(x_k\right)-P_i\left(x_j\right)}},\mathrm{for}k=1,2,\dots ,n,$$

(101)

with roots $x_j$ satisfying $P_n\left(x_j\right)=0$ and $j=1,2,\dots ,n$. This allows to write the interpolation of $f\left(x\right)$ as

$$f\left(x\right)=p_{n-1}\left(x\right)=\sum _{k=1}^{n}f\left(x_k\right)B_k\left(x\right)=\mathit{V}\left(f\right)\mathit{V}\left(B\right).$$

(102)

This put us into position to introduce vectors for numeric representation as follows

$$\mathit{V}\left(f\right)={\left(f\left(x_1\right),\dots ,f\left(x_n\right)\right)}^T,$$

(103)

and

$$\mathit{V}\left(B\right)=\left(B_1\left(x\right),\dots ,B_n\left(x\right)\right),$$

(104)

delivering a direct representation of the interpolation polynomial $p_{n-1}\left(x\right)$ as an analytic representation

$$p_{n-1}\left(x\right)=\mathit{V}\left(f\right)\mathit{V}\left(B\right).$$

(105)

This formula is straight forward and can be implemented in a direct way.

Using the interpolating function (105), we may approximate an indefinite integral ${\mathcal{J}}^{\pm }\left(f\right)$ as follows:

$$\begin{array}{ccc}\hfill {\mathcal{J}}^{+}\left(f\right)\left(x\right)& =& {\int }_a^{x}f\left(x\right)dx\approx {\mathcal{J}}_n^{+}\left(f\right)\left(x\right)\hfill \\ & =& {\int }_a^x\sum _{k=1}^{n}f\left(x_k\right)B_k\left(x\right)dx\hfill \\ & \approx & \sum _{i=1}^n{B}_i\left(x\right){\int }_a^{x_i}\sum _{k=1}^{n}f\left(x_k\right)B_k\left(x\right)dx\hfill \\ & =& \sum _{i=1}^n{B}_i\left(x\right)\sum _{k=1}^{n}f\left(x_k\right){\int }_a^{x_i}{B}_k\left(x\right)dx\hfill \\ & =& \sum _{i=1}^n{B}_i\left(x\right)\sum _{k=1}^{n}f\left(x_k\right)A_{k,i}^{+}\hfill \\ & =& \sum _{i=1}^n{B}_i\left(x\right){\mathcal{I}}_i^{+},\hfill \end{array}$$

(106)

with

$$A_{k,i}^{+}={\int }_a^{x_i}{B}_k\left(x\right)dx,k,i=1,2,\dots ,n.$$

(107)

The same steps can be used to approximate ${\mathcal{J}}^{-}\left(f\right)\left(x\right)$ by:

$$\begin{array}{ccc}\hfill {\mathcal{J}}^{-}\left(f\right)\left(x\right)& =& {\int }_x^{b}f\left(x\right)dx\approx {\mathcal{J}}_n^{-}\left(f\right)\left(x\right)\hfill \\ & =& {\int }_x^b\sum _{k=1}^{n}f\left(x_k\right)B_k\left(x\right)dx\hfill \\ & \approx & \sum _{i=1}^n{B}_i\left(x\right){\int }_{x_i}^b\sum _{k=1}^{n}f\left(x_k\right)B_k\left(x\right)dx\hfill \\ & =& \sum _{i=1}^n{B}_i\left(x\right)\sum _{k=1}^{n}f\left(x_k\right){\int }_{x_i}^b{B}_k\left(x\right)dx\hfill \\ & =& \sum _{i=1}^n{B}_i\left(x\right)\sum _{k=1}^{n}f\left(x_k\right)A_{k,i}^{-}\hfill \\ & =& \sum _{i=1}^n{B}_i\left(x\right){\mathcal{I}}_i^{-},\hfill \end{array}$$

(108)

with

$$A_{k,i}^{-}={\int }_{x_i}^b{B}_k\left(x\right)dx,k,i=1,2,\dots ,n.$$

(109)

In addition we have the approximation of a quadrature rule as

$$\begin{array}{ccc}\hfill \mathcal{J}\left(f\right)\left(x\right)& =& {\mathcal{J}}^{+}\left(f\right)\left(x\right)+{\mathcal{J}}^{-}\left(f\right)\left(x\right)\hfill \\ & \approx & \sum _{i=1}^n{B}_i\left(x\right)\sum _{k=1}^{n}f\left(x_k\right)\left(A_{k,i}^{+}+A_{k,i}^{-}\right)\hfill \\ & =& \sum _{i=1}^n{B}_i\left(x\right)\sum _{k=1}^{n}f\left(x_k\right)A_{k,i},\hfill \end{array}$$

(110)

delivering

$$\mathcal{J}\left(f\right)\left(x\right)\approx {\mathcal{J}}_n\left(f\right)\left(x\right)=\sum _{i=1}^n{B}_i\left(x\right)\sum _{k=1}^{n}f\left(x_k\right)A_{k,i}.$$

(111)

Note that the selection of the orthogonal polynomial type determines the coefficient matrix $A_{k,i}$.

Now we are in position to set

$$\mathit{V}\left(f\right)={\left(f\left(x_1\right),\dots ,f\left(x_n\right)\right)}^T,$$

(112)

and

$$\mathit{V}\left(B\right)=\left(B_1\left(x\right),\dots ,B_n\left(x\right)\right),$$

(113)

so that an approximation of ${\mathcal{J}}^{\pm }$ can be defined as ${\mathcal{J}}_n^{\pm }$ by

$${\mathcal{J}}^{\pm }\left(f\right)\left(x\right)\approx \left({\mathcal{J}}_n^{\pm }f\right)\left(x\right)=\mathit{V}\left(B\right)A^{\pm }\mathit{V}\left(f\right).$$

(114)

Thus, we obtain a practicable approximation of an indefinite integral using generalized Lagrange interpolation. In addition we find

$$\begin{array}{ccc}\hfill {\mathcal{J}}^{+}\left(f\right)\left(x\right)+{\mathcal{J}}^{-}\left(f\right)\left(x\right)& =& \mathcal{J}\left(f\right)\left(x\right)\hfill \\ & \approx & \left({\mathcal{J}}_n^{+}f\right)\left(x\right)+\left({\mathcal{J}}_n^{-}f\right)\left(x\right)\hfill \\ & =& {\mathcal{J}}_n\left(f\right)\left(x\right)\hfill \\ & =& \mathit{V}\left(B\right)\left(A^{+}+A^{-}\right)\mathit{V}\left(f\right),\hfill \end{array}$$

(115)

We may now use the introduced matrices and vectors to formulate the implicit Runge Kutta technique, which provides an approximation to an initial value problem.

4. Runge Kutta Method

So far, we’ve demonstrated that indefinite integration techniques are defined by the matrices $A_{i,j}^{+}$. Because the matrices’ representations are based on distinct underlying basis systems and discretization methods, it is only appropriate to refer to various Runge-Kutta procedures.

Given the discretization $\mathit{V}\left(t\right)=\left(t_0,t_1..,t_N\right)$, the matrix $A_{i,j}^{+}$ defined as in (60), (65), (88), or (107), and weights $b_i$ based on (14) the one-step Runge-Kutta iteration (12) can be transformed to an iterative multistep formula using the following iterations:

$${\mathit{V}}^{(k+1)}\left(k\right)=f\left(t_0+t_{k}h,y_k+h{A}^{+}.{\mathit{V}}^{(k+1)}\left(k\right)\right).$$

(116)

$$y_{k+1}=y_k+h\mathit{V}\left(b\right).{\mathit{V}}^{(k+1)}\left(k\right),\mathrm{with}k=0,1,2,\dots \mathrm{N},$$

(117)

where $h=t_{k+1}-t_k$ is the individual step length. The notation $V^{\left(k\right)}(•)$ denotes the kth iteration step and $V^{\left(1\right)}(•)=V(•)$.

The resultant discrete $N+1$ solution points $\left(t_k,y\left(t_k\right)\right)=\left(t_k,y_k\right)$, with $k=0,1,\dots ,N$ may now be used in an interpolation to get a representation of order N as an analytical function of the solution to the initial value problem:

$$y\left(t\right)=y_0+\sum _{k=1}^{N}y\left(t_k\right)B_k\left(t\right)=y_0+\mathit{V}\left(y\right).\mathit{B}\left(t\right),$$

(118)

where $\mathit{B}\left(t\right)={\left(B_1\left(t\right),\dots ,B_N\left(t\right)\right)}^T$ is the set of basis functions which are either Sinc, Sinc-Exp functions, Lagrange, or generalized Lagrange polynomials, respectively.

5. Stability

In previous sections, we looked at the convergence and error behavior of the Sinc-, L-, and gL-RK algorithms. The $A^{+}$ matrix of indefinite integration serves as the basis for all three approaches. In this section, we’d like to discuss another element of these RK methods: their stability. The geometric features of the RK approaches can be exploited to do this. We shall employ the stability function $R\left(z\right)$, often known as the Padé approximation.

These findings are based on a stability study by Dahlquist [28]. A linear stability analysis assumes that any initial value problem,

$$y^{\prime }=f(t,y),\mathrm{with}y\left(0\right)=y_0,$$

(119)

may be linearized to

$$y^{\prime }=\lambda y,$$

(120)

which represents the behavior of stiff systems in the vicinity of a stationary point. The eigenvalues of the Jacobian matrix determined at the stability point are denoted as $\lambda \in \mathbb{C}$. If the real parts of the $\lambda$’s are negative we speak of stable behavior; i.e. the equation is stable in the sense of Lyapunov when $\Re \left(\lambda \right)\le 0$. It is widely known that the solutions to the differential system in (119) decline exponentially as $t\to \infty$. Applying a RK one-step approach which in general is a relation of the form $y_{n+1}={\mathsf{\Phi }}_h\left(y_n\right)$ to this equation results to

$$y_{n+1}=R\left(z\right)y_n=R{\left(z\right)}^n{y}_0,$$

(121)

where $z=\lambda h$. The RK method’s stability function is denoted as $R\left(z\right)$. This means necessarily that if $\left|R\right(z\left)\right|\le 1$, we have stability. The stability domain, defined as $\mathcal{S}=\{z\in \mathbb{C}:|R\left(z\right)|\le 1\}$, is the set of points in the complex plane where the calculated solution stays bounded after several steps of computing. To meet the stability condition, all eigenvalues $\lambda$, multiplied by step size h, must be in $\mathcal{S}$. If $z=h\lambda$, the solution of (120) will be $exp\left(z\right)$ in a single h-step operation. In the same time span, the approximate solution computed with a RK method will be multiplied by a z-specific function $R\left(z\right)$. Values of z in the left half-plane are of particular importance since the precise solution is bounded in this instance, and proper problem modelling requires the computed solution to behave similarly.

As discussed above the following tableau defines an s-order RK method

(122)

In contrast to the literature, we shall refer to it as an s-order RK method rather than an s-stage RK method because here $A^{+}$ is always of order $s\times s$. The vector $V\left(y\right)$, consisting of the s order values, satisfies

$$V\left(y\right)=\mathbf{1}{y}_0+h\lambda A^{+}V\left(y\right)=\mathbf{1}{y}_0+z{A}^{+}V\left(y\right),$$

(123)

where $y_0$ is the initial approximation. It follows that

$$V\left(y\right)={\left(I-z{A}^{+}\right)}^{-1}\mathbf{1}{y}_0.$$

(124)

Substitute this into the solution approximation determined at the end of the steps, we get

$$y_1=y_0+z\mathit{V}\left(b\right){\left(I-z{A}^{+}\right)}^{-1}\mathbf{1}{y}_0=R\left(z\right)y_{0,}$$

(125)

where

$$R\left(z\right)=1+z\mathit{V}\left(b\right){\left(I-z{A}^{+}\right)}^{-1}\mathbf{1}.$$

(126)

Continuing the iteration, we discover

$$y_{n+1}=R\left(\lambda h\right)y_n,$$

(127)

where

$$R\left(z\right)=1+z\mathit{V}\left(b\right){\left(I-z{A}^{+}\right)}^{-1}\mathbf{1}={\displaystyle \frac{det\left(I-z{A}^{+}+z\mathbf{1}\otimes \mathit{V}\left(b\right)\right)}{det\left(I-z{A}^{+}\right)}},$$

(128)

here ⊗ represents the Kronecker product. As already mentioned the function $R\left(z\right)$ is called the stability function. Thus the condition that $y_n\to 0$ as $n\to \infty$ is equivalent to $\left|R\right(z\left)\right|\le 1$. This motivates the definition of a region of absolute stability which is a set $\mathcal{S}\subset \mathbb{C}$[28]. A Runge-Kutta method is called A-stable if the region of absolute stability contains the set ${\mathbb{C}}^{-}=\{z\in \mathbb{C}:\Re \left(z\right)<0\}$. As Dahlquist stated in his work, “In most applications A-stability is not a necessary property.” [28]. This indicates that for application purposes, the word global stability has precedence [11]. The stability function $R\left(z\right)$ of a Runge-Kutta method with matrix $A^{+}$ and $\mathit{V}\left(b\right)$ is thus defined by

$$R\left(z\right)={\displaystyle \frac{det\left(I-z{A}^{+}+z\mathbf{1}\otimes \mathit{V}\left(b\right)\right)}{det\left(I-z{A}^{+}\right)}},$$

(129)

where $\mathbf{1}$ denotes the vector with all ones. According to Ehle [29], the stability functions of high-order implicit Runge-Kutta methods are Padé approximations to the exponential function $y=e^z$. Ehle’s Conjecture states that only the diagonal Padé and the first two subdiagonals are A-stable. $R\left(z\right)$ is a rational function of the form

$$R\left(z\right)={\displaystyle \frac{P_{k,l}\left(z\right)}{Q_{l,k}(-z)}}={\widehat{P}}_{[k/l]},$$

(130)

with $P_{k,l}\left(z\right)$ a kth degree polynomial $P_{k,l}\left(z\right)={\sum }_{j=0}^k\left(\begin{array}{c}k\\ j\end{array}\right){\displaystyle \frac{(k+l-j)!}{(k+l)!}}{z}^j$ and $Q_{l,k}\left(z\right)=P_{l,k}(-z)$ is a polynomial of degree l. Here ${\widehat{P}}_{[k/l]}$ denotes the Padé approximation of order $(k,l)$.

In addition, we wish to show how $R\left(z\right)$ with $z\in \mathbb{C}$ creates geometric shapes that correspond to ordered star structures. To do this, we first apply the Ehle representation of the Pade approximation to polynomials that are the same as the Kummer confluent hypergeometric function ${}_{1}F_1(a;b;z)$ in a finite representation [20]

$$P_{k,l}\left(z\right)={}_{1}F_1(-k;-k-l;z),$$

(131)

and

$$Q_{l,k}\left(z\right)={}_{1}F_1(-l;-l-k;-z).$$

(132)

The answer to Dahlquist’s initial value problem (120) is $y=exp\left(z\right)$, which is represented as a rational polynomial using Ehle’s polynomial representation. Considering Ehle’s solution, we can see that the Padé representation ${\widehat{P}}_{[k/l]}$ is dependent on the order of the polynomials k and l. This results in k zeros and l poles for the solution, depending on the order of approximation. The zeros of the divisor polynomial determine the poles. The link between the Padé approximation and the solution method’s stability function $R\left(z\right)$ allows us to investigate the interaction between them. Because the Padé approximation gives an approximate solution for the exponential function, it is beneficial to describe it as an approximation of order $[k/l]$. This can be expressed as the following relationship:

$${\widehat{P}}_{[k/l]}\left(z\right)={\displaystyle \frac{P_{k,l}\left(z\right)}{Q_{l,k}(-z)}}\simeq e^z,$$

(133)

or in a more accurate form

$${\widehat{P}}_{[k/l]}\left(z\right)-e^z=\mathcal{O}\left(z^{p+1}\right).$$

(134)

If $p\ge k+l$ then ${\widehat{P}}_{[k/l]}\left(z\right)$ is unique.

The key aspect here is that a Padé approximation or stability function is determined only by the basis $B_i\left(x\right)$ and the type of discretization used. Thus, an RK method’s stability is exclusively defined by the matrix $A^{+}$ structure and the weights $b_i$ given by (13) and (14), respectively.

Hairer et al. and Iserles et al. conducted investigations on this type of problem using order stars in the 1990s and compiled the findings in their works, which included multiple references on this subject [7,8,51].

Example 1 

(Geometry of Padé Approximation.). We use the Dahlquist solution as an example to show how this geometric analysis approach works. We use a Padé approximant ${\widehat{P}}_{[4/4]}\left(z\right)$ to compare with the solution $exp\left(z\right)$. The right-hand half determines the absolute value of the ratio $\left|{\widehat{P}}_{[4/4]}{e}^{-z}\right|$, while the contour for $\left|{\widehat{P}}_{[4/4]}{e}^{-z}\right|\le 1$ is also calculated. We color this region blue, beginning with red for value 1 and ending with 0. This coloration indicates the absolute magnitude of the ratio. The zeros of the Padé approximant are shown as dots in red (•), encircled by a color for values less than one. The stars in green (★) on $\mathbb{C}$’s right half-plane depict the Padé approximant’s poles. The hue around these poles is white, suggesting that the absolute ratio is larger than one in this area. The general distribution of zeros and poles is star-shaped see Figure 1. This is also the meaning behind the name, which means “ordered stars.”

Order stars is a simplified depiction of the context for the geometric study of the distribution of zeros and poles when approximating an entire function, such as the exponential function. The idealization is that we can establish the structure of the Padé approximants using the order of the Kummer functions, yielding an analytical representation of the solution function. The second point of simplification is that the solution function, $y=exp\left(z\right)$, is known. However, this is rarely the case when using an RK technique. Nonetheless, the situation is not hopeless, since we know that a Padé approximant exists for generic functions; for a comprehensive discussion see [52,53,54]. ▴

Example 2 

(Order Stars for gL-RK Methods.). As a fundamental illustration of the gL-RK method’s stability, we take one-parameter Gegenbauer polynomials as a basis, as well as independent discretization with Legendre and Chebyshev polynomials.

First, using the Gegenbauer $C_n^{\left(m\right)}\left(x\right)$ polynomials as basis, we have to deal with m as a influencing parameter. On the other hand, because we know that the discretization of the matrix $A^{+}$ affects the poles in the complex plane; compare with the Ehle representation (133), and that the basis affects the placement of the zeros due to relation (128), we can influence the stability by varying both driving factors. In Figure 2, we show the computations with Gegenbauer polynomials as a basis, as well as the Legendre and Chebyshev polynomials used to discretize the matrix $A^{+}$. Even when the parameter value $m>3/2$, no differences are seen as compared to the discussion in [55]. This leads us to the conclusion that the discretization of matrix $A^{+}$ is primarily responsible for stability. We reiterate that the step size is implicitly defined by the approximation order and polynomial type. The approach does not use a set step size, but rather a variable one defined by the root distribution of the polynomial type.

6. Applications

This section provides examples of the use of implicit Runge-Kutta procedures. The numerical processes have been implemented in Mathematica and are available in our book, “Applied Sinc Approximations: Polynomial and Sinc Methods for Solving Calculus Problems in Mathematica” [34]. We will not detail the implementation here, but will instead provide findings for particular scenarios. The interested reader can discover the essential functions and methods in [34].

Example 3 

(Gaussian Problem solved using Sinc-RK.). We define the initial value problem by presenting the essential information. The dependant variable y, represented here as a scalar number, satisfies the first-order linear differential equation. The parameter μ determines the coefficient for the independent variables. We investigate the initial value problem at $t_0=a=0$, where $y\left(0\right)=1$. In addition to the initial value for y, we must define the interval of interest, which runs from $a=0$ to $b=T=4$.

$$y^{\prime }=-\mu ty\left(t\right),$$

(135)

where $y\left(0\right)=y_0=1$, $\mu =1$, and $(a,b)=(0,4)$. The exact solution of the initial value problem is $y\left(t\right)=exp\left(-t^2/2\right)$. The exact solution corresponds to a Gaussian function.

The numerical solution is obtained via a Runge-Kutta-based Sinc technique using the information from the problem described above. The order of approximation s and the number of steps N in the Runge-Kutta iteration are the two most important parameters. Because we are working with Sinc approximation, we must keep in mind that the order of approximation is $\sigma =2s+1$, and the discretization of the steps takes $m_{\mathit{RK}}=2N+1$.

This information gained is depicted in Figure 3. The left panel shows the symbolic solution (solid line) and the interpolant (dashed) of discrete Sinc-RK points (red). The right panel displays the local absolute errors caused by the difference between the symbolic and the interpolant (solid line). The red dots indicate the discrete errors of the Sinc-RK approximation. The approximant error is nearly homogenous on a low level, as one would anticipate. The discrete errors of the Sinc-RK iteration occur at places on the local error curve as expected. The non-equidistant distribution of the Sinc points in the Sinc-RK iteration leads to very small errors in the initial range of the approximation. At the same time it can be seen that the error is kept almost at a homogeneous level for the end area of the iteration. This error behavior indicates that the Sinc-RK iteration is stable. We’ll go over this property in greater depth later. The global error for the setup with $s=13$ and $N=32$ is ${10}^{-5}$. This is surprising given the comparatively minimal numerical effort. At this point, the issue arises: which of the two variables, the order s or number of Sinc points $N_{\mathrm{RK}}$, affects the convergence behavior? Before delving into this subject, examine another aspect: the nonlinearity of the system of equations (11).

To investigate the behavior of the Sinc-RK approximation for nonlinear functions f on the right-hand side of the ODE, we make a minor adjustment by squaring the dependant variable. We also increase the value μ to 50. We look at the solution behavior of this ODE on the standard interval $t\in [0,1)$. The remainder of the entries stay unchanged.

$$y^{\prime }=-\mu t{y}^2,$$

(136)

where $y\left(0\right)=1$ and $\mu =50$.

To compare with the numerical solution we compute the symbolic representation of the initial value problem once more to get $y\left(t\right)=1/\left(1+25t^2\right)$. In fact, we obtained half of the Runge function from this initial value problem.

We carry out the numerical approximation as previously using the same discretization. Figure 4 displays the calculation’s results. The left panel depicts the symbolic solution (solid line) superimposed with the numerical solution (dashed). The numerical solution is an interpolation of discrete RK data (red dots). We note that typical approaches for interpolating the Runge function provide challenges, particularly with equidistant step sizes. The findings given show that an approximation (interpolation) with a modest error is attainable. The left panel of Figure 4 depicts the interpolation’s attained local errors (solid line) as well as the RK approximation based on Sinc points. As in the previous calculation, we observe an almost homogeneous local error distribution of the interpolated RK data. The magnitude of the error is as before ${10}^{-5}$ at $s=6$ and $N=32$ Sinc points. We should remark, however, that the interval has been reduced by one-quarter. However, we observed in our computations that for calculations over a four-times longer time frame, the numerical solution of the nonlinear system of equations does not converge to the appropriate extent, producing untrustworthy results. As noted earlier, this is a general problem with implicit RK methods that we exclude here. When the non-linear system of determining equations for the $k_i$ (11) can be solved, we can be confident in the approximation obtained.

We now wish to return to the issue we posed earlier: how much impact do the approximation order s and the number of Sinc points $N_{\mathrm{RK}}$ have on the accuracy of the Sinc-RK method? To have a clear understanding, we will investigate the discrete Sinc-RK data.

The discrete global errors are determined using the discrete form of the ${\mathbb{L}}_2$ norm defined as

$$E_m={\left(\sum _{k=1}^m{\left|y\left(t_k\right)-y_k\right|}^2\right)}^{1/2},k=1,2,3,\dots ,m,$$

(137)

where m is either the order $m=2s+1$ or the number of steps $m=2N+1$ utilized in the Sinc-RK computation. Here, $y_k$ represents the RK approximation at $t_k$. The discrete ${\mathbb{L}}_2$ norm compares the exact value of $y\left(t_k\right)$ to $y_k$.

To study the important effect in the Sinc-RK approximation, we must distinguish between two influencing variables: approximation order s and discretization N of the approximation. The first step is to investigate the influence on approximation order s. To do this, we maintain the discretization N constant while varying the order s from small to large values. Because we utilize a Sinc approximation to calculate the matrix $A^{+}$, which effectively defines the order, we anticipate an error law that decays after exhibiting root exponential behavior.

This sort of examination is depicted in Figure 5. In addition to altering the approximation order s, we computed the error decay for various discretization N. As anticipated, the global error satisfies the equation $E_s=\mathcal{O}\left(exp\left(-b\sqrt{s}\right)\right)$. Because the error curves for different N are so close together, the change in discretization N is minimal. The findings in Figure 5 indicate that the errors of a Sinc-RK approximation is primarily dictated by the Sinc approximation to $A^{+}$.

The results of the fitting are tabulated in Table 1. It is apparent that the decay rate is nearly the same for different discretization numbers N. This means that the Sinc-RK method follows the expected error law $E_s\sim \mathcal{O}\left(exp\left(-\beta \sqrt{s}\right)\right)$.

In a second computation, we’ll look at how the global error varies with the number of discretization points N along the RK trajectory. In order to accomplish this, we fix the order s and use (137) to calculate the global error $E_N$ because the Sinc points change along the trajectory. The results of these computations are depicted in Figure 6. We can observe that a power law $E_N\sim \mathcal{O}\left(N^{-\delta }exp\left(-\beta \sqrt{N}\right)\right)$ has a considerable impact on the error in N. The exponential dependency $\mathrm{in}\sqrt{N}$ is rather modest and may be mostly ignored (see Table 2). Under this assumption, we observe the well-known error tendency of RK techniques [8]. We can also detect from Figure 6 that s has an exponential dependency as previously discussed. This is seen by the rise in the distance between the computed values.

The results of the fitting in Figure 6 are tabulated in Table 2. It is apparent that the decay rate is almost the same for different approximation orders s. This means that the Sinc-RK method follows the error law $E_N\sim \mathcal{O}\left(N^{-\delta }exp\left(-\beta \sqrt{N}\right)\right)$ where β is a small quantity compared to δ.

In conclusion, the error laws for scalar ODEs follow the error laws of the Sinc methods in a Sinc-RK approximation, resulting in a slightly modified power law for discretization along the trajectory with Sinc points. This suggests the error level is mostly determined by the approximation’s order s.▴

Limit cycles are a feature seen in higher dimensional ODEs [31,32]. In the following example, we will compute a limit cycle for a system of first-order ODEs.

Example 4 

(Limit Cycle of a First Order System using L-RK.). Birkhoff and Rota demonstrate that coupled systems of first-order ODEs have stable limit cycles [32]. We extend their model by parameterizing the equations’ major coupling term. This parameter transforms the limit cycle from a rounded square over a circle to a more or less square shape. The limit cycle is a rounded square for $\lambda <1$, a circle for $\lambda =1$, but deformed to a diamond for $\lambda >1$. The underlying differential equations have the following form.

$$\begin{array}{ccc}\hfill u^{\prime }& =& u-v-u^3-\lambda u{v}^2,\hfill \end{array}$$

(138)

$$\begin{array}{ccc}\hfill v^{\prime }& =& u+v-v^3-\lambda u^{2}v,\hfill \end{array}$$

(139)

where λ is the coupling parameter and $u\left(0\right)=1$, $v\left(0\right)=3/2$ are the initial conditions. To assess the influence of λ, we adjust the coupling parameter at a predetermined interval. There is no exact solutions to this system of ODEs.

We use the L-RK approach to calculate the limit cycle structure. To do this, we employ the same starting conditions for varying parameters λ to ensure a comparable outcome. The parameters λ are supplied in discrete steps, and the approximants are calculated with the Lagrange Runge-Kutta approach. We set the approximation order s to $s=64$ and the number of discrete points on the trajectory to $n=176$. The greater order assures a minimal error, and the quantity of discrete points ensures computation stability. We employ orthogonal Chebyshev polynomials to determine the discretization locations.

The computation results are displayed in Figure 7 and Figure 8. Figure 7 depicts the approximants (dots) and Lagrange interpolants (solid lines) for various coupling constants λ. The approximants and interpolants match precisely, resulting in only negligible errors owing to interpolation. We employ Chebyshev polynomials as a basis since similar to Sinc discretization, convergence is better near the interval’s start and end. This is due to the root distribution being relatively near to the interval’s endpoints. We also explored Gegenbauer and Legendre polynomials, which produced tiny oscillations around the ends. Given the identical starting conditions and an increase in the coupling constant λ, the approximants exhibit oscillations after a delay time τ. When λ is raised, a transition zone is generated and the amplitude of the oscillations decrease without damping phenomena. Furthermore, we notice the creation of approximant structures as λ grows. The oscillations’ operating time, τ, rises with λ.

A quasi phase diagram $v\left(t\right)$ over $u\left(t\right)$ (Figure 8 left panel) clearly represents the reduction in amplitudes, the distortion of the oscillation, and the delay τ for the commencement of the oscillations. Given the minor errors present, it is evident that limit cycles are emerging. The limit cycles include circular structures for $\lambda =1$, rounded square structures for $\lambda <1$, and clear square structures for $\lambda >1$. We can also notice that the limit cycles for $\lambda >1$ exhibit tiny oscillations in the transient range τ. However, these do not cause instability and are governed by the order s of the approximants and the relatively large number n of discretization points. These oscillations will fade away if we keep increasing s and n. The right panel in Figure 8 depicts a limit cycle created by approximants with varying initial conditions. The coupling constant here is $\lambda =3$.

Figure 8 (right panel) shows that the center of the phase plane is unstable, whereas the limit cycle is stable.▴

Example 5 

(Troesch’s Problem as Initial Value Problem Solved using gL-RK.). First, let us define the problem so that we can assess the relevance of the strategy. Weibel first developed Troesch’s equation while investigating a plasma confinement problem in 1959 [30]. The problem was formulated as a non-linear two point boundary problem (BVP), which Troesch studied one year later [33] in the representation:

$$u″\left(x\right)=\lambda sinh\left(\lambda u\right(x\left)\right),$$

(140)

with $u\left(0\right)=0$ and $u\left(1\right)=1$, where $x\in [0,1]$. The inconspicuous problem comprises one parameter λ, which might change in principle in the range $\lambda \in (0,\infty )$. However, the numerical approximation for values $\lambda >5$ presented significant complications [33,37]. This was the starting point for several investigations into the problem, ending in the 1970s with an analytic solution given by Bulrisch, Stoer, and Bulrisch as a toy problem [35,38]. Since then, it is known that the equivalent initial value problem may be addressed using Jacobi elliptic functions. The exact solution to the related Cauchy initial value problem (IVP)

$$u″\left(x\right)=\lambda sinh\left(\lambda u\right(x\left)\right),$$

(141)

with $u\left(0\right)=0$ and $u^{\prime }\left(0\right)=s$, where s is a real value and satisfies $s>0$. The IVP possesses the exact solution

$$u\left(x\right)={\displaystyle \frac{2}{\lambda }}\mathit{Arcsinh}\left({\displaystyle \frac{s}{2}}{\displaystyle \frac{\mathit{sn}\left(\lambda x\left|k^2\right.\right)}{\mathit{cn}\left(\lambda x\left|k^2\right.\right)}}\right),\mathit{with}{k}^2=1-{\displaystyle \frac{s^2}{4}},$$

(142)

which is important for testing the numerical approximation in an accurate manner. Stoer and Bulirsch highlighted the relationship between (140) and (141) while discussing shooting approaches for solving BVP utilizing IVP methods [39]. These writers have already pointed out that the BVP has a moveable logarithmic singularity around its upper border. However, they demonstrated that, although being outside the interval $[0,1]$, this singularity is caused by cn’s critical behavior at the terminus. In (142), sn and cn are Jacobian elliptic functions of modulus k that depend on the initial slope s. If $K\left(k^2\right)$ represents the quarter period of cn, then cn has a zero at

$$x_s={\displaystyle \frac{K\left(k^2\right)}{\lambda }},$$

(143)

and therefore $u\left(x\right)$ has a logarithmic singularity there [39]. But $K\left(k^2\right)$ has the expansion

$$K\left(k^2\right)\simeq ln\left({\displaystyle \frac{4}{\sqrt{1-k^2}}}\right)+{\displaystyle \frac{1}{4}}\left(ln\left({\displaystyle \frac{4}{\sqrt{1-k^2}}}\right)-1\right)\left(1-k^2\right)+\dots ,$$

(144)

or, rewritten in terms of s,

$$K\left(k^2\right)\simeq ln\left({\displaystyle \frac{8}{\left|s\right|}}\right)+{\displaystyle \frac{s^2}{16}}\left(ln\left({\displaystyle \frac{8}{\left|s\right|}}\right)-1\right)+\dots .$$

(145)

For the solution of the actual boundary-value problem, i.e.,

$$u(1;s)=1,$$

(146)

we require the values of s that depend on the choice of λ. For different λ values, the starting value s is determined by solving (146) with the exact solution in this expression. To solve the equation for s, a high level of precision in the root finding problem is necessary to provide a valid result. All of these relationships are known since Bulrisch and Stoer investigated the topic [35,38].

However, as Troesch [33] and Ehrlich [36] noted, this two-point BVP is numerically unstable and intricate. At approximately the same time, Roberts and Shipman [40] proposed that this problem could only be solved by combining three methods: the perturbation technique, the parallel shooting method, and the continuation method. They suggest that a mix of these strategies is required since none of them alone is effective enough to provide a solution. In fact, their technique only provides valid solutions for equation (140) for $1<\lambda <5$. Jones [41] utilized a modified Newton approach to solve this problem. Miele, Aggarwal, and Tietze [42] attempted to apply a modified quasi-linearization approach with limited success. Troesch discovered its numerical solution in [37], utilizing the shooting method. In [46], this problem was solved using the decomposition technique; in [48,49], the variational iteration method was used; in [40], a combination of the multipoint shooting method with the continuation and perturbation technique was applied to the problem; in [42], the quasi linearization method had a limited success; in [43], the method of transformation groups appeared; in [44], the invariant embedding method delivered hardly convincing results; in [45], the inverse shooting was applied and many more. Overall, the methodologies utilized were fairly restrictive and did not result in a definitive numerical procedure. As a result, the problem requires more numerical efforts to be addressed effectively.

We shall apply gL-RK techniques to Troesch’s IVP, demonstrating that a suitable choice of discretization points and the employment of orthogonal polynomials may produce satisfactory numerical results. Using these features correctly can successfully solve Troesch’s problem for a wide variety of λ values. We should point out that Sinc approaches have already been used to solve this problem, albeit with inefficient discretization [47].

To solve the problem (141) numerically as an IVP, we begin with the precise solution (142), which includes the unknown parameter s. If we wish to solve the problem using an IVP, we must first determine the initial slope s. To find s, we use the boundary condition (146) $u(1;s)=1$. This boundary condition is transformed into an initial condition by calculating s numerically from the precise solution using the root finding approach for a given λ. We utilize the parameter $\lambda =7$, which is greater than the number provided in the literature.

This, in turn, yields the exact answer for the specified λ value. We will then use the precise answer as a reference for estimating errors.

To solve the initial value problem numerically, we employ the Jacobi $P_m^{(\alpha ,\beta )}\left(x\right)$ polynomials as a basis for collocation and discretization. Our observations indicate that given the order s, there is a lower constraint at which the Newton method cannot solve the implicit determining equations. Only when this lower threshold is exceeded can a stable solution be achieved. In the following computations, we use $s=12$ as the minimum approximation order. The gL-RK approach is discretized using Jacobi points of size $n=64$. We use the synchronized variant, with parameters $\alpha =1/3$ and $\beta =1/4$.

Figure 9 depicts the results of the computations with $\lambda =7$. The left panel shows the exact solution (solid line), the interpolant (dashed), and the approximant (dots). The numerical approximations are visually indistinguishable from the exact result. The right panel displays the absolute values of local errors for both the interpolant (solid line) and the approximant (dots). The approximant’s discrete error is several orders of magnitude smaller than the interpolant’s local error. However, at times, the interpolant and approximant errors coincide (which is not rectified in this graph).

We solved the approximations for Troesch’s initial value problem using a similar calculating procedure as above, but with the parameter λ variable. The different solutions are illustrated in Figure 10. The related local errors have the same order of magnitude as in Figure 9. The graphs in Figure 10 indicate that the function y rises extremely quickly as t approaches 1 for larger λ. The shift in slope grows steeper as λ increases. The gL-RK technique efficiently handles the fast expansion of the function y, even with small values for approximation order $s=12$ and RK number of support points $n=64$. We employed Jacobi polynomials for both interpolation and discretization.

The example demonstrates how the gL-RK approach may be used to efficiently compute a reasonably stiff problem such as the Troesch equation.▴

Example 6 

(The nonlinear Klein-Gordon Equation Solved using gL-RK.). In this example, we will provide a solution technique for nonlinear partial differential equations (nPDE) that is based on Lie’s method of infinitesimal transformations [50] and use the gL-RK method for numerical solution. This method enables us to convert the PDE to an ODE and construct the resulting initial value problem. We shall demonstrate this with the cubic Klein-Gordon (cKG) equation.

The cubic Klein-Gordon equation is well known and can be written as

$$u_{\mathit{tt}}-u_{\mathit{xx}}=\alpha u+\beta u^3,$$

(147)

where α and β are real parameters. The determination of infinitesimal symmetries (Lie symmetries) reveals the following transformation under which equation (147) stays invariant

$$\begin{array}{ccc}\hfill {\xi }_1(x,t,u)& =& -k_1-2k_{3}t,\hfill \\ \hfill {\xi }_2(x,t,u)& =& -k_2-2k_{3}t,\hfill \\ \hfill {\varphi }_1(x,t,u)& =& 0.\hfill \end{array}$$

(148)

Here the parameters $k_i$ with $i=1,2,3$ are group constants which can be freely chosen. The choice of these parameters selects a transformation group under which (147) will be invariant. In addition this selection allows to reduce the nPDE to a nonlinear ODE which solves equation (147).

Selecting the sub-group $k_1=c$, $k_2=1$, and $k_3=0$; i.e. translations in spatial and temporal variables, we are able to generate via a Lie reduction the following ODE

$$u″={\displaystyle \frac{c^2}{c^2-1}}\left(\alpha u\left(\zeta \right)+\beta u{\left(\zeta \right)}^3\right),$$

(149)

where $u\left(\zeta \right)=u(t-x/c)$. It is straight forward to solve this ODE as an initial value problem with $u\left(0\right)=\sqrt{2}$, $u^{\prime }\left(0\right)=0$ by selecting $\mu =\sqrt{2}$ and $\alpha =1$ and $\beta =-1$; here ${\mu }^2=c^2/\left(c^2-1\right)$. The exact solution using the given initial values is $u\left(\zeta \right)=\sqrt{2}\sqrt{1-tanh{\left(\sqrt{2}\zeta \right)}^2}$.

Since our target is the solution of the boundary value problem (147) with $u(x\to \pm \infty ,t)=0$ and for $t=0$, we will use $u(x,0)=\sqrt{2}\mathit{sech}\left(\sqrt{2}x\right)$. This allows us to approximate the natural boundary conditions at $\pm \infty$ by values at $x=-5$ and $t=0$ in the following setup of the initial value problem in the reduced variable $\zeta =x-ct$.

Note this definition is possible because we are looking for moving wave solutions of the problem. The application of the gL-RK method to the “initial” value problem delivers the approximant. We should remark, that the initial point of integration ${\zeta }_0$ is different from zero.

Figure 11 shows the gL-RK calculation results for u and $u^{\prime }$. The left panel displays the exact solution (solid line), the interpolant (dashed), and the approximant (dots). Figure 11’s right panel displays the absolute local error $E_n\left(\zeta \right)$ for u as well as the absolute local error of the approximant (dots). With a small number of steps in the RK technique and an approximation order of $s=24$, we reach a global error of ${10}^{-6}$. For practical applications, this is enough to derive the solution structure.

When ζ is substituted by the dependency in x and t; $\zeta =x-ct$, the similarity solution yields the solution to the boundary value problem. Figure 12 illustrates the wave’s form as a simple solitary peak, as well as its derivative. We see that the basic peak travels at a constant speed $c=\sqrt{2}$.

When we vary the cubic Klein-Gordon equation’s parameters α and β and adjust the initial conditions, we get an oscillating wave in space and time. The used initial conditions are $u\left(0\right)=5.4$, $u^{\prime }\left(0\right)=0$ on the interval $\zeta \in (0,24)$. The parameters are set to $v=3$, $\alpha =-1$, and $\beta =-1$.

Again we use the gL-RK approach to recalculate the approximant and interpolant for this initial value problem using $s=24$ and $n=128$ RK steps. Jacobian polynomials are used to define the basis and discretize the data.

Figure 13 shows the results. The solution’s structure is periodic in both u and $u^{\prime }$. This is clearly a nonlinear solution, as the minima and maxima of u are rather sharply curved, and the derivatives are flatter at their extrema. The right panel of Figure 13 displays a phase-space picture that clearly depicts the solution’s periodicity.

If we return from the similarity variable $\zeta =x-vt$ to the spatial and temporal coordinates (see Figure 14), the wave formation can be represented in an $(x,t)$ coordinate system. The tilting of the wave crests in this coordinate system is attributable to the finite velocity v. However, it is obvious that the wave nature occurs in the spatiotemporal continuum, but only a limited range can be represented by the finite integration with the gL-RK method.

In conclusion, the gL-RK approach, when paired with the Lie symmetry method, produces useful results for partial differential equations. This necessitates reframing a boundary-value problem as an initial value problem. We understand that this is not always practicable, but when combined, it may result in significant time and resource savings.▴

7. Conclusions

The technique of implicit Runge-Kutta procedures shows that, using the collocation method and therefore the representation of the matrix $A^{+}$ and the weights $b_i$, a well-defined algorithm exists for determining a numerical implicit Runge-Kutta process. The many representations of the matrix $A^{+}$ generate variability, allowing for tailored adaptation to a given problem. This technique works for both stiff and non-stiff initial value problems. In addition to numerical integration, comprehending the matrix $A^{+}$ enables direct comments about the method’s stability when the pole and zero distributions of the stability function are analyzed.

The method’s stability may be graphically evaluated using proven methods like order stars. Our numerical investigations demonstrated that parameter-independent orthogonal polynomials are not restricted in stability. This is an open problem for parameter-dependent orthogonal polynomials when the number of parameters is greater than one. For example, for Gegenbauer polynomials $G_n^{\left(m\right)}$, there is an upper constraint in m that is smaller than $3/2$. We discovered both a variable upper bound and a fixed lower bound at $m_c=1/2$. Additionally, there are no restrictions for orders $s<5$. The upper bound was estimated with a rational expression $r_{[6/6]}\left(s\right)$ and a Thiele continued fraction. The approximant approaches the fixed value $r_{[6/6]}(\infty )=3/2$ for $s\to \infty$. The interval for stable L-RK techniques using Gegenbauer polynomials is provided by $-1/2<m<r_{[6/6]}\left(s\right)$ for $s\ge 5$ and $m>-1/2$ for $s<5$.

The approach presented is simple to implement and highly efficient in numerical calculation. Furthermore, the strategy based on Sinc approximations has the well-known feature of exponential convergence. The zeros of orthogonal polynomials are sensitive in the classes of resultant functions, allowing for a priori estimate of the global error. It is also critical that a priori error formulae exist for calculating the global error, allowing automated computation by defining error tolerances. An a posteriori additional error analysis is obsolete.

The examples provided demonstrate that implicit RK techniques based on indefinite integrals provide an approach that is highly adaptable while requiring little implementation work. The main idea is that just the matrix $A^{+}$ and the weights $b_i$ must be known. Thus, the procedure is specified using a Butcher table. The variability derives from the discrete representation of indefinite integrals, which may be obtained using orthogonal polynomials or Sinc techniques.