# Stepwise global error control in Euler’s method using the DP853 triple and
the Taylor remainder term

J. S. C. Prentice
Senior Research Officer
Mathsophical Ltd.
Johannesburg, South Africa

###### Abstract

We report on a novel algorithm for controlling global error in a step-by-step
(stepwise) sense, in the numerical solution of a scalar, autonomous, nonstiff
or weakly stiff problem. The algorithm exploits the remainder term of a Taylor
expansion of the solution. It requires the use of the DP853 triple to solve an
auxiliary problem which, in turn, enables the remainder term to be determined.
A quenching process then allows the solution generated by Euler’s method to be
controlled. We have achieved tolerances on the relative global error as strict
as $10^{-10}.$

2010 MSC: 65L05, 65L06, 65L70

Key words and phrases: Ordinary differential equations, Local error, Global
error, Runge-Kutta, Error control

##  1 Introduction

We seek to solve numerically the real-valued scalar autonomous problem

$\displaystyle y^{\prime}\left(x\right)$ $\displaystyle=f\left(y\right)$ (1)
$\displaystyle y\left(x_{0}\right)$ $\displaystyle=y_{0}$ $\displaystyle x\in$
$\displaystyle\left[x_{0},x_{N}\right],$

using an algorithm based on Runge-Kutta methods, with stepwise (step-by-step)
global error control. Work has been done on stepwise estimation and
suppression of the global error in these types of problems, and we will
discuss this in a later section. There, we will also refer to our own earlier
work in this regard. In this brief introduction, it suffices to state that we
intend to describe a novel approach to controlling the global error in Euler’s
method applied to (1), in a stepwise sense, by exploiting the remainder term
of a low-order Taylor expansion. Our algorithm will make use of the DP853
embedded triple [9]. Remarkably enough, we have been able to achieve relative
global error tolerances as strict as $10^{-10}$.

##  2 Relevant Concepts

###  2.1 The Taylor-Lagrange function

Let $y\left(x\right)$ be a real-valued univariate function, and assume that
$y\left(x\right)$ is as differentiable as is required in this paper. Taylor’s
theorem [34] provides the following result:

$y\left(x\right)=y\left(x_{0}\right)+y^{\prime}\left(\xi_{x},\mu_{x}\right)\left(x-x_{0}\right)=y_{0}+f\left(\mu_{x}\right)\left(x-x_{0}\right),$
(2)

where $\mu_{x}=y\left(\xi_{x}\right)$ and $x_{0}<\xi_{x}<x.$ The second term
on the RHS is the remainder term, presented here in Lagrange form [1],[3]. By
differentiating (2) with respect to $x,$we find

$\displaystyle y^{\prime}\left(x\right)$
$\displaystyle=f\left(\mu_{x}\right)+\left(x-x_{0}\right)\frac{\partial
f\left(\mu_{x}\right)}{\partial\mu_{x}}\frac{d\mu_{x}}{dx}$
$\displaystyle=f\left(\mu_{x}\right)+f_{y}\left(\mu_{x}\right)\left(x-x_{0}\right)\frac{d\mu_{x}}{dx}.$
(3)

Now, the LHS of this expression exists (see (1)), and so we will assume that
the RHS also exists. Hence, we assume that $\frac{d\mu_{x}}{dx}$ exists. This
suggests that $\mu_{x}$ is a function of $x,$ and we will sometimes write
$\mu\left(x\right)$ or simply $\mu$ in place of $\mu_{x}.$ We will refer to
$\mu_{x}$ as the Taylor-Lagrange function, or simply Lagrange function, for
this particular situation. In (2) and (3), we have used the notation

$\displaystyle y^{\prime}\left(\xi_{x},\mu_{x}\right)$ $\displaystyle\equiv
y^{\prime}\left(\mu\left(x\right)\right)=f\left(\mu\left(x\right)\right)\equiv
f\left(\mu_{x}\right)$ $\displaystyle
y^{\prime\prime}\left(\xi_{x},\mu_{x}\right)$
$\displaystyle\equiv\frac{dy^{\prime}\left(\mu\left(x\right)\right)}{d\mu\left(x\right)}=\left.\frac{dy^{\prime}\left(y\right)}{dy}\right|_{y=\mu_{x}}=\left.\frac{df}{dy}\right|_{y=\mu_{x}}\equiv
f_{y}\left(\mu_{x}\right).$

Rearranging (3) and using (1) and (2) gives

$\displaystyle\frac{d\mu}{dx}$
$\displaystyle=\frac{f\left(y\right)-f\left(\mu\right)}{f_{y}\left(\mu\right)\left(x-x_{0}\right)}$
$\displaystyle=\frac{f\left(y_{0}+f\left(\mu\right)\left(x-x_{0}\right)\right)-f\left(\mu\right)}{f_{y}\left(\mu\right)\left(x-x_{0}\right)}$
(4) $\displaystyle\equiv g\left(x,\mu\right).$

This is an initial-value problem that can, in principle, be solved to yield
the Lagrange function $\mu\left(x\right)$ for a suitable initial value. Once
$\mu\left(x\right)$ is known, the remainder term
$f\left(\mu\right)\left(x-x_{0}\right)$ is easily computed. Then, adding the
remainder term to $y\left(x_{0}\right)$ yields, from (2), an approximation to
$y\left(x\right).$ For interest’s sake, we refer the reader to our recent work
in this regard [19].

###  2.2 Runge-Kutta Methods

The topic of Runge-Kutta (RK) methods is very broad, indeed. We will assume
the reader has suitable familiarity with these methods, and list Butcher [4]
as a good general reference. Here, we will simply establish notation and
terminology relevant to this paper.

To solve $\mu^{\prime}=g\left(x,\mu\right)$ numerically, we partition
$\left[x_{0},x_{N}\right]$ using the nodes $x_{0}<x_{1}<\ldots<x_{N}$ and we
use the RK method

$\mu_{i+1}=\mu_{i}+h_{i+1}\Phi\left(x_{i},\mu_{i}\right).$

Here, $\mu_{i}$ and $\mu_{i+1}$ are the approximate numerical solutions at
$x_{i}$ and $x_{i+1},$ respectively, $h_{i+1}\equiv x_{i+1}-x_{i}$ is the
stepsize, and $\Phi\left(x,\mu\right)$ is known as the increment function. The
increment function defines the RK method. We will make use of several RK
methods in this work - denoted RK1, RK4, RK7 and DP853 - and we will describe
them when appropriate. Also, the consistency property of RK methods means that

$\displaystyle\lim_{h_{i+1}\rightarrow 0}\Phi\left(x_{i},\mu_{i}\right)$
$\displaystyle=g\left(x_{i},\mu_{i}\right)$
$\displaystyle\lim_{h_{i+1}\rightarrow 0}\Phi_{\mu}\left(x_{i},\mu_{i}\right)$
$\displaystyle=g_{\mu}\left(x_{i},\mu_{i}\right)$ (5)
$\displaystyle\lim_{h_{i+1}\rightarrow
0}\Phi_{\mu\mu}\left(x_{i},\mu_{i}\right)$
$\displaystyle=g_{\mu\mu}\left(x_{i},\mu_{i}\right),$

where the subscript $\mu$ denotes partial differentiation with respect to
$\mu.$

####  2.2.1 The DP853 embedded triple

The DP853 embedded triple [9] comprises three RK methods (of orders 3, 5 and
8) that use differing combinations of the same internal stages to generate
their output. There are 12 stages in DP853, each of which has the general form

$k_{q}=h_{i+1}g\left(x_{i}+c_{q}h_{i+1},\mu_{i}^{V}+a_{q,1}k_{1}+a_{q,2}k_{2}+\ldots+a_{q,q-1}k_{q-1}\right)$

for each $q\in\left[1,2,3,\ldots,12\right].$ The methods are explicit, so that
each stage is dependent on the previous stages, and they must be computed
sequentially. We use the superscripts $L,H$ and $V$ to denote the third-,
fifth- and eighth-order solutions, respectively. Hence, the outputs of DP853
are given by

$\displaystyle\mu_{i+1}^{L}$ $\displaystyle=$
$\displaystyle\mu_{i}^{L}+\sum\limits_{q=1}^{12}b_{q}^{L}k_{q}\equiv\mu_{i}^{L}+h_{i+1}\Phi^{L}\left(x_{i},\mu_{i}^{L}\right)$
(6) $\displaystyle\mu_{i+1}^{H}$ $\displaystyle=$
$\displaystyle\mu_{i}^{H}+\sum\limits_{q=1}^{12}b_{q}^{H}k_{q}\equiv\mu_{i}^{H}+h_{i+1}\Phi^{H}\left(x_{i},\mu_{i}^{H}\right)$
$\displaystyle\mu_{i+1}^{V}$ $\displaystyle=$
$\displaystyle\mu_{i}^{V}+\sum\limits_{q=1}^{12}b_{q}^{V}k_{q}.\equiv\mu_{i}^{V}+h_{i+1}\Phi^{V}\left(x_{i},\mu_{i}^{V}\right)$

wherein we have implicitly defined coefficients $b_{q}^{\cdots}$ and
appropriate increment functions. It transpires that DP853 has $b_{q}^{H,V}=0$
for $q\in\left[2,3,4,5\right]$ and $b_{q}^{L}=0$ for
$q\in\left[2,3,4,5,6,7,8,10,11\right].$

We modify (6) slightly: we use $\mu_{i}^{H}$ as input for $\mu_{i+1}^{L},$
i.e.

$\mu_{i+1}^{L}=\mu_{i}^{H}+\sum\limits_{q=1}^{12}b_{q}^{L}k_{q}.$ (7)

We discuss the value of this approach in the next subsection.

###  2.3 Local and Global Error

The RK local error at $x_{i+1}$ is

$\delta_{i+1}\equiv\mu\left(x_{i+1}\right)-\left(\mu\left(x_{i}\right)+h_{i+1}\Phi\left(x_{i},\mu\left(x_{i}\right)\right)\right)$

and the RK global error at $x_{i+1}$ is

$\Delta_{i+1}\equiv\mu\left(x_{i+1}\right)-\mu_{i+1}=\mu\left(x_{i+1}\right)-\left(\mu_{i}+h_{i+1}\Phi\left(x_{i},\mu_{i}\right)\right),$

where $\mu\left(x_{i}\right)$ and $\mu\left(x_{i+1}\right)$ are the exact
solutions at $x_{i}$ and $x_{i+1},$ respectively. The difference between the
two concepts is subtle: the local error definition requires that
$\mu\left(x_{i}\right)$ is known.

An RK method of order $p,$ written RK$p,$ has

$\displaystyle\delta_{i+1}$ $\displaystyle\propto
h_{i+1}^{p+1}+\ldots=L_{i+1}h_{i+1}^{p+1}+\ldots$ (8)
$\displaystyle\Delta_{i+1}$ $\displaystyle\propto
h^{p}+\ldots=G_{i+1}h^{p}+\ldots$ (9)

In the latter expression, the symbol $h$ represents the mean stepsize over the
subinterval $\left[x_{0},x_{i+1}\right].$ The ellipses represent higher-order
terms.

The local and global errors are related. We have shown [20][21] that

$\Delta_{i+1}=\delta_{i+1}+\Delta_{i}+h_{i+1}\Phi_{\mu}\left(x_{i},\mu_{i}\right)\Delta_{i}+\ldots$
(10)

where $\Delta_{i}$ is the global error in the input $\mu_{i}.$ Hence, ignoring
higher order terms, (7) yields

$\Delta_{i+1}^{L}=\delta_{i+1}^{L}+\left(1+h_{i+1}\Phi_{\mu}\left(x_{i},\mu_{i}^{H}\right)\right)\Delta_{i}^{H}.$

Compare this with

