# Adaptive Piecewise Poly-Sinc Methods for Ordinary Differential Equations

Omar Khalil  Mathematics Department, German University in Cairo, New Cairo
City, Cairo 11835, Egypt  Hany El-Sharkawy  Mathematics Department, German
University in Cairo, New Cairo City, Cairo 11835, Egypt  Department of
Mathematics, Faculty of Science, Ain Shams University, 11566 Abbassia, Cairo,
Egypt  Maha Youssef  Institute of Applied Analysis and Numerical Simulation,
University of Stuttgart, Pfaffenwaldring 57, D-70569 Stuttgart, Germany  Gerd
Baumann  Mathematics Department, German University in Cairo, New Cairo City,
Cairo 11835, Egypt  Faculty of Natural Sciences, University of Ulm, Albert-
Einstein-Allee 11, D-89069 Ulm, Germany

###### Abstract

We propose a new method of adaptive piecewise approximation based on Sinc
points for ordinary differential equations. The adaptive method is a piecewise
collocation method which utilizes Poly-Sinc interpolation to reach a preset
level of accuracy for the approximation. Our work extends the adaptive
piecewise Poly-Sinc method to function approximation, for which we derived an
a priori error estimate for our adaptive method and showed its exponential
convergence in the number of iterations. In this work, we show the exponential
convergence in the number of iterations of the a priori error estimate
obtained from the piecewise collocation method, provided that a good estimate
of the exact solution of the ordinary differential equation at the Sinc points
exists. We use a statistical approach for partition refinement. The adaptive
greedy piecewise Poly-Sinc algorithm is validated on regular and stiff
ordinary differential equations.

Keywords— adaptive approximation; Poly-Sinc interpolation; Sinc methods;
Lagrange interpolation; initial value problems; boundary value problems;
exponential convergence; regular differential equations; stiff differential
equations

##  1 Introduction

Numerous phenomena in engineering, physics, and mathematics are modeled either
by initial value problems (IVPs) or by boundary value problems (BVPs)
described by ordinary differential equations (ODEs). Accordingly, the
numerical solution of IVPs for deterministic and random ODEs is a basic
problem in the sciences. For a review of the state of the art on theory and
algorithms for numerical initial value solvers, we refer to the monographs [1,
2, 3, 4, 5, 6, 7] and the references therein.

Exact solutions may not be available for some ODEs. This has led to the
development of a number of methods to estimate the a posteriori error, which
is based on the residual of the ODE [8], forming the basis for adaptive
methods for ODEs. The a posteriori error estimates have been derived for
different numerical methods, such as piecewise polynomial collocation methods
[9, 10] and Galerkin methods [11, 12, 13, 14]. An a posteriori error estimate
in connection with adjoint methods was developed in [15]. Kehlet et al. [16]
incorporated numerical round-off errors in their a posteriori estimates. An a
posteriori error estimate based on the variational principle was derived in
[17]. Convergence rates for the adaptive approximation of ODEs using a
posteriori error estimation were discussed in [18, 19]. A less common form is
the a priori error estimate [20, 8]. Hybrid a priori–a posteriori error
estimates for ODEs were developed in [21, 22]. An advantage of the a priori
error estimate over the a posteriori error estimate is that the a priori error
estimate does not require the computation of the residual of the ODE. However,
some knowledge about the exact solution of the ODE is required for the a
priori error estimate. It was shown in [23, 24] that the a priori error
estimate of the Poly-Sinc approximation is exponentially convergent in the
number of Sinc points, provided that the exact solution belongs to the set of
analytic functions.

We propose an adaptive piecewise method, in which the points in a given
partition are used as partitioning points. This piecewise property allows for
a greater flexibility of constructing the polynomials of arbitrary degree in
each partition. Recently, we developed an a priori error estimate for the
adaptive method based on piecewise Poly-Sinc interpolation for function
approximation [25]. In this work [25], we used a statistical approach for
partition refinement in which we computed the fraction of a standard deviation
[26, 27, 28] as the ratio of the mean absolute deviation to the sample
standard deviation. It was shown in [29] that the ratio approaches
$\sqrt{\frac{2}{\pi}}\approx 0.798$ for an infinite number of normal samples.
We extend the work [25] for regular and stiff ODEs. In this paper, we discuss
the adaptive piecewise Poly-Sinc method for regular and stiff ODEs, and show
that the exponentially convergent a priori error estimate for our adaptive
method differs from that for function approximation [25] by a small constant.

This paper is organized as follows. Section 2 provides an overview of the
Poly-Sinc approximation, the residual computation, the indefinite integral
approximation, and the collocation method. Section 3 discusses the piecewise
collocation method, which is the cornerstone of the adaptive piecewise Poly-
Sinc algorithm. In Section 4, we present the adaptive piecewise Poly-Sinc
algorithm for ODEs and the statistical approach for partition refinement. We
also demonstrate the exponential convergence of the a priori error estimate
for our adaptive method. We validate our adaptive Poly-Sinc method on regular
ODEs and ODEs whose exact solutions exhibit an interior layer, a boundary
layer, and a shock layer in Section 5. Finally, we present our concluding
remarks in Section 6.

##  2 Background

###  2.1 Poly-Sinc Approximation

A novel family of polynomial approximation called Poly-Sinc interpolation
which interpolate data of the form $\\{x_{k},y_{k}\\}_{k=-M}^{N}$ where
$\\{x_{k}\\}_{k=-M}^{N}$ are Sinc points, were derived in [23, 30] and
extended in [24]. The interpolation to this type of data is accurate provided
that the function $y$ with values $y_{k}=y(x_{k})$ belong to the space of
analytic functions [30, 31]. For the ease of presentation and discussion, we
assume that $M=N$. Poly-Sinc approximation was developed in order to mitigate
the poor accuracy associated with differentiating the Sinc approximation when
approximating the derivative of functions [23]. Moreover, Poly-Sinc
approximation is characterized by its ease of implementation. Theoretical
frameworks on the error analysis of function approximation, quadrature, and
the stability of the Poly-Sinc approximation were studied in [23, 24, 32, 33].
Furthermore, Poly-Sinc approximation was used to solve BVPs in ordinary and
partial differential equations [34, 35, 36, 31, 37, 38]. We start with a brief
overview of Lagrange interpolation. Then, we discuss the generation of Sinc
points using conformal mappings.

####  2.1.1 Lagrange Interpolation

Lagrange interpolation is a polynomial interpolation scheme [39], which is
constructed by Lagrange basis polynomials

$u_{k}(x)=\dfrac{g(x)}{(x-x_{k})g^{\prime}(x_{k})},\quad k=1,2,\ldots,m,$

where $\\{x_{k}\\}_{k=1}^{m}$ are the interpolation points and
$g(x)=\prod_{l=1}^{m}(x-x_{l})$. The Lagrange basis polynomials satisfy the
property

$u_{k}(x_{j})=\begin{cases*}1{,}&if $k=j$,\\\ 0{,}&if $k\neq j$.\end{cases*}$

Hence, the polynomial approximation in the Lagrange form can be written as

$y_{h}(x)=\sum_{k=1}^{m}y(x_{k})u_{k}(x),$ (1)

where $y_{h}(x)$ is a polynomial of degree $m-1$ and it interpolates the
function $f(x)$ at the interpolation points, i.e., $y_{h}(x_{k})=y(x_{k})$.
For Sinc points, the polynomial approximation $y_{h}(x)$ becomes

$y_{h}(x)=\sum_{k=-N}^{N}y(x_{k})u_{k}(x),$ (2)

where $m=2N+1$ is the number of Sinc points. If the coefficients $y(x_{k})$
are unknown, then we replace $y(x_{k})$ with $c_{k}$, and Equations (1) and
(2) become

$y_{\rm{c}}(x)=\sum_{k=1}^{m}c_{k}\,u_{k}(x)$ (3)

and

$y_{\rm{c}}(x)=\sum_{k=-N}^{N}c_{k}\,u_{k}(x),$ (4)

respectively.

####  2.1.2 Conformal Mappings and Function Space

We introduce some notations related to Sinc methods [23, 30, 24]. Let
$\varphi\ :\ D\to D_{d}$ be a conformal map that maps a simply connected
region $D\subset\mathbb{C}$ onto the strip

$D_{d}=\\{z\in\mathbb{C}:|\imaginary(z)|<d\\},$

where $d$ is a given positive number. The region $D$ has a boundary $\partial
D$, and let $a$ and $b$ be two distinct points on $\partial D$. Let
$\psi=\varphi^{-1},\psi\ :\ D_{d}\to D$ be the inverse conformal map. Let
$\Gamma$ be an arc defined by

$\Gamma=\\{z\in[a,b]\ :\ z=\psi(x),\ x\in\mathbb{R}\\},$

where $a=\psi(-\infty)$ and $b=\psi(\infty)$. For real finite numbers $a,\,b$,
and $\Gamma\subseteq{\mathbb{R}}$, $\varphi(x)={\ln}((x-a)/(b-x))$ and
$x_{k}=\psi(kh)=(a+be^{kh})/(1+e^{kh})$ are the Sinc points with spacing
$h(d,\beta_{\mathsf{s}})=\left(\dfrac{\pi
d}{\beta_{\mathsf{s}}N}\right)^{1/2},\ \beta_{\mathsf{s}}>0$ [40, 30]. Sinc
points can be also generated for semi-infinite or infinite intervals. For a
comprehensive list of conformal maps, see [24, 30].

We briefly discuss the function space for $y$. Let $\rho=e^{\varphi}$,
$\alpha_{\mathsf{s}}$ be an arbitrary positive integer number, and
$\mathbb{L}_{\alpha_{\mathsf{s}},\beta_{\mathsf{s}}}(D)$ be the family of all
functions that are analytic in $D=\varphi^{-1}(D_{d})$ such that for all $z\in
D$, we have

$|y(z)|\leq
C\dfrac{|\rho(z)|^{\alpha_{\mathsf{s}}}}{[1+|\rho(z)|]^{\alpha_{\mathsf{s}}+\beta_{\mathsf{s}}}}.$

We next set the restrictions on $\alpha_{\mathsf{s}},\beta_{\mathsf{s}}$, and
$d$ such that $0<\alpha_{\mathsf{s}}\leq 1$, $0<\beta_{\mathsf{s}}\leq 1$, and
$0<d<\pi$. Let $\mathbb{M}_{\alpha_{\mathsf{s}},\beta_{\mathsf{s}}}(D)$ be the
set of all functions $g$ defined on $D$ that have finite limits
$g(a)=\lim_{z\to a}g(z)$ and $g(b)=\lim_{z\to b}g(z)$, where the limits are
taken from within $D$, and such that
$y\in\mathbb{L}_{\alpha_{\mathsf{s}},\beta_{\mathsf{s}}}(D)$, where

$y=g-\dfrac{g(a)+\rho\,g(b)}{1+\rho}.$

The transformation guarantees that $y$ vanishes at the endpoints of $(a,b)$.
We assume that $y$ is analytic and uniformly bounded by $\mathcal{B}(y)$,
i.e., $|y(x)|\leq\mathcal{B}(y)$, in the larger region

$D_{2}=D\cup_{t\in(a,b)}B(t,r),$

where $r>0$ and $B(t,r)=\\{z\in\mathbb{C}:|z-t|<r\\}$.

###  2.2 Residual

The residual is used as a measure of the accuracy of the adaptive Poly-Sinc
method. The general form of a second-order ODE can be expressed as [41]

$F(x,y,y^{\prime},y^{\prime\prime})=0.$ (5)

An exact solution $y$ satisfies (5). If the exact solution $y$ is unknown, we
replace it with the approximation $y_{\rm{c}}$ and Equation (5) becomes

$F(x,y_{\rm{c}},y_{\rm{c}}^{\prime},y_{\rm{c}}^{\prime\prime})=R(x),$ (6)

