# A Self-Adaptive Algorithm of the Clean Numerical Simulation (CNS) for Chaos

Shijie Qin  Shijun Liao  State Key Laboratory of Ocean Engineering, Shanghai
200240, China  Center of Marine Numerical Experiment, School of Naval
Architecture, Ocean and Civil Engineering, Shanghai Jiao Tong University,
Shanghai 200240, China  School of Physics and Astronomy, Shanghai Jiao Tong
University, Shanghai 200240, China

###### Abstract

The background numerical noise $\varepsilon_{0}$ is determined by the maximum
of truncation error and round-off error. For a chaotic system, the numerical
error $\varepsilon(t)$ grows exponentially, say,
$\varepsilon(t)=\varepsilon_{0}\exp(\kappa\,t)$, where $\kappa>0$ is the so-
called noise-growing exponent. This is the reason why one can not gain a
convergent simulation of chaotic systems in a long enough interval of time by
means of traditional algorithms in double precision, since the background
numerical noise $\varepsilon_{0}$ might stop decreasing because of the use of
double precision. This restriction can be overcome by means of the clean
numerical simulation (CNS), which can decrease the background numerical noise
$\varepsilon_{0}$ to any required tiny level. A lot of successful applications
show the novelty and validity of the CNS. In this paper, we further propose
some strategies to greatly increase the computational efficiency of the CNS
algorithms for chaotic dynamical systems. It is highly suggested to keep a
balance between truncation error and round-off error and besides to
progressively enlarge the background numerical noise $\varepsilon_{0}$, since
the exponentially increasing numerical noise $\varepsilon(t)$ is much larger
than it. Some examples are given to illustrate the validity of our strategies
for the CNS.

###### keywords:

Chaos; Clean Numerical Simulation (CNS); self-adaptive algorithm;
computational efficiency.

##  1 Introduction

For a chaotic dynamical system, the sensitivity dependence on initial
conditions (SDIC) was first discovered by Poincaré [1], and then this
phenomenon was discovered once again by Lorenz [2] with a more familiar name
“butterfly-effect”. Due to the SDIC, a very weak, small-scale disturbance of
the initial condition will give rise to a huge deviation of numerical solution
of the chaotic system after a long enough temporal interval [3, 4, 5].
Furthermore, it was found that a chaotic dynamical system not only has the
sensitivity dependence on initial conditions (SDIC) but also possesses the
sensitivity dependence on numerical algorithms (SDNA), as reported by Lorenz
[6, 7]. All of these phenomena are due to the exponential increase of noise
(or uncertainty) of chaotic systems, but unfortunately artificial numerical
noises (i.e. truncation errors and round-off errors) are always inevitable for
almost all of the numerical algorithms. Thus, for a chaotic dynamical system,
calculated trajectories of computer-generated simulations obtained by means of
different numerical algorithms (with single/double precision) and different
time steps are mostly quite different. Naturally, such kind of non-
replicability/unreliability of chaotic solution has brought plenty of heated
debates on the credence of the numerical simulation of chaotic dynamical
system [8], and someone even made an extremely pessimistic conclusion that
“for chaotic systems, numerical convergence cannot be guaranteed forever” [9].
In addition, it has been recently reported that “shadowing solutions can be
almost surely nonphysical”, which “invalidates the argument that small
perturbations in a chaotic system can only have a small impact on its
statistical behavior” [10].

To gain a reproducible/reliable numerical simulation of chaotic systems, Liao
[11] proposed a brand-new numerical strategy, namely the “Clean Numerical
Simulation” (CNS) [12, 13, 14], to control the background numerical noise,
say, truncation error and round-off error, during a temporal interval
$t\in[0,T_{c}]$, where $T_{c}$ is the so-called “critical predictable time”
and this temporal interval should be long enough for calculating statistics.
In the frame of the CNS [11, 12, 13, 14, 15, 16, 17, 18], the temporal
truncation error and the spatial truncation error are able to be decreased to
a required small level via using the Taylor expansion method with a high
enough order in the temporal dimension and adopting a fine enough
discretization method in the spatial dimension (such as the high-order spatial
Fourier expansion), respectively. Significantly, all of the physical and
numerical variables/parameters should be represented by means of the multiple
precision (MP) [19] with a large enough number of significant digits, and thus
the round-off error is also able to be decreased to a required small level.
Moreover, an additional numerical simulation of the identical chaotic system
with the even smaller numerical noise is required and performed in order to
determine such a “critical predictable time” $T_{c}$, so that the numerical
noise could be negligible and thus the computer-generated solution of a
chaotic system is reproducible/reliable within the whole spatial computational
domain and in the temporal interval $[0,T_{c}]$. In this way, different from
some other general numerical algorithms, the CNS is able to give the
reproducible/reliable numerical simulation of a chaotic dynamical system
within a finite but long enough temporal interval.

Here it should be emphasized that although our CNS strategy is based on the
classical Taylor series method [20] as well as the multiple precision [19],
the scientific significance of this strategy is mainly about the “critical
predictable time” $T_{c}$: the CNS can greatly reduce the background numerical
noise, i.e. truncation error and round-off error, to any a required tiny level
so that the numerical noise is negligible compared with the “true” physical
solution, and thus the corresponding numerical result of a chaotic system is
reproducible/reliable in an interval of time $[0,T_{c}]$ that is long enough
for statistics, as described in the next section. In other words, the results
of chaotic dynamical systems given by the CNS can be regarded as a “clean”
benchmark solution [15, 16, 17, 18], which is the main purpose of proposing
this CNS strategy. By contrast, solely adopting the Taylor series method [20,
21, 22, 23] to solve a chaotic system for high precision, one usually does not
focus on the “critical predictable time” $T_{c}$ and thus obtain a mixture of
the “true” physical solution and the “false” numerical noise, which are mostly
at the same order of magnitude, since the background numerical noise of a
simulation of chaos should increase exponentially (and quickly) until to the
same level of “true” physical solution, which is not considered by traditional
numerical strategies.

For the computer-generated simulation of a chaotic dynamical system given by a
certain numerical algorithm, it is well-known that the averaged level of
numerical noise should increase exponentially within a temporal interval
$[0,T_{c}]$, i.e.

$\varepsilon(t)=\varepsilon_{0}\exp(\kappa\,t),\hskip
28.45274ptt\in[0,T_{c}],$ (1)

where the noise-growing exponent $\kappa$ is a positive constant (usually
corresponding to the maximum Lyapunov exponent for a chaotic dynamical system
with a finite degree of freedom), $T_{c}$ is the above-mentioned critical
predictable time in the frame of the CNS strategy, $\varepsilon_{0}$
represents the level of initial/background numerical noise (which is
determined by the initial truncation error and the initial round-off error),
and $\varepsilon(t)$ denotes the averaged level of evolving noise for a
computer-generated simulation, respectively. Considering that there might be
another increasing pattern of numerical noise [24, 25] that is more
meticulous, the exponential growing (1) is still suitable for a long-time
simulation. Theoretically, the critical predictable time $T_{c}$ is determined
by a given value of the critical numerical noise $\varepsilon_{c}$, i.e.
$\varepsilon_{c}=\varepsilon_{0}\exp(\kappa\,T_{c})$ that leads to

$T_{c}=\frac{1}{\kappa}\ln\left(\frac{\varepsilon_{c}}{\varepsilon_{0}}\right).$
(2)

Obviously, if the value of $\varepsilon_{c}$ is unchanged, the smaller the
level of the initial/background numerical noise $\varepsilon_{0}$, the larger
the critical predictable time $T_{c}$.

Unfortunately, it is impossible in practice to obtain the evolving noise
$\varepsilon(t)$ with high accuracy, because we do not know the true
(physical) solution of a numerical simulation of chaos. Thus, a practical
approach with satisfied numerical precision is required to calculate the
$\varepsilon(t)$. Let $\mathbf{x}\in\Omega$ represent the dimensional vector
in a chaotic dynamical system, $\phi(\mathbf{x},t)$ denote the solution of
numerical (computer-generated) simulation that is reproducible/convergent
within $t\in[0,T_{c}]$ possessing the initial/background numerical noise
$\varepsilon_{0}$, and $\phi^{\prime}(\mathbf{x},t)$ denote another solution
(using the identical initial/boundary conditions and physical parameters) that
is reliable within $t\in[0,T_{c}^{\prime}]$ possessing the initial/background
numerical noise $\varepsilon_{0}^{\prime}$ which is smaller than
$\varepsilon_{0}$. Due to the exponentially growing property (1) of numerical
noise for a chaotic dynamical system, there is $T_{c}^{\prime}>T_{c}$ and that
$\phi^{\prime}(\mathbf{x},t)$ within $t\in[0,T_{c}]$ must be superior and much
closer to the physical solution (true solution) compared with
$\phi(\mathbf{x},t)$. And thus, $\phi^{\prime}(\mathbf{x},t)$ could be seen as
a benchmark solution to help us determine the numerical noise of
$\phi(\mathbf{x},t)$ within $\mathbf{x}\in\Omega$ approximately. Therefore,
practically, the evolving noise $\varepsilon(t)$ is obtained via comparing
$\phi(\mathbf{x},t)$ (possessing the background numerical noise
$\varepsilon_{0}$) with a superior numerical solution
$\phi^{\prime}(\mathbf{x},t)$ (possessing the smaller background numerical
noise $\varepsilon_{0}^{\prime}$).