$\Delta_{i+1}^{L}=\delta_{i+1}^{L}+\left(1+h_{i+1}\Phi_{\mu}\left(x_{i},\mu_{i}^{L}\right)\right)\Delta_{i}^{L}$

which follows from (6). The second term on the RHS is the contribution to
$\Delta_{i+1}^{L}$ due to the global error of the input $\mu_{i}^{H}$ or
$\mu_{i}^{L}.$ Since we expect
$\left|\Delta_{i}^{H}\right|\ll\left|\Delta_{i}^{L}\right|,$ we see that (7)
reduces the magnitude of $\Delta_{i+1}^{L}$ compared with what it would
otherwise have been.

###  2.4 Error Control

####  2.4.1 Local error

It is possible to adjust the stepsize so as to control the local error. If
$\varepsilon_{\delta}>0$ denotes a user-defined tolerance on the local error,
then

$\left|\delta_{i+1}\right|=\left|L_{i+1}h_{i}^{p+1}\right|\leqslant\varepsilon_{\delta}\Rightarrow
h_{i+1}\leqslant\left(\frac{\varepsilon_{\delta}}{\left|L_{i+1}\right|}\right)^{\frac{1}{p+1}}$

provides an estimate of the stepsize necessary to bound
$\left|\delta_{i+1}\right|.$ In this expression we have ignored higher-order
terms present in the expansion of $\delta_{i+1}$ (see (8)). To cater for this
omission, it is common practice to use a safety factor $0<\eta<1,$ so that we
have

$h_{i+1}=\eta\left(\frac{\varepsilon_{\delta}}{\left|L_{i+1}\right|}\right)^{\frac{1}{p+1}}$

as the desired value of $h_{i+1}.$ We will use $\eta=0.85$ throughout this
study. Furthermore, it is possible to control the local error density
$\rho_{i+1},$ or local error per unit step as it is often known,

$\rho_{i+1}\equiv\frac{\left|\delta_{i+1}\right|}{h_{i}}=\left|L_{i+1}h_{i+1}^{p}\right|\leqslant\varepsilon_{\rho}\Rightarrow
h_{i+1}=\eta\left(\frac{\varepsilon_{\rho}}{\left|L_{i+1}\right|}\right)^{\frac{1}{p}},$
(11)

where $\varepsilon_{\rho}=\varepsilon_{\delta}/h_{i+1}.$

There are two virtues in controlling the error density. Firstly, the exponent
is larger $\left(1>\frac{1}{p}>\frac{1}{p+1}\right)$ and, anticipating that
$\frac{\varepsilon_{\rho}}{\left|L_{i+1}\right|}\ll 1,$ we will find that
$h_{i+1}$ will be smaller, meaning that
$\left|\delta_{i+1}\right|\leqslant\varepsilon_{\delta}$ is satisfied, as
well.

Secondly, if $h_{i+1}$ is suitably small we find, from (10), the recursion

$\Delta_{i+1}=\delta_{i+1}+\Delta_{i}\Rightarrow\Delta_{i+1}\leqslant\left|\delta_{i+1}\right|+\Delta_{i}.$
(12)

With $\Delta_{0}=0$ (because the initial value is known exactly) and assuming
$\left|\delta_{i+1}\right|\leqslant\varepsilon_{\rho}h_{i+1},$ we find

$\Delta_{N}\leqslant\varepsilon_{\rho}\sum\limits_{i=0}^{N-1}h_{i+1}=\varepsilon_{\rho}\left(x_{N}-x_{0}\right).$

The recursion (12) is monotonically increasing, so that
$\varepsilon_{\rho}\left(x_{N}-x_{0}\right)$ is a rough bound on the global
error that can be made a priori. Hence, for a desired value of $\Delta_{N},$
we could estimate the necessary local tolerance $\varepsilon_{\rho}$ via

$\varepsilon_{\rho}\leqslant\frac{\Delta_{N}}{x_{N}-x_{0}}.$ (13)

We prefer to control relative error, as in

$\frac{\rho_{i+1}}{\left|\mu_{i+1}\right|}\equiv\frac{\left|\delta_{i+1}\right|}{\left|\mu_{i+1}\right|h_{i+1}}\leqslant\varepsilon_{\rho}\Rightarrow
h_{i+1}=\eta\left(\frac{\varepsilon_{\rho}\left|\mu_{i+1}\right|}{\left|L_{i}\right|}\right)^{\frac{1}{p}},$
(14)

when $\left|\mu_{i+1}\right|>1,$ and absolute error, as in (11), when
$\left|\mu_{i+1}\right|\leqslant 1.$ The difference between the two cases is
in the numerator on the RHS of (14).

For relative error control, we have
$\rho_{i+1}\leqslant\varepsilon_{\rho}\left|\mu_{i+1}\right|,$ and for
absolute error control we have $\rho_{i+1}\leqslant\varepsilon_{\rho}.$ These
two cases can be conveniently expressed as

$\rho_{i+1}\leqslant\varepsilon_{\rho}\max\left(1,\left|\mu_{i+1}\right|\right).$

The RHS is a hinge function (of $\left|\mu_{i+1}\right|),$ and is somewhat
different to the often-used linear expression

$\rho_{i+1}\leqslant\varepsilon_{\rho}\left(1+\left|\mu_{i+1}\right|\right).$

We see that the linear expression gives
$\rho_{i+1}\leqslant\varepsilon_{\rho}$ only when $\left|\mu_{i+1}\right|=0.$
For all other values of $\left|\mu_{i+1}\right|,$ the bound imposed on
$\rho_{i+1}$ is larger than $\varepsilon_{\rho}$ or, since
$\rho_{i+1}/\left|\mu_{i+1}\right|\leqslant\varepsilon_{\rho}\left(1/\left|\mu_{i+1}\right|+1\right),$
the bound imposed on $\rho_{i+1}/\left|\mu_{i+1}\right|$ is larger than
$\varepsilon_{\rho}.$ Hence, we prefer to use the hinge function, since it
respects the user-defined bound $\varepsilon_{\rho}.$

####  2.4.2 Global error

Global error control is a more difficult task, particularly in a stepwise
sense. The easiest implementable form of global error control requires an
estimation of the error after the integration has been done, and then solving
the problem again using a new, smaller stepsize. This is known as
reintegration, and is an a posteriori form of control. Much work has been done
regarding global error estimation (see, for example,
[2][5][6][7][8][10][11][12][16][17][32][33][36]). However, there are cases to
be made for controlling the global error as the integration proceeds
[23][26][27][28][29][30][31]. In this regard, much less work has been done.
The quasi-consistent peer methods (see [13][14][15][35][36] and many
references therein) have shown very good performance, particularly when
dealing with strongly stiff systems, but these methods are not quite stepwise
control methods. They typically require a ”reboot” of the method when a
certain tolerance has been breached, which results in the method having to
backtrack several, perhaps many, nodes, resulting in what is effectively a
partial reintegration. In this paper, and in our other work
[22][24][25][26][27], we insist on a backtrack of no more than a single node,
when required, using a process we refer to as quenching.

##  3 The Algorithm

We intend to construct an algorithm that will allow for local control of the
global error in a RK1 method, by determining the remainder term in (2). This
is a more subtle approach than the direct extrapolation technique we have
previously used [22][24][25]. We believe the best way to describe our
algorithm is in a step-by-step manner that follows the sequential
implementation of the algorithm itself.

  1. 1.

Using $f\left(y\right)$ and (1), we construct $g\left(x,\mu\right).$

  2. 2.

We choose a point $x_{\mu}$ very close to $x_{0},$ and use a high-order RK
method - such as RK4 or even RK7 - to find a solution $y_{\mu}$ of (1) at
$x_{\mu}.$ We assume $y_{\mu}$ to be sufficiently accurate so as to be
regarded as practically exact. From (2), we then have

$y_{\mu}-y_{0}-f\left(\mu\right)\left(x_{\mu}-x_{0}\right)=0,$

which can be solved using a nonlinear solver to find $\mu$ (see the Appendix;
there is an important subtlety to consider). We label this value $\mu_{1}$ and
we relabel $x_{\mu}$ as $x_{1}.$ Moreover, we observe that (2) is trivially
satisfied when $x=x_{0},$ so that $\mu$ is arbitrary at $x_{0}.$ There is no
harm, then, in setting $\mu_{0}=0.$ This labelling is merely a logistical
convenience, and it provides the ordered pairs
$\left\\{\left(x_{0},\mu_{0}\right),\left(x_{1},\mu_{1}\right)\right\\}$ as
the first two values in our numerical solution for $\mu.$ NB: We also relabel
$y_{\mu}$ as $y_{1}.$

  3. 3.

We now use a pair of RK methods to determine $\mu$ at $x_{2}=x_{1}+h_{2},$
with $\left(x_{1},\mu_{1}\right)$ as input. The stepsize $h_{2}$ is user-
defined (more on this point later). These RK methods are the third- and fifth-
order components of the DP853 embedded triple. We need two values of $\mu$ at
$x_{2},$ obtained with methods of differing order, to be able to implement
local error control. We will denote these two solutions $\mu_{2}^{L}$ and
$\mu_{2}^{H}$. Using DP853 sets the order parameter $p=3$ in the steps that
follow. However, other users may choose to use RK methods of different orders,
so we will persist in using the symbol $p$ for the sake of generality.

  4. 4.

The next step is to implement local error control (LEC) in $\mu$ at $x_{2}.$
We assume that $\mu_{2}^{H}$ is sufficiently accurate, relative to
$\mu_{2}^{L},$ such that

$\mu_{2}^{H}-\mu_{2}^{L}=L_{2}h_{2}^{p+1}\Rightarrow
L_{2}=\frac{\mu_{2}^{H}-\mu_{2}^{L}}{h_{2}^{p+1}}$

is feasible. We now set $\mu_{2}=\mu_{2}^{H}$ and apply the process

$\displaystyle\left|\mu_{2}\right|$ $\displaystyle>1\text{ and
}\left|L_{2}h_{2}^{p}\right|>\varepsilon_{\rho}\left|\mu_{2}\right|\Rightarrow
h_{2}=\eta\left(\frac{\varepsilon_{\rho}\left|\mu_{2}\right|}{\left|L_{2}\right|}\right)^{\frac{1}{p}}$
$\displaystyle\left|\mu_{2}\right|$ $\displaystyle\leqslant 1\text{ and
}\left|L_{2}h_{2}^{p}\right|>\varepsilon_{\rho}\Rightarrow
h_{2}=\eta\left(\frac{\varepsilon_{\rho}}{\left|L_{2}\right|}\right)^{\frac{1}{p}}$

to compute a new stepsize, if necessary. We refer to this error control as
primary LEC.

There is a potentially valuable refinement that we can implement at this
stage. Recall that we intend to determine the remainder term
$f\left(\mu\right)\left(x-x_{0}\right)$. If we assume
$\mu_{2}\approx\mu\left(x_{2}\right)-L_{2}h_{2}^{p+1}$ then we have

$\displaystyle f\left(\mu_{2}\right)\left(x_{2}-x_{0}\right)$
$\displaystyle=f\left(\mu\left(x_{2}\right)-L_{2}h_{2}^{p+1}\right)\left(x_{2}-x_{0}\right)$
$\displaystyle=f\left(\mu\left(x_{2}\right)\right)\left(x_{2}-x_{0}\right)-L_{2}h_{2}^{p+1}f_{y}\left(\mu\left(x_{2}\right)\right)\left(x_{2}-x_{0}\right).$

The term
$L_{2}h_{2}^{p+1}f_{y}\left(\mu\left(x_{2}\right)\right)\left(x_{2}-x_{0}\right)\approx
L_{2}h_{2}^{p+1}f_{y}\left(\mu_{2}\right)\left(x_{2}-x_{0}\right)$ can be
thought of as a local error in
$y_{0}+f\left(\mu_{2}\right)\left(x_{2}-x_{0}\right).$ We can apply the same
reasoning as above to compute a new stepsize for a given tolerance on this
particular local error, as in