where $R(x)$ is the residual. The residual in (6) for the $i-$th iteration
becomes

$F\left(x,y_{\rm{c}}^{(i)},\left(y_{\rm{c}}^{(i)}\right)^{\prime},\left(y_{\rm{c}}^{(i)}\right)^{\prime\prime}\right)=R^{(i)}(x),\quad
i=1,2,\ldots,\kappa,$

where $\kappa$ is the number of iterations.

We will denote the residual for integral and differential equations with
$R_{\mathrm{I}}$ and $R_{\mathrm{D}}$, respectively. The residual is used as
an indicator for partition refinement as discussed in Algorithm 4 (see Section
4).

###  2.3 Error Analysis

We briefly discuss the error analysis for Poly-Sinc approximation over the
global interval $[a,b]$. At the end of this section, we will discuss the error
analysis of Poly-Sinc approximation for IVPs and BVPs.

For the Poly-Sinc approximation on a finite interval [23, 24, 32, 33, 38], it
was shown that

${\max_{x\in[a,b]}}\
|y(x)-y_{h}(x)|\leq\dfrac{A}{r^{m}}\sqrt{N}e^{-\beta\sqrt{N}},$

where $y(x)$ is the exact solution and $y_{h}(x)$ is its Poly-Sinc
approximation, $A$ is a constant independent of $N$, $m=2N+1$ is the number of
Sinc points in the interval, $r$ is the radius of the ball containing the $m$
Sinc points, and $\beta>0$ is the convergence rate parameter. On a finite
interval $[a,b]$, it was shown that [24, 38]

$\max_{x\in[a,b]}\ |y(x)-y_{h}(x)|\leq
A\left|\dfrac{b-a}{2r}\right|^{m}N^{3/2}\tanh[-4](\dfrac{\eta}{4\sqrt{N}})\exp(-\dfrac{\pi^{2}\sqrt{N}}{2\eta}){,}$
(7)

where $\eta$ is a positive constant. Inequality (7) can be written as

$\max_{x\in[a,b]}\ |y(x)-y_{h}(x)|\leq
A\left|\dfrac{b-a}{2r}\right|^{m}N^{\alpha}\exp(-\gamma N^{\beta}).$

Next, we discuss the collocation method for IVPs and BVPs.

###  2.4 Collocation Method

A collocation method [42, 43] is a technique in which a system of algebraic
equations is constructed from the ODE via the use of collocation points. Here,
we adopt the Poly-Sinc collocation method [44, 36], in which the collocation
points are the Sinc points and the basis functions are the Lagrange
polynomials with Sinc points.

####  2.4.1 Initial Value Problem

The IVP is transformed into an integral equation. We briefly discuss the
approximation of indefinite integrals using Poly-Sinc methods ([45] § 9.3).
Define

$(\mathcal{J}^{+}w\,y)(x)=\int_{a}^{x}y(t)w(t)\differential{t},$ (8)

where the weight function $w(x)$ is positive on the interval $(a,b)$ and has
the property that the moments $\int_{a}^{b}x^{j}w(x)\differential{x}$ do not
vanish for $j=0,1,2,\ldots$. Let $A^{+}$ be an $m\times m$ matrix whose
entries are

$[A^{+}]_{k\,j}=\int_{a}^{x_{k}}u_{j}(x)w(x)\differential{x},$

where ${u_{j}(x)},j=-N,\ldots,N$, are the Lagrange basis polynomials stacked
in a vector $L(x)=({u_{-N}(x)},\ldots,{u_{N}(x)})^{\top},$ and
$(\cdot)^{\top}$ is the transpose operator. The interpolation points
$\\{x_{j}\\}_{j=-N}^{N}$ are the Sinc points generated in the interval $[a,b]$
as discussed in Section 2.1.2.

Then, the indefinite integral (8) can be approximated as

$\displaystyle(\mathcal{J}_{m}^{+}w\,y)(x)$
$\displaystyle=\sum_{j=-N}^{N}({\mathcal{J}_{m}^{+}}w\,y)(x_{j})u_{j}(x)$
$\displaystyle\approx\sum_{j=-N}^{N}\left[\sum_{k=-N}^{N}y(x_{k})[A^{+}]_{j\,k}\right]u_{j}(x)$
$\displaystyle=\sum_{j=-N}^{N}u_{j}(x)\sum_{k=-N}^{N}y(x_{k})[A^{+}]_{j\,k}$
$\displaystyle={L(x)^{\top}A^{+}Vy,}$

where $Vy=(y(x_{-N}),\ldots,y(x_{N}))^{\top}$. We state the following theorem
for IVPs [30].

######  Theorem 1 (Initial Value Problem ([30] §1.5.8)).

If $y\in\mathbb{M}_{\alpha_{\mathsf{s}},\beta_{\mathsf{s}}}(D)$, then, for all
$N>1$

$\lVert\mathcal{J}^{+}w\,y-\mathcal{J}_{m}^{+}w\,y\rVert={O}(\epsilon_{N}),$

where $\epsilon_{N}=\sqrt{N}e^{-\beta\sqrt{N}}$.

####  2.4.2 Boundary Value Problem

For a BVP, the collocation method solves for the unknown coefficients $c_{k}$
in (3) or (4) by setting

$R_{\mathrm{D}}(x_{k})=0,\quad k=-N,\ldots,N.$

However, we replace the two equations corresponding to $x_{-N}$ and $x_{N}$
with the boundary conditions $y(a)=y_{a}$ and $y(b)=y_{b}$, respectively. We
state the following theorem for BVPs.

######  Theorem 2 (Boundary Value Problem ([30] §1.5.6)).

If $y\in\mathbb{M}_{\alpha_{\mathsf{s}},\beta_{\mathsf{s}}}(D)$ and
$\bm{c}=(c_{-N},\ldots,c_{N})^{\top}$ is a complex vector of order $m$, such
that for some $\delta>0$,

$\left(\sum_{j=-N}^{N}\left|y(x_{j})-c_{j}\right|^{2}\right)^{1/2}<\delta,$

then,

$\lVert y-L(x)^{\top}\bm{c}\rVert<C\epsilon_{N}+\delta.$

##  3 Piecewise Collocation Method

We discuss the piecewise collocation method, in which the domain $I=[a,b]$ is
discretized into $K\in\mathbb{N}$ non-overlapping partitions
$I_{n}=[x^{n-1},x^{n}),n=1,2,\ldots,K-1$, $x^{0}=a,\ x^{K}=b$ with
$I_{K}=[x^{K-1},b]$ and $\cup_{n=1}^{K}I_{n}=I=[a,b]$. The space of piecewise
discontinuous polynomials can be defined as

$\mathfrak{D}_{k}(I)=\\{v:v|_{I_{n}}\in\mathbb{P}_{k}(I_{n}),\ \
n=1,{2,}\ldots,K\\},$

where $\mathbb{P}_{k}(I_{n})$ denotes the space of polynomials of degree at
most $k$ on $I_{n}$. The piecewise collocation method solves the collocation
method in Section 2.4 over partitions. The approximate solution in the global
partition $[a,b]$ can be written as

$y_{h}(x)=\sum_{k=1}^{K}y_{h,\,k}(x)\mathbbm{1}_{x\in
I_{k}}=\sum_{k=1}^{K}\mathbbm{1}_{x\in
I_{k}}\sum_{j=1}^{m_{k}}y_{j,\,k}u_{j,\,k}(x),$ (9)

where $y_{h,\,k}(x)=\sum_{j=1}^{m_{k}}y_{j,\,k}u_{j,\,k}(x)$ is the Lagrange
interpolation in the $k-$th partition. The basis functions

$u_{j,\,k}(x)=\dfrac{g_{k}(x)}{(x-x_{j,\,k})g_{k}^{\prime}(x_{j,\,k})},\quad
j=1,2,\ldots,m_{k},\ k=1,2\ldots,K,$

where $\\{x_{j,\,k}\\}_{j=1}^{m_{k}}$ are the interpolation points in the
$k-$th partition, $g_{k}(x)=\prod_{l=1}^{m_{k}}(x-x_{l,\,k})$, and $m_{k}$ is
the number of points in the $k-$th partition. The function
$\mathbbm{1}_{\mathcal{C}}$ is an indicator function which outputs 1 if the
condition $\mathcal{C}$ is satisfied and otherwise 0. If the coefficients
$y_{j,\,k}$ are unknown, then we replace $y_{j,\,k}$ with $c_{j,\,k}$, and
Equation (9) becomes

$y_{\rm{c}}(x)=\sum_{k=1}^{K}y_{{\rm{c}},\,k}(x)\mathbbm{1}_{x\in
I_{k}}=\sum_{k=1}^{K}\mathbbm{1}_{x\in
I_{k}}\sum_{j=1}^{m_{k}}c_{j,\,k}u_{j,\,k}(x).$ (10)

The residual for the $k-$th partition can be written as

$F\left(x,y_{{\rm{c}},\,k},\left(y_{{\rm{c}},\,k}\right)^{\prime},\left(y_{{\rm{c}},\,k}\right)^{\prime\prime}\right)=R_{k}(x),\quad
x\in I_{k},\ k=1,{2,}\ldots,K.$

The collocation method solves for the unknowns $c_{j,\,k}$ by setting
$R_{k}(x_{j,\,k})=0,$ $\,j=1,{2,}\ldots,m_{k},k=1,{2,}\ldots,K$, which we
discuss next for IVPs and BVPs.

###  3.1 Initial Value Problem

In this section, we provide examples for first-order and second-order IVPs.

#### Relaxation Problem

We discuss the piecewise collocation method for a first-order IVP in integral
form. Consider the following relaxation or decay Equation [46] on the interval
$[a,b]$

$\dfrac{\differential{y(x)}}{\differential{x}}=-\alpha y(x),\quad y(a)=y_{a},$
(11)

where $\alpha>0$ is the relaxation parameter. The exact solution is
$y(x)=y_{a}\exp(-\alpha(x-a))$. We transform the IVP (11) into an integral
form

$y(x)=y_{a}-\alpha\int_{a}^{x}y(t)\differential{t}.$

The residual becomes

$R_{\mathrm{I}}(x)=y_{\rm{c}}(x)-y_{a}+\alpha\int_{a}^{x}y(t)\differential{t}{.}$

We approximate the indefinite integral as discussed in Section 2.4.1 and the
approximate residual becomes

$\tilde{R}_{\mathrm{I}}(x,a,y_{a},y_{\rm{c}}(x))=y_{\rm{c}}(x)-y_{a}+\alpha(\mathcal{J}_{m}^{+}y_{\rm{c}})(x).$

The domain $[a,b]$ is partitioned as discussed in Section 3. For the $k-$th
partition, $k=1,{2,}\ldots,K$, we replace $a$ with $x^{k-1}$ and $y_{a}$ with
$y_{{\rm{c}},\,k-1}(x^{k-1})$. The approximate residual becomes

$\tilde{R}_{\mathrm{I},\,k}(x,x^{k-1},{y_{{\rm{c}},\,k-1}(x^{k-1})},{y_{{\rm{c}},\,k}(x))}={y_{{\rm{c}},\,k}(x)}-{y_{{\rm{c}},\,k-1}(x^{k-1})}+\alpha(\mathcal{J}_{m}^{+}y_{\rm{c}})(x),\qquad
x\in I_{k}.$

We remove the equation corresponding to the leftmost Sinc point in each
partition, and replace them with the conditions

$\displaystyle y_{\rm{c}}(x^{0})$ $\displaystyle=y_{a},$ (12a) $\displaystyle
y_{{\rm{c}},\,k}(x^{k-1})$ $\displaystyle=y_{{\rm{c}},\,k-1}(x^{k-1}),\quad
k=2,{3,}\ldots,K,$ (12b)

and the set of equations