Up to now, the above-mentioned CNS strategy has been applied to many chaotic
dynamical systems successfully with the corresponding computer-generated
simulations being reproducible and of course reliable. For example, via using
some general numerical algorithms with double precision, one always obtains
the reproducible numerical solutions of the well-known Lorenz system only in a
short temporal interval, i.e. $t\in[0,32]$ approximately [11]. By contrast,
via using the CNS strategy, a reproducible/convergent numerical simulation of
the same chaotic Lorenz system within quite a long temporal interval, i.e.
$t\in[0,10000]$, was obtained for the first time by Liao and Wang [14].
Besides, Liao and Li [26] studied the evolution of the microscopic physical
uncertainty of initial condition for the famous three-body system (which is
chaotic) by means of the CNS, and they found that the uncertainty finally
becomes macroscopical, which leads to the random escape of the three-body
system as well as the behavior of symmetry breaking. It indicates that the
uncertainty of microscopic physics could be the origin of large-scale
randomness for the well-known three-body system. Furthermore, the numerical
noise of the CNS strategy can be controlled to be much smaller than the
uncertainty of microscopic physics, and thus via using the CNS, Lin et al.
[27] theoretically provided rigorous evidence to demonstrate that the
microscopic thermal fluctuation should be the origin of large-scale randomness
of the two-dimensional turbulent Rayleigh-Bénard convection. Significantly,
with the help of China’s national supercomputer, the CNS strategy was applied
to investigate the periodic orbits of the famous three-body problem, and more
than $2000$ brand-new families of periodic orbits were discovered successfully
by Li et al. [28, 29, 30]. Those newly found periodic orbits were reported
twice in the popular magazine New Scientist [31, 32], because, for the three-
body problem, there are only three families of chaotic periodic orbits that
had ever been found since Newton mentioned this famous problem three hundred
years ago! It is also worth noting that, according to a known periodic orbit
as well as three equal masses, and integrating the governing equations by
means of the CNS, Li et al. [33] obtained $135445$ brand-new periodic orbits
of arbitrarily unequal masses of the three-body system, including $13315$
stable ones. In addition, using the CNS in quite a long temporal interval, Xu
et al. [17] obtained the reliable/reproducible trajectories of a free-fall
disk that is chaotic under some certain physical parameters, and the CNS
strategy is able to help him accurately forecast the position and posture of
the chaotic free-fall disk near the bifurcation point.

As for spatiotemporal chaos, Hu & Liao [15] and Qin & Liao [16] proposed an
efficient CNS strategy utilized in physical space to numerically solve the 1D
complex Ginzburg-Landau equation (CGLE) and the damped driven sine-Gordon
equation (SGE), respectively, which further demonstrates the effectiveness of
the CNS strategy that can exactly maintain both the statistical properties and
symmetric features of the spatiotemporal chaotic systems in which general
numerical algorithms with double precision always fail. Recently, taking the
CNS strategy as a tool, it has been found that the statistical features (such
as the probability density function) of some chaotic dynamical systems are
extremely sensitive to the tiny noise/disturbance, and thus this kind of chaos
is called ultra-chaos by Liao & Qin [18]. As a brand-new concept, the ultra-
chaos might deepen and enrich our understandings about chaos and turbulence.
Furthermore, with the help of CNS, Qin & Liao [34] provide rigorous evidence
that numerical noises as a kind of tiny artificial stochastic disturbances
have quantitatively and qualitatively large-scale influences on a sustained
turbulence. In a word, the above-mentioned investigations demonstrate the
effectiveness and potential of the CNS for complex chaotic dynamical systems.

Although the CNS is able to be applied to obtain the reproducible/convergent
numerical simulation of a chaotic dynamical system within a long enough
temporal interval, it is more time-consuming compared with some other general
numerical algorithms with double precision [14, 27]. In this paper, according
to the exponentially growing property of noise $\varepsilon(t)$ in (1), we
propose a modified strategy of the CNS, called the “self-adaptive CNS”, to
significantly increase the computational efficiency of the CNS algorithm. To
illustrate its validity, we apply the CNS with the self-adaptive precision to
some chaotic systems, such as the Lorenz equation, the hyper-chaotic Rössler
system, the three-body problem, and the damped driven sine-Gordon equation.

##  2 Basic ideas of the self-adaptive CNS

In this section, let us use the Lorenz equations [2]