$\displaystyle\frac{\left|L_{2}f_{y}\left(\mu_{2}\right)\left(x-x_{0}\right)\right|h_{2}^{p+1}}{h_{2}\max(1,\left|y_{0}+f\left(\mu_{2}\right)\left(x_{2}-x_{0}\right)\right|)}\leqslant\varepsilon_{\rho}$
$\displaystyle\Rightarrow h_{2}$
$\displaystyle=\eta\left(\frac{\varepsilon_{\rho}\max(1,\left|y_{0}+f\left(\mu_{2}\right)\left(x_{2}-x_{0}\right)\right|)}{\left|L_{2}f_{y}\left(\mu_{2}\right)\left(x_{2}-x_{0}\right)\right|}\right)^{\frac{1}{p}}.$
(15)

Note the hinge function in the numerator. The division by $h_{2}$ corresponds
to error per unit step control. If this value of $h_{2}$ is smaller than the
value obtained from the earlier process, then it is chosen as the new
stepsize. We refer to this error control as secondary LEC. Now, given that we
propagate $\mu_{1}^{H}$ in DP3 (see (7)) and that, for strict tolerances,
stepsizes may generally be small, $L_{2}h_{2}^{p+1}$ may be a good
approximation to the global error $\Delta_{2}.$ Consequently, (15) constitutes
a form of global error control, or at least global error suppression, in
$y_{0}+f\left(\mu_{2}\right)\left(x_{2}-x_{0}\right).$ This is because the
quantity
$\Delta_{2}f_{y}\left(\mu\left(x_{2}\right)\right)\left(x_{2}-x_{0}\right)$ is
the global error in $y_{0}+f\left(\mu_{2}\right)\left(x_{2}-x_{0}\right).$ The
nett result is that secondary LEC will likely lead to small global errors in
$y_{0}+f\left(\mu_{2}\right)\left(x_{2}-x_{0}\right),$ a favourable and useful
outcome. Secondary LEC should be used with some caution, however. In general,
the denominator contains the factor $\left(x_{i}-x_{0}\right)$ which could
become large for large intervals. This might lead to excessively small values
for $h_{i}.$ Perhaps a lower limit on the stepsize should be considered to
counteract this. Perhaps the secondary LEC stepsize should not be allowed to
be smaller than one-tenth of the primary LEC stepsize, say. We did not need to
use such a lower limit in our numerical calculations, but it is something that
is worth considering.

Effectively, what we have done here is to attempt to control the relative
local error per step in both $\mu_{2}$ and
$y_{0}+f\left(\mu_{2}\right)\left(x_{2}-x_{0}\right).$ We then compute a new
value for $\mu_{2}$ using the new stepsize (i.e. we repeat Step 3 using the
new stepsize), and then use the new values for $\mu_{2}^{L}$ and $\mu_{2}^{H}$
so obtained to find a new value for $L_{2}.$ In the event that no stepsize
adjustment is necessary, we simply proceed directly to step 5 (note: in the
event of repeating Step 3, we do not then repeat Step 4, as well - we move
directly to step 5).

  5. 5.

The global error in $\mu_{2}$ can now be estimated using

$\Delta_{2}=\mu_{2}^{V}-\mu_{2}^{L}$

where $\mu_{2}^{V}$ is the eighth-order solution from DP853. We prefer to
overestimate $\Delta_{2},$ so we do not use
$\Delta_{2}=\mu_{2}^{V}-\mu_{2}^{H},$ even though we have set
$\mu_{2}=\mu_{2}^{H}.$ Of course, it is the relative global error

$\Delta_{2}^{rel}\equiv\frac{\Delta_{2}}{\max\left(1,\left|\mu_{2}\right|\right)}$

that actually interests us.

  6. 6.

We use Euler’s Method (RK1) to find $y_{2},$ as in

$y_{2}=y_{1}+h_{2}f\left(y_{1}\right).$

It is, of course, the global error in $y_{2}$ that we ultimately wish to
control.

  7. 7.

We find the remainder term

$R_{2}=f\left(\mu_{2}\right)\left(x_{2}-x_{0}\right).$

  8. 8.

We estimate the relative global error in the remainder term $\Delta R_{2}$ by

$\displaystyle f\left(\mu_{2}\right)\left(x_{2}-x_{0}\right)$
$\displaystyle=f\left(\mu\left(x_{2}\right)-\Delta_{2}\right)\left(x_{2}-x_{0}\right)$
$\displaystyle\approx
f\left(\mu\left(x_{2}\right)\right)\left(x_{2}-x_{0}\right)-f_{y}\left(\mu_{2}\right)\left(x_{2}-x_{0}\right)\Delta_{2}$
$\displaystyle\text{ \ \ \
}+\frac{f_{yy}\left(\mu_{2}\right)\left(x_{2}-x_{0}\right)\Delta_{2}^{2}}{2}$
$\displaystyle\Rightarrow\Delta R_{2}$
$\displaystyle=\frac{f_{yy}\left(\mu_{2}\right)\left(x_{2}-x_{0}\right)\Delta_{2}^{2}-2f_{y}\left(\mu_{2}\right)\left(x_{2}-x_{0}\right)\Delta_{2}}{2\max(1,\left|R_{2}\right|)}.$

  9. 9.

We now have the Euler solution $y_{2}$ and the Taylor approximation

$\displaystyle y_{2}^{T}$ $\displaystyle=$ $\displaystyle
y_{0}+f\left(\mu_{2}\right)\left(x_{2}-x_{0}\right)$ $\displaystyle\approx$
$\displaystyle\underset{y\left(x_{2}\right)}{\underbrace{\left[y_{0}+f\left(\mu\left(x_{2}\right)\right)\left(x_{2}-x_{0}\right)\right]}}+\Delta
R_{2}.$

The term in square brackets is the exact solution $y\left(x_{2}\right).$
Therefore, if $\Delta R_{2}$ is suitably small,

$\Delta y_{2}\approx y_{2}^{T}-y_{2},$

where $\Delta y_{2}$ denotes the global error in $y_{2}.$ In fact, we measure
the relative global error

$\Delta y_{2}^{rel}\equiv\frac{\Delta
y_{2}}{\max(1,\left|y_{2}^{T}\right|)}\approx\frac{y_{2}^{T}-y_{2}}{\max(1,\left|y_{2}^{T}\right|)}$

and we test the inequality

$\left|\Delta y_{2}^{rel}\right|>\left|\varepsilon_{g}-\left|\Delta
y_{2}^{T}\right|\right|,$ (16)

where $\varepsilon_{g}>0$ is a user-defined tolerance on the relative global
error in the Euler solution, and

$\Delta
y_{2}^{T}\equiv\frac{f_{yy}\left(\mu_{2}\right)\left(x_{2}-x_{0}\right)\Delta_{2}^{2}-2f_{y}\left(\mu_{2}\right)\left(x_{2}-x_{0}\right)\Delta_{2}}{2\max(1,\left|y_{2}^{T}\right|)}.$

If (16) is true we quench the Euler solution $y_{2}$ by setting

$y_{2}=y_{0}+f\left(\mu_{2}\right)\left(x_{2}-x_{0}\right),$

confident that

$\left|\Delta
y_{2}^{rel}\right|\leqslant\max\left(\varepsilon_{g},\left|\Delta
y_{2}^{T}\right|\right).$

Of course, we prefer that $\left|\Delta y_{2}\right|$ be bounded by the
imposed tolerance $\varepsilon_{g},$ and it is likely that it will be if
strict local error control has been carried out earlier (small
$\varepsilon_{\rho})$. But even if $\left|\Delta
y_{2}^{T}\right|>\varepsilon_{g}$ so that $\left|\Delta y_{2}^{T}\right|$ is
effectively the bound on $\left|\Delta y_{2}^{rel}\right|,$ at least we know
the value of $\left|\Delta y_{2}^{T}\right|.$ In such a case, if $\left|\Delta
y_{2}^{T}\right|$ is acceptably small, we would still be satisfied with the
result.

NB: It is imperative to track the growth of $\left|\Delta
y_{\cdots}^{T}\right|$ and, if $\left|\Delta y_{\cdots}^{T}\right|$ becomes
too large, we backtrack to a node $x_{k}$ where $\left|\Delta
y_{k}^{T}\right|$ is acceptably small, and then repeat the process from Step
1. That is to say, we treat $\left(x_{k},y_{k}^{T}\right)$ as a new initial
value $\left(x_{0},y_{0}\right)$ and ”reboot” the algorithm from this new
initial value. This reboot process is executed only if it is deemed necessary.

Obviously, if the test (16) is false, we do not quench the solution, and
simply proceed to the next step. An argument can be made that one should
always quench while $\left|\Delta y_{\cdots}^{T}\right|$ is acceptably small.
We have chosen not to do this in our numerical examples, in order to make
explicit the effects of quenching and not quenching. Nevertheless, persistent
quenching does have its merits.

  10. 10.

We now have a solution $y_{2}$ that satisfies a tolerance on its global
accuracy. Furthermore, this error control has been implemented in a stepwise
(local) manner. We are ready to compute $y_{3}.$ We return to Step 3, update
the indices by one and continue until we eventually reach the endpoint
$x_{N}$.

##  4 Numerical Examples

We consider the following examples for the purpose of numerical experiments:

Table 1: Details of numerical examples |
---|---
# | ODE | Interval | $y_{0}$ | Solution
1 | $y^{\prime}=y$ | $\left[0,5\right]$ | $2$ | $y\left(x\right)=2e^{x}$
2 | $y^{\prime}=y^{2}$ | $\left[-10,-3\right]$ | $0.1$ | $y\left(x\right)=-\frac{1}{x}$
3 | $y^{\prime}=\frac{y}{4}\left(1-\frac{y}{20}\right)$ | $\left[0,20\right]$ | $1$ | $y\left(x\right)=\frac{20}{1+19e^{-x/4}}$
4 | $y^{\prime}=\frac{1}{y}$ | $\left[5,25\right]$ | $1$ | $y\left(x\right)=\sqrt{2x-9}$
5 | $y^{\prime}=\cos y$ | $\left[a,b\right]$ | $-1$ | $x=\ln\left(\sec y+\tan y\right)$
6 | $y^{\prime}=-y$ | $\left[0,10\right]$ | $1$ | $y\left(x\right)=e^{-x}$

where, in #5,

$\displaystyle a$ $\displaystyle=$
$\displaystyle-1.2261911708835170708130609674719$ $\displaystyle b$
$\displaystyle=$ $\displaystyle 1.2261911708835170708130609674719.$

Some of these problems are weakly stiff, and the solution in #5 is given
implicitly.

Our results are summarized in Tables $2-8$. In these tables we have used the
notation

$\text{ \ \ \ }\begin{array}[]{ccl}\max\Delta&\text{
}&\begin{array}[]{l}\text{largest value of }\left|\Delta
y_{i}^{rel}\right|\text{ on the relevant interval.}\\\ \end{array}\\\
N&&\begin{array}[]{l}\text{Number of nodes used.}\\\ \end{array}\\\
\text{Q}&&\begin{array}[]{l}\text{Number of nodes where quenching was
needed.}\\\ \end{array}\end{array}$ $\text{ \ \ \
}\begin{array}[]{ccl}\text{P}&&\begin{array}[]{l}\text{Number of nodes where
local error control in }\mu_{i}\text{ }\\\ \text{was primary.}\end{array}\\\
\text{S}&&\begin{array}[]{l}\text{Number of nodes where local error control in
}\mu_{i}\text{ }\\\ \text{was secondary.}\end{array}\\\
\bigstar&&\begin{array}[]{l}\text{Number of steps where stepsize control was
due to stability }\\\ \text{considerations (see the later section {Stiffness
and Stability}).}\end{array}\\\ h_{2}&&\begin{array}[]{l}\text{size of
}h_{2}=x_{2}-x_{1},\text{rounded to first significant digit.}\\\
\lx@intercol\hfil\hfil\lx@intercol\end{array}\\\
\max\Delta_{E}&&\begin{array}[]{l}\text{largest value of }\left|\Delta
y_{i}^{rel}\right|\text{ using the }N\text{ nodes in}\\\ \text{RK1 without any
error control at all.}\end{array}\end{array}$