$\tilde{R}_{\mathrm{I},\,k}(x_{j,\,k},x^{k-1},y_{{\rm{c}},\,k-1}(x^{k-1}),y_{{\rm{c}},\,k}(x_{j,\,k}))=0,\
j=2,{3,}\ldots,m_{k},k=1,{2,}\ldots,K.$ (13)

The set of equations (12b) is known as the continuity equations at the
interior boundaries [47]. The collocation algorithm for the IVP (11) is
outlined in Algorithm 1.

1

input : $K:$ number of partitions

$m_{k}:$ number of Sinc points in the $k-$th partition

output : $y_{\rm{c}}(x)$: approximate solution

Replace $y(x)$ with the global approximate solution (10).

Solve for the $m_{k}\,K$ unknowns $\\{c_{j,k}\\}_{j=1,\,k=1}^{m_{k},\,K}$
using the initial condition (12a), continuity Equation (12b), and the set of
equations for the residual (13).

Algorithm 1 Piecewise Poly-Sinc Algorithm (IVP (11)).

###  3.2 Hanging Bar Problem

We discuss the piecewise collocation method for a second-order IVP in integral
form. Considering the following IVP on the interval $[a,b]$

$\displaystyle\forall
x\in\left[a,b\right],\quad-(\tilde{K}(x)y^{\prime}(x))^{\prime}$
$\displaystyle=f(x),$ (14) $\displaystyle y(a)$ $\displaystyle=y_{a},$
$\displaystyle y^{\prime}(a)$ $\displaystyle=\tilde{y}_{a}.$

where $y(x)$ is the sought-for solution. In the context of the hanging bar
problem [48], $y(x)$ and $\tilde{K}(x)$ are the displacement and the material
property of the bar at the position $x$, respectively. For simplicity, we set
$\tilde{K}(x)=1$. Equation (14) can be written as a system of first-order
equations

$\displaystyle y^{\prime}=q,\quad q^{\prime}=-f,\quad y(a)=y_{a},\
q(a)=\tilde{y}_{a}.$ (15)

The integral form of (15) is

$\displaystyle y(x)$ $\displaystyle=y_{a}+\int_{a}^{x}q(t)\differential{t},$
(16) $\displaystyle q(x)$
$\displaystyle=\tilde{y}_{a}-\int_{a}^{x}f(t)\differential{t}.$ (17)

Plugging (17) into (16), we obtain

$\begin{split}y(x)&=y_{a}+\int_{a}^{x}\left[\tilde{y}_{a}-\int_{a}^{{t}}f(s)\differential{s}\right]\differential{t}\\\
&=y_{a}+\tilde{y}_{a}(x-a)-\int_{a}^{x}\int_{a}^{t}f(s)\differential{s}\differential{t}.\end{split}$

Using integration by parts $\int u\differential{v}=uv-\int v\differential{u}$
[49] with $u(t)=\int_{a}^{t}f(s)\differential{s}$,

$\displaystyle\int_{a}^{x}\int_{a}^{t}f(s)\differential{s}\differential{t}$
$\displaystyle=t\int_{a}^{t}f(s)\differential{s}\bigg{|}_{t=a}^{t=x}-\int_{a}^{x}tf(t)\differential{t}$
$\displaystyle=x\int_{a}^{x}f(s)\differential{s}-\int_{a}^{x}tf(t)\differential{t}$
$\displaystyle=x\int_{a}^{x}f(s)\differential{s}-\int_{a}^{x}sf(s)\differential{s}$
$\displaystyle=\int_{a}^{x}(x-s)f(s)\differential{s},$

where we set $t=s$. Thus, the integral form of the solution to (14) becomes

$y(x)=y_{a}+(x-a)\tilde{y}_{a}-x\int_{a}^{x}f(s)\differential{s}+\int_{a}^{x}sf(s)\differential{s}.$

The residual can be written as

$R_{\mathrm{I}}(x,a,y_{a},\tilde{y}_{a},y_{\rm{c}}(x))=y_{\rm{c}}(x)-y_{a}-(x-a)\tilde{y}_{a}+x\int_{a}^{x}f(s)\differential{s}-\int_{a}^{x}sf(s)\differential{s}.$
(18)

Approximating the indefinite integral in (18), the approximate residual
becomes

$\tilde{R}_{\mathrm{I}}(x,a,y_{a},\tilde{y}_{a},y_{\rm{c}}(x))=y_{\rm{c}}(x)-y_{a}-(x-a)\tilde{y}_{a}+x(\mathcal{J}_{m}^{+}f)(x)-(\mathcal{J}_{m}^{+}x\,f)(x).$
(19)

For the $k-$th partition, $k=1,2,\ldots,K$, we replace $a$ with $x^{k-1}$ ,
$y_{a}$ with $y_{{\rm{c}},\,k-1}(x^{k-1})$, and $\tilde{y}_{a}$ with
$y_{{\rm{c}},\,k-1}^{\prime}(x^{k-1})$. The approximate residual becomes

$\tilde{R}_{\mathrm{I},\,k}(x,x^{k-1},{y_{{\rm{c}},\,k-1}(x^{k-1})},{y_{{\rm{c}},\,k-1}^{\prime}(x^{k-1})},{y_{{\rm{c}},\,k}(x)})={y_{{\rm{c}},\,k}(x)}-{y_{{\rm{c}},\,k-1}(x^{k-1})}\\\
-(x-{x^{k-1}}){y_{{\rm{c}},\,k-1}^{\prime}(x^{k-1})}+x(\mathcal{J}_{m}^{+}f)(x)-(\mathcal{J}_{m}^{+}x\,f)(x),\
x\in I_{k}.$

We remove the equations corresponding to the leftmost and rightmost Sinc
points in each partition, and replace them with the conditions

$\displaystyle{y_{\rm{c}}(x^{0})}$ $\displaystyle=y_{a},$ (20a)
$\displaystyle{y_{\rm{c}}^{\prime}(x^{0})}$ $\displaystyle=\tilde{y}_{a},$
(20b) $\displaystyle y_{{\rm{c}},\,k}(x^{k-1})$
$\displaystyle=y_{{\rm{c}},\,k-1}(x^{k-1}),\quad k=2,\ldots,K,$ (20c)
$\displaystyle y_{{\rm{c}},\,k}^{\prime}(x^{k-1})$
$\displaystyle=y_{{\rm{c}},\,k-1}^{\prime}(x^{k-1}),\quad k=2,\ldots,K,$ (20d)

and the set of equations

$\tilde{R}_{\mathrm{I},\,k}(x_{j,\,k},x^{k-1},y_{{\rm{c}},\,k-1}(x^{k-1}),y_{{\rm{c}},\,k-1}^{\prime}(x^{k-1}),y_{{\rm{c}},\,k}(x_{j,\,k}))=0,$
(21)

$j=2,\ldots,m_{k}-1,k=1,\ldots,K$. Equations (20c)–(20d) are known as the
continuity equations at the interior boundaries [47]. The collocation
algorithm for the IVP (14) is outlined in Algorithm 2.

1

input : $K:$ number of partitions

$m_{k}:$ number of Sinc points in the $k-$th partition

output : $y_{\rm{c}}(x)$: approximate solution

Replace $y(x)$ with the global approximate solution (10).

Solve for the $m_{k}\,K$ unknowns $\\{c_{j,k}\\}_{j=1,\,k=1}^{m,\,K}$ using
initial conditions (20a)–(20b), continuity Equations (20c)–(20d), and the set
of equations for the residual (21).

Algorithm 2 Piecewise Poly-Sinc Algorithm (IVP (14)).

###  3.3 Boundary Value Problem

The collocation method for the BVP is similar to that of the IVP in Section
3.2, except that we replace the set of equations (20) with

$\displaystyle y_{\rm{c}}(x^{0})$ $\displaystyle=y_{a}{,}$ (22a)
$\displaystyle y_{\rm{c}}(x^{K})$ $\displaystyle=y_{b}{,}$ (22b)
$\displaystyle y_{{\rm{c}},\,k}(x^{k-1})$
$\displaystyle=y_{{\rm{c}},\,k-1}(x^{k-1}),\quad k=2,{3,}\ldots,K{,}$ (22c)
$\displaystyle y_{{\rm{c}},\,k}^{\prime}(x^{k-1})$
$\displaystyle=y_{{\rm{c}},\,k-1}^{\prime}(x^{k-1}),\quad k=2,{3,}\ldots,K,$
(22d)

and the set of equations of the residual for the BVP becomes

$R_{\mathrm{D},\,k}(x_{j,\,k})=0,\ j=2,{3,}\ldots,m_{k}-1,k=1,{2,}\ldots,K,$
(23)

where $R_{\mathrm{D},\,k}$ is the residual of the differential equation in the
$k-$th partition. The piecewise Poly-Sinc collocation algorithm for the BVP is
outlined in Algorithm 3.

input : $K:$ number of partitions

$m_{k}:$ number of Sinc points in the $k-$th partition

output : $y_{\rm{c}}(x)$: approximate solution

Replace $y(x)$ with the global approximate solution (10).

Solve for the $m_{k}\,K$ unknowns $\\{c_{j,k}\\}_{j=1,\,k=1}^{m,\,K}$ using
boundary conditions (22a)–(22b), continuity Equations (22c)–(22d), and the set
of equations for the residual (23).

Algorithm 3 Piecewise Poly-Sinc Algorithm (BVP).

##  4 Adaptive Piecewise Poly-Sinc Algorithm

This section introduces the greedy algorithmic approach used in adaptive
piecewise Poly-Sinc methods. The core feature used is the non-overlapping
properties of Sinc points and the uniform exponential convergence on each
partition of the approximation interval. Greedy algorithms seek the “best”
candidate of possible solutions at a given step [50]. Greedy algorithms have
been applied to model order reduction for parametrized partial differential
equations [51, 52]. The adaptive piecewise Poly-Sinc algorithm is greedy in
the sense that it makes a choice that aims to find the “best” approximation
for the solution of the ODE in the current step [50]. The algorithm takes an
iterative form in which it computes the $L^{2}$ norm values of the residual
for all partitions constituting the global interval $I=[a,b]$. At the $i-$th
step, the algorithm refines the partitions for which the $L^{2}$ norm values
of the residual are relatively large. By refining the partitions as discussed
above, it is expected that the mean value of the $L^{2}$ norm values over all
partitions decreases in each step. As the iteration proceeds, the algorithm
expects to find the “best” polynomial approximation for the solution of the
ODE.

###  4.1 Algorithm Description

We discuss the adaptive algorithm for the piecewise Poly-Sinc approximation.
The following steps of the adaptive algorithm are performed in an iterative
loop [53]

$\mathsf{SOLVE}\to\mathsf{ESTIMATE}\to\mathsf{MARK}\to\mathsf{REFINE}.$

The adaptive piecewise Poly-Sinc algorithm is outlined in Algorithm 4. The
refinement strategy is performed as follows. For the $i-$th iteration, we
compute the set of $L^{2}$ norm values
$\left\\{\left\|R_{{k}}^{(i)}(x)\right\|_{L^{2}(I_{k}^{(i)})}\right\\}_{k=1}^{K_{i}}$
over the $K_{i}$ partitions, from which the sample mean
$\overline{R_{i}}=\frac{1}{K_{i}}\sum_{j=1}^{K_{i}}\left\|R_{j}^{(i)}(x)\right\|_{L^{2}(I_{j}^{(i)})}$
and the sample standard deviation [54]

$s_{i}=\frac{1}{\sqrt{K_{i}-1}}\sqrt{\sum_{j=1}^{K_{i}}\left(\left\|R_{j}^{(i)}(x)\right\|_{L^{2}(I_{j}^{(i)})}-\overline{R_{i}}\right)^{2}}$