$\left\\{\begin{array}[]{l}\dot{x}(t)=\sigma\,[y(t)-x(t)],\\\
\dot{y}(t)=R\,x(t)-y(t)-x(t)\,z(t),\\\ \dot{z}(t)=x(t)\,y(t)+b\,z(t),\\\
\end{array}\right.$ (3)

in the case of

$\sigma=10,\;\;R=28,\;\;b=-8/3,$ (4)

under the initial condition

$x(0)=-15.8,\;\;y(0)=-17.48,\;\;z(0)=35.64,$ (5)

as one of the most famous chaotic systems (with one positive Lyapunov
exponent) to briefly describe the basic ideas of the self-adaptive CNS.

The CNS algorithm for the Lorenz system (3)-(5) is mainly based on the $M$th-
order Taylor series in the temporal interval $[t,t+\Delta t]$:

$x(t+\Delta t)\approx x(t)+\sum_{m=1}^{M}x^{[m]}(t)\,(\Delta t)^{m},$ (6)
$y(t+\Delta t)\approx y(t)+\sum_{m=1}^{M}y^{[m]}(t)\,(\Delta t)^{m},$ (7)
$z(t+\Delta t)\approx z(t)+\sum_{m=1}^{M}z^{[m]}(t)\,(\Delta t)^{m},$ (8)

where $\Delta t$ is the time step and

$x^{[m]}(t)=\frac{1}{m!}\frac{d^{m}x(t)}{dt^{m}},\;\;y^{[m]}(t)=\frac{1}{m!}\frac{d^{m}y(t)}{dt^{m}},\;\;z^{[m]}(t)=\frac{1}{m!}\frac{d^{m}z(t)}{dt^{m}}$
(9)

are the high-order temporal derivatives. Differentiating both sides of Eqs.
(3) $(m-1)$ times with respect to $t$ and then dividing them by $m!$, we
obtain the iterative formulae

$x^{[m]}(t)=\frac{\sigma}{m}\,\left[y^{[m-1]}(t)-x^{[m-1]}(t)\right],$ (10)
$y^{[m]}(t)=\frac{1}{m}\,\left[R\,x^{[m-1]}(t)-y^{[m-1]}(t)-\sum_{i=0}^{m-1}x^{[i]}(t)\,z^{[m-1-i]}(t)\right],$
(11)
$z^{[m]}(t)=\frac{1}{m}\,\left[\sum_{i=0}^{m-1}x^{[i]}(t)\,y^{[m-1-i]}(t)+b\,z^{[m-1]}(t)\right],$
(12)

for arbitrary $m\geq 1$. Note that parallel technology can be applied to
calculate the sum terms in (6)-(8), (11) and (12).

According to (2), the background numerical noises $\varepsilon_{0}$ must be
small enough if one needs a reliable (reproducible) chaotic solution within a
large temporal interval $t\in[0,T_{c}]$. It is worth noting that the
background numerical noises $\varepsilon_{0}$ in (2) is a constant, which is
determined by the maximum of the spatio-temporal truncation error (resulting
from the truncation of an infinite number of series) and the round-off error
(resulting from a limited number of significant digits of data). The basic
idea of the above-mentioned CNS algorithm is to greatly decrease both the
temporal truncation error and round-off error so that the background numerical
noise $\varepsilon_{0}$ is small enough for a numerical simulation to be
reproducible/reliable within a given temporal interval $t\in[0,T_{c}]$.

Obviously, if the temporal Taylor expansion (6)-(8) is given a large enough
order $M$, the temporal truncation error is able to be decreased under a
required small level. More importantly, different from other traditional
algorithms, we express all of the physical and numerical variables/parameters
by means of the multiple precision (MP) via choosing the significant digits
with a large enough number $N_{s}$, and thus the round-off error is also able
to be decreased under a required small level. In this way, both the temporal
truncation error and round-off error are able to be decreased under a required
small level via the CNS.

Note that computer-generated solutions of chaotic Lorenz system (3)-(5) given
by some general numerical algorithms with double precision are
reproducible/convergent in quite a short temporal interval $t\in[0,32]$. In
$2014$, Liao & Wang [14] obtained a reproducible/convergent numerical
simulation $(x,y,z)$ of the above-mentioned Lorenz system in quite a long
temporal interval $t\in[0,10000]$ (Lorenz unit time) by means of a parallel
algorithm of the CNS using the $3500\hskip 0.56905pt$th-order Taylor expansion
($M=3500$) with the constant time step $\Delta t=0.01$ and $4180$-digit
multiple precision ($N_{s}=4180$) for all physical and numerical
variables/parameters, whose reproducibility/reliability (from the mathematical
viewpoint) was confirmed by means of another simulation
$(x^{\prime},y^{\prime},z^{\prime})$ given by the CNS with the smaller
background numerical noise using the $3600\hskip 0.56905pt$th-order Taylor
expansion ($M=3600$) with the time step $\Delta t=0.01$ and $4515$-digit
multiple precision ($N_{s}=4515$). For simplicity, define the relative error

$\varepsilon(t)=\frac{\big{|}x^{\prime}(t)-x(t)\big{|}+\big{|}y^{\prime}(t)-y(t)\big{|}+\big{|}z^{\prime}(t)-z(t)\big{|}}{\big{|}x^{\prime}(t)\big{|}+\big{|}y^{\prime}(t)\big{|}+\big{|}z^{\prime}(t)\big{|}},$
(13)

where $(x,y,z)$ and $(x^{\prime},y^{\prime},z^{\prime})$ are the two CNS
results mentioned above. It was found that (1) indeed holds with the noise-
growing exponent $\kappa\approx 0.91$ (which corresponds to the maximum
Lyapunov exponent of this Lorenz system), say, $\kappa/\ln 10\approx 0.40$,
indicating that the background numerical noise $\varepsilon_{0}$ will be
enlarged nearly $10^{4000}$ times at $t=10000$. This is the reason why Liao &
Wang [14] had to use the $4180$-digit multiple precision ($N_{s}=4180$) and
the $3500\hskip 0.56905pt$th-order Taylor expansion ($M=3500$) in their CNS
algorithm so as to greatly decrease the background numerical noise
$\varepsilon_{0}$ to a very tiny level!

Frankly speaking, one hardly uses 3500th-order Taylor expansion and data in a
multiple precision with 4180 significant digits in practice. However, from a
theoretical viewpoint, the reproducible/convergent chaotic simulation of the
Lorenz equations in such a long interval of time is very important, since it
gives, for the first time, direct evidence that one can indeed gain a
reproducible/convergent trajectory of chaotic systems in a long enough
interval of time. It invalidates the argument that “for chaotic systems,
numerical convergence cannot be guaranteed forever” [9], although a large
number of calculations are required: it took 220.9 hours (i.e. about 9 days
and 5 hours) using 1200 CPUs of the National Supercomputer TH-1A at Tianjian,
China [14]. This kind of convergent simulation in $t\in[0,10000]$ can be used
as a benchmark solution of the chaotic Lorenz system (3)-(5) to verify the
modified CNS algorithms, as described below.

How to increase the computational efficiency of the CNS?

###  2.1 Keeping a balance between truncation error and round-off error

Note that the background numerical noises $\varepsilon_{0}$ in (1) is
determined by the maximum of the truncation error and round-off error. So, for
solving this problem, it is the optimum that the temporal truncation error is
at the same level as the round-off error. So, we should keep a balance between
the temporal truncation error and the round-off error. Unlike Liao and Wang
[14] who used a constant time step, the variable stepsize (VS) strategy [35]
is able to be applied to the above-mentioned CNS algorithm by means of an
allowed tolerance $tol$ (whose value is given) of the governing equations.
Referring to Barrio et al. [35], the optimal time stepsize is given by

$\Delta t=min\hskip
0.56905pt\left(\frac{tol^{\frac{1}{M}}}{\|x_{i}^{[M-1]}(t)\|_{\infty}^{\frac{1}{M-1}}},\frac{tol^{\frac{1}{M+1}}}{\|x_{i}^{[M]}(t)\|_{\infty}^{\frac{1}{M}}}\right),$
(14)

where $M$ denotes the order of Taylor expansion, $tol$ denotes the allowed
tolerance, $\|~{}\|_{\infty}$ is the infinite norm for the variable $x_{i}$
($i=1,2,3$), and $x_{1}(t)$, $x_{2}(t)$, $x_{3}(t)$ correspond to $x(t)$,
$y(t)$, $z(t)$, respectively. Considering that the parallel technology is
applied to calculate the sum terms in (6)-(8), (11) and (12), here we use the
empirical formula

$M=\left\lceil-1.5\log_{10}(tol)\right\rceil,$ (15)

to determine a proper order of Taylor expansion for the high calculating
efficiency. Furthermore, considering the round-off error of data should be
controlled at the same level of the temporal truncation error, we choose

$tol=10^{-N_{s}},$ (16)

where $N_{s}$ denotes the number of significant digits chosen by means of the
multiple precision. In this way, we can control the background numerical noise
by means of choosing the number $N_{s}$ for multiple precision and keeping a
balance between the temporal truncation error and round-off error via (15) and
(16) with an optimal value of the time step via (14).

The above-mentioned strategy can greatly increase the computational efficiency
of the CNS algorithm. For the chaotic Lorenz system (3)-(5), according to (1),
here $\varepsilon_{0}=10^{-N_{s}}$ and $\kappa\approx 0.91$ (i.e. $\kappa/\ln
10\approx 0.40$), we should choose $N_{s}=4020$ so as to guarantee that the
numerical noise $\varepsilon(t)$ is nearly at the level of $10^{-20}$ at
$t=10000$. In fact, using $N_{s}=4020$ and the corresponding
$tol=10^{-N_{s}}=10^{-4020}$ and
$M=\left\lceil-1.5\log_{10}(tol)\right\rceil=\left\lceil 1.5\hskip
0.85358ptN_{s}\right\rceil=6030$, we obtain a convergent simulation by means
of a parallel CNS algorithm using 1200 CPUs of the National Supercomputer TH-2
at Guangzhou, China, which agrees in the accuracy of more than 20 significant
digits in the whole interval of time $t\in[0,10000]$ with the benchmark
solution given by Liao & Wang [14]. Note that it took only 96.8 hours (i.e.
about 4 days and 1 hours), just about 44% of the CPU time required by Liao &
Wang [14] in a supercomputer. This illustrates that the computational
efficiency of the CNS algorithm can be indeed greatly increased by using an
optimal time step (14) and keeping a balance between the truncation error and
round-off error, as mentioned above.

###  2.2 Using self-adaptive multiple-precision

The background numerical noise $\varepsilon_{0}$ is determined by the maximum
of the truncation error and round-off error. According to (1), one had to use
very small background numerical noise $\varepsilon_{0}$ so as to gain a
convergent chaotic simulation in a long interval of time. This is indeed true.
For example, to gain the convergent benchmark solution of the chaotic Lorenz
system in $t\in[0,10000]$, Liao & Wang [14] used the $3500\hskip 0.56905pt$th-
order Taylor expansion ($M=3500$) with the time step $\Delta t=0.01$ in the
$4180$-digit multiple precision ($N_{s}=4180$). The corresponding background
numerical noise is indeed rather small. But, unfortunately, it is rather time-
consuming.

Note that a key point of the CNS is to determine the critical predictable time
$T_{c}$. Since there exists a balance between the truncation error and round-
off error, for the chaotic Lorenz system (3)-(5) in the last section, it is
reasonable to assume that the background numerical noise should be equal to
the round-off error, say, $\varepsilon_{0}=10^{-N_{s}}$, where $N_{s}$ is the
initial significant digit number of the multiple-precision (MP). Then,
according to (1) we have

$\varepsilon(t)=\varepsilon_{0}\exp(\kappa\,t)=10^{-(N_{s}-\kappa\,t/\ln
10)},\hskip 28.45274ptt\in[0,T_{c}],$ (17)

where the noise-growing exponent $\kappa\approx 0.91$ is known and is
generally equal to the leading Lyapunov exponent of a temporal chaos, and
further

$\varepsilon_{c}=10^{-N_{s}+\kappa\,T_{c}/\ln 10},$ (18)

which gives the relationship between the initial significant digit number
$N_{s}$ of the multiple-precision (MP) and the critical predictable time
$T_{c}$:

$N_{s}=\left\lceil\frac{\gamma\,\kappa\,T_{c}}{\ln
10}-\log_{10}\varepsilon_{c}\right\rceil,$ (19)

where $\varepsilon_{c}$ denotes the critical numerical noise that is close to
the order of magnitude of the true physical solution, $\lceil\cdot\rceil$
stands for the ceiling function, $\gamma\geq 1$ is a constant used here as a
kind of safety factor, respectively.

According to (17), the numerical noise $\varepsilon(t)$ increases
exponentially. Thus, after a short time such as at $t=t^{\prime}$,
$\varepsilon(t^{\prime})$ becomes much larger than the background numerical
noise $\varepsilon_{0}$. So, it is unnecessary to keep the background
numerical noise $\varepsilon_{0}$ being the same in the whole interval
$t\in[0,T_{c}]$. In theory, according to (1), using a larger background
numerical noise $\varepsilon^{\prime}_{0}$ does not influence the numerical
result in $t\geq t^{\prime}$, as long as
$\varepsilon_{0}<\varepsilon^{\prime}_{0}\leq\varepsilon(t)$. Note that, for
the above-mentioned CNS algorithm of the Lorenz system (3)-(5), the larger
background numerical noise corresponds to the multiple-precision with a
smaller number $N_{s}$ of significant digits and a larger allowed tolerance
$tol=10^{-N_{s}}$ that further leads to a smaller order
$M=\left\lceil-1.5\log_{10}(tol)\right\rceil=\left\lceil 1.5\hskip
0.85358ptN_{s}\right\rceil$ of Taylor expansion.

Thus, according to (19), after integrating a temporal interval
$t\in[0,t^{\prime})$, where $t^{\prime}<T_{c}$, it is sufficient to use a
smaller number of significant digits

$N_{s}=\left\lceil\frac{\gamma\,\kappa\,(T_{c}-t^{\prime})}{\ln
10}-\log_{10}\varepsilon_{c}\right\rceil,$ (20)

and a larger allowed tolerance $tol=10^{-N_{s}}$ that further leads to a
smaller order $M=\left\lceil-1.5\log_{10}(tol)\right\rceil=\left\lceil
1.5\hskip 0.85358ptN_{s}\right\rceil$ of Taylor expansion, respectively, to
gain the CNS result for $t\geq t^{\prime}$. In practice, it is unnecessary to
change $N_{s}$, $tol$ and the corresponding $M$ at each time step but at some
given times such as $t^{\prime}=100$, $500$, $1000$ and so on.

Table 1: CPU times of the self-adaptive CNS algorithm with the adjustable
multiple-precision (MP) and an optimal variable time step for the chaotic
Lorenz system (3)-(5) in $t\in[0,10000]$, i.e. $T_{c}=10000$, where the number
$N_{s}$ of significant digits is determined by (21) with the allowed tolerance
$tol=10^{-N_{s}}$ that further leads to the order
$M=\left\lceil-1.5\log_{10}(tol)\right\rceil=\left\lceil 1.5\hskip
0.85358ptN_{s}\right\rceil$ of Taylor expansion with the optimal time step
(14), taking the safety factor $\gamma=1.1$ and using different values of
$\Delta T$.

$\Delta T$ CPU time (hours) 25 37.7 50 37.2 100 37.4 500 40.3 1000 42.6

Substituting $\kappa=0.91$ into (20) and choosing $\varepsilon_{c}=10^{-2}$,
we have the following relationship

$N_{s}=\left\lceil\frac{\gamma\,(T_{c}-t^{*})}{2.53}+2\right\rceil\approx\left\lceil\frac{\gamma\,(T_{c}-t^{*})}{2.53}\right\rceil,$
(21)

where $t^{*}=n\,\Delta T$ with $n=0,1,2,...$ and $\Delta T$ being a constant
such as $\Delta T=25$, $50$, $100$, $500$, $1000$ and so on. In practice,
there is $T_{c}-t^{*}>500$ for the high enough remaining precision, say, the
value of $N_{s}$ is stopped decreasing when $t>9500$ for the long time
simulation with $t\in[0,10000]$ in this section. Taking the safety factor
$\gamma=1.1$ and using different values of $\Delta T$, the corresponding CPU
times of this self-adaptive CNS algorithm with the adjustable multiple-
precision (MP) and an optimal variable time step, for the chaotic Lorenz
system (3)-(5) in $t\in[0,10000]$, i.e. $T_{c}=10000$, are listed in Table 1.
It indicates that the required CPU time of the above-mentioned self-adaptive
CNS algorithm is not very sensitive to the value of $\Delta T$, and thus we
can choose $\Delta T=0.5\%\hskip 0.85358ptT_{c}=50$ for the relatively higher
computational efficiency.

Table 2: Convergent result of the chaotic Lorenz system (3)-(5) in
$t\in[0,10000]$, i.e. $T_{c}=10000$, given by the self-adaptive CNS algorithm
using $tol=10^{-N_{s}}$ and
$M=\left\lceil-1.5\log_{10}(tol)\right\rceil=\left\lceil 1.5\hskip
0.85358ptN_{s}\right\rceil$ with the optimal time step (14), where $N_{s}$ is
determined by (21) with taking $\gamma=1.1$ and using $\Delta T=100$.

$t$ $x$ $y$ $z$ 1000 13.881997000862393623 19.918303160406394373
26.901943308376105536 2000 -6.8738836932050180481 -1.4848348276698421977
31.349521074674276721 3000 1.6932902170011335241 3.6003418650451083164
21.410875101298497293 4000 -7.6926663606916323997 -13.499590676622338604
14.199428882538458225 5000 -6.0844510954990075032 -10.813737089458431017
12.739116756422288312 6000 0.21673563458354078642 2.1042785739999677006
22.124608735478140521 7000 -11.394859731998057561 -16.575389386215504779
23.681268415272744261 8000 -1.2658734776208739301 -2.3362702560379947755
17.495968339114928401 9000 13.479653230046728502 17.282101858684218362
29.238196888213967777 10000 -15.817277998368267071 -17.366868329556944701
35.558386165882592794

The CNS algorithm described in § 2.1, together with the above-mentioned self-
adaptive strategy, can greatly increase the computational efficiency of the
CNS. For example, for the chaotic Lorenz system (3)-(5) in $t\in[0,10000]$,
using $tol=10^{-N_{s}}$ and
$M=\left\lceil-1.5\log_{10}(tol)\right\rceil=\left\lceil 1.5\hskip
0.85358ptN_{s}\right\rceil$ with the optimal time step (14), where $N_{s}$ is
determined by (21) with taking $\gamma=1.1$, $T_{c}=10000$ and $\Delta T=50$,
we successfully obtain a reproducible/convergent numerical simulation by means
of the parallel CNS algorithm using 1200 CPUs of the National Supercomputer
TH-2 at Guangzhou, China, which agrees with the benchmark solution given by
Liao & Wang [14] in the accuracy of at least 20 significant digits in the
whole interval of time $t\in[0,10000]$, as shown in Table 2. Note that it took
only 37.2 hours (i.e. about 1 day and 13 hours), just about 17% of the CPU
time of the previous CNS algorithm applied by Liao & Wang [14] who likewise
used a supercomputer. Thus, the self-adaptive CNS algorithm mentioned above
has indeed much higher computational efficiency than the previous CNS with the
constant background numerical noise.

---
Figure 1: Evolving noises $\varepsilon(t)$ of the CNS results of the Lorenz
system (3)-(5) in the whole interval of time $t\in[0,10000]$, given by the CNS
algorithm combined with the variable stepsize strategy (marked by CNS-VS, red
solid line) using $N_{s}=4020$, $tol=10^{-N_{s}}=10^{-4020}$,
$M=\left\lceil-1.5\log_{10}(tol)\right\rceil=\left\lceil 1.5\hskip
0.85358ptN_{s}\right\rceil=6030$ with the optimal time step (14), and given by
the self-adaptive CNS algorithm (marked by CNS-SA, blue solid line) using
$tol=10^{-N_{s}}$, $M=\left\lceil-1.5\log_{10}(tol)\right\rceil=\left\lceil
1.5\hskip 0.85358ptN_{s}\right\rceil$ with the optimal time step (14), where
$N_{s}$ is determined by (21) with taking $\gamma=1.1$, $T_{c}=10000$ and
$\Delta T=100$. These evolving noises $\varepsilon(t)$ are obtained via the
comparison with the previous CNS algorithm applied by Liao & Wang [14] that
has the constant stepsize (marked by CNS-CS). Black dashed line:
$\log_{10}(\varepsilon)=0.40\hskip 0.56905ptt-4020$.

In addition, the numerical precision of the above-mentioned self-adaptive CNS
algorithm can be guaranteed: as shown in Fig. 1, the results given by this
strategy (marked by CNS-SA) and the CNS algorithm combined with the variable
stepsize strategy described in § 2.1 (marked by CNS-VS) are in the accuracy of
at least 20 significant digits in the whole interval of time $t\in[0,10000]$,
which is obtained via the comparison with the previous CNS algorithm applied
by Liao & Wang [14] that has the constant stepsize and the constant background
numerical noise (marked by CNS-CS).

##  3 Some examples

###  3.1 Self-adaptive CNS for hyper-chaotic Rössler system

The chaotic four-dimensional Rössler system [36]

$\left\\{\begin{array}[]{l}\dot{x}(t)=-\,y(t)-z(t),\\\
\dot{y}(t)=x(t)+a\,y(t)+w(t),\\\ \dot{z}(t)=b+x(t)\,z(t),\\\
\dot{w}(t)=-\,c\,z(t)+d\,w(t),\\\ \end{array}\right.$

in the case of

$a=0.25,\;\;b=3,\;\;c=0.5,\;\;d=0.05,$ (22)

under the initial condition

$x(0)=-20,\;\;y(0)=z(0)=0,\;\;w(0)=15,$ (23)

has aroused wide concern as a typical hyper-chaotic system [36, 37, 38], since
it has two positive Lyapunov exponents. How can we gain a convergent chaotic
simulation of Rössler system (3.1)-(23) in the accuracy of 20 significant
digits in $t\in[0,10000]$?

The CNS algorithm for the hyper-chaotic Rössler system (3.1)-(23) is also
based on the $M$th-order truncated Taylor series in the temporal interval
$[t,t+\Delta t]$:

$x(t+\Delta t)\approx x(t)+\sum_{m=1}^{M}x^{[m]}(t)\,(\Delta t)^{m},$ (24)
$y(t+\Delta t)\approx y(t)+\sum_{m=1}^{M}y^{[m]}(t)\,(\Delta t)^{m},$ (25)
$z(t+\Delta t)\approx z(t)+\sum_{m=1}^{M}z^{[m]}(t)\,(\Delta t)^{m},$ (26)
$w(t+\Delta t)\approx w(t)+\sum_{m=1}^{M}w^{[m]}(t)\,(\Delta t)^{m},$ (27)

where the high-order derivatives are governed by

$x^{[m]}(t)=\frac{1}{m}\,\left[-\,y^{[m-1]}(t)-z^{[m-1]}(t)\right],$ (28)
$y^{[m]}(t)=\frac{1}{m}\,\left[x^{[m-1]}(t)+a\,y^{[m-1]}(t)+w^{[m-1]}(t)\right],$
(29)
$z^{[m]}(t)=\frac{1}{m}\,\left[B_{m}+\sum_{i=0}^{m-1}x^{[i]}(t)\,z^{[m-1-i]}(t)\right],$
(30)
$w^{[m]}(t)=\frac{1}{m}\,\left[-\,c\,z^{[m-1]}(t)+d\,w^{[m-1]}(t)\right],$
(31)

for arbitrary $m\geq 1$ and

$B_{m}=\left\\{\begin{array}[]{lr}b,\hskip 28.45274ptm=1,\\\ 0,\hskip
28.45274ptm>1.\end{array}\right.$ (32)

Note that the parallel technology can be applied to calculate the sum terms in
(24)-(27) and (30).

Table 3: Convergent result of the hyper-chaotic Rössler system (3.1)-(23) in
$t\in[0,10000]$ given by the CNS parallel algorithm.

$t$ $x$ $y$ $z$ $w$ 1000 -33.992602 -5.5093173 0.087878252 20.503330 2000
-13.578396 9.9097524 0.23536319 13.972979 3000 -76.551626 32.540905
0.039413405 40.392583 4000 -27.968165 -19.878204 0.10477112 24.712276 5000
-21.983158 20.457146 0.14308917 23.787920 6000 -11.968879 20.979201 0.32526220
26.481453 7000 -5.9175355 11.379490 1.0511723 17.005655 8000 -18.606119
-8.2632834 0.15776751 17.709076 9000 -16.668563 13.935457 0.18978718 30.452723
10000 -56.166749 27.803911 0.053903377 29.797559

It is easy for us to know that the maximum Lyapunov exponent of the hyper-
chaotic Rössler system (3.1)-(23) is about 0.11, which gives us the
corresponding noise-growing exponent $\kappa\approx 0.11$ in (1). In this
case, if $\varepsilon_{0}=10^{-N_{s}}$, we have the numerical noise evolution

$\varepsilon(t)\approx\varepsilon_{0}\exp(0.11\hskip 0.85358ptt)\approx
10^{-(N_{s}-0.048\hskip 0.85358ptt)}.$

If our CNS simulation should be in the accuracy of at least 8 significant
digits in the whole interval of $t\in[0,10000]$, we have $T_{c}=10000$ and

$-(N_{s}-0.048\hskip 0.85358ptT_{c})\leq-8,$

which gives $N_{s}\geq 488$, indicating that we should choose $N_{s}=488$.

Similarly, in the frame of the CNS, the background numerical noise (i.e.
truncation error and round-off error) of this system can be decreased under a
required tiny level by means of choosing a large enough order $M$ of the
Taylor expansion (24)-(27) and a large enough number $N_{s}$ of significant
digits for multiple-precision. First, following Liao & Wang [14] who used a
constant time step, we obtain a reproducible/convergent simulation of the
hyper-chaotic Rössler system (3.1)-(23) in $t\in[0,10000]$ by means of a
parallel CNS algorithm using the 415th-order Taylor expansion ($M=415$) with a
fixed time step $\Delta t=0.01$ in the multiple precision of 488 significant
digits ($N_{s}=488$), as listed in Table 3. In fact, this convergent
simulation result agrees in the accuracy of more than 8 significant digits in
the whole interval of time $t\in[0,10000]$ compared with the benchmark
solution given by another CNS using the 500th-order Taylor expansion ($M=500$)
with a fixed time step $\Delta t=0.01$ in the multiple precision of 550
significant digits ($N_{s}=550$). It takes 5804 seconds (i.e. about 1 hours
and 37 minutes) using 50 Intel’s CPUs (Xeon Silver 4114) on our local cluster.

Then, we apply the strategy of keeping a balance between truncation error and
round-off error (mentioned in § 2.1) to increase the computational efficiency
of the CNS algorithm. The variable stepsize (VS) scheme is applied with an
optimal time step determined by (14), where $i=1,2,3,4$ is for this hyper-
chaotic Rössler system and thus $x_{1}(t)$, $x_{2}(t)$, $x_{3}(t)$, $x_{4}(t)$
correspond to $x(t)$, $y(t)$, $z(t)$, $w(t)$, respectively. Considering that
the parallel technology is applied to calculate the sum terms in (24)-(27) and
(30), here we adopt the empirical formula

$M=\left\lceil-1.5\log_{10}(tol)\right\rceil,$ (33)

to choose a proper order of Taylor expansion for the high calculating
efficiency. Besides, (16) is used to keep the balance between the round-off
error and the truncation error. In this way, we can control the background
numerical noise $\varepsilon_{0}$ by choosing the number $N_{s}$ of
significant digits for multiple-precision, say, $\varepsilon_{0}$ is at the
level of $10^{-N_{s}}$.

By means of the CNS algorithm described in § 2.1 using a fixed value of
$N_{s}=488$ for multiple-precision, $tol=10^{-488}$ for the allowed tolerance,
$M=\left\lceil-1.5\log_{10}(tol)\right\rceil=\left\lceil 1.5\hskip
0.85358ptN_{s}\right\rceil=732$ for the order of Taylor expansion, and an
optimal time step given by (14), we obtain a reproducible/convergent numerical
simulation of the hyper-chaotic Rössler system (3.1)-(23) in $t\in[0,10000]$,
which gives exactly the same result as those listed in Table 3. And it takes
608 seconds (i.e. about 10 minutes) using 50 Intel’s CPUs (Xeon Silver 4114)
on our local cluster, say, only 10.5% CPU time of the previous CNS algorithm
(i.e. 5804 seconds) with a fixed time step. The convergence of this CNS result
is confirmed by comparing it with another CNS result with the even smaller
background numerical noise, given by a fixed value of $N_{s}=550$ for
multiple-precision, $tol=10^{-550}$ for the allowed tolerance,
$M=\left\lceil-1.5\log_{10}(tol)\right\rceil=\left\lceil 1.5\hskip
0.85358ptN_{s}\right\rceil=825$ for the order of Taylor expansion, and the
optimal time step via (14).

---
Figure 2: CPU times of the self-adaptive CNS algorithm for the hyper-chaotic
Rössler system (3.1)-(23) in $t\in[0,10000]$, given by different values of
$\Delta T$, using $tol=10^{-N_{s}}$ and
$M=\left\lceil-1.5\log_{10}(tol)\right\rceil=\left\lceil 1.5\hskip
0.85358ptN_{s}\right\rceil$ with the optimal time step (14), where $N_{s}$ is
determined by (34) with taking $\gamma=1.1$ and $T_{c}=10000$. Red circle:
computed results; Black dashed line: CPU time (seconds) equals to $0.04\hskip
0.85358pt\Delta T+248$.

In addition, the computational efficiency can be further increased by means of
the self-adaptive CNS algorithm described in § 2.2. Substituting $\kappa=0.11$
into (20) and choosing $\varepsilon_{c}=10^{-2}$, we have the following
relationship

$N_{s}=\left\lceil\frac{\gamma\,(T_{c}-t^{*})}{20.9}+2\right\rceil\approx\left\lceil\frac{\gamma\,(T_{c}-t^{*})}{20.9}\right\rceil,$
(34)

where $t^{*}=n\,\Delta T$ with the non-negative integer $n$. In practice,
there is $T_{c}-t^{*}>1000$ for the high enough remaining precision, say, the
value of $N_{s}$ is stopped decreasing when $t>9000$ for the long time
simulation with $t\in[0,10000]$ in this section. Using different values of
$\Delta T$, the corresponding CPU times of this self-adaptive CNS algorithm
with the adjustable multiple-precision (MP) and an optimal variable time step,
for the hyper-chaotic Rössler system (3.1)-(23) in $t\in[0,10000]$, i.e.
$T_{c}=10000$, are illustrated in Fig. 2. It indicates that there is an
approximate linear relationship that the required CPU time (seconds) equals to
$0.04\hskip 0.85358pt\Delta T+248$. Since the corresponding slope 0.04 is
rather small, it reconfirms the conclusion that the required CPU time of the
above-mentioned self-adaptive CNS algorithm is not very sensitive to the value
of $\Delta T$. Thus, in practice we can choose $\Delta T=0.5\%\hskip
0.85358ptT_{c}=50$ for the relatively higher computational efficiency.

For the hyper-chaotic Rössler system (3.1)-(23) in $t\in[0,10000]$, using
$tol=10^{-N_{s}}$ and $M=\left\lceil-1.5\log_{10}(tol)\right\rceil=\left\lceil
1.5\hskip 0.85358ptN_{s}\right\rceil$ (according to (16) and (33),
respectively) with the optimal time step (14), where $N_{s}$ is determined by
(34) with taking $\gamma=1.1$, $T_{c}=10000$ and $\Delta T=50$, we
successfully obtain the same reproducible/convergent numerical simulation by
means of the parallel CNS algorithm together with the above-mentioned self-
adaptive strategy using 50 Intel’s CPUs (Xeon Silver 4114) on our local
cluster, as listed in Table 3, which agrees with the benchmark solution given
by another CNS, using the 500th-order Taylor expansion ($M=500$) with a fixed
time step $\Delta t=0.01$ in the multiple precision of 550 significant digits
($N_{s}=550$) mentioned above, in the accuracy of at least 8 significant
digits in the whole interval of time $t\in[0,10000]$. Especially, it takes
only 250 seconds (i.e. about 4 minutes), say, only 4.3% of the CPU time (i.e.
5804 seconds) of the previous CNS algorithm with a fixed time step and a fixed
value of $N_{s}$ for multiple-precision. This further verifies the high
computational efficiency of the self-adaptive CNS algorithm mentioned in § 2.

###  3.2 Self-adaptive CNS for three-body problem

Here let us consider the well-known three-body problem [1, 39, 40, 41], i.e.
the motion of three celestial objects/bodies under their mutual gravitation.
Let $x_{1}$, $x_{2}$, $x_{3}$ denote three Cartesian coordinates and
$\mathbf{r}_{i}=(x_{1,i},x_{2,i},x_{2,i})$ denotes the corresponding position
vector of the $i$th body. Considering Newton’s Law of Gravitation, the motion
of three bodies is governed by the following non-dimensional equation

$\ddot{x}_{k,i}=\sum_{j=1,j\neq
i}^{3}\rho_{j}\frac{x_{k,j}-x_{k,i}}{R^{3}_{i,j}},\hskip 28.45274ptk=1,2,3,$
(35)

where

$R_{i,j}=\left[\sum_{k=1}^{3}(x_{k,j}-x_{k,i})^{2}\right]^{\frac{1}{2}}$ (36)

and

$\rho_{i}=\frac{m_{i}}{m_{1}},\hskip 28.45274pti=1,2,3$ (37)

denotes the ratio of mass, in which $m_{i}$ denotes the mass of the $i$th
body.

Similarly, in the frame of the CNS, the background numerical noise (i.e.
truncation error and round-off error) of solving the three-body problem (35)
can be decreased under a required tiny level by means of choosing a large
enough order $M$ of the Taylor expansion and a large enough number $N_{s}$ of
significant digits for multiple-precision. For more details, please refer to
Liao [42].

Without loss of generality, in this paper we follow Liao [42] to consider the
motion of three bodies with the initial positions

$\mathbf{r}_{1}=(0,0,-1)+d\mathbf{r}_{1},\hskip
28.45274pt\mathbf{r}_{2}=(0,0,0),\hskip
28.45274pt\mathbf{r}_{3}=-(\mathbf{r}_{1}+\mathbf{r}_{2}),$ (38)

as well as the initial velocities

$\dot{\mathbf{r}}_{1}=(0,-1,0),\hskip
28.45274pt\dot{\mathbf{r}}_{2}=(1,1,0),\hskip
28.45274pt\dot{\mathbf{r}}_{3}=-(\dot{\mathbf{r}}_{1}+\dot{\mathbf{r}}_{2}),$
(39)

where $d\mathbf{r}_{1}=\delta\,(1,0,0)$ denotes the micro-level physical
uncertainty with $\delta=10^{-60}$. For simplicity, we consider the case of
equal masses, say, $\rho_{j}=1$ with $j=1,2,3$.

It is easy for us to know that the maximum Lyapunov exponent of the above-
mentioned three-body problem is about 0.168, which gives us the corresponding
noise-growing exponent $\kappa\approx 0.168$ in (1). In this case, if
$\varepsilon_{0}=10^{-N_{s}}$, we have the numerical noise evolution

$\varepsilon(t)\approx\varepsilon_{0}\exp(0.168\hskip 0.85358ptt)\approx
10^{-(N_{s}-0.073\hskip 0.85358ptt)}.$

If our CNS simulation should be in the accuracy of at least 11 significant
digits in the whole interval of $t\in[0,1000]$, we have $T_{c}=1000$ and

$-(N_{s}-0.073\hskip 0.85358ptT_{c})\leq-11,$

which gives $N_{s}\geq 84$, indicating that we should choose $N_{s}=84$.

First, following Liao [42] who used a constant time step, we obtain a
reproducible/convergent simulation of the three-body problem (35)-(39) in
$t\in[0,1000]$ by means of a CNS algorithm using the 45th-order Taylor
expansion ($M=45$) with a fixed time step $\Delta t=0.01$ in the multiple
precision of 84 significant digits ($N_{s}=84$). In fact, this convergent
simulation result agrees in the accuracy of at least 11 significant digits in
the whole interval of time $t\in[0,1000]$ compared with the benchmark solution
(with the even smaller background numerical noise) given by another CNS using
the 60th-order Taylor expansion ($M=60$) with a fixed time step $\Delta
t=0.01$ in the multiple precision of 100 significant digits ($N_{s}=100$). It
takes 1327 seconds (i.e. about 22 minutes) using Intel’s CPU (Xeon Silver
4114) on our local cluster.

Then, we apply the strategy of keeping a balance between truncation error and
round-off error (mentioned in § 2.1) to increase the computational efficiency
of the CNS algorithm. The variable stepsize (VS) scheme is applied with an
optimal time step determined by (14), where $x_{i}$ is replaced by $x_{k,i}$
for this three-body problem. Considering that (16) is used to keep the balance
between the round-off error and the truncation error, and there is no parallel
technology applied in the CNS algorithm for solving the three-body problem
(35)-(39), here we adopt the optimal order of Taylor expansion [43]

$M=\left\lceil 1.15\hskip 0.85358ptN_{s}+1\right\rceil.$ (40)

In this way, we can control the background numerical noise $\varepsilon_{0}$
by choosing the number $N_{s}$ of significant digits for multiple-precision,
say, $\varepsilon_{0}$ is at the level of $10^{-N_{s}}$.

By means of the CNS algorithm described in § 2.1 using a fixed value of
$N_{s}=84$ for multiple-precision, $tol=10^{-84}$ for the allowed tolerance,
$M=\left\lceil 1.15\hskip 0.85358ptN_{s}+1\right\rceil=98$ for the order of
Taylor expansion, and an optimal time step given by (14), we obtain a
reproducible/convergent numerical simulation of the three-body problem
(35)-(39) in $t\in[0,1000]$, and it takes 370 seconds (i.e. about 6 minutes)
using Intel’s CPU (Xeon Silver 4114) on our local cluster, say, only 28% CPU
time of the previous CNS algorithm (i.e. 1327 seconds) with a fixed time step.
Furthermore, this CNS result is in the accuracy of more than 11 significant
digits in the whole interval of time $t\in[0,1000]$, compared with the above-
mentioned CNS benchmark solution.

---
Figure 3: CPU times of the self-adaptive CNS algorithm for the three-body
problem (35)-(39) in $t\in[0,1000]$, given by different values of $\Delta T$,
using $tol=10^{-N_{s}}$ and $M=\left\lceil 1.15\hskip
0.85358ptN_{s}+1\right\rceil$ with the optimal time step (14), where $N_{s}$
is determined by (41) with taking $T_{c}=1000$. Red circle: computed results;
Black dashed line: CPU time (seconds) equals to $0.16\hskip 0.85358pt\Delta
T+83$.

In addition, the computational efficiency can be further increased by means of
the self-adaptive CNS algorithm described in § 2.2. Substituting
$\kappa=0.168$ into (20) and choosing $\gamma=1.1$, $\varepsilon_{c}=10^{-2}$,
we have the following relationship

$N_{s}=\left\lceil 0.08\hskip 0.85358pt(T_{c}-t^{*})+2\right\rceil,$ (41)

where $t^{*}=n\,\Delta T$ with the non-negative integer $n$. In practice,
there is $T_{c}-t^{*}>200$ for the high enough remaining precision, say, the
value of $N_{s}$ is stopped decreasing when $t>800$ for the long time
simulation with $t\in[0,1000]$ in this section. Using different values of
$\Delta T$, the corresponding CPU times of this self-adaptive CNS algorithm
with the adjustable multiple-precision (MP) and an optimal variable time step,
for the three-body problem (35)-(39) in $t\in[0,1000]$, i.e. $T_{c}=1000$, are
illustrated in Fig. 3. It indicates that there is an approximate linear
relationship that the required CPU time (seconds) equals to $0.16\hskip
0.85358pt\Delta T+83$. Since the slope 0.16 of this linear relationship is
rather small, it once again confirms the conclusion that the required CPU time
of the above-mentioned self-adaptive CNS algorithm is not very sensitive to
the value of $\Delta T$. Thus, in practice we can choose $\Delta T=0.5\%\hskip
0.85358ptT_{c}=5$ for the relatively higher computational efficiency.

For the three-body problem (35)-(39) in $t\in[0,1000]$, using
$tol=10^{-N_{s}}$ and $M=\left\lceil 1.15\hskip 0.85358ptN_{s}+1\right\rceil$
(according to (16) and (40), respectively) with the optimal time step (14),
where $N_{s}$ is determined by (41) with taking $T_{c}=1000$ and $\Delta T=5$,
we successfully obtain the same reproducible/convergent numerical simulation
by means of the above-mentioned self-adaptive CNS algorithm using Intel’s CPU
(Xeon Silver 4114) on our local cluster, which agrees with the benchmark
solution given by another CNS (with the even smaller background numerical
noise), using the 60th-order Taylor expansion ($M=60$) with a fixed time step
$\Delta t=0.01$ in the multiple precision of 100 significant digits
($N_{s}=100$) mentioned above, in the accuracy of at least 11 significant
digits in the whole interval of time $t\in[0,1000]$. Especially, it takes only
84 seconds (i.e. less than 2 minutes), say, only 6.3% of the CPU time (i.e.
1327 seconds) of the previous CNS algorithm with a fixed time step and a fixed
value of $N_{s}$ for multiple-precision. This result also verifies the high
computational efficiency of the self-adaptive CNS algorithm mentioned in § 2.

###  3.3 Self-adaptive CNS for spatiotemporal chaos

Let us consider here a spatiotemporal chaos, i.e. a chain of pendulums coupled
through the elastic restoring force, governed by the damped driven sine-Gordon
equation [44, 45, 46]:

$u_{tt}(x,t)=u_{xx}(x,t)-sin[u(x,t)]-\alpha\hskip
0.56905ptu_{t}(x,t)+\Gamma\hskip 0.56905ptsin(\omega\hskip
0.56905ptt-\lambda\hskip 0.56905ptx),$ (42)

subject to a periodic boundary condition

$u(x+l,t)=u(x,t),$ (43)

where the subscript represents the spatial/temporal partial derivative, $x$
and $t$ denote the variables in the spatial and temporal dimensions, $u(x,t)$
represents the angle of a pendulum, $\alpha$ denotes a constant related to the
damped friction, $\Gamma$ denotes a constant related to the external force
field, $\omega$ is the temporal frequency and $\lambda=2\pi/l$ is the spatial
frequency, $l$ denotes the total calculating length of the system,
respectively. Without loss of generality, we follow Chacón et al. [44] to
consider the following case

$\omega=\frac{3}{5},\;\;\alpha=\frac{1}{10},\;\;\Gamma=\frac{461}{500},\;\;l=500,\;\;\lambda=\frac{2\pi}{l}=\frac{\pi}{250},$
(44)

with the initial condition

$u(x,0)=0,\;\;u_{t}(x,0)=0.$ (45)

As reported by Qin & Liao [16], the above-mentioned model corresponds to a
spatiotemporal chaos, whose statistics are extremely sensitive to a small
disturbance: such kind of chaos belongs to the so-called ultra-chaos, which is
in a higher level of disorders than a normal-chaos, as reported by Liao & Qin
[18].

Similarly, the CNS algorithm for the sine-Gordon equation (42)-(45) is also
based on a high enough order of Taylor expansion in the temporal dimension for
decreasing the temporal truncation error under a required tiny level, but
combined with a high enough order of the spatial Fourier expansion for a fine
enough spatial discretization for decreasing the spatial truncation error
under a required tiny level. First, the spatial interval $x\in[0,l\hskip
0.56905pt)$ is discretized uniformly by $N$ equidistant points, say,
$x_{k}=l\,k/N$, where $k=0,1,2,...,N-1$ and $x_{k}$ denotes the $k$th discrete
points in the physical space. Note that the parallel technology is applied in
the spatial discretization. Then, in the temporal dimension, the $M$th-order
Taylor expansion method is used in the temporal interval $[t,t+\Delta t]$,
say,

$u(x_{k},t+\Delta t)\approx u(x_{k},t)+\sum_{m=1}^{M}u^{[m]}(x_{k},t)\,(\Delta
t)^{m},\hskip 28.45274pt0\leq k\leq N,$ (46)

where $\Delta t$ is the time step and

$u^{[m]}(x_{k},t)=\frac{1}{m!}\frac{\partial^{m}u(x_{k},t)}{\partial t^{m}}.$
(47)

From (46), we have the first-order derivative with respect to time:

$u_{t}(x_{k},t+\Delta t)=u^{[1]}(x_{k},t+\Delta t)\approx
u^{[1]}(x_{k},t)+\sum_{m=1}^{M}\,(m+1)\,u^{[m+1]}(x_{k},t)\,\left(\Delta
t\right)^{m}.$ (48)

The high-order temporal derivatives in (46) and (48) can be obtained via
differentiating both sides of Eq. (42) with respect to $t$, while the
corresponding spatial derivatives are approximated by means of the $N$th-order
Fourier spectral expression, i.e.

$u^{[m]}(x,t)\approx\frac{1}{2}\,a_{m,0}(t)+\sum_{n=1}^{\frac{N}{2}-1}\left[a_{m,n}(t)\,cos(\lambda\hskip
0.56905ptnx)+b_{m,n}(t)\,sin(\lambda\hskip
0.56905ptnx)\right]+a_{m,\frac{N}{2}}(t)\,cos\hskip
0.28453pt\left(\frac{\lambda\hskip 0.56905ptNx}{2}\right),$ (49)

where the Fast Fourier Transform (FFT) algorithm as well as parallel
technology can be adopted. For more details, please refer to Qin & Liao [16].
Besides, all physical and numerical variables/parameters are in the multiple
precision with a large enough number $N_{s}$ of significant digits so as to
decrease the round-off error under a required tiny level. As reported by Hu &
Liao [15] and Qin & Liao [16], this kind of parallel CNS algorithm in physical
space for spatio-temporal chaos has much higher computational efficiency than
the previous CNS algorithm in spectrum space [27]. To further increase the
computational efficiency, Qin & Liao [16] applied the VS scheme in the
temporal dimension with a given allowed tolerance $tol$ for solving the
governing equation, using an optimal time step determined by:

$\Delta t=min\hskip
0.56905pt\Bigg{(}\frac{tol^{\frac{1}{M}}}{\|u^{[M-1]}(x_{k},t)\|_{\infty}^{\frac{1}{M-1}}},\frac{tol^{\frac{1}{M+1}}}{\|u^{[M]}(x_{k},t)\|_{\infty}^{\frac{1}{M}}}\Bigg{)},$
(50)

where $\|~{}\|_{\infty}$ is the infinite norm for the variable $x_{k}$. We
adopted the empirical formula $M=\left\lceil-\log_{10}(tol)-10\right\rceil$ to
determine a proper order of Taylor expansion for the high calculating
efficiency [16]. In addition, we use (16) to balance the round-off error at
the same level of the temporal truncation error. In this way, one can control
the background numerical noise $\varepsilon_{0}$ only by means of choosing the
number $N_{s}$ of significant digits for multiple precision.

In this way, Liao & Qin [18] obtained a convergent chaotic simulation $u(x,t)$
of the damped driven sine-Gordon equation (42)-(45) in $t\in[0,3600]$ by means
of a parallel algorithm of the CNS using $N=2^{16}=65536$ and $N_{s}=230$,
corresponding to $tol=10^{-230}$ and $M=220$ according to (16) and (54),
respectively, with the optimal time step via (50). It took 202.6 hours (about
8 days and 11 hours) using 256 Intel’s CPUs (Xeon Silver 4114) on our local
cluster. To confirm its convergence in the whole interval $t\in[0,3600]$, say,
$T_{c}=3600$, Liao & Qin [18] obtained another CNS result $u^{\prime}(x,t)$
with the even smaller background numerical noise using the same $N=65536$ but
$N_{s}=240$, corresponding to $tol=10^{-240}$ and $M=230$ according to (16)
and (54), respectively, with the optimal time step via (50). The deviation of
$u(x,t)$ from $u^{\prime}(x,t)$ is given by

$\varepsilon(t)=\frac{\sum\limits_{n=0}^{\frac{N}{2}}\left|(c_{n}^{\prime})^{2}-(c_{n})^{2}\right|}{\sum\limits_{n=0}^{\frac{N}{2}}|c^{\prime}_{n}|^{2}},$
(51)

where $c_{n}$ and $c^{\prime}_{n}$ are the the complex coefficients of the
spatial Fourier expansion of $u(x,t)$ and $u^{\prime}(x,t)$ at a given time
$t$, respectively. It was found that the deviation evolves in a power law

$\varepsilon(t)\approx\varepsilon_{0}\exp\hskip 0.56905pt(\kappa\,t),$ (52)

where $\varepsilon_{0}=10^{-N_{s}}$ is the background numerical noise and
$\kappa\approx 0.14$ is the noise-growing exponent. For more details, please
refer to Liao & Qin [18].

Similarly, to further increase the computational efficiency, here we adopt the
self-adaptive CNS algorithm described in § 2.2. Substituting $\kappa=0.14$
into (20) and choosing $\varepsilon_{c}=10^{-2}$, we have the following
relationship

$N_{s}=\left\lceil\frac{\gamma\,(T_{c}-t^{*})}{16.4}+2\right\rceil,$ (53)

where $t^{*}=n\,\Delta T$ with the non-negative integer $n$. In practice,
$\Delta T=0.5\%\hskip 0.85358ptT_{c}=18$ is chosen for the higher
computational efficiency and there is $T_{c}-t^{*}>600$ for the high enough
remaining precision, say, the value of $N_{s}$ is stopped decreasing when
$t>3000$ for the long time simulation with $t\in[0,3600]$ in this section.
Considering that (16) is used to keep the balance between the round-off error
and the truncation error, and there is no parallel technology applied in the
CNS algorithm for solving the sine-Gordon equation (42)-(45), here we also
adopt the optimal order of Taylor expansion [43]

$M=\left\lceil 1.15\hskip 0.85358ptN_{s}+1\right\rceil.$ (54)

We adopt the self-adaptive CNS algorithm mentioned above for solving the sine-
Gordon equation (42)-(45), using $N=65536$ for the spatial discretization, the
allowed tolerance $tol=10^{-N_{s}}$ of governing equations via (16), and the
order $M=\left\lceil 1.15\hskip 0.85358ptN_{s}+1\right\rceil$ of Taylor
expansion via (54) with the optimal time step (50), respectively, where the
self-adaptive number of $N_{s}$ for multiple-precision is determined by (53)
with taking $\gamma=1.2$, $T_{c}=3600$ and $\Delta T=18$. Here it should be
emphasized that, according to the definition (51), the power law (52) means an
averaged evolution of deviation, and the real deviations at different discrete
points in the physical space have fluctuations compared with this averaged
deviation. Thus we choose a relatively large safety factor $\gamma=1.2$ for
enough precision.

Note that the same reproducible/convergent numerical result in the averaged
accuracy of 5 significant digits in the whole interval of time $t\in[0,3600]$
(compared with the benchmark solution given by another CNS using the same
$N=65536$ but a fixed value of $N_{s}=240$, corresponding to $tol=10^{-240}$
and $M=230$ according to (16) and (54), respectively) is obtained, by means of
the self-adaptive CNS algorithm mentioned above, which takes 76.9 hours (i.e.
about 3 days and 5 hours) using 256 Intel’s CPUs (Xeon Silver 4114) on our
local cluster, say, only about 38% of the CPU time (i.e. 202.6 hours) of the
previous CNS algorithm. This illustrates that the self-adaptive CNS algorithm
described in § 2 can indeed greatly increase the computational efficiency for
a spatiotemporal chaos.

##  4 Conclusion

The background numerical noise $\varepsilon_{0}$ is determined by the maximum
of truncation error and round-off error. For a chaotic dynamical system, the
numerical error $\varepsilon(t)$ grows exponentially, say,
$\varepsilon(t)=\varepsilon_{0}\exp(\kappa\,t)$, where $\kappa>0$ is the so-
called noise-growing exponent. This is the reason why one can not gain a
convergent simulation of chaotic systems in a long enough interval of time by
means of traditional algorithms in double precision, since the background
numerical noise $\varepsilon_{0}$ might stop decreasing because of the use of
double precision. This restriction can be overcome by means of the clean
numerical simulation (CNS) [11, 12, 13, 14], which can decrease the background
numerical noise $\varepsilon_{0}$ to any a required tiny level. A lot of
successful applications show the novelty and validity of the CNS [15, 16, 17,
18, 26, 27, 28, 29, 30, 33, 34]. In this paper, we propose some strategies to
greatly increase the computational efficiency of the CNS algorithms for
chaotic dynamical systems. It is highly suggested to keep a balance between
truncation error and round-off error and besides to progressively enlarge the
background numerical noise $\varepsilon_{0}$, since the exponentially
increasing numerical noise $\varepsilon(t)$ is much larger than it. To
illustrate the validity of our strategies, we apply the CNS algorithm combined
with the self-adaptive precision to some chaotic dynamical systems, such as
the Lorenz system, the hyper-chaotic Rössler system, the three-body problem,
and a spatiotemporal chaos governed by the damped driven sine-Gordon equation.
All of our results indicate that the self-adaptive CNS algorithm can indeed
greatly increase the computational efficiency for chaotic systems.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or
personal relationships that could have appeared to influence the work reported
in this paper.

## Acknowledgments

This work is partly supported by National Natural Science Foundation of China
(No. 12272230) and Shanghai Pilot Program for Basic Research - Shanghai Jiao
Tong University (No. 21TQ1400202). The parallel algorithms for the Lorenz
system in this paper were performed on TH-2 at National Supercomputer Centre
in Guangzhou, China.

## References

## References

  * [1] H. Poincaré, Sur le problème des trois corps et les équations de la dynamique, Acta Math. 13 (1) (1890) A3–A270.
  * [2] E. N. Lorenz, Deterministic nonperiodic flow, J. Atmos. Sci. 20 (2) (1963) 130–141.
  * [3] J. C. Sprott, J. C. Sprott, Chaos and time-series analysis, Vol. 69, Citeseer, 2003\.
  * [4] Q. Wang, Forward and adjoint sensitivity computation of chaotic dynamical systems, J. Comput. Phys. 235 (2013) 1–13.
  * [5] N. Kuznetsov, G. Leonov, T. Mokaev, A. Prasad, M. Shrimali, Finite-time Lyapunov dimension and hidden attractor of the Rabinovich system, Nonlinear Dyn. 92 (2) (2018) 267–285.
  * [6] E. N. Lorenz, Computational chaos-a prelude to computational instability, Physica D 35 (3) (1989) 299–317.
  * [7] E. N. Lorenz, Computational periodicity as observed in a simple system, Tellus Ser. A-Dyn. Meteorol. Oceanol. 58 (5) (2006) 549–557.
  * [8] T. E. Nazaré, E. G. Nepomuceno, S. A. Martins, D. N. Butusov, A note on the reproducibility of chaos simulation, Entropy 22 (9) (2020) 953.
  * [9] J. Teixeira, C. A. Reynolds, K. Judd, Time step sensitivity of nonlinear atmospheric models: numerical convergence, truncation error growth, and ensemble design, J. Atmos. Sci. 64 (1) (2007) 175–189.
  * [10] N. Chandramoorthy, Q. Wang, On the probability of finding nonphysical solutions through shadowing, J. Comput. Phys. 440 (2021) 110389.
  * [11] S. Liao, On the reliability of computed chaotic solutions of non-linear differential equations, Tellus Ser. A-Dyn. Meteorol. Oceanol. 61 (4) (2009) 550–564.
  * [12] S. Liao, On the numerical simulation of propagation of micro-level inherent uncertainty for chaotic dynamic systems, Chaos Solitons Fractals 47 (2013) 1–12.
  * [13] S. Liao, Physical limit of prediction for chaotic motion of three-body problem, Commun. Nonlinear Sci. Numer. Simul. 19 (3) (2014) 601–616.
  * [14] S. Liao, P. Wang, On the mathematically reliable long-term simulation of chaotic solutions of lorenz equation in the interval [0, 10000], Sci. China-Phys. Mech. Astron. 57 (2) (2014) 330–335.
  * [15] T. Hu, S. Liao, On the risks of using double precision in numerical simulations of spatio-temporal chaos, J. Comput. Phys. 418 (2020) 109629.
  * [16] S. Qin, S. Liao, Influence of numerical noises on computer-generated simulation of spatio-temporal chaos, Chaos Solitons Fractals 136 (2020) 109790.
  * [17] T. Xu, J. Li, Z. Li, S. Liao, Accurate predictions of chaotic motion of a free fall disk, Phys. Fluids 33 (3) (2021) 037111.
  * [18] S. Liao, S. Qin, Ultra-chaos: an insurmountable objective obstacle of reproducibility and replication, Adv. Appl. Math. Mech. 14 (4) (2022) 799–815.
  * [19] P. Oyanarte, MP-a multiple precision package, Comput. Phys. Commun. 59 (2) (1990) 345–358.
  * [20] G. Corliss, Y. F. Chang, Solving ordinary differential equations using taylor series, ACM Trans. Math. Softw. 8 (2) (1982) 114–144.
  * [21] Y. F. Chang, G. Corliss, ATOMFT: solving ODEs and DAEs using Taylor series, Comput. Math. Appl. 28 (10-12) (1994) 209–233.
  * [22] R. Barrio, Performance of the Taylor series method for ODEs/DAEs, Appl. Math. Comput. 163 (2) (2005) 525–545.
  * [23] N. S. Nedialkov, J. D. Pryce, Solving differential-algebraic equations by Taylor series (III): The DAETS code, J. Numer. Anal. Ind. Appl. Math. 3 (1-2) (2008) 61–80.
  * [24] Y. C. Li, Linear hydrodynamic stability, Not. Am. Math. Soc. 65 (10) (2018) 1255–1259.
  * [25] Y. C. Li, R. D. J. G. Ho, A. Berera, Z. C. Feng, Superfast amplification and superfast nonlinear saturation of perturbations as a mechanism of turbulence, J. Fluid Mech. 904 (2020).
  * [26] S. Liao, X. Li, On the inherent self-excited macroscopic randomness of chaotic three-body systems, Int. J. Bifurcation Chaos 25 (09) (2015) 1530023.
  * [27] Z. Lin, L. Wang, S. Liao, On the origin of intrinsic randomness of Rayleigh-Bénard turbulence, Sci. China-Phys. Mech. Astron. 60 (1) (2017) 1–13.
  * [28] X. Li, S. Liao, More than six hundred new families of Newtonian periodic planar collisionless three-body orbits, Sci. China-Phys. Mech. Astron. 60 (12) (2017) 129511.
  * [29] X. Li, Y. Jing, S. Liao, Over a thousand new periodic orbits of a planar three-body system with unequal masses, Publ. Astron. Soc. Jpn. 70 (4) (2018) 64\.
  * [30] X. Li, S. Liao, Collisionless periodic orbits in the free-fall three-body problem, New Astron. 70 (2019) 22–26.
  * [31] L. Crane, Infamous three-body problem has over a thousand new solutions, New Scientist(https://www.newscientist.com/article/2148074-infamous-three-body-problem-has-over-a-thousand-new-solutions/) (2017).
  * [32] C. Whyte, Watch the weird new solutions to the baffling three-body problem, New Scientist(https://www.newscientist.com/article/2170161-watch-the-weird-new-solutions-to-the-baffling-three-body-problem/) (2018).
  * [33] X. Li, X. Li, S. Liao, One family of 13315 stable periodic orbits of non-hierarchical unequal-mass triple systems, Sci. China-Phys. Mech. Astron. 64 (1) (2021) 1–6.
  * [34] S. Qin, S. Liao, Large-scale influence of numerical noises as artificial stochastic disturbances on a sustained turbulence, J. Fluid Mech. 948 (2022) A7.
  * [35] R. Barrio, F. Blesa, M. Lara, VSVO formulation of the Taylor method for the numerical solution of ODEs, Comput. Math. Appl. 50 (1) (2005) 93–111.
  * [36] O. Rossler, An equation for hyperchaos, Phys. Lett. A 71 (2-3) (1979) 155–157.
  * [37] H. Zhang, X. Ma, M. Li, J. Zou, Controlling and tracking hyperchaotic Rössler system via active backstepping design, Chaos Solitons Fractals 26 (2) (2005) 353–361.
  * [38] N. Stankevich, A. Kazakov, S. Gonchenko, Scenarios of hyperchaos occurrence in 4D Rössler system, Chaos 30 (12) (2020) 123129.
  * [39] F. Diacu, P. Holmes, Celestial encounters: the origins of chaos and stability, Princeton university press, 1999.
  * [40] M. Hénon, C. Heiles, The applicability of the third integral of motion: some numerical experiments, Astron. J. 69 (1964) 73.
  * [41] M. Valtonen, H. Karttunen, The three-body problem, Cambridge University Press, 2006\.
  * [42] S. Liao, Physical limit of prediction for chaotic motion of three-body problem, Commun. Nonlinear Sci. Numer. Simul. 19 (3) (2014) 601–616.
  * [43] À. Jorba, M. Zou, A software package for the numerical integration of ODEs by means of high-order Taylor methods, Exp. Math. 14 (1) (2005) 99–117.
  * [44] R. Chacón, A. Bellorín, L. Guerrero, J. González, Spatiotemporal chaos in sine-Gordon systems subjected to wave fields: Onset and suppression, Phys. Rev. E 77 (4) (2008) 046212.
  * [45] J. B. Keller, D. W. McLaughlin, G. C. Papanicolaou, Surveys in Applied Mathematics, Springer, 2013.
  * [46] M. A. Ferré, M. G. Clerc, S. Coulibally, R. G. Rojas, M. Tlidi, Localized structures and spatiotemporal chaos: comparison between the driven damped sine-Gordon and the Lugiato-Lefever model, Eur. Phys. J. D 71 (6) (2017) 1–8.