Table 2: Numerical results for $\varepsilon_{g}=10^{-2},\varepsilon_{\rho}=10^{-4}$ |
---|---
# | $\max\Delta$ | $N$ | Q | P | S | $\bigstar$ | $h_{2}$ | $\max\Delta_{E}$
1 | $9.3\times 10^{-3}$ | $71$ | $16$ | $0$ | $0$ | $1$ | $1.4\times 10^{-3}$ | $2.0\times 10^{-1}$
2 | $8.8\times 10^{-3}$ | $91$ | $1$ | $0$ | $0$ | $1$ | $1.4\times 10^{-3}$ | $1.3\times 10^{-2}$
3 | $9.4\times 10^{-3}$ | $221$ | $1$ | $0$ | $0$ | $1$ | $1.4\times 10^{-3}$ | $1.2\times 10^{-2}$
4 | $6.3\times 10^{-3}$ | $221$ | $0$ | $0$ | $0$ | $1$ | $1.4\times 10^{-3}$ | $6.3\times 10^{-3}$
5 | $8.9\times 10^{-3}$ | $46$ | $2$ | $0$ | $0$ | $1$ | $1.4\times 10^{-3}$ | $2.1\times 10^{-2}$
6 | $9.2\times 10^{-3}$ | $121$ | $2$ | $0$ | $0$ | $1$ | $1.4\times 10^{-3}$ | $1.5\times 10^{-2}$

Table 3: Numerical results for $\varepsilon_{g}=10^{-4},\varepsilon_{\rho}=10^{-6}$ |
---|---
# | $\max\Delta$ | $N$ | Q | P | S | $\bigstar$ | $h_{2}$ | $\max\Delta_{E}$
1 | $8.0\times 10^{-5}$ | $88$ | $75$ | $63$ | $0$ | $1$ | $1.4\times 10^{-3}$ | $1.6\times 10^{-1}$
2 | $9.2\times 10^{-5}$ | $91$ | $32$ | $0$ | $0$ | $1$ | $1.4\times 10^{-3}$ | $1.3\times 10^{-2}$
3 | $9.4\times 10^{-5}$ | $221$ | $115$ | $0$ | $0$ | $1$ | $1.4\times 10^{-3}$ | $1.2\times 10^{-2}$
4 | $9.5\times 10^{-5}$ | $223$ | $64$ | $11$ | $0$ | $1$ | $1.4\times 10^{-3}$ | $5.6\times 10^{-3}$
5 | $7.9\times 10^{-5}$ | $63$ | $47$ | $30$ | $8$ | $1$ | $1.4\times 10^{-3}$ | $1.5\times 10^{-2}$
6 | $9.4\times 10^{-5}$ | $121$ | $53$ | $63$ | $0$ | $1$ | $1.4\times 10^{-3}$ | $1.5\times 10^{-2}$

Table 4: Numerical results for $\varepsilon_{g}=10^{-6},\varepsilon_{\rho}=10^{-8}$ |
---|---
# | $\max\Delta$ | $N$ | Q | P | S | $\bigstar$ | $h_{2}$ | $\max\Delta_{E}$
1 | $7.2\times 10^{-7}$ | $327$ | $324$ | $310$ | $0$ | $1$ | $1.2\times 10^{-3}$ | $3.9\times 10^{-2}$
2 | $7.2\times 10^{-7}$ | $108$ | $91$ | $0$ | $41$ | $1$ | $1.4\times 10^{-3}$ | $1.0\times 10^{-2}$
3 | $9.1\times 10^{-7}$ | $221$ | $211$ | $0$ | $0$ | $1$ | $1.4\times 10^{-3}$ | $1.2\times 10^{-2}$
4 | $5.4\times 10^{-7}$ | $311$ | $308$ | $142$ | $0$ | $1$ | $1.0\times 10^{-3}$ | $2.1\times 10^{-3}$
5 | $8.9\times 10^{-7}$ | $214$ | $209$ | $158$ | $38$ | $1$ | $1.4\times 10^{-3}$ | $3.9\times 10^{-3}$
6 | $9.2\times 10^{-7}$ | $276$ | $251$ | $40$ | $219$ | $1$ | $1.2\times 10^{-3}$ | $4.0\times 10^{-3}$

Table 5: Numerical results for $\varepsilon_{g}=10^{-8},\varepsilon_{\rho}=10^{-10}$ |
---|---
# | $\max\Delta$ | $N$ | Q | P | S | $\bigstar$ | $h_{2}$ | $\max\Delta_{E}$
1 | $8.3\times 10^{-15}$ | $1474$ | $1472$ | $1459$ | $0$ | $1$ | $2.7\times 10^{-4}$ | $8.6\times 10^{-3}$
2 | $7.7\times 10^{-9}$ | $330$ | $323$ | $86$ | $219$ | $0$ | $1.0\times 10^{-3}$ | $3.4\times 10^{-3}$
3 | $7.9\times 10^{-9}$ | $695$ | $692$ | $663$ | $0$ | $0$ | $5.4\times 10^{-4}$ | $4.8\times 10^{-3}$
4 | $1.1\times 10^{-15}$ | $951$ | $949$ | $890$ | $0$ | $0$ | $3.3\times 10^{-4}$ | $6.5\times 10^{-4}$
5 | $8.9\times 10^{-10}$ | $949$ | $946$ | $754$ | $176$ | $1$ | $3.5\times 10^{-4}$ | $9.0\times 10^{-4}$
6 | $9.5\times 10^{-9}$ | $1233$ | $1223$ | $206$ | $1012$ | $1$ | $2.6\times 10^{-4}$ | $8.9\times 10^{-4}$

Table 6: Numerical results for $\varepsilon_{g}=10^{-10},\varepsilon_{\rho}=10^{-12}$ |
---|---
# | $\max\Delta$ | $N$ | Q | P | S | $\bigstar$ | $h_{2}$ | $\max\Delta_{E}$
1 | $1.5\times 10^{-14}$ | $6776$ | $6774$ | $6766$ | $0$ | $0$ | $1.1\times 10^{-4}$ | $1.9\times 10^{-3}$
2 | $4.7\times 10^{-11}$ | $1460$ | $1457$ | $425$ | $1011$ | $0$ | $2.2\times 10^{-4}$ | $7.6\times 10^{-4}$
3 | $1.8\times 10^{-15}$ | $3139$ | $3137$ | $3115$ | $0$ | $0$ | $1.2\times 10^{-4}$ | $1.1\times 10^{-3}$
4 | $1.2\times 10^{-15}$ | $4361$ | $4359$ | $4278$ | $0$ | $0$ | $6.1\times 10^{-5}$ | $1.5\times 10^{-4}$
5 | $8.6\times 10^{-11}$ | $4354$ | $4350$ | $3518$ | $820$ | $0$ | $1.4\times 10^{-4}$ | $2.0\times 10^{-4}$
6 | $2.2\times 10^{-15}$ | $5676$ | $5674$ | $975$ | $4691$ | $0$ | $1.0\times 10^{-4}$ | $1.9\times 10^{-4}$

Table 7: Numerical results for $\varepsilon_{g}=10^{-2},\varepsilon_{\rho}=10^{-3}$ |
---|---
# | $\max\Delta$ | $N$ | Q | P | S | $\bigstar$ | $h_{2}$ | $\max\Delta_{E}$
1 | $9.3\times 10^{-3}$ | $71$ | $16$ | $0$ | $0$ | $1$ | $1.4\times 10^{-3}$ | $2.0\times 10^{-1}$
2 | $8.8\times 10^{-3}$ | $91$ | $1$ | $0$ | $0$ | $1$ | $1.4\times 10^{-3}$ | $1.3\times 10^{-2}$
3 | $9.4\times 10^{-3}$ | $221$ | $1$ | $0$ | $0$ | $1$ | $1.4\times 10^{-3}$ | $1.2\times 10^{-2}$
4 | $6.3\times 10^{-3}$ | $221$ | $0$ | $0$ | $0$ | $1$ | $1.4\times 10^{-3}$ | $6.3\times 10^{-3}$
5 | $8.9\times 10^{-3}$ | $46$ | $2$ | $0$ | $0$ | $1$ | $1.4\times 10^{-3}$ | $2.1\times 10^{-2}$
6 | $9.2\times 10^{-3}$ | $121$ | $2$ | $0$ | $0$ | $1$ | $1.4\times 10^{-3}$ | $1.5\times 10^{-2}$

Table 8: Numerical results for $\varepsilon_{g}=10^{-6},\varepsilon_{\rho}=10^{-7}$ |
---|---
# | $\max\Delta$ | $N$ | Q | P | S | $\bigstar$ | $h_{2}$ | $\max\Delta_{E}$
1 | $9.5\times 10^{-7}$ | $162$ | $159$ | $141$ | $0$ | $1$ | $1.2\times 10^{-3}$ | $7.9\times 10^{-2}$
2 | $7.2\times 10^{-7}$ | $92$ | $75$ | $0$ | $7$ | $1$ | $1.4\times 10^{-3}$ | $1.2\times 10^{-2}$
3 | $9.1\times 10^{-7}$ | $221$ | $211$ | $0$ | $0$ | $1$ | $1.4\times 10^{-3}$ | $1.2\times 10^{-2}$
4 | $9.4\times 10^{-7}$ | $243$ | $240$ | $47$ | $0$ | $1$ | $1.4\times 10^{-3}$ | $3.6\times 10^{-3}$
5 | $8.9\times 10^{-7}$ | $110$ | $106$ | $69$ | $17$ | $1$ | $1.4\times 10^{-3}$ | $7.9\times 10^{-3}$
6 | $9.4\times 10^{-7}$ | $150$ | $137$ | $14$ | $71$ | $1$ | $1.2\times 10^{-3}$ | $8.3\times 10^{-3}$

In every single case, $\max\Delta$ is less than $\varepsilon_{g}.$ For any
given example, the value of $N$ tends to increase as $\varepsilon_{g}$ is
decreased. The number of steps where stepsize control was due to stability
considerations was never more than one for any of the examples. The value of
$1.4\times 10^{-3}$ for $h_{2}$ occurs often; this is explained in the next
section. The number of nodes where quenching was needed tends to increases as
$\varepsilon_{g}$ is decreased. LEC was not needed at all for
$\varepsilon_{g}=10^{-2},$ but for all other cases LEC was necessary and, in
some cases, at almost every node. Primary LEC occurred more often than
secondary LEC, but both types were definitely present. There were, however,
some cases where there was no secondary LEC even though primary LEC was
needed. We included $\max\Delta_{E}$ purely for interests’ sake - clearly
$\max\Delta$ is less than $\max\Delta_{E}$ for every case, usually by many
orders of magnitude (except #4 in Table 2, where no quenching or LEC was
required). In Tables $2-6,$ $\varepsilon_{\rho}=\varepsilon_{g}/100,$ and in
Tables 7 and 8, $\varepsilon_{\rho}=\varepsilon_{g}/10,$ suggesting that a
looser bound on LEC still yields acceptable results (see comments in the
Appendix regarding initial stepsize).