are computed [26, 27, 28]. The residual $R_{j}^{(i)}(x)$ for the $j-$th
partition and the $i-$th iteration is discussed in Section 2.2. The partitions
with the indices
$\mathfrak{I}_{i}=\left\\{j:\left\|R_{j}^{(i)}(x)\right\|_{L^{2}(I_{j}^{(i)})}-\overline{R_{i}}\geq\omega_{i}\,s_{i}\right\\}$
are marked for refinement, where the statistic [29]

$\omega_{i}=\dfrac{\frac{1}{K_{i}}\sum_{j=1}^{K_{i}}\left|\left\|R_{j}^{(i)}(x)\right\|_{L^{2}(I_{j}^{(i)})}-\overline{R_{i}}\right|}{s_{i}}.$

Using Hölder’s inequality for sums with $p=q=2$ ([55] § 3.2.8), one can show
that $\omega_{i}\leq\sqrt{\frac{K_{i}-1}{K_{i}}}<1$. We restrict to second-
order moments only. The points in the partitions with the indices
$\mathfrak{I}_{i}$ are used as partitioning points and $m=2N+1$ Sinc points
are inserted in the newly created partitions. The algorithm terminates when
the stopping criterion is satisfied. The approximate solution
$y_{\rm{c}}^{(i)}(x),i=1,\ldots,\kappa$, for the $i-$th iteration is computed
using the collocation method outlined in Algorithms 1 and 2 for IVPs and
Algorithm 3 for BVPs. We note that, for partition refinement, the residual is
computed in its differential form $R_{\rm{D}}(x)$.

The definite integral in the $L^{2}$ norm [56] is numerically computed using a
Sinc quadrature [30], i.e.,

$\lVert
f(x)\rVert_{L^{2}([a,b])}^{2}=\int_{a}^{b}|f(x)|^{2}\differential{x}\approx
h\sum_{k=-N}^{N}\dfrac{1}{\varphi^{\prime}(x_{k})}f^{2}(x_{k}),$

where $\\{x_{k}\\}_{k=-N}^{N}\in[a,b]$ are the quadrature points, which are
also Sinc points, and $\varphi(x)$ is the conformal mapping in Section 2.1.2.
The supremum norm on an interval $I=[a,b]$ is approximated as ([57] Table 2.1)

$\lVert f(x)\rVert_{I}\approx\max\ \\{|f(x_{k})|\\}_{k=-N}^{N},$

where $\\{x_{k}\\}_{k=-N}^{N}$ are the Sinc points on $I$, whose generation is
discussed in Section 2.1.2.

###  4.2 Error Analysis

We state below the main theorem.

######  Theorem 3 (Estimate of Upper Bound [25]).

Let $y$ be in $\mathbb{M}_{\alpha_{\mathsf{s}},\beta_{\mathsf{s}}}(\varphi)$,
analytic and bounded in $D_{2}$, and let $y_{h}^{(i)}(x)$ be the piecewise
Poly-Sinc approximation in the $i$-th iteration. Let
$\xi_{i}=\operatorname*{arg\,max}_{k}\left|\mathsf{Len}(I_{k}^{(i)})\right|$
be the index of the largest partition in the $i-$th iteration and
$\mathsf{Len}(I_{k}^{(i)})$ be the length of the $k-$th partition in the
$i-$th iteration. Let $K_{i}$ be the number of partitions in the $i-$th
iteration. Then, there exists a constant $A$, independent of the $i-$th
iteration, such that

$\max_{x\in[a,b]}\left|y(x)-y_{h}^{(i)}(x)\right|\leq K_{i}\,E_{i},$

where
$E_{i}=\dfrac{A}{(2r_{\xi_{i}})^{m_{\xi_{i}}}}\,\lambda^{m_{\xi_{i}}(i-1)}(b-a)^{m_{\xi_{i}}}$.

For fitting purposes, we compute the mean value of the error estimate, i.e.,

$\max_{x\in[a,b]}\dfrac{1}{K_{i}}\left|y(x)-y_{h}^{(i)}(x)\right|\leq E_{i}.$

We state the following theorem on collocation.

input : Global partition $[a,b]$, threshold $\varepsilon_{\mathrm{stop}}$, $N$

output : $\mathsf{S}$: set of points

$\kappa$: number of iterations

$\mathsf{R}=\\{\overline{R_{i}}\\}_{i=1}^{\kappa}$: set of mean values of
$\lVert
R_{k}^{(i)}(x)\rVert_{L^{2}(I_{k}^{(i)})},\,k=1,\ldots,K_{i},\,i=1,\ldots,\kappa$

$y_{\rm{c}}^{(\kappa)}(x)$: approximate solution

1

_2_ init:  _ $\mathsf{S}=\\{\\}$, $\mathsf{P}=\\{x^{0}=a,x^{1}=b\\}$,
$\mathsf{R}=\\{\\}$_

3 __ Start with $2N+1$ points in the global interval $[a,b]$. Append the
points to $\mathsf{S}$ and $\mathsf{P}.$

(Solve). Compute $y_{\rm{c}}^{(i)}(x)$ for IVP (Sections 3.1 and 3.2) or BVP
(Section 3.3).

4 (Estimate). Compute $\\{\|R^{(1)}(x)\|\\}_{k=1}^{K_{1}=1}$ over the global
interval $[a,b]$. Append the mean value $\overline{R_{1}}$ in $\mathsf{R}$.

5 (Mark). Set $\mathcal{I}_{1}=1$.

6 (Refine).  Use the points in the partition with index $\mathcal{I}_{1}=1$ as
partitioning points. Insert $2N+1$ Sinc points in each of the newly created
partitions. Update $\mathsf{S}$ and $\mathsf{P}$.

7 Set $i=2$.

8

9

10while  _$\overline{R_{i-1}} >\varepsilon_{\mathrm{stop}}$_ do

11 (Solve). Compute $y_{\rm{c}}^{(i)}(x)$ for IVP (Sections 3.1 and 3.2) or
BVP (Section 3.3).

12 (Estimate). Compute
$\left\\{\left\|R_{k}^{(i)}(x)\right\|_{L^{2}(I_{k}^{(i)})}\right\\}_{k=1}^{K_{i}}$
over $K_{i}$ partitions. Compute the sample mean value $\overline{R_{i}}$ and
the sample standard deviation $s_{i}$. Append $\overline{R_{i}}$ in
$\mathsf{R}$.

13 (Mark). Identify the partitions with indices
$\mathfrak{I}_{i}=\left\\{j:\left\|R_{j}^{(i)}(x)\right\|_{L^{2}(I_{j}^{(i)})}-\overline{R_{i}}\geq\omega_{i}\,s_{i}\right\\},\
j=1,\ldots,K_{i}$.

14 (Refine). Use the $2N+1$ points in the partitions with indices
$\mathfrak{J}_{i}$ as partitioning points. Insert $2N+1$ points in each of the
newly created partitions. Update $\mathsf{S}$ and $\mathsf{P}.$

15 $i\leftarrow i+1$.

16

17 end while

18$\kappa\leftarrow i-1$.

19

Algorithm 4 Adaptive Piecewise Poly-Sinc Algorithm.

######  Theorem 4.

Let $y$ be in
${\mathbb{M}_{\alpha_{\mathsf{s}},\beta_{\mathsf{s}}}(\varphi)}$, analytic and
bounded in $D_{2}$, and let $y_{\rm{c}}^{(i)}(x)$ be the piecewise Poly-Sinc
approximation in the $i$-th iteration with the estimated coefficients
$c_{j,\,k}^{(i)}$ using the piecewise collocation method in Section 3. Let
$m_{k}^{(i)}$ be the number of points in the $k-$th partition for the $i-$th
iteration and $n_{i}=\sum_{k=1}^{K_{i}}m_{k}^{(i)}$ be the total number of
points in the $i-$th iteration. Let
$\bm{c}_{i}=(c_{1,1}^{(i)},\ldots,c_{1,m_{1}^{(i)}}^{(i)},\ldots,c_{1,K_{i}}^{(i)},\ldots,c_{1,m_{K_{i}}^{(i)}}^{(i)})$
be the $n_{i}\times 1$ vector of estimated coefficients in
$y_{\rm{c}}^{(i)}(x)$ and $\bm{y}_{i}=y(x_{j,\,k}^{(i)})$ be the corresponding
vector of the exact values of $y(x)$ at the Sinc points
$\\{x_{j,\,k}^{(i)}\\}$ with

$\|\bm{c}_{i}-\bm{y}_{i}\|_{1}<\delta,\qquad i=1,2,\ldots,\kappa.$

Then,

$\max_{x\in[a,b]}|y(x)-y_{\rm{c}}^{(i)}(x)|\leq
K_{i}\,E_{i}+\delta\left(\frac{1}{\pi}\ln(\overline{\overline{m}})+1.07618\right),$

where $\overline{\overline{m}}=\max_{k,\,i}\\{m_{k}^{(i)}\\}$ and
$\|\cdot\|_{1}$ denotes the $\ell_{1}$ norm ([58] Ch. 5).

###### Proof.

The derivation follows that of ([30] § 1.5.6).

$\displaystyle\max_{x\in[a,b]}|y(x)-y_{\rm{c}}^{(i)}(x)|$
$\displaystyle=\|y(x)-y_{\rm{c}}^{(i)}(x)\|_{I}$
$\displaystyle=\|y(x)-y_{\rm{c}}^{(i)}(x)+y_{h}^{(i)}(x)-y_{h}^{(i)}(x)\|_{I}$
$\displaystyle\leq\|y(x)-y_{h}^{(i)}(x)\|_{I}+\|y_{\rm{c}}^{(i)}(x)-y_{h}^{(i)}(x)\|_{I}$
$\displaystyle\begin{split}&\leq\bigg{\|}y(x)-{y_{h}^{(i)}(x)}\bigg{\|}_{I}+\\\
&\quad\bigg{\|}\sum_{k=1}^{K_{i}}\mathbbm{1}_{x\in{I_{k}^{(i)}}}\sum_{j=1}^{m_{k}^{(i)}}(c_{j,\,k}^{(i)}-y_{j,\,k}^{(i)})u_{j,\,k}^{(i)}(x)\bigg{\|}_{I}\\\
\end{split}$ $\displaystyle\leq
K_{i}\,E_{i}+\sum_{k=1}^{K_{i}}\bigg{\|}\sum_{j=1}^{m_{k}^{(i)}}(c_{j,\,k}^{(i)}-y_{j,\,k}^{(i)})u_{j,\,k}^{(i)}(x)\bigg{\|}_{I_{k}^{(i)}}$
$\displaystyle\leq
K_{i}\,E_{i}+\sum_{k=1}^{K_{i}}\bigg{\|}\sum_{j=1}^{m_{k}^{(i)}}|c_{j,\,k}^{(i)}-y_{j,\,k}^{(i)}|\cdot|u_{j,\,k}^{(i)}(x)|\bigg{\|}_{I_{k}^{(i)}}$
$\displaystyle\leq
K_{i}\,E_{i}+\frac{\delta}{K_{i}}\sum_{k=1}^{K_{i}}\bigg{\|}\sum_{j=1}^{m_{k}^{(i)}}|u_{j,\,k}^{(i)}(x)|\bigg{\|}_{I_{k}^{(i)}}$
$\displaystyle\leq
K_{i}\,E_{i}+\frac{\delta}{K_{i}}\sum_{k=1}^{K_{i}}\left(\frac{1}{\pi}\ln(m_{k}^{(i)})+1.07618\right)$
$\displaystyle\leq
K_{i}\,E_{i}+\frac{\delta}{K_{i}}K_{i}\left(\frac{1}{\pi}\ln(\overline{\overline{m}})+1.07618\right)$
$\displaystyle=K_{i}\,E_{i}+\delta\left(\frac{1}{\pi}\ln(\overline{\overline{m}})+1.07618\right),$

where
$\big{\|}\sum_{j=1}^{m_{k}^{(i)}}|u_{j,\,k}^{(i)}(x)|\big{\|}_{I_{k}^{(i)}}\approx\frac{1}{\pi}\ln(m_{k}^{(i)})+1.07618$
is the Lebesgue constant for Poly-Sinc approximation [33, 59, 31]. On average,
the term
$|c_{j,\,k}^{(i)}-y_{j,\,k}^{(i)}|<\frac{\delta}{n_{i}}<\frac{\delta}{K_{i}\
\min_{k}\\{m_{k}^{(i)}\\}}<\frac{\delta}{K_{i}}$. ∎

For fitting purposes, we compute the mean value of the error estimate, i.e.,

$\displaystyle{\max_{x\in[a,b]}}\
\dfrac{1}{K_{i}}\left|y(x)-y_{\rm{c}}^{(i)}(x)\right|$ $\displaystyle\leq
E_{i}+\frac{\delta}{K_{i}}\left(\frac{1}{\pi}\ln(\overline{\overline{m}})+1.07618\right)$
(24)
$\displaystyle<E_{i}+\delta\left(\frac{1}{\pi}\ln(\overline{\overline{m}})+1.07618\right).$

##  5 Results

The results in this section were computed using Mathematica [60]. We tested
our adaptive algorithm on regular and stiff ODEs. The Sinc spacing is
$h=\frac{\pi}{\sqrt{e\,N}}$. For all examples, we set $e=1/2$ and the number
of points per partition to be constant, i.e., $m_{k}^{(i)}=m=2N+1$. A
precision of $200$ digits is used.

###  5.1 Norms

The supremum norm has a theoretical advantage. However, its computation is
slower than that of the $L^{2}$ norm [25]. Hence, we use the $L^{2}$ norm in
our computations.

###  5.2 Initial Value Problem

We test our adaptive piecewise Poly-Sinc algorithm on regular first-order and
second-order IVPs.

######  Example 1 (Relaxation Problem).

We start with the relaxation problem in Section 3.1. We set $a=0$ and the
exact solution becomes $y(x)=\exp(-\alpha\,x)$. We set the exponential decay
parameter $\alpha=20$ and confine the domain of the solution to the interval
$[0,1]$. The approximate solution $y_{\rm{c}}(x)$ is computed as discussed in
Section 3.1.

We set the number of Sinc points to be inserted in all partitions as
$m=2N+1=5$. The stopping criterion $\varepsilon_{\mathrm{stop}}=10^{-6}$ was
used. The algorithm terminates after $\kappa=7$ iterations and the number of
points $|{\mathsf{S}}|=530$.

Figure 1a shows the approximate solution $y_{\rm{c}}^{(7)}(x)$. A proper
subset of the set of points $\mathsf{S}$ is shown as red dots, which are
projected onto the approximate solution $y_{\rm{c}}^{(7)}(x)$. This proper
subset is used to observe the approximate solution $y_{\rm{c}}^{(7)}(x)$. We
plot the statistic $\omega_{i}$ as a function of the iteration index
$i,\,i=2,{3,}\ldots,7$, in Figure 1b. The oscillations are decaying and the
statistic $\omega_{i}$ is converging to an asymptotic value. The mean value
$\overline{\omega_{i}}\approx 0.6$ is denoted by a horizontal line.

We perform the least-squares fitting of the logarithm of the set $\mathsf{R}$
to the logarithm of the upper bound (24). Figure 2a shows the least-squares
fitted model (24) to the set $\mathsf{R}$. The dots represent the set
$\mathsf{R}$ and the solid line represents the least-squares fitted model
(24). Figure 2b shows the residual, absolute local approximation error, and
the mean value for the last iteration. The mean value $\overline{R_{7}}$ is
below the threshold value $10^{-6}$. The $L^{2}$ norm of the approximation
error $\|y(x)-y_{\rm{c}}^{(7)}(x)\|\approx 1.5\times 10^{-7}$.

######  Example 2 (Hanging Bar Problem).

We apply the collocation method on the hanging bar problem (14) and
$y(x)=e^{x}(x-1)^{2}$ [61]. The approximate solution $y_{\rm{c}}(x)$ is
computed as discussed in Section 3.1.

We set the number of Sinc points to be inserted in all partitions as
$m=2N+1=7$. The stopping criterion $\varepsilon_{\mathrm{stop}}=10^{-6}$ was
used. The algorithm terminates after $\kappa=3$ iterations and the number of
points $|{\mathsf{S}}|=350$.

Figure 3 shows the approximate solution $y_{\rm{c}}^{(3)}(x)$. A proper subset
of the set of points $\mathsf{S}$ is shown as red dots, which are projected
onto the approximate solution $y_{\rm{c}}^{(3)}(x)$.

We perform the least-squares fitting of the logarithm of the set $\mathsf{R}$
to the logarithm of the upper bound (24). Figure 4a shows the least-squares
fitted model (24) to the set $\mathsf{R}$. Figure 4b shows the residual,
absolute local approximation error, and the mean value for the last iteration.
The mean value $\overline{R_{3}}$ is below the threshold value $10^{-6}$. The
$L^{2}$ norm of the approximation error $\|y(x)-y_{\rm{c}}^{(3)}(x)\|\approx
5.82\times 10^{-9}$.

a

b

Figure 1: (a) The approximating polynomial $y_{\rm{c}}^{(7)}(x)$. A proper
subset of the set of points $\mathsf{S}$ is shown. (b) Plot of the statistic
$\omega_{i}$. ( ) $\overline{\omega_{i}}\approx 0.6$.

a

b

Figure 2:  (a) Fitting the upper bound (24) with the set $\mathsf{R}$. (b)
Visualization of the residual, absolute local approximation error, and the
mean value.  Figure 3: The approximating polynomial $y_{\rm{c}}^{(3)}(x)$. A
proper subset of the set of points $\mathsf{S}$ is shown.

a

b

Figure 4: (a) Fitting the upper bound (24) with the set $\mathsf{R}$. (b)
Visualization of the residual, absolute local approximation error and the mean
value.

###  5.3 Boundary Value Problem

We discuss a number of stiff BVPs [8] based on the general linear second-order
BVP

$-(a(x)\,y^{\prime})^{\prime}+b(x)y^{\prime}+c(x)\,y=f(x),\quad x\in[0,1],$
(25)

where $a(x)>0,\,b(x),\,c(x)$ are the coefficients, and $f(x)$ is the source
term.

######  Example 3.

We study the BVP (25) with $a(x)=x+0.01,\,b(x)=c(x)=0$, and $f(x)=1$. The
exact solution is

$y(x)=\mathrm{c}_{1}{\ln}(1+100x)+\mathrm{c}_{2}-x,$

where $\mathrm{c}_{1}=1/{\ln}(101)$ and $\mathrm{c}_{2}=0$ are obtained from
the boundary conditions $y(0)=y(1)=0$. The exact solution experiences a
boundary layer near $x=0$.

We set the number of Sinc points to be inserted in all partitions as
$m=2N+1=5$. The stopping criterion $\varepsilon_{\mathrm{stop}}=10^{-6}$ was
used. The algorithm terminates after $\kappa=10$ iterations and the number of
points $|{\mathsf{S}}|=2055$.

Figure 5a shows the approximate solution $y_{\rm{c}}^{(10)}(x)$. A proper
subset of the set of points $\mathsf{S}$ is shown as red dots, which are
projected onto the approximate solution $y_{\rm{c}}^{(10)}(x)$. We plot the
statistic $\omega_{i}$ as a function of the iteration index
$i,\,i=2,{3,}\ldots,10$, in Figure 5b. It is observed that the oscillations
are decaying and the statistic $\omega_{i}$ is converging to an asymptotic
value. The mean value $\overline{\omega_{i}}\approx 0.64$.

a

b

Figure 5: (a) The approximating polynomial $y_{\rm{c}}^{(10)}(x)$. A proper
subset of the set of points $\mathsf{S}$ is shown. (b) Plot of the statistic
$\omega_{i}$. ( ) $\overline{\omega_{i}}\approx 0.64$.

We perform the least-squares fitting of the logarithm of the set $\mathsf{R}$
to the logarithm of the upper bound (24). Figure 6a shows the least-squares
fitted model (24) to the set $\mathsf{R}$. Figure 6b shows the residual,
absolute local approximation error, and the mean value for the last iteration.
The plot of the $L^{2}$ norm values of the residual over the partitions
demonstrates that fine partitions are formed near $x=0$ due to the presence of
the boundary layer. The mean value $\overline{R_{10}}$ is below the threshold
value $10^{-6}$. The $L^{2}$ norm of the approximation error
$\|y(x)-y_{\rm{c}}^{(10)}(x)\|\approx 1.12\times 10^{-8}$. The threshold value
in [8] is $0.05$.

a

b

Figure 6: (a) Fitting the upper bound (24) with the set $\mathsf{R}$. (b)
Visualization of the residual, absolute local approximation error, and the
mean value.

######  Example 4.

We study the BVP (25) with $a(x)=0.01,\,b(x)=0,\,c(x)=1$, and $f(x)=1/x$.
Using the variation of parameters method [62], the exact solution of this
problem is

$y(x)=-5\,\mathrm{Ei}(-10x)\exp(10x)+5\,\mathrm{Ei}(10x)\exp(-10x)+\mathrm{c}_{1}\exp(10x)+\mathrm{c}_{2}\exp(-10x),$

where $\mathrm{Ei}(x)=\int_{-\infty}^{x}\frac{\exp(t)}{t}\differential{t}$ is
the exponential integral function [55] and

$\mathrm{c}_{1}=-\mathrm{c}_{2}=\dfrac{-5\,\mathrm{Ei}(-10)\exp(10)+5\,\mathrm{Ei}(10)\exp(-10)}{\exp(-10)-\exp(10)}.$

The exact solution $y(x)$ experiences a boundary layer near $x=0$ and a slope
change at approximately $x=0.5$. Equation (25) is multiplied by the factor $x$
so that the residual $R(x)$ does not contain a singularity at $x=0$.

We set the number of Sinc points to be inserted in all partitions as
$m=2N+1=5$. The stopping criterion $\varepsilon_{\mathrm{stop}}=10^{-6}$ was
used. The algorithm terminates after $\kappa=9$ iterations and the number of
points $|{\mathsf{S}}|=1630$.

We perform the least-squares fitting of the logarithm of the set $\mathsf{R}$
to the logarithm of the upper bound (24). Figure 7a shows the least-squares
fitted model (24) to the set $\mathsf{R}$. Figure 7b shows the residual,
absolute local approximation error, and the mean value for the last iteration.
The plot of the $L^{2}$ norm values of the residual over the partitions shows
that fine partitions are formed near $x=0$ due to the presence of the boundary
layer. The mean value $\overline{R_{9}}$ is below the threshold value
$10^{-6}$. The $L^{2}$ norm of the approximation error
$\|y(x)-y_{\rm{c}}^{(9)}(x)\|\approx 1.6\times 10^{-6}$. The threshold value
in [8] is $0.01$.

The approximating polynomial $y_{\rm{c}}^{(9)}(x)$ and a proper subset of the
set of points $\mathsf{S}$ are shown in Figure 8a. The corresponding plot for
the statistic $\omega_{i}$ a is shown in Figure 8c. The oscillations are
decaying and the mean value $\overline{\omega_{i}}\approx 0.66$. It was
mentioned that Equation (25) was multiplied by the factor $x$ so that the
residual $R(x)$ does not contain a singularity at $x=0$. We replace the
residual $R(x)$ with the quantity $y(x)-y_{\rm{c}}(x)$ in Algorithm 4 and the
BVP (25) contains the term $1/x$. We set the number of Sinc points to be
inserted in all partitions as $m=2N+1=5$. The stopping criterion
$\varepsilon_{\mathrm{stop}}=10^{-6}$ was used. The algorithm terminates after
$\kappa=8$ iterations and the number of points $|{\mathsf{S}}|=730$, which is
smaller than the one obtained from multiplying the residual by $x$. This is
expected since the exact solution $y(x)$ is used. The approximating polynomial
$y_{\rm{c}}^{(8)}(x)$ and a proper subset of the set of points $\mathsf{S}$
are shown in Figure 8b. The corresponding plot for the statistic $\omega_{i}$
is shown in Figure 8d. The mean value $\overline{\omega_{i}}\approx 0.61$. It
is observed that the statistic $\omega_{i}$ oscillates around the mean
$\overline{\omega_{i}}\approx 0.61$.