In Figures $1-10$ (for convenience, all figures are collected at the end of
the paper, after the Appendix), we show some error curves for two of the
examples (#3 and #5) for all five global tolerances considered here. In all of
these figures, we adopt the following conventions: horizontal red line is
$\varepsilon_{g};$ black dots indicate $\left|\Delta y_{i}^{rel}\right|,$ the
relative global error estimate in the RK1 solution; green curve indicates
$\left|\Delta y_{i}^{T}\right|,$ the estimate of the relative global error in
$y_{i}^{T}$; magenta circles indicate the true global error in $y_{i}$ at
those nodes where quenching was not needed; blue circles indicate the true
global error in $y_{i}$ at those nodes where quenching was implemented.
Clearly, quenching does not occur when $\left|\Delta y_{i}^{rel}\right|$ dips
below $\varepsilon_{g},$ as expected. In all cases, the true global error in
$y_{i}$ (either magenta circles or blue circles) lies below $\varepsilon_{g},$
as desired. Note that $\left|\Delta y_{i}^{T}\right|$ lies above the quenched
error, but below $\varepsilon_{g},$ thus serving the purpose of a good
estimator. The exception occurs when $\varepsilon_{g}=10^{-10},$ where
$\left|\Delta y_{i}^{T}\right|$ is similar to the quenched error, and both are
of the order of machine precision. This probably exposes the limit of the
algorithm’s capabilities, at least on the computational platform used here
[18], perhaps indicating that a tolerance this strict is not practical. We
also see that quenching becomes more prevalent as $\varepsilon_{g}$ decreases,
because the RK1 solution cannot satisfy such strict tolerances.

Although we only show results for two examples, the other examples exhibit
similar behaviour. Of course, details for all the examples are to be found in
the tables.

##  5 Reboot

In Step 9 above, we briefly mentioned the reboot process. In Figure 11, we
show the outcome of this process for example #3, for two reboot tolerances, on
the interval $\left[0,50\right].$ The reboot tolerance $\varepsilon_{rb}$ is
the maximum allowed magnitude of $\left|\Delta y_{i}^{T}\right|,$ and the
reboot process is designed to prevent $\left|\Delta y_{i}^{T}\right|$ from
becoming too significant with respect to $\varepsilon_{g}.$ In the first of
these figures, we use $\varepsilon_{rb}=10^{-10},$ and in the second we use
$\varepsilon_{rb}=10^{-7}.$ It is clear that no reboots were required in the
second case, but in the first case nine reboots were needed (the reboots occur
whenever $\left|\Delta y_{i}^{T}\right|$ breaches
$\varepsilon_{rb}=10^{-10})$. Admittedly, the value
$\varepsilon_{rb}=10^{-10}$ is strict, and was chosen for the purpose of
demonstration. Ordinarily, we would probably be satisfied with
$\varepsilon_{rb}=\varepsilon_{g}/1000.$ For the record, the case
$\varepsilon_{rb}=10^{-7}$ has $N=521,Q=125,P=0,S=0$ and $\bigstar=1,$ and the
case $\varepsilon_{rb}=10^{-10}$ has $N=709,Q=111,P=27$ and $S=0$. Each reboot
also required a single stepsize adjustment based on stability/stiffness (see
next section), giving $\bigstar=10.$

The reboot process is, effectively, a type of quench. The idea is to restart
the integration process at the reboot node $x_{rb}$ with the value of $y$
given by

$y_{rb}=y_{0}+f\left(\mu_{rb}^{V}\right)\left(x_{rb}-x_{0}\right),$

following which we set new values for the initial point

$\displaystyle x_{0}$ $\displaystyle=$ $\displaystyle x_{rb}$ $\displaystyle
y_{0}$ $\displaystyle=$ $\displaystyle y_{rb}$

and we proceed from Step 1 described earlier. The reboot node $x_{rb}$ is the
node that precedes the node at which $\left|\Delta y_{i}^{T}\right|$ breaches
$\varepsilon_{rb}.$ Hence, our reboot process requires a backtrack of only one
node. Note that, since we have changed the initial point, we need to determine
a new $g\left(x,\mu\right).$ We necessarily assume that using $\mu_{rb}^{V}$
in $y_{rb}$ will yield the most accurate value of $y_{rb}$ that we can
achieve.

##  6 Stiffness and Stability

We have

$\displaystyle g\left(x,\mu\right)$ $\displaystyle=$
$\displaystyle\frac{f\left(y_{0}+f\left(\mu\right)\left(x-x_{0}\right)\right)-f\left(\mu\right)}{f_{y}\left(\mu\right)\left(x-x_{0}\right)}$
$\displaystyle=$
$\displaystyle\frac{f\left(y_{0}+f\left(\mu\right)\left(x-x_{0}\right)\right)-f\left(y_{0}+\mu-
y_{0}\right)}{f_{y}\left(\mu\right)\left(x-x_{0}\right)}$
$\displaystyle\approx$
$\displaystyle\frac{f\left(y_{0}+f\left(\mu\right)\left(x-x_{0}\right)\right)-f\left(y_{0}\right)-f_{y}\left(y_{0}\right)\left(\mu-
y_{0}\right)}{\left(\frac{f_{y}\left(\mu\right)}{f\left(\mu\right)}\right)f\left(\mu\right)\left(x-x_{0}\right)}$
$\displaystyle=$
$\displaystyle\left(\frac{f\left(\mu\right)}{f_{y}\left(\mu\right)}\right)\left(\frac{f\left(y_{0}+f\left(\mu\right)\left(x-x_{0}\right)\right)-f\left(y_{0}\right)}{f\left(\mu\right)\left(x-x_{0}\right)}\right)$
$\displaystyle+\left(\frac{f\left(\mu\right)}{f_{y}\left(\mu\right)}\right)\left(\frac{f_{y}\left(y_{0}\right)\left(y_{0}-\mu\right)}{f\left(\mu\right)\left(x-x_{0}\right)}\right)$
$\displaystyle\approx$
$\displaystyle\frac{f\left(\mu\right)f_{y}\left(y_{0}\right)}{f_{y}\left(\mu\right)}+\frac{f_{y}\left(y_{0}\right)\left(y_{0}-\mu\right)}{f_{y}\left(\mu\right)\left(x-x_{0}\right)}$
$\displaystyle=$
$\displaystyle\frac{f_{y}\left(y_{0}\right)}{f_{y}\left(\mu\right)}\left(f\left(\mu\right)+\frac{y_{0}-\mu}{x-x_{0}}\right).$

Hence,

$g_{\mu}\left(x,\mu\right)=-\frac{f_{y}\left(y_{0}\right)f_{yy}\left(\mu\right)}{f_{y}^{2}\left(\mu\right)}\left(f\left(\mu\right)+\frac{y_{0}-\mu}{x-x_{0}}\right)+\frac{f_{y}\left(y_{0}\right)}{f_{y}\left(\mu\right)}\left(f_{y}\left(\mu\right)-\frac{1}{x-x_{0}}\right).$

For $\mu_{1}\approx y_{0}$ we find

$\displaystyle g_{\mu}\left(x_{1},\mu_{1}\right)$ $\displaystyle\approx$
$\displaystyle-\frac{f_{y}\left(y_{0}\right)f_{yy}\left(y_{0}\right)}{f_{y}^{2}\left(y_{0}\right)}\left(f\left(y_{0}\right)+\frac{y_{0}-y_{0}}{x_{1}-x_{0}}\right)$
$\displaystyle+\frac{f_{y}\left(y_{0}\right)}{f_{y}\left(y_{0}\right)}\left(f_{y}\left(y_{0}\right)-\frac{1}{x_{1}-x_{0}}\right)$
$\displaystyle=$
$\displaystyle-\frac{f_{yy}\left(y_{0}\right)f\left(y_{0}\right)}{f_{y}\left(y_{0}\right)}+f_{y}\left(y_{0}\right)-\frac{1}{x_{1}-x_{0}}$

and if $x_{1}\approx x_{0}$ and $y^{\prime}=f\left(y\right)$ is nonstiff or
weakly stiff, we would likely find that $\frac{-1}{x_{1}-x_{0}}$ is the
dominant term. So,

$g_{\mu}\left(x_{1},\mu_{1}\right)\approx-\left(\frac{1}{x_{1}-x_{0}}\right)=-1000$

for the values used in this study. We note that if $f_{y}\left(y_{0}\right)$
is small and positive, the first term on the RHS could contribute to stiffness
in $\mu^{\prime}=g\left(x,\mu\right);$ if $f_{y}\left(y_{0}\right)$ is
negative the second term contributes to stiffness in
$\mu^{\prime}=g\left(x,\mu\right);$ if $f_{y}\left(y_{0}\right)$ is negative
and has suitably small magnitude, the first term could cancel somewhat with
the third term, thus reducing the stiffness in
$\mu^{\prime}=g\left(x,\mu\right).$

For the examples considered here, we found

Table 9: $\mu_{1}$ and $g_{\mu}\left(x_{1},\mu_{1}\right)$ |
---|---
# | $y_{0}$ | $\mu_{1}$ | $g_{\mu}\left(x_{1},\mu_{1}\right)$ |
1 | $2$ | $2.0010$ | $-999.00$ |
2 | $0.1$ | $0.1000$ | $-999.85$ |
3 | $1$ | $1.0001$ | $-999.76$ |
4 | $1$ | $1.0005$ | $-999.99$ |
5 | $-1$ | $-0.9997$ | $-998.98$ |
6 | $1$ | $0.8885$ | $-1001.00$ |

consistent with our assumption $\mu_{1}\approx y_{0}$ and our result
$g_{\mu}\left(x_{1},\mu_{1}\right)\approx-1000.$

For stability, we demand that, when solving the Dahlquist equation

$\mu^{\prime}=-\lambda\mu,\text{ }\lambda>0,\text{
}\mu\left(0\right)=\mu_{0}=1$

using any of the RK methods in DP853, we find a solution at $x_{1}$ that has
the same properties as the true solution, i.e.

$0<e^{-\lambda x_{1}}<1.$

In other words, we demand that if

$\mu_{1}=\mu_{0}+h\Phi\left(x_{0},\mu_{0}\right),$

then

$0<\mu_{1}<1.$ (17)

Usually, stability with regard to stiffness concerns only the upper bound in
this expression, but we have observed that the methods in DP853 yield negative
values of $\mu_{1}$ for some values of $h.$ Hence, we have determined the
range of values of $h$ such that both bounds in (17) are satisfied by all
three methods in DP853. This has yielded

$-1.3764<-\lambda h<0.$

For our situation involving $\mu^{\prime}=g\left(x,\mu\right),$ the stiffness
constant $\lambda$ for the step $\left[x_{1},x_{2}\right]$ is estimated as

$\lambda=\left|g_{\mu}\left(x_{1},\mu_{1}\right)\right|$

so that

$h_{2}=\frac{1.3764}{\left|g_{\mu}\left(x_{1},\mu_{1}\right)\right|}\approx
1.3764\times 10^{-3}$ (18)

is the largest value of $h_{2}=x_{2}-x_{1}$ that we should be willing to use.
This value (rounded to $1.4\times 10^{-3})$ appears numerous times in the
tables of results presented earlier, as the value of $h_{2}$ that was used.
This restriction on the magnitude of the stepsize is performed at every step
prior to the numerical solution on that step being computed. Any stepsize
adjustment via LEC that is subsequently required is, of course, performed
after the numerical solution on that step has been found.

Note that we do not need to determine the stability regions of DP853 in the
complex plane, because we are dealing only with univariate problems in this
work. The stability interval $\left[-1.3764,0\right]$ on the negative real
axis is all that we require. For systems of differential equations, the
stiffness constant would be estimated using the eigenvalues of the Jacobian of
the system, which could be complex numbers. In such a case, the entire
stability region of each method in DP853 in the complex plane would need to be
known.

##  7 Future Research

We regard this paper as the first in a larger project, and so we have
restricted our study to that of scalar, autonomous, nonstiff or weakly stiff
problems. We anticipate that future research will address the following:

  1. 1.

Non-autonomous problems $y^{\prime}=f\left(x,y\right)$, which will require two
unknowns, $\xi$ and $\mu,$ to be determined for the remainder term. This means
that we will need to solve the system

$\left[\begin{array}[]{c}\xi^{\prime}\\\
\mu^{\prime}\end{array}\right]=\left[\begin{array}[]{c}g_{1}\left(x,\xi,\mu\right)\\\
g_{2}\left(x,\xi,\mu\right)\end{array}\right].$

  2. 2.

Non-autonomous systems

$\left[\begin{array}[]{c}y_{1}^{\prime}\\\ \vdots\\\
y_{m}^{\prime}\end{array}\right]=\left[\begin{array}[]{c}f_{1}\left(x,y_{1},\ldots,y_{m}\right)\\\
\vdots\\\ f_{m}\left(x,y_{1},\ldots,y_{m}\right)\end{array}\right],$

which we anticipate will yield $m+1$ equations for each component, giving
$m^{2}+m$ equations overall that will need to be solved to find the remainder
term.

  3. 3.

The use of tandem RK methods instead of DP853. For example, we may use RK2,
RK4, and RK7 to generate low-, high- and very high-order solutions for $\xi$
and $\mu,$ even if these RK methods are not embedded. If this proves
successful, it means that we are not bound to DP853, and that might provide
flexibility for the algorithm.

  4. 4.

Strongly stiff systems must be considered, wherein we might have to utilize
the implicit form of RK1 and suitably stable methods in place of DP853.
Indeed, we believe our algorithm has a modular character, so that the RK
methods we have used could be replaced with other methods that may have
properties more suitable to the problem at hand, whilst preserving the
structure of the algorithm as a whole. As a matter of record, we have solved
the stiff problem $y^{\prime}=-1000y$ quite easily using our algorithm,
although very small stepsizes were needed $(h\sim 10^{-3})$ due to the
explicit nature of DP853. Replacing DP853 with A-stable methods might be a
useful line of enquiry.

  5. 5.

The use of low-order RK methods other than Euler’s method.

  6. 6.

Can a stricter estimate of the global error $\left|\Delta y_{i}^{T}\right|$ be
made, and is it necessary to do so?

##  8 Concluding Comments

We have described a novel algorithm for controlling global error in a stepwise
sense. The algorithm exploits the remainder term of a Taylor expansion of the
solution, and requires the use of the DP853 triple to solve an auxiliary
problem which enables this remainder term to be determined. A quenching
process then allows the solution generated by Euler’s method to be controlled.
We have achieved tolerances on the relative global error as strict as
$10^{-10}.$ Admittedly, we have restricted the calculations in this paper to
scalar, autonomous, nonstiff or weakly stiff problems. Naturally, we hope that
future research will expand the scope of our algorithm.

## References

  * [1] T. Apostol, Calculus, Wiley, 1967.
  * [2] R. Aïd, L. Levacher, Numerical investigations on global error estimation for ordinary differential equations, J. Comput. Appl. Math., 82 (1997), 21-39
  * [3] E. Boman and R. Rogers, Lagrange’s Form of the Remainder, Chapter 5 in Real Analysis, LibreTexts${}^{\text{TM}}$ Mathematics (https://math.libretexts.org).
  * [4] J.C. Butcher, Numerical Methods for Ordinary Differential Equations, Wiley, Chichester, 2003.
  * [5] M. Calvo, D.J. Higham, J.I. Montijano, L. Rández, Global error estimation with adaptive explicit Runge–Kutta methods, IMA J. Numer. Anal., 16 (1996), 47-63
  * [6] M. Calvo, S. González-Pinto, J.I. Montijano, Global error estimation based on the tolerance proportionality for some adaptive Runge–Kutta codes, J. Comput. Appl. Math., 218 (2008), 329-341
  * [7] J.R. Dormand, R.R. Duckers, P.J. Prince, Global error estimation with Runge–Kutta methods, IMA J. Numer. Anal., 4 (1984), 169-184
  * [8] J.R. Dormand, M.A. Lockyer, N.E. McGorrigan, P.J. Prince, Global error estimation with Runge–Kutta triples, Comput. Math. Appl., 18 (1989), 835-846
  * [9] E. Hairer, S.P. Norsett and G. Wanner, Solving Ordinary Differential Equations I - Nonstiff Problems, Springer, Berlin, 2000\.
  * [10] G.Yu. Kulikov, R. Weiner, Doubly quasi-consistent parallel explicit peer methods with built-in global error estimation, J. Comput. Appl. Math., 233 (2010), 2351-2364
  * [11] G.Yu. Kulikov, R. Weiner, Global error estimation and control in linearly-implicit parallel two-step peer W-methods, J. Comput. Appl. Math., 236 (2011), 1226-1239
  * [12] G.Yu. Kulikov, Cheap global error estimation in some Runge–Kutta pairs, IMA J. Numer. Anal., 33 (2013), 136-163
  * [13] G.Yu. Kulikov, R. Weiner, Variable-stepsize interpolating explicit parallel peer methods with inherent global error control, SIAM J. Sci. Comput., 32 (2010), 1695-1723
  * [14] G.Yu. Kulikov, Global error control in adaptive Nordsieck methods, SIAM J. Sci. Comput., 34 (2012), A839-A860
  * [15] G.Yu. Kulikov, Adaptive Nordsieck formulas with advanced global error control mechanisms, Russian J. Numer. Anal. Math. Modelling, 28 (2013), 321-352
  * [16] J. Lang, J.G. Verwer, On global error estimation and control for initial value problems, SIAM J. Sci. Comput., 29 (2007), 1460-1475
  * [17] T. Macdougall, J.H. Verner, Global error estimators for order 7, 8 Runge–Kutta pairs, Numer. Algorithms, 31 (2002), 215-231
  * [18] Intel i7 9700K CPU, Strix Z390-F mainboard, 48GB DDR4 RAM, Windows 10 64-bit, Matlab R2022b.
  * [19] J.S.C. Prentice, Enhancing the accuracy of the Taylor polynomial by determining the remainder term, arXiv.org (Cornell University Library), 2023, 14p, [arXiv:2302.09216]
  * [20] J.S.C. Prentice, General error propagation in the RKGL method, J. Comput. Appl. Math., 228 (2009), 344–354, DOI: 10.1016/j.cam.2008.09.032
  * [21] J.S.C. Prentice, Error in Runge–Kutta methods, International Journal of Mathematical Education in Science and Technology, 44:3 (2013) 434-442, DOI: 10.1080/0020739X.2012.714492
  * [22] J.S.C. Prentice, Stepwise global error control in an explicit Runge-Kutta method using local extrapolation with high-order selective quenching, Journal of Mathematics Research, 3, 2 (2011) 126-136
  * [23] J.S.C. Prentice, Runge-Kutta Methods: Local error control does not imply global error control, Journal of Pure and Applied Mathematics: Advances and Applications, 6, 1 (2011) 71-84; also at [arXiv:1110.5728]
  * [24] J.S.C. Prentice, Relative Global Error Control in the RKQ Algorithm for Systems of Ordinary Differential Equations, Journal of Mathematics Research, 3, 4 (2011) 59-66
  * [25] J.S.C. Prentice, Nyström Methods in the RKQ Algorithm for Initial-value Problems, arXiv.org (Cornell University Library), 2011, 7p, [arXiv:1110.6749]
  * [26] J.S.C. Prentice, Global Error Control in the Runge-Kutta Solution of a Hamiltonian System, arXiv.org (Cornell University Library), 2011, 12p, [arXiv:1111.6996]
  * [27] J.S.C. Prentice, Global Error Control in a Lotka-Volterra System using the RKQ Algorithm, Journal of Mathematical and Computational Science, 2, 2 (2012) 413-424
  * [28] J.S.C. Prentice, Error-induced Stiffness in the Runge-Kutta Solution of a Lotka-Volterra System, Applied Mathematical Sciences, 6, 52 (2012) 2599-2604
  * [29] J.S.C. Prentice, Event Location and Global Error in Runge-Kutta Methods, Applied Mathematical Sciences, 6, 64 (2012) 3163-3167
  * [30] J.S.C. Prentice, The Relative Efficiency of Quenching and Reintegration for Global Error Control in Runge-Kutta Methods, Journal of Mathematical and Computational Science, 2, 4 (2012) 880-888
  * [31] J.S.C. Prentice, False convergence in the nonlinear shooting method, Journal of Mathematical and Computational Science, 5, 4 (2015) 454-461
  * [32] L.F. Shampine, L.S. Baca, Global error estimates for ODEs based on extrapolation methods, SIAM J. Sci. Stat. Comput., 6 (1985), 1-14
  * [33] R.D. Skeel, Thirteen ways to estimate global error, Numer. Math., 48 (1986), 1-20
  * [34] M. Spivak, Calculus 3rd ed., Publish or Perish, Houston, 1994\.
  * [35] R. Weiner, G.Yu. Kulikov, H. Podhaisky, Variable-stepsize doubly quasi-consistent parallel explicit peer methods with global error control, Appl. Numer. Math., 62 (2012), 1591-1603
  * [36] R. Weiner, G.Yu. Kulikov, Local and global error estimation and control within explicit two-step peer triples, J. Comp. Appl. Math., 262, (2014), 261-270

##  9 Appendix

###  9.1 Newton’s Method for $\mu_{1}$

We must solve

$y_{\mu}-y_{0}-f\left(\mu_{1}\right)\left(x_{\mu}-x_{0}\right)=0$

for $\mu_{1}.$ It is not enough to merely solve this equation in the form
presented here. Recall from (2) that we actually have

$y\left(x_{\mu}\right)=y\left(x_{0}\right)+y^{\prime}\left(\xi_{1},\mu_{1}\right)\left(x_{\mu}-x_{0}\right),$

where

$\mu_{1}\equiv y\left(\xi_{1}\right).$

Strictly speaking, then, we must solve

$F_{\xi}\left(\xi_{1}\right)\equiv
y_{\mu}-y_{0}-y^{\prime}\left(\xi_{1},y\left(\xi_{1}\right)\right)\left(x_{\mu}-x_{0}\right)=0$

to find $\xi_{1}.$ We use Newton’s Method in the following manner:

$\displaystyle\text{Use RK7 on }\left[x_{0},\xi_{1}^{j}\right]\text{ to
estimate }y\left(\xi_{1}^{j}\right)$ $\displaystyle\text{Use RK7 on
}\left[x_{0},\xi_{1}^{j}+\delta\xi\right]\text{ to estimate
}y\left(\xi_{1}^{j}+\delta\xi\right)$ $\displaystyle
F_{\xi}^{\prime}\left(\xi_{1}^{j}\right)=\frac{F_{\xi}\left(\xi_{1}^{j}+\delta\xi\right)-F_{\xi}\left(\xi_{1}^{j}\right)}{\delta\xi}$
$\displaystyle\xi_{1}^{j+1}=\xi_{1}^{j}-\frac{F_{\xi}\left(\xi_{1}^{j}\right)}{F_{\xi}^{\prime}\left(\xi_{1}^{j}\right)},$

where we use $\xi_{1}^{0}=\left(x_{\mu}+x_{0}\right)/2$ as the starting value
for the Newton iteration. We have used $\delta\xi=10^{-5}$ in all our
calculations, and we found that no more than three Newton iterations were ever
needed for all the examples that we have considered. The stopping criterion we
used for the iteration scheme was simply

$\left|F_{\xi}\left(\xi_{1}^{j}\right)\right|<10^{-14}.$

Once we have found an acceptable value for $\xi_{1},$ we use RK7 on
$\left[x_{0},\xi_{1}\right]$ to estimate $y\left(\xi_{1}\right),$ and we take
this estimate to be $\mu_{1}.$ In the implementations of RK7 here, we used
five steps on the relevant intervals. We note that if the problem is strongly
stiff, we would probably prefer to use an A-stable RK method instead of RK7.

###  9.2 Suggested local tolerance and initial stepsize

####  9.2.1 Local tolerance

Obtain a ”quick and dirty” solution for $\mu^{\prime}=g\left(x,\mu\right)$
using an RK method of moderate order, such as RK4 (which has only four
internal stages, and so is reasonably efficient for this purpose). Denote this
solution $\left(x_{j},\mu_{j}\right).$ Assume that the resulting relative
global error in $y_{j}=y_{0}+f\left(\mu_{j}\right)\left(x_{j}-x_{0}\right)$ is
given by

$\Delta_{j}^{rel}=\frac{f_{y}\left(\mu_{j}\right)\left(x_{j}-x_{0}\right)\delta\mu_{j}}{\max\nolimits_{j}\left(1,\left|y_{0}+R\left(x_{j},\mu_{j}\right)\right|\right)}=\frac{f_{y}\left(\mu_{j}\right)\left(x_{j}-x_{0}\right)\delta\mu_{j}}{\max\nolimits_{j}\left(1,\left|y_{0}+f\left(\mu_{j}\right)\left(x_{j}-x_{0}\right)\right|\right)}$

where $\delta\mu_{j}$ denotes the global error in $\mu_{j}.$ We demand that
this error be equal in magnitude, at most, to the specified global tolerance
$\varepsilon_{g}.$ This gives

$\displaystyle\left|\frac{f_{y}\left(\mu_{j}\right)\left(x_{j}-x_{0}\right)\delta\mu_{j}}{\max\nolimits_{j}\left(1,\left|y_{0}+f\left(\mu_{j}\right)\left(x_{j}-x_{0}\right)\right|\right)}\right|\leqslant\varepsilon_{g}$
$\displaystyle\Rightarrow$
$\displaystyle\left|\delta\mu_{j}\right|\leqslant\frac{\varepsilon_{g}\max\nolimits_{j}\left(1,\left|y_{0}+f\left(\mu_{j}\right)\left(x_{j}-x_{0}\right)\right|\right)}{\left|f_{y}\left(\mu_{j}\right)\left(x_{j}-x_{0}\right)\right|}.$

Choose

$\delta\mu_{m}\equiv\max\nolimits_{j}\left|\frac{\varepsilon_{g}\max\nolimits_{j}\left(1,\left|y_{0}+f\left(\mu_{j}\right)\left(x_{j}-x_{0}\right)\right|\right)}{\left|f_{y}\left(\mu_{j}\right)\left(x_{j}-x_{0}\right)\right|}\right|.$

Assume

$\delta\mu_{m}=\left(x_{N}-x_{0}\right)\varepsilon_{\rho}$

and so, using (13),

$\varepsilon_{\rho}=\frac{\delta\mu_{m}}{\left(x_{N}-x_{0}\right)}$

is an estimate for the value of $\varepsilon_{\rho}$ consistent with the
specified global tolerance $\varepsilon_{g}.$ We also consider a user-defined
default value for $\varepsilon_{\rho},$ denoted $\varepsilon_{\rho}^{D},$ and
we choose the local tolerance via

$\varepsilon_{\rho}=\min\left(\varepsilon_{\rho}^{D},\frac{\delta\mu_{m}}{\left(x_{N}-x_{0}\right)}\right).$

In this study we have generally used
$\varepsilon_{\rho}^{D}=\varepsilon_{g}/100,$ but we have done some
calculations using $\varepsilon_{\rho}^{D}=\varepsilon_{g}/10$ (see Tables 7
and 8).

Lastly, we note that if the problem is strongly stiff, we would probably
prefer to use an A-stable RK method instead of RK4.

####  9.2.2 Initial stepsize

Assuming that DP3 is formally equivalent to a Taylor expansion up to third-
order, we treat the fourth-order term

$\frac{\mu^{\left(iv\right)}\left(x_{1},\mu_{1}\right)h_{2}^{4}}{24}$

as the numerical error at $x_{2}$ obtained using DP3. The relevant stepsize
$h_{2}$ is estimated from

$\left|\frac{\mu^{\left(iv\right)}\left(x_{1},\mu_{1}\right)h_{2}^{4}}{24h_{2}\max(1,\left|\mu_{1}\right|)}\right|=\left|\frac{\mu^{\left(iv\right)}\left(x_{1},\mu_{1}\right)h_{2}^{3}}{24\max(1,\left|\mu_{1}\right|)}\right|=\varepsilon_{\rho}.$

Using

$\displaystyle\mu^{\left(iv\right)}=g^{\prime\prime\prime}$
$\displaystyle=\left(\frac{\partial}{\partial
x}+\frac{d\mu}{dx}\frac{\partial}{\partial\mu}\right)\left(\frac{\partial}{\partial
x}+\frac{d\mu}{dx}\frac{\partial}{\partial\mu}\right)\left(\frac{\partial}{\partial
x}+\frac{d\mu}{dx}\frac{\partial}{\partial\mu}\right)g$
$\displaystyle=\left(\frac{\partial}{\partial
x}+g\frac{\partial}{\partial\mu}\right)\left(\frac{\partial}{\partial
x}+g\frac{\partial}{\partial\mu}\right)\left(\frac{\partial}{\partial
x}+g\frac{\partial}{\partial\mu}\right)g$

gives

$h_{2}=\left|\frac{24\max(1,\left|\mu_{1}\right|)}{g^{\prime\prime\prime}\left(x_{1},\mu_{1}\right)}\right|^{\frac{1}{3}}.$
(19)

Naturally, we also specify a default value $h_{2}^{D},$ and we invoke the
stability consideration (18), so that

$h_{2}=\min\left(h_{2}^{D},\frac{1.3764}{\left|g_{\mu}\left(x_{1},\mu_{1}\right)\right|},\left|\frac{24\max(1,\left|\mu_{1}\right|)}{g^{\prime\prime\prime}\left(x_{1},\mu_{1}\right)}\right|^{\frac{1}{3}}\right).$

In our calculations we used $h_{2}^{D}=0.1,$ and found that the initial
stepsize was always set either by (18) or (19).


# Non Existence of Manifold Structures on Topological constructions

Mohammad Al Attar

(March 2022)

Abstract: We show the non existence of manifold structure on certain
topological constructions. The main results of this paper are: (1) The wedge
of manifold with itself is not a manifold. (2) The cone of a manifold is not a
manifold. The main tools for the main results are homological algebra and
algebraic topology.

Topological manifolds are spaces that resemble euclidean space locally. That
is, at a sufficiently small region, the manifold looks like (a region of)
$\mathbb{R}^{n}$. For instance, a line is a one dimensional manifold. A
surface is a two dimensional manifold and so on. A piece of paper is a real
world model of a two dimensional manifold. The earth, locally, looks like a 2
dimensional region, and is in some sense locally flat. However, to discuss the
notion of flatness from the point of view of curvature, one requires
additional structure on a topological manifold. In particular, a topological
manifold is not equipped with data the takes into account curvature. However,
manifolds that take into account curvature are called Riemannian manifolds,
which are, informally, manifolds, such that the notions of differentiation and
integration are well defined; in other words are smooth manifolds, with the
additional structure that you can define the notions of length on a linear
approximation of the manifold at a given region of any point. It is well known
that there are only two ’reasonable’ 1-dimensional manifolds from a
topological sense: The real line and the circle. It is clear that these two
spaces are topologically distinct. For 2-manifolds, the situation gets a bit
more complicated. In particular ’reasonable’ 2-dimensional manifolds are
characterised by the number of holes they have. The number of holes is called
the genus. For instance, the two sphere has no holes and therefore has genus
0. However, one may obtain a genus one surface by attaching a handle to the
two sphere; obtaining a torus. This process is an inductive process, yielding
an orientable manifold at each stage. However, one may also obtain non
orientable manifolds by attaching the boundary of the mobius band. In higher
dimensions, the situation gets much more complicated and is an active area of
research. In this paper, we mainly discuss topological manifolds of arbitrary
dimensions and show that certain topological constructions do not yield a
topological manifold structure. The background assumes a solid background in
algebraic topology. Hatcher’s book [1] and [3] are very good sources for
algebraic topology. $[2]$ is a nice introduction to topological manifolds.
$[4]$ Is the common reference for point set topology. The structure of this
paper is as follows. Chapter 2 discusses some backgrounds needed for and used
in chapter 3. The third chapter contains the main results. The author has
embarked on the goal to write complete proofs of the results discussed in that
section.

##  1 Background

In this section, we discuss some of the background necessary for the second
section. In pursuit of understanding the second section better, we define a
few key notions. We illustrate some of the notions defined with a few
examples. The examples in this section are standard and can be found in $[1]$
and $[2]$ along with proofs. Our notation follows as that in $[1]$ and $[3]$

######  Definition 1.1.

A topological manifold of dimension $n$ is a topological space $M$ that is
second countable, Hausdorff and locally euclidean of dimension $n$: For each
$p\in M$ there exists a pair $(U,\phi)$, where $U$ is an open subset of $M$
containing $p$ and $\phi:U\rightarrow\mathbb{R}^{n}$ is a homeomorphism. The
pair $(U,\phi)$ is called a chart on $M$. $U$ is called a coordinate domain ,
and $\phi$ is called a coordinate map. The collection of charts is called an
atlas on $M$. Letting $r^{1},...,r^{n}$ denote the standard coordinates on
$\mathbb{R}^{n}$, we have $\phi=(x^{1},...,x^{n})$ , where
$x^{i}=r^{i}\circ\phi$. Thus, one calls the data $(U,x^{1},...,x^{n})$ or more
simply $(x^{1},...,x^{n})$, a local coordinate system on $M$.

For a point $p$ in a manifold $M$, one typically identifies $p$ with its image
under a coordinate map. To see this in action, let $(U,x^{1}...,x^{n})$ be a
chart on $M$ . Then one may blur the distinction between $p$ and
$(x^{1}(p),...,x^{n}(p))$. The reason for such a temporary identification is
as follows: Given a point $p\in M$, a chart about $p$ identifies, atleast
topologically, a neighborhood of $p$ with an open region of $\mathbb{R}^{n}$.
Thus, we would like the coordinate system on the region surrounding the point
to be induced from the standard coordinates on euclidean space. This is indeed
the case, under our identification. More precisely, identifying $U$ with its
image under the coordinate map, we obtain $x^{i}=r^{i}|_{U}$. Note that this
identification is not canonical, since obviously different charts may give
rise to different coordinates of the point.

######  Example 1.1.

$\mathbb{R}^{n}$ for each $n$ with the standard topology is a manifold of
dimension $n$. For each $n$, $\mathbb{S}^{n}$ is a compact topological
manifold of dimension $n$. $\mathbb{R}^{n}$ has a global coordinate system.
Namely, $(\mathbb{R}^{n},r^{1},...,r^{n})$. $\mathbb{S}^{n}$ cannot have a
global coordinate system, for otherwise we would have a homeomorphism
$\mathbb{R}^{n}\approx\mathbb{S}^{n}$.

######  Definition 1.2.

Let $(C_{*},\partial)$ be a chain complex of abelian groups. Then,
$\partial^{2}=0$. We thus define for each n n-th homology group of $C_{*}$:
$H_{n}(C_{*})=\frac{im\partial_{n}}{ker\partial_{n+1}}$.

######  Example 1.2.

Let $X$ be a space and $A\subseteq X$ a subspace. Then we have a pair $(X,A)$.
For each $n$, $S_{n}(A)$ is the free abelian group generated by all the
singular n simplices in $A$. Since an singular n simplex in $A$ is an n
simplex in $X$, $S_{n}(A)\subseteq S_{n}(X)$. We have a short exact sequence
of chain complexes: $0\rightarrow S_{*}(A)\hookrightarrow S_{*}(X)\rightarrow
S_{*}(X,A)\rightarrow 0$.Where we define $S_{n}(X,A)$ to be the group
$S_{n}(X)$ modulo $S_{n}(A)$. We note that we may realize $S_{n}(X,A)$ as the
free abelian group generated by the $n$ simplices in $X$ with image not
contained in $A$.

######  Theorem 1.1.

(Excision) Let $X$ be a space and $Z\subseteq A\subseteq X$ with
$\overline{Z}\subseteq intA$. Then the inclusion
$(X-Z,A-Z)\hookrightarrow(X,A)$ induces an isomorphism $H_{n}(X-Z,A-Z)\cong
H_{n}(X,A)$ for each $n$.

######  Example 1.3.

Suppose $X$ is a $T_{1}$ space and $x\in X$. Then for any open neighborhood
$U$ of $x$,

  1. 1.

$H_{n}(X,U)\cong H_{n}(X-\\{x\\},U-\\{x\\})$.

  2. 2.

$H_{n}(X,X-\\{x\\})\cong H_{n}(U,U-\\{x\\})$.

First one is a direct consequence of excision. Let us establish the second
one: Note that $\\{x\\}\subseteq U$ and thus $X\backslash U\subseteq
X\backslash\\{x\\}$. Therefore, $H_{n}(X\backslash X\backslash
U,X\backslash\\{x\\}\backslash X\backslash U)=H_{n}(U,U\backslash\\{x\\})\cong
H_{n}(X,X\backslash\\{x\\})$.

######  Definition 1.3.

If $X$ is a space we define the $\textbf{cone }of$ $X$, denoted $CX$, to be
the quotient space $\frac{X\times[0,1]}{X\times 1}$. If $X$ and $Y$ are spaces
then choosing $x\in X$ and $y\in Y$, we define the wedge sum of $X$ and $Y$ to
be the quotient space $\frac{X\coprod Y}{x\sim y}$.

######  Example 1.4.

If $X$ and $Y$ are homeomorphic spaces, then $CX$ is homeomorphic to $CY$;
$CX\approx CY$. The reason for this is the functoriality of the cone map: That
is, if $f:X\rightarrow Y$ is continuous, then the cone map is given by
$Cf:CX\rightarrow CY$, $(Cf)([x,t]=[f(x),t]$. The cone map respects
composition and preserves the identity . The respecting of said properties
allows us to deduce that homeomorphic spaces have homeomorphic cones.

######  Theorem 1.2.

(Relative Mayer-Vietoris) Suppose we have a pair of space $(X,Y)=(intA\cup
intB,intC\cup intD)$ , where $C\subseteq A$ and $D\subseteq B$. Then, we have
the following long exact sequence in homology:

$...\rightarrow H_{n}(A\cap B,C\cap D)\rightarrow H_{n}(A,C)\oplus
H_{n}(B,D)\rightarrow H_{n}(X,Y)\rightarrow...$

##  2 Main Results

The proof technique of the next result follows from trying to get control of
the situation by using the Relative Mayer Vietoris sequence from algebraic
topology. The reason for doing so is to reduce problems about manifolds to
problems about local homology groups $H_{*}(X,X\backslash\\{x\\})$ for $x\in
X$.

######  Theorem 2.1.

Suppose $M$ is a manifold of dimension $n\geq 2$. Then, $M\vee M$ is not a
topological manifold.

###### Proof.

We first make the claim that if $M$ is an $n$ dimensional manifold and
$p_{1},...,p_{m}\in M$. Then

$H_{n}(M,M\backslash\\{p_{1},...,p_{m}\\})\cong\mathbb{Z}^{m}$

Indeed, it suffices by inducting on $m$ to establish the $m=2$ case. To that
end, choose open sets $U_{1}$, $U_{2}$ about $p_{1}$ and $p_{2}$ respectively
for which $U_{1}\cap U_{2}=\varnothing$. By excision,
$H_{n}(M,M\backslash\\{p_{1},p_{2}\\})$ $\cong H_{n}(U_{1}\cup
U_{2},U_{1}\backslash\\{p_{1}\\}\cup U_{2}\backslash\\{p_{2}\\})\cong
H_{n}(U_{1},U_{1}-p_{1})\oplus
H_{n}(U_{2},U_{2}-p_{2})\cong\mathbb{Z}\oplus\mathbb{Z}$. The second to last
isomorphism follows from the relative version of Mayer Vietoris. Now to prove
the theorem, suppose $M\vee M$ is a topological manifold of dimension $n$.
Then,

$H_{n}(M\vee M,M\vee M\backslash\\{\star\\})\cong\mathbb{Z}$.

Let $p_{1}$ and $p_{2}$ denote the base points. Now we claim:

$H_{n}(M\vee M,(M\vee M)\backslash\\{\star\\})\cong\mathbb{Z}\oplus\mathbb{Z}$

To prove the result above, we use the Relative Mayer Vietoris sequence.
Clearly it suffices to show:

$H_{n}(M\vee M,(M\vee M)\backslash\\{\star\\})\cong
H_{n}(M,M\backslash\\{p_{1}\\})\oplus H_{n}(M,M\backslash\\{p_{1}\\})$

To show the above, we invoke the Relative Mayer-Vietoris sequence, note that
$M\vee M=M_{1}\cup M_{2}$, where we have homeomorphisms $M_{1}\approx M$ and
$M_{2}\approx M$ that restrict to homeomorphisms
$M_{i}\backslash\\{\star\\}\approx M\backslash\\{p_{i}\\}$ for each $i$. Thus
consider $(M\vee M,M\vee M\backslash\\{\star\\})=(M_{1}\cup
M_{2},M_{1}\backslash\\{\star\\})\cup M_{2}\backslash\\{\star\\})$. Hence we
have by Relative Mayer-Vietoris

$...\rightarrow H_{n}(\star)\rightarrow
H_{n}(M_{1},M_{1}\backslash\\{\star\\})\oplus
H_{n}(M_{2},M_{2}\backslash\\{\star\\})\rightarrow H_{n}(M\vee M,M\vee
M\backslash\\{\star\\})\rightarrow H_{n-1}(\star)\rightarrow...$

Since $n\geq 2$, $H_{n}(\star)=H_{n-1}(\star)=0$ . Hence the long exact
sequence above gives the isomorphism we seek.

∎

The motivation behind the second result is the case when when we are dealing
with $\mathbb{S}^{1}$. We have $C\mathbb{S}^{1}=D^{2}$ . The cone of a point
$x$, $Cx=I$ . None of these cones admit a manifold structure (without
boundary!). Thus, one can naturally ask if $CM$ is ever a manifold when $M$ is
a manifold. The intuition of the proof technique is to get the discussion
closer to a discussion of local homology groups. This is because, as is known
(see [1] and the sixth example in the preceding section), that the dimension
of a manifold is characterized by its local homology groups
$H_{n}(X,X\backslash x)$ for each $x\in X$. Thus, to get the closer discussion
to local homology groups, one requires to go through a degree of technical
details. Here is another motivating example

######  Theorem 2.2.

(Motivating) If $M$ is a compact 2-manifold then $CM$ is not a 3-manifold.

###### Proof.

Suppose for a contradiction that $CM$ is a 3 manifold. Since $M$ is compact,
so is its cone. Since the cone of any space is contractible, it must be simply
connected. Therefore by the topological Poincare Conjecture,
$CM\approx\mathbb{S}^{3}$. However, $H_{3}(\mathbb{S}^{3})=\mathbb{Z}$,
whereas since cones are homotopy equivalent to a point, we have $H_{3}(CM)=0$.

∎

Squeezing more out of the proof above, we obtain the following example:

######  Example 2.1.

$C\mathbb{R}P^{2}$ is not a topological manifold of dimension 3.

######  Theorem 2.3.

If $M$ is a manifold then $CM$ is not a topological manifold.

###### Proof.

We recall that $CM$ is obtained by forming a closed cylinder of finite height,
via attaching a copy of $M$ to each $t$, where $t\in[0,1]$ then identifying
the points of the in the copy of $M$ at $t=1$; obtaining a ”vertex” point. Let
$\star$ denote the vertex point. If $CM$ is a manifold, then it must be a
manifold of dimension $n+1$. Hence, $CM\backslash\\{\star\\}$ is also a
manifold of dimension $n+1$. Now note that
$CM\backslash\\{\star\\}=M\times[0,1)$. Since $0$ is a deformation retract of
$[0,1)$, $\tilde{H}_{m}(CM\backslash\\{\star\\})\cong\tilde{H}_{m}(M)$ for
each $m$. Now note that if $p\in M\times 0$ , then $p=(x,0)$. Thus,

$\tilde{H}_{m}(CM\backslash\\{\star\\},CM\backslash\\{\star,p\\})\cong\tilde{H}_{m}(M\times[0,1),(M\times[0,1))\backslash\\{p\\})$

To establish the result we seek, it suffices to establish the following two
claims.

Claim 1:
$\tilde{H}_{m}(M\times[0,1),(M\times[0,1))\backslash\\{p\\})\cong\tilde{H_{m}}(M\times
0,(M\times 0)\backslash\\{p\\})$

To prove the claim it suffices to show that the pair
$(M\times[0,1),(M\times[0,1)\backslash\\{p\\}))$ is homotopy equivalent to
$(M\times 0,(M\times 0)\backslash\\{p\\})$. To see this, first note that $0$
is a deformation retract of $[0,1)$ via the map

$H:[0,1)\times I\rightarrow[0,1)$, given by $H(s,t)=(1-t)s$

Now $H$ induces a continuous map $H^{M}:(M\times[0,1))\times I\rightarrow
M\times[0,1)$ given by:

$H^{M}((m,s),t)=(m,H(s,t))$

It is clear that $H^{M}$ deformation retracts $M\times[0,1)$ onto $M\times 0$.
Now, the retraction associated with the homotopy above is the map

$r:M\times[0,1)\rightarrow M\times 0$ , given by $r(x,s)=(x,0)$

Thus, $r((M\times[0,1))\backslash\\{p\\})=(M\times 0)\backslash\\{x\\}$ . As a
consequence, we get that $H$ induces a homotopy equivalence
$(M\times[0,1),M\times[0,1)\backslash\\{p\\})\simeq(M\times 0,(M\times
0)\backslash\\{p\\})$.

Claim 2: $\tilde{H_{m}}(M\times 0,(M\times
0)\backslash\\{p\\})\cong\tilde{H}_{m}(M,M\backslash\\{x\\})$

This follows because we have a canonical homeomorphism $M\rightarrow
M\times\\{0\\}$ which induces a homeomorphism of pairs
$(M,M\backslash\\{x\\})\approx(M\times 0,(M\times 0)\backslash\\{p\\})$.

However, since the top local homology groups must differ for dimensional
reasons, we obtain a contradiction.

∎

##  3 Conclusion

In this paper, we first provided a proof of showing that wedging a manifold
with itself (of dimension atleast 2) is not a manifold. Then, we have written
a proof to show that the cone of a manifold is not a topological manifold The
result can be extended in various directions. For instance, $C\mathbb{S}^{1}$
is a topological manifold with boundary. Thus, one can investigate the
conditions conditions for which the cone of a manifold (as a manifold with
boundary) admits a smooth structure. Moreover, one can seek to find an
analogous result to the results presented in this paper using cohomology.

##  4 References

[1] Pi.math.cornell.edu. 2022. Allen Hatcher’s Homepage. [online] Available
at: ¡https://pi.math.cornell.edu/ hatcher/¿

[2] Lee, J., 2011. Introduction to topological manifolds. New York: Springer.

[3] Rotman, J., 2013. An Introduction to Algebraic Topology. New York, NY:
Springer.

[4] Munkres, J., 2013. Topology. Harlow: Pearson.