a

b

Figure 7: (a) Fitting the upper bound (24) with the set $\mathsf{R}$. (b)
Visualization of the residual, absolute local approximation error, and the
mean value.

a

b

c

d

Figure 8: Visualization of the approximating polynomial and statistic
$\omega_{i}$. (a) Approximating polynomial $y_{\rm{c}}^{(9)}(x)$ obtained by
multiplying Equation (25) by $x$. (b) Approximating polynomial
$y_{\rm{c}}^{(8)}(x)$ obtained by $y(x)-y_{\rm{c}}(x)$. (c) Plot of
$\omega_{i},\,i=2,{3,}\ldots,9$. ( ) $\overline{\omega_{i}}\approx 0.66$. (d)
Plot of $\omega_{i},\,i=2,{3,}\ldots,8$. ( ) $\overline{\omega_{i}}\approx
0.61$.

######  Example 5.

We study a variation of Example 4, in which the source term
$f(x)=\frac{1}{\sqrt{x}}$. Even though the source term has a singularity at
$x=0$, its definite integral over the domain $[0,1]$ is finite. The exact
solution [62] is

$y(x)={\rm{c}_{1}}\exp(-10x)+{\rm{c}_{2}}\exp(10x)+\exp(10x)\left(-5\left(\sqrt{\frac{\pi}{10}}\text{erf}\left(\sqrt{10x}\right)\right)\right)\\\
-\exp(-10x)\left(-5\left(\sqrt{\frac{\pi}{10}}\text{erfi}\left(\sqrt{10x}\right)\right)\right),$

where

${\rm{c}}_{1}=-{\rm{c}}_{2}=-\frac{\sqrt{\frac{5\pi}{2}}\left(e^{20}\text{erf}\left(\sqrt{10}\right)-\text{erfi}\left(\sqrt{10}\right)\right)}{e^{20}-1},$

$\mathrm{\erf}(x)=\frac{2}{\sqrt{\pi}}\int_{0}^{x}e^{-t^{2}}\differential{t}$
([55] § 7.1.1) is the error function,
$\mathrm{erfi}(x)\equiv-\imath\erf(\imath\,x)=\frac{2}{\sqrt{\pi}}\int_{0}^{x}e^{t^{2}}\differential{t}$
[63], and $\imath^{2}=-1$. The solution $y(x)$ has a boundary layer near
$x=0$. Equation (25) is multiplied by the factor $\sqrt{x}$ so that the
residual $R(x)$ does not contain a singularity at $x=0$.

We set the number of Sinc points to be inserted in all partitions as
$m=2N+1=5$. The stopping criterion $\varepsilon_{\mathrm{stop}}=10^{-6}$ was
used. The algorithm terminates after $\kappa=7$ iterations and the number of
points $|{\mathsf{S}}|=1183$.

We performed least-squares fitting of the logarithm of the set $\mathsf{R}$ to
the logarithm of the upper bound (24). Figure 9a shows the least-squares
fitted model (24) to the set $\mathsf{R}$. Figure 9b shows the residual,
absolute local approximation error, and the mean value for the last iteration.
Fine partitions are formed near $x=0$ due to the presence of the boundary
layer, as seen in the plot of the $L^{2}$ norm values of the residual over the
partitions. The mean value $\overline{R_{7}}$ is below the threshold value
$10^{-6}$. The $L^{2}$ norm of the approximation error
$\|y(x)-y_{\rm{c}}^{(7)}(x)\|\approx 2.18\times 10^{-7}$.

a

b

Figure 9:  (a) Fitting the upper bound (24) with the set $\mathsf{R}$. (b)
Visualization of the residual, absolute local approximation error, and the
mean value.

The approximating polynomial ${y_{\rm{c}}^{(7)}(x)}$ and a proper subset of
the set of points $\mathsf{S}$ are shown in Figure 10a. The corresponding plot
for the statistic $\omega_{i}$ a is shown in Figure 10c. The statistic
$\omega_{i}$ oscillates around the median value ${\tilde{\omega}_{i}}\approx
0.49$. It was mentioned that Equation (25) was multiplied by the factor
$\sqrt{x}$ so that the residual $R(x)$ does not contain a singularity at
$x=0$. We replace the residual $R(x)$ with the quantity $y(x)-y_{\rm{c}}(x)$
in Algorithm 4 and the BVP (25) contains the term $1/\sqrt{x}$. We set the
number of Sinc points to be inserted in all partitions as $m=2N+1=5$. The
stopping criterion $\varepsilon_{\mathrm{stop}}=10^{-6}$ was used. The
algorithm terminates after $\kappa=6$ iterations and the number of points
$|{\mathsf{S}}|=595$, which is smaller than the one obtained by multiplying
the residual by $\sqrt{x}$. This is expected since the exact solution $y(x)$
is used. The approximating polynomial $y_{\rm{c}}^{(6)}(x)$ and a proper
subset of the set of points $\mathsf{S}$ are shown in Figure 10b. The
corresponding plot for the statistic $\omega_{i}$ is shown in Figure 10d. The
oscillations are decaying and the statistic $\omega_{i}$ is converging to an
asymptotic value. The mean value $\overline{\omega_{i}}$ is $0.58$.

######  Example 6.

We study a variation of Example 4, in which the source term
$f(x)=\frac{\exp(x)-1}{x}$. The source term has a removable singularity at
$x=0$, since $\lim\limits_{x\to 0}f(x)=1$. The algorithm can directly solve
the BVP, even though the residual contains the term $\frac{\exp(x)-1}{x}$. The
exact solution [62] is

$y(x)={\rm{c}_{1}}\exp(-10x)+{\rm{c}_{2}}\exp(10x)-5\exp(10x)(\text{Ei}(-9x)-\text{Ei}(-10x))\\\
+5\exp(-10x)(\text{Ei}(11x)-\text{Ei}(10x)),$

where

${\rm{c}}_{1}=\frac{5\left(e^{20}\text{Ei}(-10)-e^{20}\text{Ei}(-9)-\text{Ei}(10)+\text{Ei}(11)-e^{20}{\ln}\left(\frac{10}{9}\right)-e^{20}{\ln}\left(\frac{11}{10}\right)\right)}{e^{20}-1}$

and

${\rm{c}_{2}}=\frac{5\left(-e^{20}\text{Ei}(-10)+e^{20}\text{Ei}(-9)+\text{Ei}(10)-\text{Ei}(11)+{\ln}\left(\frac{10}{9}\right)+{\ln}\left(\frac{11}{10}\right)\right)}{e^{20}-1}.$

We set the number of Sinc points to be inserted in all partitions as
$m=2N+1=5$. The stopping criterion $\varepsilon_{\mathrm{stop}}=10^{-6}$ was
used. The algorithm terminates after $\kappa=8$ iterations and the number of
points $|{\mathsf{S}}|=605$.

Figure 11a shows the approximate solution $y_{\rm{c}}^{(8)}(x)$. A proper
subset of the set of points $\mathsf{S}$ is shown as red dots, which are
projected onto the approximate solution $y_{\rm{c}}^{(8)}(x)$. We plot the
statistic $\omega_{i}$ as a function of the iteration index
$i,\,i=2,{3,}\ldots,8$, in Figure 11b. The oscillations are decaying and the
statistic $\omega_{i}$ is converging to an asymptotic value. The mean value
$\overline{\omega_{i}}\approx 0.65$.

a

b

c

d

Figure 10: Visualization of the approximating polynomial and the statistic
$\omega_{i}$. (a) Approximating polynomial $y_{\rm{c}}^{(7)}(x)$ obtained from
multiplying Equation (25) by $x$. (b) Approximating polynomial
$y_{\rm{c}}^{(6)}(x)$ obtained from $y(x)-y_{\rm{c}}(x)$. (c) Plot of
$\omega_{i},\,i=2,{3,}\ldots,7$. ( ) ${\tilde{\omega}_{i}}\approx 0.49$. (d)
Plot of $\omega_{i},\,i=2,{3,}\ldots,6$. ( ) $\overline{\omega_{i}}\approx
0.58$.

a

b

Figure 11:  (a) The approximating polynomial $y_{\rm{c}}^{(8)}(x)$. A proper
subset of the set of points $\mathsf{S}$ is shown. (b) Plot of the statistic
$\omega_{i}$. ( ) $\overline{\omega_{i}}\approx 0.65$.

We perform the least-squares fitting of the logarithm of the set $\mathsf{R}$
to the logarithm of the upper bound (24). Figure 12a shows the least-squares
fitted model (24) to the set $\mathsf{R}$. Figure 12b shows the residual,
absolute local approximation error, and the mean value for the last iteration.
The mean value $\overline{R_{8}}$ is below the threshold value $10^{-6}$. The
$L^{2}$ norm of the approximation error $\|y(x)-y_{\rm{c}}^{(10)}(x)\|\approx
3.1\times 10^{-7}$.

a

b

Figure 12: (a) Fitting the upper bound (24) with the set $\mathsf{R}$. (b)
Visualization of the residual, absolute local approximation error, and the
mean value.

######  Example 7.

We study the BVP (25) with $a(x)=0.02,\,b(x)=1,\,c(x)=0$, and $f(x)=1$. The
exact solution is given by

$y(x)=\mathrm{c}_{1}+\mathrm{c}_{2}\exp(50x)+x,$

where $\mathrm{c}_{1}=-\mathrm{c}_{2}=\frac{1}{-1+\exp(50)}$. The solution
$y(x)$ has a boundary layer near $x=1$.

We set the number of Sinc points to be inserted in all partitions as
$m=2N+1=5$. The stopping criterion $\varepsilon_{\mathrm{stop}}=10^{-6}$ was
used. The algorithm terminates after $\kappa=9$ iterations and the number of
points $|{\mathsf{S}}|=1055$.

Figure 13a shows the approximate solution $y_{\rm{c}}^{(9)}(x)$. A proper
subset of the set of points $\mathsf{S}$ is shown as red dots, which are
projected onto the approximate solution $y_{\rm{c}}^{(9)}(x)$. We plot the
statistic $\omega_{i}$ as a function of the iteration index
$i,\,i=2,{3,}\ldots,9$, in Figure 13b. The oscillations are decaying and the
statistic $\omega_{i}$ is converging to an asymptotic value. The mean value
$\overline{\omega_{i}}$ is $0.62$.

a

b

Figure 13: (a) The approximating polynomial $y_{\rm{c}}^{(9)}(x)$. A proper
subset of the set of points $\mathsf{S}$ is shown. (b) Plot of the statistic
$\omega_{i}$. ( ) $\overline{\omega_{i}}\approx 0.62$.

We perform the least-squares fitting of the logarithm of the set $\mathsf{R}$
to the logarithm of the upper bound (24). Figure 14a shows the least-squares
fitted model (24) to the set $\mathsf{R}$. Figure 14b shows the residual,
absolute local approximation error, and the mean value for the last iteration.
The plot of the $L^{2}$ norm values of the residual over the partitions shows
that fine partitions are formed near $x=1$ due to the presence of the boundary
layer. One observation is that the mean value $\overline{R_{9}}$ is below the
threshold value $10^{-6}$. The $L^{2}$ norm of the approximation error
$\|y(x)-y_{\rm{c}}^{(9)}(x)\|\approx 2.36\times 10^{-8}$ and the threshold
value in [8] is $0.02$.

a

b

Figure 14: (a) Fitting the upper bound (24) with the set $\mathsf{R}$. (b)
Visualization of the residual, absolute local approximation error, and the
mean value.

In this example, we increase the number of points per partition to $m=2N+1=7$
to examine the effect of the increase on the convergence of the algorithm. The
algorithm terminates after $\kappa=5$ iterations and the number of points
$|{\mathsf{S}}|=350$. Figure 15 shows the set $\mathsf{R}$ for $m=2N+1=5$ Sinc
points and $m=2N+1=7$ Sinc points. It is observed that increasing the number
of Sinc points per partition leads to faster convergence and a fewer number of
iterations.

Figure 15: Plotting the set $\mathsf{R}$ for $m=2N+1=5$ Sinc points ( ) and
$m=2N+1=7$ Sinc points ( ).

######  Example 8.

We study the following BVP [64, 28, 65]

$-(\upsilon(x)y^{\prime})^{\prime}=2\left[1+\alpha(x-\bar{x})\left(\arctan(\alpha(x-\bar{x}))+\arctan(\alpha\bar{x})\right)\right]$
(26)

with boundary conditions $y(0)=y(1)=0$, where $\alpha>0$ and

$\upsilon(x)=\dfrac{1}{\alpha}+\alpha(x-\bar{x})^{2}.$

For large values of $\alpha$, the BVP (26) has an interior layer close to
$\bar{x}$ [28]. The exact solution is given by

$y(x)=(1-x)\left[\arctan(\alpha(x-\bar{x}))+\arctan(\alpha\bar{x})\right].$

We use the values reported in [64], i.e., $\alpha=100$ and $\bar{x}=0.36388$.
This value of $\bar{x}$ was chosen so that
$\lim\limits_{\alpha\to\infty}y(\bar{x}^{+})\approx 2$ [64].

We set the number of Sinc points to be inserted in all partitions as
$m=2N+1=7$. The stopping criterion $\varepsilon_{\mathrm{stop}}=10^{-12}$ was
used. The algorithm terminates after $\kappa=15$ iterations and the number of
points $|{\mathsf{S}}|=$ 21,469.

Figure 16a shows the approximate solution $y_{\rm{c}}^{(15)}(x)$. A proper
subset of the set of points $\mathsf{S}$ is shown as red dots, which are
projected onto the approximate solution $y_{\rm{c}}^{(15)}(x)$. We plot the
statistic $\omega_{i}$ as a function of the iteration index
$i,\,i=2,{3,}\ldots,15$, in Figure 16b. One finding is that the oscillations
are decaying and the statistic $\omega_{i}$ is converging to an asymptotic
value. The mean value $\overline{\omega_{i}}$ is $0.46$.

a

b

Figure 16: (a) The approximating polynomial $y_{\rm{c}}^{(15)}(x)$. A proper
subset of the set of points $\mathsf{S}$ is shown. (b) Plot of the statistic
$\omega_{i}$. ( ) $\overline{\omega_{i}}\approx 0.46$.

We perform the least-squares fitting of the logarithm of the set $\mathsf{R}$
to the logarithm of the upper bound (24), where the parameter $\delta$ is
multiplied by $10^{-12}$ and $10^{-12}\delta=O(10^{-19})$. Figure 17a shows
the least-squares fitted model (24) to the set $\mathsf{R}$. The residual,
absolute local approximation error, and the mean value for the last iteration
are shown in Figure 17b. Fine partitions are formed near $x=\overline{x}$ due
to the presence of the interior layer, as shown in the plot of the $L^{2}$
norm values of the residual over the partitions. The mean value
$\overline{R_{15}}$ is below the threshold
$\varepsilon_{\mathsf{stop}}=10^{-12}$.

We compare the $L^{2}$ norm of the approximation error of our adaptive
piecewise Poly-Sinc method with other methods in Table 1. Method [28] requires
a parameter for the construction of refinement intervals. The $L^{2}$ norm
value of the approximation error is smaller than those reported in [65, 28].

Table 1: Comparison of the $L^{2}$ norm of the approximation error for different methods. Method |  Adaptive PW PS |  [28] |  ([65] Table 13(g))
---|---|---|---
$\|y(x)-y_{\mathrm{method}}(x)\|$ |  $1.104\times 10^{-14}$ |  $3\times 10^{-5}$ |  $1.9\times 10^{-12}$

a

b

Figure 17: (a) Fitting the upper bound (24) with the set $\mathsf{R}$. (b)
Visualization of the residual, absolute local approximation error, and the
mean value.

######  Example 9.

We consider the BVP [66, 67]

$-\epsilon\,y^{\prime\prime}-x\,y^{\prime}=\epsilon\pi^{2}\cos(\pi x)+\pi
x\,\sin(\pi x),\quad x\in[-1,1],$

with boundary conditions $y(-1)=-2,\ y(1)=0$ and $\epsilon>0$ is a parameter.
The exact solution follows as

$y(x)=\cos(\pi x)+\dfrac{\erf(x/\sqrt{2\epsilon})}{\erf(1/\sqrt{2\epsilon})}.$

The exact solution has a shock layer near $x=0$ [66]. We set
$\epsilon=10^{-6}$ [67].

We set the number of Sinc points to be inserted in all partitions as
$m=2N+1=5$. The stopping criterion $\varepsilon_{\mathrm{stop}}=10^{-11}$ was
used. The algorithm terminates after $\kappa=16$ iterations and the number of
points $|{\mathsf{S}}|=18530$.

Figure 18a shows the approximate solution $y_{\rm{c}}^{(16)}(x)$. A proper
subset of the set of points $\mathsf{S}$ is shown as red dots, which are
projected onto the approximate solution $y_{\rm{c}}^{(16)}(x)$. We plot the
statistic $\omega_{i}$ as a function of the iteration index
$i,\,i=2,{3,}\ldots,16$, in Figure 18b. The oscillations are decaying and the
statistic $\omega_{i}$ is converging to an asymptotic value. The mean value
$\overline{\omega_{i}}\approx 0.55$.

We perform the least-squares fitting of the logarithm of the set $\mathsf{R}$
to the logarithm of the upper bound (24), where the parameter $\delta$ is
multiplied by the factor $10^{-9}$ and $10^{-9}\delta=O(10^{-12})$. Figure 19a
shows the least-squares fitted model (24) to the set $\mathsf{R}$. The
residual, absolute local approximation error, and the mean value for the last
iteration are shown in Figure 19b. The plot of the $L^{2}$ norm values of the
residual over the partitions show that fine partitions are formed near $x=0$
due to the presence of the shock layer. The mean value $\overline{R_{16}}$ is
below the threshold $\varepsilon_{\mathsf{stop}}=10^{-11}$.

We compare the supremum norm of the approximation error of our adaptive
piecewise Poly-Sinc method with other methods in Table 2. $B-$splines were
used as basis functions [67]. The supremum norm result is in the same order as
that of [67]. Our adaptive method can reach a smaller value if we set
$\varepsilon_{\mathsf{stop}}<10^{-11}$.

a

b

Figure 18: (a) The approximating polynomial $y_{\rm{c}}^{(16)}(x)$. A proper subset of the set of points $\mathsf{S}$ is shown. (b) Plot of the statistic $\omega_{i}$. Table 2: Comparison of the supremum norm of the approximation error for different methods. Method |  Adaptive PW PS |  [67]
---|---|---
$\|y(x)-y_{\mathrm{method}}(x)\|_{\infty}$ |  $1.215\times 10^{-10}$ |  $4.3\times 10^{-10}$

We observe that the processing times differ among the many problems we look
at. This is due to two factors: first, the fact that we are aiming for a very
exact outcome; and second, the fact that there are many sorts of challenges.
Therefore, listing the computing time in seconds is meaningless since it would
only indicate the computational power used on our machine, which will vary for
various users. Overall, it is evident that using adaptive techniques takes
longer than using a simple collocation approach or finite element methods with
a lower accuracy. In our examples, we used a stopping criterion of $10^{-6}$
or less. As a result, our goal is to complete the computation in the most
accurate way feasible rather than the quickest way possible. This will
inevitably lengthen the processing time for some problems, such as stiff or
layer problems.

a

b

Figure 19: (a) Fitting the upper bound (24) with the set $\mathsf{R}$. (b)
Visualization of the residual, absolute local approximation error, and the
mean value.

##  6 Conclusions

In this paper, we developed an adaptive piecewise collocation method based on
Poly-Sinc interpolation for the approximation of solutions to ODEs. We showed
the exponential convergence in the number of iterations of the a priori error
estimate obtained from the piecewise collocation method, and provided that a
good estimate of the exact solution $y(x)$ at the Sinc points exists. We used
a statistical approach for partition refinement, in which we computed the
fraction of a standard deviation as the ratio of the mean absolute deviation
to the sample standard deviation. We demonstrated by several examples that an
exponential error decay is observed for regular ODEs and ODEs whose exact
solutions exhibit an interior layer, a boundary layer, and a shock layer. We
showed that our adaptive algorithm can deliver results with high accuracy at
the expense of the slower computation for stiff ODEs.

Author Contributions: Conceptualization, M.Y. and G.B.; methodology, O.K.,
M.Y. and G.B.; software, O.K. and G.B.; validation, O.K., H.E.-S., M.Y. and
G.B.; formal analysis, O.K. and H.E.-S.; investigation, H.E.-S. and G.B.;
resources, not applicable; data curation, not applicable; writing—original
draft preparation, O.K.; writing—review and editing, G.B. and O.K.;
visualization, O.K.; supervision, H.E.-S., M.Y. and G.B.; project
administration, H.E.-S.; funding acquisition, there is no funding of this
article. All authors have read and agreed to the published version of the
manuscript.
Funding: M.Y. is funded by BMBF under the contract 05M20VSA.
Institutional Review Board Statement: Not applicable.
Informed Consent Statement:  Not applicable.
Data Availability Statement: Not applicable.
Acknowledgments: The authors thank Frank Stenger for his insightful
conversations and continuous guidance.
Conflicts of Interest: The authors declare no conflict of interest.

## References

  * [1] E. Hairer, S. P. Nørsett, and G. Wanner, _Solving Ordinary Differential Equations I: Nonstiff Problems_. Berlin, Heidelberg: Springer Berlin Heidelberg, 1993.
  * [2] E. Hairer and G. Wanner, _Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems_. Berlin, Heidelberg: Springer Berlin Heidelberg, 1996.
  * [3] E. Hairer, G. Wanner, and C. Lubich, _Geometric Numerical Integration: Structure-Preserving Algorithms for Ordinary Differential Equations_. Berlin, Heidelberg: Springer Berlin Heidelberg, 2006.
  * [4] U. M. Ascher, R. M. M. Mattheij, and R. D. Russell, _Numerical Solution of Boundary Value Problems for Ordinary Differential Equations_. Philadelphia, PA: SIAM, 1995.
  * [5] I. Stakgold, _Boundary Value Problems of Mathematical Physics: Volume I & II_. Philadelphia, PA: SIAM, 2000.
  * [6] O. Axelsson and V. A. Barker, _Finite Element Solution of Boundary Value Problems: Theory and Computation_. Philadelphia, PA: SIAM, 2001.
  * [7] H. B. Keller, _Numerical Methods for Two-Point Boundary-Value Problems_. Mineola, NY: Dover Publications, Inc., 2018.
  * [8] K. Eriksson, D. Estep, P. Hansbo, and C. Johnson, “Introduction to adaptive methods for differential equations,” _Acta numerica_ , vol. 4, pp. 105–158, 1995.
  * [9] K. Wright, “Adaptive methods for piecewise polynomial collocation for ordinary differential equations,” _BIT Numerical Mathematics_ , vol. 47, no. 1, pp. 197–212, 2007.
  * [10] Z. Tao, Y. Jiang, and Y. Cheng, “An adaptive high-order piecewise polynomial based sparse grid collocation method with applications,” _Journal of Computational Physics_ , vol. 433, p. 109770, 2021.
  * [11] A. Logg, “Multi-adaptive Galerkin methods for ODEs. I,” _SIAM J. Sci. Comput._ , vol. 24, no. 6, pp. 1879–1902, 2003.
  * [12] ——, “Multi-adaptive Galerkin methods for ODEs. II. Implementation and applications,” _SIAM J. Sci. Comput._ , vol. 25, no. 4, pp. 1119–1141, 2003/04.
  * [13] M. Baccouch, “Analysis of a posteriori error estimates of the discontinuous Galerkin method for nonlinear ordinary differential equations,” _Appl. Numer. Math._ , vol. 106, pp. 129–153, 2016.
  * [14] ——, “A posteriori error estimates and adaptivity for the discontinuous Galerkin solutions of nonlinear second-order initial-value problems,” _Appl. Numer. Math._ , vol. 121, pp. 18–37, 2017.
  * [15] Y. Cao and L. Petzold, “A posteriori error estimation and global error control for ordinary differential equations by the adjoint method,” _SIAM J. Sci. Comput._ , vol. 26, no. 2, pp. 359–374, 2004.
  * [16] B. Kehlet and A. Logg, “A posteriori error analysis of round-off errors in the numerical solution of ordinary differential equations,” _Numer. Algorithms_ , vol. 76, no. 1, pp. 191–210, 2017.
  * [17] K.-S. Moon, A. Szepessy, R. Tempone, and G. E. Zouraris, “A variational principle for adaptive approximation of ordinary differential equations,” _Numer. Math._ , vol. 96, no. 1, pp. 131–152, 2003.
  * [18] ——, “Convergence rates for adaptive approximation of ordinary differential equations,” _Numer. Math._ , vol. 96, no. 1, pp. 99–129, 2003.
  * [19] K.-S. Moon, E. von Schwerin, A. Szepessy, and R. Tempone, “An adaptive algorithm for ordinary, stochastic and partial differential equations,” in _Recent advances in adaptive computation_ , ser. Contemporary Mathematics. American Mathematical Society, Providence, RI, 2005, vol. 383, pp. 325–343.
  * [20] C. Johnson, “Error estimates and adaptive time-step control for a class of one-step methods for stiff ordinary differential equations,” _SIAM J. Numer. Anal._ , vol. 25, no. 4, pp. 908–926, 1988.
  * [21] D. Estep and D. French, “Global error control for the continuous Galerkin finite element method for ordinary differential equations,” _RAIRO Modél. Math. Anal. Numér._ , vol. 28, no. 7, pp. 815–852, 1994.
  * [22] D. Estep, V. Ginting, and S. Tavener, “A posteriori analysis of a multirate numerical method for ordinary differential equations,” _Comput. Methods Appl. Mech. Engrg._ , vol. 223/224, pp. 10–27, 2012.
  * [23] F. Stenger, “Polynomial function and derivative approximation of Sinc data,” _Journal of Complexity_ , vol. 25, no. 3, pp. 292–302, 2009.
  * [24] F. Stenger, M. Youssef, and J. Niebsch, _Improved Approximation via Use of Transformations_. New York, NY: Springer New York, 2013, pp. 25–49.
  * [25] O. A. Khalil, H. A. El-Sharkawy, M. Youssef, and G. Baumann, “Adaptive piecewise Poly-Sinc methods for function approximation,” 2022, submitted to Applied Numerical Mathematics.
  * [26] G. F. Carey and D. L. Humphrey, “Finite element mesh refinement algorithm using element residuals,” in _Codes for Boundary-Value Problems in Ordinary Differential Equations_ , B. Childs, M. Scott, J. W. Daniel, E. Denman, and P. Nelson, Eds. Berlin, Heidelberg: Springer Berlin Heidelberg, 1979, pp. 243–249.
  * [27] G. F. Carey, “Adaptive refinement and nonlinear fluid problems,” _Computer Methods in Applied Mechanics and Engineering_ , vol. 17-18, pp. 541–560, 1979.
  * [28] G. F. Carey and D. L. Humphrey, “Mesh refinement and iterative solution methods for finite element computations,” _International Journal for Numerical Methods in Engineering_ , vol. 17, no. 11, pp. 1717–1734, 1981.
  * [29] R. C. Geary, “The ratio of the mean deviation to the standard deviation as a test of normality,” _Biometrika_ , vol. 27, no. 3-4, pp. 310–332, 10 1935\.
  * [30] F. Stenger, _Handbook of Sinc Numerical Methods_. Boca Raton, FL: CRC Press, 2011.
  * [31] M. Youssef and R. Pulch, “Poly-Sinc solution of stochastic elliptic differential equations,” _Journal of Scientific Computing_ , vol. 87, no. 3, pp. 1–19, 2021, article no. 82.
  * [32] F. Stenger, H. A. M. El-Sharkawy, and G. Baumann, _The Lebesgue Constant for Sinc Approximations_. Cham, Switzerland: Springer International Publishing, 2014, pp. 319–335.
  * [33] M. Youssef, H. A. El-Sharkawy, and G. Baumann, “Lebesgue constant using Sinc points,” _Advances in Numerical Analysis_ , vol. 2016, 2016.
  * [34] M. Youssef and G. Baumann, “Collocation method to solve elliptic equations, bivariate Poly-Sinc approximation,” _J. Progress. Res. Math_ , vol. 7, no. 3, pp. 1079–1091, 2016.
  * [35] M. Youssef, “Poly-Sinc Approximation Methods,” Ph.D. dissertation, Mathematics Department, German University in Cairo, 2017.
  * [36] M. Youssef and G. Baumann, “Troesch’s problem solved by Sinc methods,” _Mathematics and Computers in Simulation_ , vol. 162, pp. 31–44, 2019.
  * [37] O. A. Khalil and G. Baumann, “Discontinuous Galerkin methods using poly-sinc approximation,” _Mathematics and Computers in Simulation_ , vol. 179, pp. 96–110, 2021.
  * [38] ——, “Convergence rate estimation of poly-Sinc-based discontinuous Galerkin methods,” _Applied Numerical Mathematics_ , vol. 165, pp. 527–552, 2021.
  * [39] J. Stoer and R. Bulirsch, _Introduction to Numerical Analysis_. New York, NY: Springer-Verlag, 2002.
  * [40] G. Baumann and F. Stenger, “Sinc-approximations of fractional operators: a computing approach,” _Mathematics_ , vol. 3, no. 2, pp. 444–480, 2015.
  * [41] E. A. Coddington, _An introduction to Ordinary Differential Equations_. Dover Publications, Inc., 1989\.
  * [42] J. Lund and K. L. Bowers, _Sinc Methods for Quadrature and Differential Equations_. Philadelphia: SIAM, 1992, vol. 32.
  * [43] F. Stenger, _Numerical Methods Based on Sinc and Analytic Functions_ , ser. Springer Series in Computational Mathematics. New York, NY: Springer-Verlag, 1993, vol. 20.
  * [44] M. Youssef and G. Baumann, “Solution of nonlinear singular boundary value problems using polynomial-Sinc approximation,” _Commun. Fac. Sci. Univ. Ank. Sér. A1 Math. Stat._ , vol. 63, no. 2, pp. 41–58, 2014.
  * [45] G. Baumann, Ed., _New sinc methods of numerical analysis_ , ser. Trends in Mathematics. Cham, Switzerland: Birkhäuser/Springer, 2021.
  * [46] F. J. Vesely, _Computational Physics : An Introduction_. New York, NY: Kluwer Academic/Plenum Publishers, 2001\.
  * [47] G. Carey and B. A. Finlayson, “Orthogonal collocation on finite elements,” _Chemical Engineering Science_ , vol. 30, no. 5, pp. 587–596, 1975.
  * [48] G. Strang, _Computational Science and Engineering_. Wellesley-Cambridge Press Wellesley, 2007.
  * [49] A.-M. Wazwaz, _Linear and Nonlinear Integral Equations: Methods and Applications_. Higher Education Press, Beijing, and Springer-Verlag Berlin Heidelberg, 2011.
  * [50] T. H. Cormen, C. E. Leiserson, R. L. Rivest, and C. Stein, _Introduction to Algorithms_. Cambridge, MA, USA: MIT press, 2009.
  * [51] B. Haasdonk and M. Ohlberger, “Reduced basis method for finite volume approximations of parametrized linear evolution equations,” _M2AN Math. Model. Numer. Anal._ , vol. 42, no. 2, pp. 277–302, 2008.
  * [52] M. A. Grepl, “Model order reduction of parametrized nonlinear reaction–diffusion systems,” _Computers & Chemical Engineering_, vol. 43, pp. 33–44, 2012.
  * [53] R. H. Nochetto, K. G. Siebert, and A. Veeser, “Theory of adaptive finite element methods: An introduction,” in _Multiscale, Nonlinear and Adaptive Approximation_ , R. DeVore and A. Kunoth, Eds. Berlin, Heidelberg: Springer Berlin Heidelberg, 2009, pp. 409–542.
  * [54] R. E. Walpole, R. H. Myers, S. L. Myers, and K. Ye, _Probability and Statistics for Engineers and Scientists_ , 9th ed. Boston, MA: Pearson Education, Inc., 2012.
  * [55] M. Abramowitz and I. A. Stegun, _Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables_ , 10th ed. New York, NY: Dover Publications, Inc., 1972.
  * [56] D. V. Cruz-Uribe and A. Fiorenza, _Variable Lebesgue Spaces: Foundations and Harmonic Analysis_. Springer Basel, 2013.
  * [57] W. Gautschi, _Numerical Analysis_. Boston, MA: Birkhäuser, 2012.
  * [58] R. A. Horn and C. R. Johnson, _Matrix Analysis_ , 2nd ed. Cambridge University Press, 2013.
  * [59] M. Youssef, H. A. El-Sharkawy, and G. Baumann, “Multivariate Lagrange interpolation at Sinc points: Error estimation and Lebesgue constant,” _Journal of Mathematics Research_ , vol. 8, no. 4, 2016.
  * [60] W. R. Inc., “Mathematica, Version 13.0.0,” 2021, Champaign, IL, USA. [Online]. Available: https://www.wolfram.com/mathematica
  * [61] B. Rivière, _Discontinuous Galerkin Methods for Solving Elliptic and Parabolic Equations: Theory and Implementation_. Philadelphia: SIAM, 2008.
  * [62] W. E. Boyce and R. C. DiPrima, _Elementary Differential Equations and Boundary Value Problems_. Hoboken, NJ: John Wiley & Sons, Inc., 2012.
  * [63] E. W. Ng and M. Geller, “A table of integrals of the error functions,” _Journal of Research of the National Bureau of Standards B_ , vol. 73, no. 1, pp. 1–20, 1969.
  * [64] H. H. Rachford and M. F. Wheeler, “An H-1-Galerkin Procedure for the Two-Point Boundary Value Problem,” in _Mathematical Aspects of Finite Elements in Partial Differential Equations_ , C. de Boor, Ed. Cambridge, MA, USA: Academic Press, 1974, pp. 353–382.
  * [65] P. Keast, G. Fairweather, and J. Diaz, “A computational study of finite element methods for second order linear two-point boundary value problems,” _Mathematics of Computation_ , vol. 40, no. 162, pp. 499–518, 1983.
  * [66] P. W. Hemker, _A Numerical Study of Stiff Two-Point Boundary Problems_ , ser. Mathematical Centre Tracts 80. Amsterdam: Mathematisch Centrum, 1977.
  * [67] U. Ascher, J. Christiansen, and R. D. Russell, “A collocation solver for mixed order systems of boundary value problem,” _Mathematics of Computation_ , vol. 33, no. 146, pp. 659–679, 1979.