# The complexity of computing in continuous time: space complexity is
precision

Manon Blanc
<EMAIL_ADDRESS>
École Polytechnique  LIX  Palaiseau  France  Université Paris-Saclay  LISN
Orsay  France  Olivier Bournez
<EMAIL_ADDRESS>
École Polytechnique  LIX  Palaiseau  France

###### Abstract

Models of computations over the integers are equivalent from a computability
and complexity theory point of view by the Church-Turing thesis. It is not
possible to unify discrete-time models over the reals. The situation is
unclear but simpler for continuous-time models, as there is a unifying
mathematical model, provided by ordinary differential equations (ODEs). Each
model corresponds to a particular class of ODEs. For example, the General
Purpose Analog Computer model of Claude Shannon, introduced as a mathematical
model of analogue machines (Differential Analyzers), is known to correspond to
polynomial ODEs. However, the question of a robust complexity theory for such
models and its relations to classical (discrete) computation theory is an old
problem. There was some recent significant progress: it has been proved that
(classical) time complexity corresponds to the length of the involved curves,
i.e. to the length of the solutions of the corresponding polynomial ODEs. The
question of whether there is a simple and robust way to measure space
complexity remains. We argue that space complexity corresponds to precision
and conversely.

Concretely, we propose and prove an algebraic characterisation of
$\mathbf{FPSPACE}$, using _continuous_ ODEs. Recent papers proposed algebraic
characterisations of polynomial-time and -space complexity classes over the
reals, but with a discrete-time: those algebras rely on discrete ODE schemes.
Here, we use classical (continuous) ODEs, with the classic definition of
derivation and hence with the more natural context of continuous-time
associated with ODEs. We characterise both the case of polynomial space
functions over the integers and the reals. This is done by proving two
inclusions. The first is obtained using some original polynomial space method
for solving ODEs. For the other, we prove that Turing machines, with a proper
representation of real numbers, can be simulated by continuous ODEs and not
just discrete ODEs. A major consequence is that the associated space
complexity is provably related to the numerical stability of involved schemas
and the associated required precision. We obtain that a problem can be solved
in polynomial space if and only if it can be simulated by some numerically
stable ODE, using a polynomial precision.

##  1 Introduction

Recently, there has been a renewed interest in models of computations over the
reals and their associated complexity classes. The fact that these models
appear in complexity issues of deep learning models (a.k.a. neural networks)
partially explains it. For example, various problems, such as the training of
fully connected neural networks, have been proved to be a
$\exists\mathbb{R}$-complete problem [2]. Complexity classes like
$\mathbf{FIXP}$ were introduced to discuss the complexity of continuous
functions’ fixed points in various contexts, such as game theory [31]. These
classes and statements are related to discrete-time models of computation over
the reals.

For discrete-time models of computations over the reals, the most famous
approaches are computable analysis, based on the Turing machine model in [58]
and [62] and algebraic models such as the Blum Shub Smale (BSS) model of
computation [8, 7]. The class $\exists\mathbb{R}$ corresponds to the
(constant-free, equivalently uniform) non-deterministic time of the BSS model
of computation. An extensive list of decision problems was proved recently to
be in this class. Both models were tailored for very different applications
and it is well-known we cannot unify existing models with the equivalent of a
Church-Turing thesis. For example, computable functions in a computable
analysis model need to be continuous, while the BSS model intends to consider
functions and problems over the polynomials that are not. It is also explained
by the fact that some models have not been introduced with the idea of
corresponding to actual physical machines but also to discuss abstract
complexity (lower and upper bounds) for associated problems.

Among models of computation over the reals, we can also distinguish
continuous-time models. This includes models of old, first-ever built
computers, such as the Differential Analysers [60]. A famous mathematical
model of such machines is the General Purpose Analog Computer model of Claude
Shannon [54]. It covers many historical machines and today’s analogue devices
[59, 61] too. It also includes various recent approaches and models from deep
learning such as Neural ODEs [23, 44] with many variants. In the context of
continuous-time, the situation is clearer than with discrete-time models, as
there is a unifying way to describe these models, provided by ODEs. Each model
corresponds to a particular class of ODEs. For example, the GPAC corresponds
to polynomial ODEs [36], and Neural ODEs are made by selecting the best
solution among a parameterised class of ODEs: see, e.g. [44].

Even if particular classes of ODES can describe such models, defining a robust
and well-defined computation theory for continuous-time computations is not an
easy problem: see [18] for the most recent survey. In short, the problem with
time complexity is that considering the time variable as a measure of time is
not robust: a curve can always be re-parameterised using a change of variable.
The problem with space complexity is similar: reparameterisation corresponds
to a change of time variable, but also of space-variable, introducing space
and time contractions: See e.g. [18, 16]. Furthermore, many problems for
simple dynamical systems are known to be undecidable, hence forbid
$\mathbf{PSPACE}$-completeness see [37] and [38].

There was a recent breakthrough in [15, 16], where the authors relate time
with the length of the solution curve of an ODE. It holds for polynomial ODEs,
and their projections, known to cover a very wide class of functions,
including all common functions or functions that can be built from them [34].
As the length of a curve is an invariant, this solves the issue of a possible
change of variable. Furthermore, the authors prove that for polynomial ODEs,
this is polynomially related to the time required to solve an ODE, hence
providing a robust notion of time for ODEs. These statements and underlying
constructions, which allow the simulation of Turing machines, led to solving
several open problems: the existence of a universal ODE [17], the proof of the
Turing-completeness of chemical reactions [32], or statements about the
hardness of several dynamical systems problems [38].

The question of whether we can give a simple equivalent defining _space-
complexity_ remains. We argue here that space complexity is polynomially
related and conversely to the numerical stability of ODEs and their associated
precision. We prove that a problem can be solved in polynomial space iff it
can be simulated by some numerically stable ODE, using a polynomial precision.
We prove this holds both for classical complexity over the discrete (functions
over the integers) and also for space complexity for real functions in the
model of computable analysis.

######  NB 1.

In the literature, there are two possible definitions for $\mathbf{FPSPACE}$,
according to whether functions with non-polynomial size values are allowed or
not. In this article, when we talk about $\mathbf{FPSPACE}$, we always assume
the outputs remain of polynomial size. Otherwise, the class is not closed by
composition: the issue is about not counting the output as part of the total
space used. Given $f$ computable in polynomial space and $g$ in logarithmic
space, $f\circ g$ (and $g\circ f$) is computable in polynomial space. But, if
exponential size output is allowed, this is not true: if we assumed only $f$
and $g$ to be computable in polynomial space since the first might give an
output of exponential size.

These questions of providing characterisations of classical complexity using
ODEs can also be seen from the so-called “implicit complexity” point of view.
Having "simple" characterisations of computability and complexity classes is
useful for various fundamental and applied science fields. We are interested
here in "algebraic" characterisations of those classes: we want to define them
as the smallest set $[f_{1},\dots,f_{k};o_{1},\dots,o_{l}]$ where the
$f_{i}$’s are functions, closed under the operators $o_{i}$’s. For example,
the set of computable functions over the integers is well-known to be:

$[0,1,\pi^{i}_{k};\textit{composition, minimisation, primitive recursion}].$

Implicit complexity aims at giving similar algebras for classes of complexity
theory: a reference survey is [24, 25]. The main benefit is to avoid the use
of the framework of Turing machines, which is rather heavy and not necessarily
well-known outside fundamental computer science. Several characterisations for
$\mathbf{PTIME}$ over the integers were proposed. The first is due to Cobham
in [26], but relies on explicit ah-hoc bounds. Other approaches have then been
proposed, see surveys [24, 25]. Recently, Bournez and Durand in [10] suggested
an algebra using the so-called "linear-length" discrete ODEs. Instead of
having explicit bounds, the linearity of the involved discrete ODE guarantees
polynomial time complexity.

Using a similar approach, Blanc and Bournez in [3] and in [4] extended the
constructions to a characterisation of $\mathbf{PTIME}$ for function over the
reals. The latter extended the result to $\mathbf{PSPACE}$, defining _robust_
ODEs as in Definition 4. However, those models rely on discrete ODEs, which
are discrete-time and less natural than continuous ODEs. We review all those
results in Section 3.2.

This paper can be related to [15, 16]: the authors of these articles provided
a characterisation of $\mathbf{PTIME}$ with _continuous_ ODEs, establishing
that time complexity corresponds to the length of the involved curve, i.e. the
motto time complexity = length. Here, we get a motto of the form space
complexity = precision. Some of our constructions have similarities with
statements in [6]. In the later paper, the authors introduce various
robustness concepts and prove they lead to tractability. See the references in
[6] for similar robustness statements. Robustness can also be associated with
a dual motivation: the authors of [39] introduced a concept of robust
undecidability, while here we want a concept of robustness leading to
tractability.

This is not the first time $\mathbf{FPSPACE}$ is characterised using
continuous ODES. However, the existing characterisation [33, 11] is obtained
with complicated conditions on ODEs, while we have a simpler statement,
linking complexity to precision in a simple manner. Notice that the latter
approach dealt with polynomial ODEs, while we do not restrict to polynomial
ODEs. We obtain our statements by revisiting the approach of the latter papers
but working over a compact domain and dealing with error correction more
finely.

Intuitively, this can also be read as being in $\mathbf{PSPACE}$ for an ODE is
consistent with having an attractor easily discretisable when there is one. We
can also define the notion of robustness, as the insensitivity to “small"
perturbations.

While discussing all these issues, we propose an algebraic characterisation of
$\mathbf{PSPACE}$, using _continuous_ ODEs with the algebra ($\mathbb{RCD}$ is
for Continuous Robust Differential). Schema robust ODE is formally defined in
Definition 1:

$\mathbb{RCD}=[0,1,\pi_{i}^{k},+,-,\times,\tanh,\cos,\pi,\frac{x}{2},\frac{x}{3};\textit{composition},\textit{robust
ODE}]$

For a function
$\mathbf{f}\mathrel{\mathop{\ordinarycolon}}\mathbb{R}^{d}\to\mathbb{R}^{d^{\prime}}$
sending every integer ${\mathbf{n}}\in\mathbb{N}^{d}$ to the vicinity of some
integer of $\mathbb{N}^{d}$, say at distance less than $1/4$, we write
$\operatorname{DP}(f)$ for its discrete part: this is the function from
$\mathbb{N}^{d}\to\mathbb{N}^{d^{\prime}}$ mapping
${\mathbf{n}}\in\mathbb{N}^{d}$ to the integer rounding of $\mathbf{f}(n)$.
For a class $\mathcal{C}$ of such functions, we write
$\operatorname{DP}(\mathcal{C})$ for the class of the discrete parts of the
functions of $\mathcal{C}$.

######  Theorem 1.

$\operatorname{DP}(\mathbb{RCD})=\mathbf{FPSPACE}\cap\mathbb{N}^{\mathbb{N}}.$

We also provide a characterisation of functions over the reals computable in
polynomial space. Inspired by [4], this is obtained by adding a limit schema
$ELim$ to $\mathbb{RCD}$. If we consider
$\overline{\mathbb{RCD}}=[0,1,\pi_{i}^{k},+,-,\times,\tanh,\cos,\pi,\frac{x}{2},\frac{x}{3};\textit{composition},\textit{robust
ODE},ELim]$ then:

######  Theorem 2 (Generic functions over the reals).

$\overline{\mathbb{RCD}}\cap\mathbb{R}^{\mathbb{R}}=\mathbf{FPSPACE}\cap\mathbb{R}^{\mathbb{R}}$
More generally:
$\overline{\mathbb{RCD}}\cap\mathbb{R}^{\mathbb{N}^{d}\times\mathbb{R}^{d^{\prime}}}=\mathbf{FPSPACE}\cap\mathbb{R}^{\mathbb{N}^{d}\times\mathbb{R}^{d}}$.

This article is organised as follows. In Section 2, we recall the concept of
dynamical systems and discuss some associated complexity issues. We introduce
the concept of robust ODE and prove that a robust ODE can be solved in
polynomial time (Theorem 4). This is obtained, using an original method for
solving ODE, optimising space, inspired by Savitch’s theorem. This provides
one direction of all the above theorems. The other direction is the object of
the following sections, starting from Section 3. We first recall some previous
results on discrete ODEs in Section 3. Using extensions of constructions from
[4], we then prove that we can simulate a Turing machine using robust
continuous ODEs in Section 4. This is obtained by simulating some discrete
ODEs using continuous ODEs, dealing with error corrections, and using the fact
that the functions are robust to a controlled error. The main result of
Section 3 is Theorem 10. It states we can simulate Turing machines robustly
with continuous ODEs when space remains polynomial. This theorem leads to the
proof of Theorem 1 in Section 4. In Section 5, we prove Theorem 2. In Section
6, we conclude and discuss future works.

#### Some basic concepts

When we say that a function over the real is computable this is always in the
sense of computable analysis: see e.g. [62, 46, 20]. A reference book for
issues related to complexity theory in computable analysis is [46]. See also
appendix 7.1.

##  2 Dynamical systems and associated complexity issues

###  2.1 Dynamical systems

Dynamical Systems are often used to model natural phenomena and in many
applied fields. The last decades have seen an impressive use of computers in
studying and analysing dynamical systems, with several visible theoretical
breakthroughs. A famous, notable example is the discovery of strange
attractors in models such as Lorenz attractors through numerical simulations
[47], with only 40 years later the mathematical proof of their existence [57].
The mathematical proof was obtained by checking some quantitative invariant
holds through computer-certified computations.

From a mathematical point of view, a _discrete-time dynamical system_ is given
by a set $D$, called _domain_ and some (possibly partial) function
$\mathbf{u}$ from $D$ to $D$. A trajectory the system is a sequence
$\mathbf{f}(t)$ evolving according to $\mathbf{u}$: that is
$\mathbf{f}(t+1)=\mathbf{u}\left(\mathbf{f}(t)\right)$ for all $t$. A
dynamical system can equivalently be described by its flow $\Phi$: by
definition, $\Phi(\mathbf{f}_{0},t)$ gives the position of the dynamics at
time $t$, for an initial position $\mathbf{f}(0)=\mathbf{f}_{0}$. It satisfies
the flow property

$\Phi(\mathbf{f}_{0},0)=\mathbf{f}_{0}\quad\Phi(\mathbf{f}_{0},t+t^{\prime})=\Phi(\Phi(\mathbf{f}_{0},t),t^{\prime})$
(1)

for all $t,t^{\prime}$. The transition dynamic $\mathbf{u}$ can be recovered
from the flow function since $\mathbf{u}(.)=\Phi(.,1)$. Hence, describing a
dynamical system by its dynamic $\mathbf{u}$ or by a flow function is
equivalent. All of this can be parametrised by some parameter $\mathbf{x}$:
$\mathbf{u}$ is also some function of $\mathbf{x}$, and
$\mathbf{f}(t+1)=\mathbf{u}\left(\mathbf{f}(t),\mathbf{x}\right)$ for all $t$,
and the flow function is $\Phi_{\mathbf{x}}(\mathbf{f}_{0},t)$.

Up to that point, we were considering discrete-time systems, but we can also
consider continuous-time dynamical systems: a _continuous-time dynamical
system_ is given by a set $D\subseteq\mathbb{R}^{d}$ and some ODE of the form
($\mathbf{f}^{\prime}$ denotes the derivative with respect to time variable
$t$)

$\mathbf{f}^{\prime}=\mathbf{u}({\mathbf{f}}(t))$ (2)

on $D$. A trajectory starting from $\mathbf{f}_{0}$ is a solution of the
associated Initial Value Problem (IVP), with $\mathbf{f}(0)=\mathbf{f}_{0}$.
If we consider $\Phi(\mathbf{f}_{0},t)$ as giving the value of the solution at
time $t$, starting from $\mathbf{f}_{0}$, it still satisfies the flow property
(1). Assuming sufficient regularity on $\Phi$ (namely that it is
differentiable), $\Phi$ satisfies some ODE, and hence giving a dynamical
system is equivalent to giving its flow: $\mathbf{u}(.)$ can be recovered by
$\mathbf{u}^{\prime}=\Phi^{\prime}(.,0)$. Here, we can still consider that all
of this is parameterised by some $\mathbf{x}$.

Hence, in a very general view, a dynamical system is given by some function
$\Phi$ satisfying the flow property. The fact that it is continuous time or
discrete time, is related to the nature of its time (i.e. second) variable,
that belongs to $\mathbb{N}$ or $\mathbb{Z}$ in the latter case, and to
$\mathbb{R}$ in the former. See [40] for a monography on the theory of
dynamical systems from a mathematical point of view.

In all previous discussions, we considered homogeneous dynamics in the sense
that $\mathbf{u}$ was only a function of $\mathbf{f}(t)$ and not also of $t$:
but, for example, for discrete-time, to cover the case
$\mathbf{f}(t+1)=\mathbf{u}\left(t,f(t)\right)$ for all $t$, it is sufficient
to consider $\overline{\mathbf{f}}(t)=(\mathbf{f}(t),t)$, and we come back to
the previous settings, as we can write
$\overline{\mathbf{f}}(t+1)=\overline{\mathbf{u}}(\overline{\mathbf{f}}(t))$,
with $\overline{\mathbf{u}}(\mathbf{f},t)=(\mathbf{u}(\mathbf{f},t),t+1)$,
that is a homogeneous dynamic.

As this is well-known, dynamical systems exhibit a very rich class of possible
behaviours. In particular, many dynamical systems are chaotic. We do not
intend here to recall how chaoticity can be defined in mathematical terms (see
[30, 40]), but this includes at least high sensitivity to initial conditions.
From our point of view, we just need to say this leads, in practice, to high
unpredictability in the long run. As observed in [52], while countless papers
exist in the literature about computations in dynamical systems, only a small
fraction of them address the problem rigorously; i.e., how far is the sought
actual quantity from the computed one? And can a computer perform such
computation up to a very small pre-specified error? An important example where
these questions are of interest is given by the reachability problem: given a
(finite) description of a dynamical system, a description of its initial
state, and the description of some “unsafe” states, the question is to tell
whether a trajectory can reach an unsafe state.

In the long run, dynamical systems may exhibit attractors. Although there is a
clear agreement about this intuitive concept, corresponding to the set of
points to which most points evolve, it is a hard task to provide a
mathematical definition covering all cases, and we refer to [49] for
discussions of many possible ways of defining this concept, and their
relations. In the general case, the attractors of dynamical systems, even if
the system is very simple, can be very rich. We refer to [52] for a
characterisation of the hardness of computing attractors from a computable
analysis point of view. We will somehow restrict to a very robust one
(numerically stable ones) for which the problem is tractable. Somehow, our
results state that the uncomputability discussed in [52] is intrinsically due
to the non-numerical stability of the considered dynamical systems there.

######  Example 1.

The van der Pol equation gives a very famous example of a dynamical system
admitting an attractor
$y^{\prime\prime}-\mu\left(1-y^{2}\right)x^{\prime}+y=0$ where $\mu$ is a
parameter. Introducing $z=y-y^{3}/3-\dot{y}/\mu$, it can be also described by
$y^{\prime}=\mu\left(y-\frac{1}{3}y^{3}-z\right)$,
$z^{\prime}=\frac{1}{\mu}y$. It is a classical mathematical exercise to prove
that the system has a limit cycle: this is usually done as an application of
Liénard’s theorem, which is established by studying the flow of the dynamics,
using qualitative arguments based on some formal mathematical statements: see,
e.g. [40].

Formally, a _discrete-time dynamical system_ is given by a set $D$, called
_domain_ and some (possibly partial) function $\mathbf{u}$ from $D$ to $D$. A
trajectory the system is a sequence $\mathbf{f}(t)$ evolving according to
$\mathbf{u}$: that is $\mathbf{f}(t+1)=\mathbf{u}\left(\mathbf{f}(t)\right)$
for all $t$. A _continuous-time dynamical system_ is given by a set
$D\subseteq\mathbb{R}^{d}$ and some ODE of the form

$\mathbf{f}^{\prime}=\mathbf{u}({\mathbf{f}}(t))$ (3)

on $D$. A trajectory starting from $\mathbf{f}_{0}$ is a solution of the
associated Initial Value Problem (IVP), with $\mathbf{f}(0)=\mathbf{f}_{0}$. A
dynamical system can equivalently be described by its flow:
$\Phi(\mathbf{f}_{0},t)$ gives the position of the dynamics at time $t$, for
an initial position $\mathbf{f}(0)=\mathbf{f}_{0}$. It satisfies the flow
property

$\Phi(\mathbf{f}_{0},0)=\mathbf{f}_{0}\quad\Phi(\mathbf{f}_{0},t+t^{\prime})=\Phi(\Phi(\mathbf{f}_{0},t),t^{\prime}).$
(4)

The dynamics or the flow can be parametrised by some parameter $\mathbf{x}$:
$\mathbf{u}$ is also some function of $\mathbf{x}$ and the flow function is
$\Phi_{\mathbf{x}}(\mathbf{f}_{0},t)$.

In the long run, dynamical systems may exhibit attractors. We refer to [49]
for discussions of many possible ways of defining this concept, and to [52]
for a characterisation of the hardness of computing attractors from a
computable analysis point of view. Somehow, our coming results state that the
uncomputability discussed in [52] is intrinsically due to the non-numerical
stability of the considered dynamical systems there.

###  2.2 Some complexity results on graphs

We need to discuss the hardness of solving IVP, or equivalently of computing
$\Phi(y,t)$. For pedagogical reasons, we first discuss the case of a simple
setting, namely the case of a (deterministic) directed graph. Indeed, observe
that a discrete-time dynamical system $(D,u)$ can also be seen as a particular
(deterministic) directed graph $G=(V,\rightarrow)$, where, in the general
case, $V$ is not necessarily finite: $G$ corresponds to $V=D$ and
$\rightarrow$ to the graph of the function $u$, i.e.
$\mathbf{x}_{t}\to\mathbf{x}_{t+1}$ iff
${\mathbf{x}}_{t+1}=\mathbf{u}\left({\mathbf{x}}_{t}\right)$. The obtained
graph is deterministic because any vertex has an outdegree $1$. Starting from
some point $\mathbf{x}_{0}$, there is at most one possible path, and
consequently, for a given time $T$, we can talk about its position at time
$t$, i.e. $\Phi(\mathbf{x}_{0},T)$ is $T$th element of this path: (as usual in
complexity theory, the length of some integer $x$ is the length of its binary
representation, denoted by $\ell(x)$).

######  Proposition 1 (The case of finite graphs).

Let $s(n)\geq\log(n)$ be space-constructible. Assume the vertices of
$G=(V,\rightarrow)$ can be encoded in binary using words of length $s(n)$.
Assume the relation $\rightarrow$ is decidable using a space polynomial in
$s(n)$. Then,

  * •

given the encoding of $\mathbf{u}\in V$ and of $\mathbf{v}\in V$, we can
decide whether there is some path from $\mathbf{u}$ to $\mathbf{v}$, in a
space polynomial in $s(n)$.

  * •

given the encoding of $\mathbf{u}\in V$, and integer $T$ in binary, we can
compute $\Phi(\mathbf{u},T)$, in a space polynomial in $s(n)$ and the length
of $T$.

The second item is even a characterisation of the complexity of the problem.
Indeed, the converse is true: If, given the encoding of $\mathbf{u}\in V$, and
integer $T$ in binary, we can compute $\Phi(\mathbf{u},T)$, in a space
polynomial in $s(n)$ and the length of $T$, then as $\rightarrow$ is given by
$\Phi(.,1)$, then $\rightarrow$ is decidable using a space polynomial in
$s(n)$.

###### Proof.

It is well-known that for finite graphs, given a directed graph
$G=(V,\rightarrow)$ and some vertices $\mathbf{u},\mathbf{v}\in V$, determine
whether there is some path between $\mathbf{u}$ and $\mathbf{v}$ in $G$,
denoted by $\mathbf{u}\stackrel{{\scriptstyle*}}{{\rightarrow}}\mathbf{v}$ is
in $\operatorname{NLOGSPACE}$: the rough idea is to guess non-
deterministically the intermediate nodes. The formal proof is detailed in
[55]. The same algorithm, working over representations of vertices, when
vertices are encoded using words of length $s(n)$ will work in
$\operatorname{NSPACE}(s(n))$ (with the addition of the binary encoding of $T$
if for the second item, if it bigger than $s(n)$). We then observe that
$\operatorname{NSPACE}(s(n))=\operatorname{SPACE}(s(n))$ from Savitch’s
theorem, recalled below. ∎

######  Theorem 3 (Savitch’s theorem, [55, Theorem 8.5]).

For any space-constructible111As proved in [55], this hypothesis can be
avoided, at the price of a slightly more complicated proof. function
$s\mathrel{\mathop{\ordinarycolon}}\mathbb{N}\to\mathbb{N}$ with $s(n)\geq\log
n$, we have
$\operatorname{NSPACE}(s(n))\subseteq\operatorname{SPACE}(s^{2}(n))$.

Recall that the key argument of the proof of Theorem 3 is to express the
question as a recursive procedure (expressing reachability in less than
$2^{t}$ steps, called ${CANYIELD(\mathbf{x}},\mathbf{y},t)$ in [55])
guaranteeing the required space complexity: we write that relation
$CANYIELD(\mathbf{x},\mathbf{y},t)$ is relation
$\mathbf{x}\rightarrow\mathbf{y}$ when $t=1$, and is relation
$\exists\mathbf{z}$ such that $CANYIELD(\mathbf{x},\mathbf{z},t/2)$ and
$CANYIELD(\mathbf{z},\mathbf{y},t/2)$ otherwise. If one prefers, this can also
be understood as “guessing” some intermediate node $\mathbf{z}$.

######  NB 2 (Attractor point of view).

We presented the above statement in terms of computing the flow
$\Phi(\mathbf{x},T)$. This could alternatively be interpreted in terms of
attractors. Indeed, when the above hypothesis holds, then dynamics is captured
by a graph. In the long run, in particular if $T$ is greater than the number
of vertices, any trajectory loops (i.e. reaches an attractor). The above
statement could then also be read as the fact that such an attractor is then
polynomial space computable.

###  2.3 Solving efficiently ODEs: what is known

This idea leads to an original method for solving ODEs. At least, this is
original for the numerical analysis literature, as far as we know.

We review what is known about the complexity of ODE solving. A more complete
survey is [38]. First, it is important to distinguish the case where we want
to solve the ODE on a bounded (hence a compact) domain, from the case of the
full domain $\mathbb{R}$: in the latter case, we might ask questions about the
evolution of the system on the long run, which is harder. Over a compact
domain, it is known that there exists some polynomial-time computable function
$u\mathrel{\mathop{\ordinarycolon}}[-1,1]\times[0,1]\rightarrow\mathbb{R}$
such that $f^{\prime}=u(f,t)$ has no computable solution, even over
$[0,\delta]$, for any $\delta>0$: see [45], extending [51, 1]. The involved
ODE has no unique solution. It is known over compact or non-compact domains
that if unicity holds, then its solution is computable [27, 28, 53]. However,
the complexity can be arbitrarily high [46, 48].

If we want to get to tractability, then some regularity hypotheses must be
assumed. A classical hypothesis is to assume the ODE to be Lipschitz.

Over a compact domain, it has been observed in several references (see e.g.
[46]) that a careful analysis of Euler’s method proves that, if
$u\mathrel{\mathop{\ordinarycolon}}B(0,1)\times[0,1]\rightarrow\mathbb{R}^{n}$,
with $B(0,1)\subseteq\mathbb{R}^{n}$, is a polynomial time computable
(right-)Lipschitz function then any solution
$f\mathrel{\mathop{\ordinarycolon}}[0,1]\rightarrow B(0,1)$ of
$f^{\prime}=u(f,t)$ must be polynomial-space computable: see the discussions
around Theorem 3.2 in [38] with the several references. Kawamura has proved in
[41] that there exists a polynomial-time computable function
$u\mathrel{\mathop{\ordinarycolon}}[-1,1]\times[0,1]\rightarrow\mathbb{R}$,
which satisfies a Lipschitz condition, such that the unique solution
$f\mathrel{\mathop{\ordinarycolon}}[0,1]\rightarrow\mathbb{R}$ takes values in
$[-1,1]$ and computing it leads to a $\mathbf{PSPACE}$-complete problem.
Hence, the question of solving ODEs over a compact domain in polynomial time
is exactly the question $\mathbf{PTIME}=\mathbf{PSPACE}$ [41], even for
$\mathcal{C}^{\infty}$-functions [42].

However, all these results are over compact domains, and dealing with non-
compact domains, i.e. in the long run, is harder. $\mathbf{PSPACE}$ membership
is not true, without stronger hypotheses. The difficulty comes from the
possibility of simulating any Turing machine by some finite-dimensional
polynomial ODE [37] over a non-compact domain. This leads to many
undecidability results for analytic, and even very simple ODEs. For example,
it follows that there is an analytic and computable function
$u\mathrel{\mathop{\ordinarycolon}}\mathbb{R}\rightarrow\mathbb{R}$ such that
the unique solution of the associated homogeneous ODE is defined on a non-
computable maximal interval of existence [37]. Futhermore, if we consider
$f_{1}(t)=e^{t}$, and $f_{i+1}(t)=e^{y_{i}(t)-1}$ then $f_{d}(t)$ is
$e^{e^{.e^{e^{t}}-1}}-1$, while all these functions are solutions of a simple
polynomial ODE over $\mathbb{R}^{d}$, namely $f^{\prime}_{1}(t)=f_{1}(t)$ and
$f^{\prime}_{d}(t)=f_{1}\dots f_{d}(t)$: i.e. a solution can grow faster than
a tower of exponentials in the description of the ODE, and hence is
necessarily intractable for time or space: see the discussion in [38, Section
3.2]. A possible way to analyse efficiency is then to analyse the complexity
of the solution assuming a bound on the growth of the function (i.e. using
parameterised complexity). It was proved in [14] that one can solve a
polynomial ODE in polynomial time assuming a bound on
$\mathbf{Y}(T)=\max_{0\leq t\leq T}\|\mathbf{f}(t)\|$. This result can be
extended to non-polynomial ODEs assuming polynomial-time computability of the
higher derivatives of $\mathbf{f}$ and an appropriate (polynomial) bound on
the growth of those derivatives [14].

The result for polynomial ODEs was later improved in [50], where it is proved
that the time $T$ and parameter $\mathbf{Y}$ can be replaced by a single
parameter, namely the length of the curve for polynomial ODEs. Furthermore,
this parameter does not need to be given as input to the algorithm. This is a
key argument for one direction of the moto “time complexity = length” we
mentioned several times. This is obtained using a non-classical method, from
the point of view of classical numerical analysis: this is not a fixed-order
numerical method, but somehow a method whose order is chosen as a function of
the inputs. A similar method was independently proposed in [56, 43],
considering parametrised complexity for analytic functions.

To get polynomial-time complexity over a non-compact domain, it is mandatory
not to use most classical methods from numerical analysis. Indeed, from the
general theory of numerical methods, presented for example in [29], every such
numerical method come with a given fixed order $k$, and from the general
theory of step-based methods, interval $[0,T]$ is divided into $N$ steps, each
of width $h_{n}\leq h$, and the error $\theta_{n}$ at step $n$ can be bounded
by $\theta_{n+1}\leq(1+\Lambda h_{n})\theta_{n}+|\epsilon_{n}|$: there is a
multiplicative error due to the Lipschitz constant $\Lambda$ of $u$, and some
numerical additive error $\epsilon_{n}$ due to the used precision. Using
Discrete Grönwall Lemma (see [29, page 213], or lemma 12 in appendix) the
final error satisfies between computed approximation $\widetilde{f}_{n}$ and
exact solution ${f}_{n}$ at step $n$ satisfies $\max_{0\leq n\leq
N}\left|\widetilde{f}_{n}-f_{n}\right|\leq e^{\Lambda
T}\left|\widetilde{f}_{0}-f_{0}\right|+\frac{e^{\Lambda
T}-1}{\Lambda}\left|\max\frac{\varepsilon_{n}}{h}\right|$. If $N=T/h$, we have
$\max_{0\leq n\leq N}\left|\widetilde{f}_{n}-f_{n}\right|\leq e^{\Lambda
T}\left|\widetilde{f}_{0}-f_{0}\right|+\frac{e^{\Lambda
T}-1}{\Lambda}\left|\max\frac{\varepsilon_{n}}{h}\right|$. To go to zero, the
last factor is chosen to be of the form $\exp(-\mathcal{O}(T))$.

The same problem happens when discussing space complexity: a non-classical
method is required to guarantee polynomial space complexity in the long run
(i.e. on the non-compact domain). As far as we know, no such method has yet
been proposed, and this is the purpose of the coming subsection. Actually, for
space complexity, in addition to all the problems mentioned, in all the above
space or time analyses, the problem is that the complexity is (possibly
implicitly) dependent on the Lipschitz constant or the length of the solution.
In a system as simple as linear dynamics, the state at time $T$ depends in
Lipchitz way from the state at time $0$, and the number or additional bits
required to guarantee some precision $2^{-n}$ growth linearly with $T$. But
the problem is that in a space polynomial in the input size, $T$ has no reason
to remain polynomial (consider, for example, a system simulating a Turing
machine, as we will consider soon). Hence, the required precision is possibly
exponential in the input size.

The above comment can be interpreted informally as the fact that “most” (this
could be “generic” in the sense of [52], i.e. (effective) descriptive theory)
dynamical systems are intrinsically unstable, and an error method introduced
at some step can make the method unavoidably incorrect in the long run unless
we have a means to “guess” what will happen.

######  NB 3 (Attractor point of view).

We presented the above statement in terms of computing the flow
$\Phi(\mathbf{f}_{0},t)$. But, this could alternatively be interpreted in
terms of attractors. The point is that computing the attractors of a given
dynamical system is hard in general, as this involves long-run behaviours.
This explains all the undecidability results obtained in [52], even for very
simple dynamics. However, as we will see, this is also explained by the fact
that the latter paper is discussing numerically unstable systems.

###  2.4 Solving efficiently ODEs: a space efficient method

This leads to an alternative approach to optimize space complexity: this can
be seen as either using a non-deterministic algorithm that “guesses” the
correct intermediate positions of the dynamics or, from the proof of Savitch’s
theorem approach, as an original recursive method to solve ODEs. As far as we
know, we have never seen such a method discussed in the literature for solving
ODEs.

Concretely: from the flow property, a strategy to compute
$\Phi(\mathbf{f}_{0},T)$ is either to use a particular numerical method if $T$
is small, says smaller than $\Delta>0$. Otherwise, we know that
$\Phi(\mathbf{f}_{0},T)=\Phi(\mathbf{z},T/2)$, where
$\mathbf{z}=\Phi(\mathbf{f}_{0},T/2)$. This always holds, so if we can compute
both quantities, we will solve the problem. The difficulty is that we cannot
precisely compute $\mathbf{z}$ in practice, but some numerical approximation
$\widetilde{\mathbf{z}}$. If the system is numerically stable, we may assume
this strategy works. The case when this strategy will not work is if the
trajectory starting from $\widetilde{\mathbf{z}}$, for the second half of the
work from time $T/2$ to $T$, has a behaviour different from the one starting
in $\mathbf{z}$: in other words, if there is a high instability somewhere,
namely in $\mathbf{z}$.

This leads to the following concept: we write $\mathbf{a}=_{n}\mathbf{b}$ for
$\|\mathbf{a}-\mathbf{b}\|\leq 2^{-n}$ for conciseness.

######  Definition 1 (Robust (continuous) ODE).

A function
$\mathbf{f}\mathrel{\mathop{\ordinarycolon}}\mathbb{R}\rightarrow\mathbb{R}$
is robustly ODE definable (from initial condition $\mathbf{g}$, and dynamic
$\mathbf{u}$) if

  1. 1.

it corresponds to the solution of the following continuous ODE:

$\mathbf{f}(0,\mathbf{x})=\mathbf{g}(\mathbf{x})\quad
and\quad\frac{\partial\mathbf{f}(t,\mathbf{x})}{\partial
t}=\mathbf{u}(\mathbf{f}(t,\mathbf{x}),\mathbf{h}(t,\mathbf{x}),t,\mathbf{x}),$
(5)

  2. 2.

and there is some rational $\Delta>0$, and some polynomial $p$ such that the
schema 5 is (polynomially) numerically stable on $[0,\Delta]$: for all integer
$n$, considering $\epsilon(n)=p(n+\ell(\mathbf{x}))$ we can compute
$\mathbf{f}(t,\mathbf{x})$ by working a precision $\epsilon(n)$: if you
consider any solution of $\tilde{\mathbf{x}}=_{\epsilon(n)}\mathbf{x}\quad
and\quad\tilde{\mathbf{h}}(t,\tilde{\mathbf{x}})=_{\epsilon(n)}{\mathbf{h}}(t,\tilde{\mathbf{x}})$,
and
$\tilde{\mathbf{f}}(0,\tilde{\mathbf{x}})=_{\epsilon(n)}\mathbf{g}(\mathbf{x})$
and $\frac{\partial\tilde{\mathbf{f}}(t,\tilde{\mathbf{x}})}{\partial
t}=_{\epsilon(n)}\mathbf{u}(\tilde{\mathbf{f}}(t,\tilde{\mathbf{x}}),\tilde{\mathbf{h}}(t,\tilde{\mathbf{x}}),t,\tilde{\mathbf{x}})$
then
$\tilde{\mathbf{f}}(t,\tilde{\mathbf{x}})=_{\epsilon(n)}\mathbf{f}(t,\mathbf{x})$
when $0\leq t\leq\Delta$.

  3. 3.

For $t\geq\Delta$, we can compute $\mathbf{f}(t,\mathbf{x})$ by computing some
approximation $\widetilde{\mathbf{f}(t/2,\mathbf{x})}$ of
$\mathbf{f}(t/2,\mathbf{x})$ at precision $\epsilon(n)$, i.e. of
$\Phi(\mathbf{g}(x),t/2)$, and then some approximation of
$\Phi(\widetilde{\mathbf{f}(t/2,\mathbf{y})},t/2)$, working at precision
$\epsilon(n)$.

######  Theorem 4.

Consider an IVP as in the previous definition. If
$\mathbf{g},\leavevmode\nobreak\ \mathbf{h}\leavevmode\nobreak\
\text{and}\leavevmode\nobreak\ \mathbf{u}$ are computable in polynomial space,
then the solution $\mathbf{f}$ can be computed in polynomial space.

###### Proof.

From definitions and above arguments, all bits of $\Phi(\mathbf{y},t)$ can be
computed non-deterministically with precision $n$ (i.e. at $2^{-n}$) using
computations with precision $\epsilon(n)$, hence is in
$\mathbf{NPSPACE}=\mathbf{PSPACE}$. From the argument of the proof of
Savitch’s theorem, this can also be turned into a deterministic polynomial
space recursive algorithm. ∎

The above theorem is the key argument to obtain one direction of our main
theorems. We now go in the reverse direction. This requires talking about
discrete ODEs, and some previous constructions.

##  3 Discrete ODEs: some previous results and constructions

###  3.1 Preliminary

We will use the concept of discrete ODE defined as follows (notice that we
will write $\frac{\delta\mathbf{f}}{\delta n}$ for discrete derivation, by
opposition of the classical $\frac{\partial\mathbf{f}}{\partial n}$ to help to
distinguish discrete vs continuous ODEs. )

######  Definition 2 (Discrete derivation, notation $\delta$).

For
$\mathbf{f}\mathrel{\mathop{\ordinarycolon}}\mathbb{N}\to\mathbb{R}^{d}\rightarrow\mathbb{R}^{d^{\prime}}$,
the _discrete derivation_ of $\mathbf{f}$ is $\frac{\delta\mathbf{f}}{\delta
n}(n,\mathbf{x})=\mathbf{f}(n+1,\mathbf{x})-\mathbf{f}(n,\mathbf{x})$.

###  3.2 Algebraic characterisation with discrete ODEs: state of the art

In this subsection, we review some of the results already obtained using
discrete ODEs.

######  NB 4.

Notice that we do not need any of these statements directly, even if we will
sometimes reuse some of their constructions (and some of their ideas).

Characterising $\mathbf{PTIME}$ over the integers:
The concept of derivation along the length was introduced in [9]. A
characterisation of $\mathbf{FPTIME}$ for functions over the integers has then
been obtained in [9]:

######  Theorem 5 (Functions over the integers [9]).

$\mathbb{LDL}\cap\mathbb{N}^{\mathbb{N}}=\mathbf{FPTIME}\cap\mathbb{N}^{\mathbb{N}}$,
for $\mathbb{LDL}=$
$[\mathbf{0},\mathbf{1},\mathbf{\pi}_{i}^{k},\mathrm{\ell}(x),\mathbf{+},\mathbf{-},\times,\textrm{sg}(x)\
;\textit{composition},\textit{linear length ODE}],$ with
$\mathbf{\pi}_{i}^{k}$ the projection function, and $\textrm{sg}(x)$ is 0 for
$x<0$ and 1 for $x>0$.

Toward the real numbers: characterising real sequences:
Later, the authors of [3] introduced

######  Definition 3 (Operation $ELim$).

Given
$\tilde{\mathbf{f}}\mathrel{\mathop{\ordinarycolon}}\mathbb{R}^{d}\times\mathbb{N}\to\mathbb{R}^{d^{\prime}}\in\mathbb{LDL}^{\bullet}$
such that for all $\mathbf{x}\in\mathbb{R}^{d}$, $n\in\mathbb{N}$,
$\|\tilde{\mathbf{f}}(\mathbf{x},2^{n})-\mathbf{f}(\mathbf{x})\|\leq 2^{-n}$
for some function $\mathbf{f}$, then $ELim(\tilde{\mathbf{f}})$ is the
(uniquely defined) corresponding function
$\mathbf{f}\mathrel{\mathop{\ordinarycolon}}\mathbb{R}^{d}\to\mathbb{R}^{d^{\prime}}$.

and then considered the class

$\overline{\mathbb{LDL}^{\bullet}}=[\mathbf{0},\mathbf{1},\mathbf{\pi}_{i}^{k},\mathrm{\ell}(x),\mathbf{+},\mathbf{-},\times,\overline{\textrm{cond}}(x),\frac{x}{2};{\textit{composition},\textit{linear
length ODE},ELim}],$

with $\overline{\textrm{cond}}(x)$ a sigmoid valuing 0 when $x<\frac{1}{4}$
and 1 when $x>\frac{3}{4}$. They proved this provides a characterisation of
functions from $\mathbb{N}$ to $\mathbb{R}$ computable in polynomial time.

######  Theorem 6 (Sequences of reals [3]).

$\overline{\mathbb{LDL}^{\bullet}}=\mathbf{FPTIME}\cap\mathbb{R}^{\mathbb{N}}$.

Characterisation of $\mathbf{PTIME}$ and $\mathbf{PSPACE}$ for functions over
the real with discrete ODEs:
The same authors later succeeded in obtaining a characterisation of functions
over the real computable in polynomial time and even space.

######  Theorem 7 ($\mathbf{FPTIME}$, , Generic functions over the reals [4]).

$\overline{\mathbb{LDL}^{\circ}}\cap\mathbb{R}^{\mathbb{N}^{d}\times\mathbb{R}^{d^{\prime}}}=\mathbf{FPTIME}\cap\mathbb{R}^{\mathbb{N}^{d}\times\mathbb{R}^{d}}$,
with

$\overline{\mathbb{LDL}^{\circ}}=[\mathbf{0},\mathbf{1},\mathbf{\pi}_{i}^{k},\mathrm{\ell}(x),\mathbf{+},\mathbf{-},\tanh,\frac{x}{2},\frac{x}{3};{\textit{composition},\textit{linear
length ODE},ELim}].$

Consider the following schema:

######  Definition 4 (Robust Discrete ODE [4]).

A bounded function $\mathbf{f}$ is robustly ODE definable if:

  1. 1.

it corresponds to the solution of the following discrete ODE:

$\mathbf{f}(0,\mathbf{x})=\mathbf{g}(\mathbf{x})\quad
and\quad\frac{\delta\mathbf{f}(t,\mathbf{x})}{\delta
t}=\mathbf{u}(\mathbf{f}(t,\mathbf{x}),\mathbf{h}(t,\mathbf{x}),t,\mathbf{x}),$
(6)

  2. 2.

where the schema (6) is (polynomially) numerically stable: there exists some
polynomial $p$ such that, for all integer $n$, writing
$\epsilon(n)=p(n+\ell(\mathbf{y}))$, if you consider any solution of
$\tilde{\mathbf{y}}=_{\epsilon(n)}\mathbf{y}\quad
and\quad\tilde{\mathbf{h}}(x,\tilde{\mathbf{y}})=_{\epsilon(n)}{\mathbf{h}}(x,\tilde{\mathbf{y}})$,
and
$\tilde{\mathbf{f}}(0,\tilde{\mathbf{y}})=_{\epsilon(n)}\mathbf{g}(\mathbf{y})$
and $\frac{\partial\tilde{\mathbf{f}}(x,\tilde{\mathbf{y}})}{\partial
x}=_{\epsilon(n)}\mathbf{u}(\tilde{\mathbf{f}}(x,\tilde{\mathbf{y}}),\tilde{\mathbf{h}}(x,\tilde{\mathbf{y}}),x,\tilde{\mathbf{y}})$
then
$\tilde{\mathbf{f}}(x,\tilde{\mathbf{y}})=_{\epsilon(n)}\mathbf{f}(x,\mathbf{y})$.

A robust discrete ODE is said to be linear if $\mathbf{u}$ is essentially
linear in $\mathbf{f}$ and $\mathbf{h}$.

Consider

$\overline{\mathbb{RLD}^{\circ}}=[\mathbf{0},\mathbf{1},\mathbf{\pi}_{i}^{k},\mathrm{\ell}(x),\mathbf{+},\mathbf{-},\tanh,\frac{x}{2},\frac{x}{3};{composition,robust\leavevmode\nobreak\
linear\leavevmode\nobreak\ ODE,ELim}].$

######  Theorem 8 ($\mathbf{FPSPACE}$, Generic functions over the reals [4]).

$\overline{\mathbb{RLD}^{\circ}}\cap\mathbb{R}^{\mathbb{R}}=\mathbf{FPSPACE}\cap\mathbb{R}^{\mathbb{R}}$
More generally:
$\overline{\mathbb{RLD}^{\circ}}\cap\mathbb{R}^{\mathbb{N}^{d}\times\mathbb{R}^{d^{\prime}}}=\mathbf{FPSPACE}\cap\mathbb{R}^{\mathbb{N}^{d}\times\mathbb{R}^{d}}$.

Notice that previous classes mix functions with integer and real arguments.
Furthermore, they all involve some various types of discrete ODEs. We need to
avoid all these issues, as we consider only continuous ODEs.

###  3.3 Simulating a discrete ODE using a continuous ODE

We first prove that it is possible to simulate a discrete ODE with a
continuous ODE. The underlying idea can be attributed to [19], and has been
improved in many ways by several authors. We present here the basic ideas,
reformulated in our context. A more precise analysis will come (Proposition
2).

######  Definition 5 ("Ideal iteration trick", [19]).

Consider the following initial value problem for a discrete ODE, given by
functions $\mathbf{g}$ and $\mathbf{u}$:

$\left\\{\begin{aligned} \mathbf{f}(0,\mathbf{x})&=\mathbf{g}(\mathbf{x})\\\
\frac{\delta\mathbf{f}}{\delta
t}(t,\mathbf{x})&=\mathbf{u}(\mathbf{f}(t,\mathbf{x}),t,\mathbf{x})\\\
\end{aligned}\right.$ (7)

Then, let
$\mathbf{G}(\mathbf{v},t,\mathbf{x})=\mathbf{u}(\mathbf{v},t,\mathbf{x})+\mathbf{v}$,
and consider the (continuous) IVP:

$\left\\{\begin{aligned}
&\mathbf{y}_{1}(0,\mathbf{x})=\mathbf{y}_{2}(0,\mathbf{x})=\mathbf{g}(\mathbf{x})\\\
&\mathbf{y}_{1}^{\prime}=c(\mathbf{G}(r(\mathbf{y}_{2}),r(t),\mathbf{x})-\mathbf{y}_{1})^{3}\theta(\sin(2\pi
t))\\\
&\mathbf{y}_{2}^{\prime}=c(r(\mathbf{y}_{1})-\mathbf{y}_{2})^{3}\theta(-\sin(2\pi
t))\\\ \end{aligned}\right.$ (8)

where $c$ a constant, $\theta(x)=0$ if $x\leq 0$ and $\theta(x)>0$ if $x>0$.
We abusively write $r(\mathbf{y})$ for the application of function
$r\mathrel{\mathop{\ordinarycolon}}\mathbb{R}\to\mathbb{R}$ componentwise on
vector $\mathbf{y}$. Here, $r$ is a rounding function: we mean, by
construction, $G$ preserves the integers, and $r$ is a function that maps a
real value close to some integer to this integer: assume, say, that for
$z\in[n-\frac{1}{4},n+\frac{1}{2}]$, $r(z)=n$, for any integer
$n\in\mathbb{Z}$.

Then, the solution of continuous ODE (8) simulates in a continuous way the
discrete ODE (7): Indeed, $\mathbf{y}_{1}$ corresponds to the actual
computation of the iterates of $\mathbf{G}$ (and hence computes the successive
values of $\mathbf{f}$) and $\mathbf{y}_{2}$ acts as a “memory" equation. Let
us detail how it works. We denote by $t=_{\epsilon}z$ the fact that
$\left|t-z\right|\leq\epsilon$.

######  NB 5.

We describe here an “ideal" computation, as $\theta(x)$ is exacly $0$ when
$x\leq 0$, and $r(z)$ is exactly some integer on suitable domains. Later in
the paper, we will deal with a not-so-ideal $\theta$ and $r$.

Initially,
$\mathbf{f}(0,\mathbf{x})=\mathbf{y}_{1}(0,\mathbf{x})=\mathbf{y}_{2}(0,\mathbf{x})=\mathbf{g}(\mathbf{x})$.
For $t\in[0,1/2]$, we have $\theta(-\sin(2\pi t))=0$, and hence
$\mathbf{y}_{2}^{\prime}=0$, so $\mathbf{y}_{2}$ is fixed and kept at value
$\mathbf{g}(\mathbf{x})$ for $t\in[0,\frac{1}{2}]$. Consequently, for
$t\in[0,1/2]$, $r(\mathbf{y}_{2})$ is also fixed and kept at value
$\mathbf{g}(\mathbf{x})$, and $r(t)$ is also fixed and kept at value $0$.
Consequently, on this interval, if we write $C(t)=c\theta(\sin(2\pi t))$, then
the dynamics of $\mathbf{y}_{1}$ is given by

$\mathbf{y}_{1}^{\prime}=C(t)(\mathbf{G}(\mathbf{g}(\mathbf{x}),0,\mathbf{x})-\mathbf{y}_{1})^{3}$
(9)

######  Lemma 1 (Analysis of ODE (9)).

The solution $\mathbf{y}_{1}(t,\mathbf{x})$ of ODE (9) is converging to
$G(\mathbf{g}(\mathbf{x}),0,\mathbf{x})$ for any initial condition.
Furthermore, for any initial condition $\mathbf{y}_{1}(0,\mathbf{x})\neq
G(\mathbf{g}(\mathbf{x}),0,\mathbf{x})$, we have
$\left\|\mathbf{y}_{1}(\frac{1}{2},\mathbf{x})-\mathbf{G}(\mathbf{g}(x),0,\mathbf{x})\right\|\leq\frac{\sqrt{2}}{2\sqrt{\int_{0}^{\frac{1}{2}}C(z)dz}}$.
In particular, for any $m\in\mathbb{N}$, we can select constant $c$ such that
for any initial condition $\mathbf{y}_{1}(0,\mathbf{x})$,

$\left\|\mathbf{y}_{1}(\frac{1}{2},\mathbf{x})-\mathbf{G}(\mathbf{g}(x),0,\mathbf{x})\right\|\leq
2^{-m}.$

###### Proof of Lemma 1, Adapted from [19].

If initially, or at any instant
$\mathbf{y}_{1}(0,\mathbf{x})=G(\mathbf{g}(x),0,\mathbf{x})$ then the result
holds, as $\mathbf{y}^{\prime}_{1}=0$, and $\mathbf{y}_{1}$ remains constant.
Otherwise, we have

$\frac{\mathbf{y}_{1}^{\prime}}{\mathbf{G}(\mathbf{g}(\mathbf{x}),0,\mathbf{x})-\mathbf{y}_{1})^{3}}=C(t).$

Integrating this equality between $0$ and $t$, we obtain

$\frac{\mathbf{1}}{2(\mathbf{G}(\mathbf{g}(\mathbf{x}),0,\mathbf{x})-\mathbf{y}_{1}(t))^{2}}-\frac{\mathbf{1}}{2(\mathbf{G}(\mathbf{g}(\mathbf{x}),0,\mathbf{x})-\mathbf{y}_{1}(0))^{2}}=\int_{0}^{t}C(z)dz,$

hence

$\frac{\mathbf{1}}{2(\mathbf{G}(\mathbf{g}(\mathbf{x}),0,\mathbf{x})-\mathbf{y}_{1}(t))^{2}}\geq\int_{0}^{t}C(z)dz.$

This yields the property. ∎

Consequently, $\mathbf{y}_{1}(t,\mathbf{x})$ will approach
$\mathbf{G}(\mathbf{g}(\mathbf{x}),0,\mathbf{x})=\mathbf{f}(1,\mathbf{x})$ on
this interval. Thus,
$\mathbf{y}_{1}(\frac{1}{2},\mathbf{x})=_{\epsilon}\mathbf{f}(1,\mathbf{x})$
and $\mathbf{y}_{2}(\frac{1}{2},\mathbf{x})=\mathbf{g}(\mathbf{x})$, for some
$\epsilon>0$, that we can consider less than $\frac{1}{4}=2^{-2}$, by
selecting a big enough constant $c$ (just taking $m=2$ above). At
$t=\frac{1}{2}$, $\mathbf{y}_{1}$ will hence have simulated one step of
discrete ODE (7).

Now, for $t\in[\frac{1}{2},1]$ the roles of $\mathbf{y}_{1}$ and
$\mathbf{y}_{2}$ are exchanged : $\mathbf{y}_{1}^{\prime}(t,\mathbf{x})=0$, so
$\mathbf{y}_{1}$ is kept fixed, $\mathbf{y}_{2}$ approaches
$r(\mathbf{y}_{1})=\mathbf{f}(1,\mathbf{x})$, thus
$\mathbf{y}_{1}(1,\mathbf{x})=_{\epsilon}\mathbf{y}_{2}(1,x)=_{\epsilon}\mathbf{f}(1,\mathbf{x})$.

By induction, from the same reasoning, we obtain that, for all
$n\in\mathbb{N}$,
$\mathbf{y}_{1}(n,\mathbf{x})=_{\epsilon}\mathbf{y}_{2}(n,\mathbf{x})=_{\epsilon}\mathbf{f}(n,\mathbf{x})$,
and actually, we also have
$\mathbf{y}_{1}(t+\frac{1}{2},\mathbf{x})=_{\epsilon}\mathbf{y}_{2}(t,\mathbf{x})=_{\epsilon}\mathbf{f}(n,\mathbf{x})$
for all $t\in[n,n+\frac{1}{2}]$, for any integer $n$.

To implement such an ODE, we have to fix a function $\theta(x)$ with the above
property. Taking $\operatorname{ReLU}(x)=\max(0,x)$ would satisfy it, but it
is not a derivable function, and hence would not lead to a (classical) ODE. We
could then take $\theta(x)=0$ for $x\leq 0$, and $\exp(-1/x)$ for $x>0$. The
point is that such a function is not real analytic. The base functions we
consider in our class $\mathbb{RCD}$ are all real analytic, and real analytic
functions are preserved by composition, so we cannot get such a function by
compositions from our base functions. Futhermore, it is known that a real
analytic function that is constant on some interval (we assumed it is $0$ for
$x\leq 0$!) is constant. Hence the above-considered function $\theta(x)$
cannot be real analytic. So, implementing this trick cannot be done directly
using our base functions, using only compositions.

In Proposition 2, we will do a similar construction, but dealing with errors
and not exact function $\theta(z)$ and $r(x)$. Furthermore, here the purpose
of function $r$ was to correct errors around integers, i.e. around
$\mathbb{N}$: this will be possibly around other $\mathbb{N}\delta$ for some
$\delta>0$.

###  3.4 Encoding of Turing machines configurations

Our proofs rely on some constructions from [4]. Concretely, we need to
simulate the execution of a Turing machine (TM) $\mathcal{M}$ by some
dynamical system over the reals. This requires to encode the configurations of
a Turing machine into some real numbers. We recall some of the definitions and
constructions from [4].

Consider a Turing machine defined by $\mathcal{M}=(\Sigma,Q,I,F,\delta)$, with
$\Sigma$ the working alphabet, $Q$ the set of states, $I,F\subseteq Q$
respectively the sets of initial and final states,
$\delta\mathrel{\mathop{\ordinarycolon}}Q\times\Sigma\rightarrow
Q\times\Sigma\times\\{\leftarrow,\rightarrow\\}$ the transition function. For
some practical reasons, similar to the ones in [4], we assume that the working
alphabet is made of the symbols $1$ and $3$, and that the blank symbol is
symbol $0$.

We explicit the encoding we will use. We assume $Q=\\{0,1,\dots,|Q|-1\\}$. Let

$\dots l_{-k}l_{-k+1}\dots l_{-1}l_{0}r_{0}r_{1}\dots r_{n}.\dots$

denote the content of the tape of the Turing machine $M$. In this
representation, the head is in front of symbol $r_{0}$, and
$l_{i},r_{i}\in\\{0,1,3\\}$ for all $i$. Furthermore, we assume that there are
no non-blank symbols between two blank symbols, i.e. that blank symbols, i.e.
symbol $0$, can only be eventually on the right, or eventually on the left.
Such a configuration $C$ can be denoted by $C=(q,l,r)$, where
$l,r\in\Sigma^{\omega}$ are words over alphabet $\Sigma=\\{0,1,3\\}$ and $q\in
Q$ denotes the internal state of $M$.

Now, write:
$\gamma_{word}\mathrel{\mathop{\ordinarycolon}}\Sigma^{\omega}\to\mathbb{R}$
for the function that maps a word $w=w_{0}w_{1}w_{2}\dots$ to the dyadic
(hence real) number $\gamma_{word}(w)=\sum_{n\geq 0}w_{n}4^{-(n+1)}$.

The idea is that configuration $C$ can also be encoded by some element
$\overline{C}=(q,\overline{l},\overline{r})\in\mathbb{N}\times\mathbb{R}^{2}$,
by considering $\overline{r}=\gamma_{word}(r)$ and
$\overline{l}=\gamma_{word}(l)$. In other words, we encode the configuration
of a bi-infinite tape Turing machine $M$ by real numbers using their radix 4
encoding, but using only digits $1$,$3$. Notice that this lives in
$Q\times[0,1]^{2}$. Denoting the image of
$\gamma_{word}\mathrel{\mathop{\ordinarycolon}}\Sigma^{\omega}\to\mathbb{R}$
by $\mathcal{I}$, this even lives in $Q\times\mathcal{I}^{2}$.

In other words, we consider the following encodings:
$\gamma_{config}(C)=(q,\overline{l},\overline{r})$ with
$\overline{l}=l_{0}4^{-1}+l_{-1}4^{-2}+\dots+l_{-k}4^{-(k+1)}+\dots$ and
$\overline{r}=r_{0}4^{-1}+r_{1}4^{-2}+\dots+l_{n}4^{-(n+1)}+\dots$.

###  3.5 Revisiting some previous constructions

We denote by $\mathbb{RCD_{*}}$ the algebra
$[0,1,\pi_{i}^{k},+,-,\times,\tanh,\cos,\pi,\frac{x}{2},\frac{x}{3};\textit{composition}]$.
This is close to the class

$\mathbb{LDL}^{\circ}=[\mathbf{0},\mathbf{1},\mathbf{\pi}_{i}^{k},\mathrm{\ell}(x),\mathbf{+},\mathbf{-},\tanh,\frac{x}{2},\frac{x}{3};{composition,linear\leavevmode\nobreak\
length\leavevmode\nobreak\ ODE}],$

considered in [5, 4], but without the function $\mathrm{\ell}(x)$, and wihtout
the possiblity of defining functions using linear length ODE (and with
multiplication added).

We will reuse some of the construction from [4] (some corrections and more
details can be found in [5]) but avoid systematically any use of linear length
ODE and the length function $\ell(x)$. Furthermore, the class considered in
[4] is mixing functions from the integers to the reals, and from the reals to
the reals, and we need to keep only functions over the reals.

The following was stated in [4, Lemma 19].

######  Lemma 2.

We denote by $Y(x,2^{m+2})$ the function $\frac{1+\tanh(2^{m+2}x)}{2}$. For
all integer $m$, for all $x\in\mathbb{R}$,
$|\operatorname{ReLU}(x)-xY(x,2^{m+2})|\leq 2^{-m}$, where
$\operatorname{ReLU}(x)=\max(0,x)$.

First, we observe that considering $Y(x,z)=\frac{1+\tanh(4xz)}{2}$ would yield
a function in $\mathbb{RCD_{*}}$ with the same property: we avoid the
computation of $2^{m}$ by a substitution of a variable, and using a
multiplication. We then write $\operatorname{ReLU\mathfrak{-s}}(Y,x)$ for
$xY(x,z)$: we have
$|\operatorname{ReLU\mathfrak{-s}}(2^{m},x)-\operatorname{ReLU}(x)|\leq
2^{-m}$.

In particular, this was used to prove we can uniformly approximate the
continuous sigmoid functions (when $1/(b-a)$ is in $\mathbb{LDL}^{\circ}$)
defined as: $\operatorname{\mathfrak{s}}(a,b,x)=0$ whenever $w\leq a$,
$\frac{x-a}{b-a}$ whenever $a\leq x\leq b$, and $1$ whenever $b\leq x$. The
above trick provides a new version of [4, Lemma 20].

######  Lemma 3 (Uniform approximation of any piecewise continuous sigmoid).

Assume $a,b,\frac{1}{b-a}$ is in $\mathbb{RCD_{*}}$. Then there is some
function $\operatorname{\mathfrak{\mathcal{C}-s}}(z,a,b,x)\in\mathbb{RCD_{*}}$
such that for all integer $m$,
$|\operatorname{\mathfrak{\mathcal{C}-s}}(2^{m},a,b,x)-\operatorname{\mathfrak{s}}(a,b,x)|\leq
2^{-m}$.

###### Proof.

Take
$\operatorname{\mathfrak{\mathcal{C}-s}}(z,a,b,x)=\frac{(x-a)Y(x-a,z2^{1+c})-(x-b)Y(x-b,z2^{1+c})}{b-a}$.
observing that
$(b-a)\operatorname{\mathfrak{s}}(a,b,x)={\operatorname{ReLU}(x-a)-\operatorname{ReLU}(x-b)}$.
From triangle inequality, it will hold, choosing $c$ with $\frac{1}{b-a}\leq
2^{c}$. ∎

The authors of [4] proved the existence of some function corresponding to a
continuous (controlled) approximation of the fractional part function:

######  Theorem 9 ([4, Lemma 28]).

There exists some function
$\xi\mathrel{\mathop{\ordinarycolon}}\mathbb{N}^{2}\to\mathbb{R}$ in
$\mathbb{LDL}^{\circ}$ such that for all $n,m\in\mathbb{N}$ and
$x\in[-2^{n},2^{n}]$, whenever $x\in[\lfloor x\rfloor+\frac{1}{8},\lfloor
x\rfloor+\frac{7}{8}]$,

$\left|\xi(2^{m},2^{n},x)-\\{x-\frac{1}{8}\\}\right|\leq 2^{-m}.$

We say that some real function is a real extension of a function over the
integers if they coincide for integer arguments It is not clear that we have a
real extension of $\xi$ in our algebra $\mathbb{RCD_{*}}$, but if we add a
real extension of such a function, from the proof of [4, Corollary 22], we
obtain the bestiary of functions considered in [4, Corollary 22]: we write
$\mathbb{RCD_{*}}+\xi$ for the algebra where some real extension of function
$\xi$ is added as a base function.

######  Corollary 1 (A bestiary of functions).

There exist

  1. 1.

$\xi_{1},\xi_{2}\mathrel{\mathop{\ordinarycolon}}\mathbb{N}^{2}\times\mathbb{R}\mapsto\mathbb{R}\leavevmode\nobreak\
\in\mathbb{RCD_{*}}+\xi$ such that, for all $n,m\in\mathbb{N}$, $\lfloor
x\rfloor\in[-2^{n}+1,2^{n}]$, whenever $x\in[\lfloor
x\rfloor-\frac{1}{2},\lfloor x\rfloor+\frac{1}{4}]$ ,
$|\xi_{1}(2^{m},2^{n},x)-\\{x\\}|\leq 2^{-m}$, and whenever $x\in[\lfloor
x\rfloor,\lfloor x\rfloor+\frac{3}{4}]$ ,
$|\xi_{2}(2^{m},2^{n},x)-\\{x\\}|\leq 2^{-m}$.

  2. 2.

$\sigma_{1},\sigma_{2}\mathrel{\mathop{\ordinarycolon}}\mathbb{N}^{2}\times\mathbb{R}\mapsto\mathbb{R}\leavevmode\nobreak\
\in\mathbb{RCD_{*}}+\xi$ such that, for all $n,m\in\mathbb{N}$, $\lfloor
x\rfloor\in[-2^{n}+1,2^{n}]$, whenever $x\in[\lfloor
x\rfloor-\frac{1}{2},\lfloor x\rfloor+\frac{1}{4}]$,
$|\sigma_{1}(2^{m},2^{n},x)-\lfloor x\rfloor|\leq 2^{-m}$, and whenever $x\in
I_{2}=[\lfloor x\rfloor,\lfloor x\rfloor+\frac{3}{4}]$,
$|\sigma_{2}(2^{m},2^{n},x)-\lfloor x\rfloor|\leq 2^{-m}$.

  3. 3.

$\lambda\mathrel{\mathop{\ordinarycolon}}\mathbb{N}^{2}\times\mathbb{R}\mapsto[0,1]\in\mathbb{RCD_{*}}+\xi$
such that for all $m,n\in\mathbb{N}$, $\lfloor x\rfloor\in[-2^{n}+1,2^{n}]$,
whenever $x\in[\lfloor x\rfloor+\frac{1}{4},\lfloor x\rfloor+\frac{1}{2}]$,
$|\lambda(2^{m},2^{n},x)-0|~\leq 2^{-m}$, and whenever $x\in[\lfloor
x\rfloor+\frac{3}{4},\lfloor x\rfloor+1]$, $|\lambda(2^{m},2^{n},x)-1|\leq
2^{-m}$.

  4. 4.

$\mod_{2}\mathrel{\mathop{\ordinarycolon}}\mathbb{N}^{2}\times\mathbb{R}\mapsto[0,1]\in\mathbb{RCD_{*}}+\xi$
such that for all $m,n\in\mathbb{N}$, $\lfloor x\rfloor\in[-2^{n}+1,2^{n}]$,
whenever $x\in[\lfloor x\rfloor-\frac{1}{4},\lfloor x\rfloor+\frac{1}{4}]$,
$|\mod_{2}(2^{m},2^{n},x)$-$\lfloor x\rfloor\mod 2|\leq 2^{-m}$.

  5. 5.

$\div_{2}\mathrel{\mathop{\ordinarycolon}}\mathbb{N}^{2}\times\mathbb{R}\mapsto[0,1]\in\mathbb{RCD_{*}}+\xi$
such that for all $m,n\in\mathbb{N}$, $\lfloor x\rfloor\in[-2^{n}+1,2^{n}]$,
whenever $x\in[\lfloor x\rfloor-\frac{1}{4},\lfloor x\rfloor+\frac{1}{4}]$,
$|\div_{2}(2^{m},2^{n},x)-\lfloor x\rfloor//2|\leq 2^{-m}$, with $//$ the
integer division.

###### Proof.

There were given by $\xi_{1}(M,N,x)=\xi(M,N,x-\frac{3}{8})-\frac{1}{2}$,
$\xi_{2}(M,N,x)=\xi(N,x-\frac{7}{8})$, $\sigma_{i}(M,N,x)=x-\xi_{i}(M,N,x)$,
$\lambda(M,N,x)=\operatorname{\mathfrak{\mathcal{C}-s}}(2M,1/4,1/2,\xi(2M,N,x-9/8))$,
$\mod_{2}(M,N,x)=1-\lambda(M,N/2,\frac{1}{2}x+\frac{7}{8})$,
$\div{2}(M,N,x)=\frac{1}{2}(\sigma_{1}(M,N,x)-\mod_{2}(M,N,x))$. ∎

Similarly, the equivalent of [4, Lemmas 23,24 and 25] still hold in
$\mathbb{RCD_{*}}+\xi$. Namely:

######  Lemma 4.

There exists $\operatorname{\mathcal{C}-if}\in\mathbb{RCD_{*}}+\xi$ such that,
$l\in[0,1]$, if we take $|d^{\prime}-0|\leq 1/4$, then
$|\operatorname{\mathcal{C}-if}(2^{m},d^{\prime},l)-0|\leq 2^{-m}$, and if we
take $|d^{\prime}-1|\leq 1/4$, then
$|\operatorname{\mathcal{C}-if}(2^{m},d^{\prime},l)-l|\leq 2^{-m}$.

######  Lemma 5.

Let $\alpha_{1},\alpha_{2},\dots,$ $\alpha_{n}$ be some integers, and
$V_{1},V_{2},\dots,V_{n}$ some constants. We write
$\operatorname{send}({\alpha_{i}\mapsto V_{i}})_{i\in\\{1,\dots,n\\}}$ for the
function that maps any $x\in[\alpha_{i}-1/4,\alpha_{i}+1/4]$ to $V_{i}$, for
all $i\in\\{1,\dots,n\\}$.

There is some function in $\mathbb{RCD_{*}}+\xi$, that we write
$\operatorname{\mathcal{C}-send}(2^{m},{\alpha_{i}\mapsto
V_{i}})_{i\in\\{1,\dots,n\\}}$, that maps any
$x\in[\alpha_{i}-1/4,\alpha_{i}+1/4]$ to a real at distance at most $2^{-m}$
of $V_{i}$, for all $i\in\\{1,\dots,n\\}$.

######  Lemma 6.

Let $N$ be some integer. Let $\alpha_{1},\alpha_{2},\dots,\alpha_{n}$ be some
integers, and $V_{i,j}$ for $1\leq i\leq n$ some constants, with $0\leq j<N$.
We write $\operatorname{send}({(\alpha_{i},j)\mapsto
V_{i,j}})_{i\in\\{1,\dots,n\\},j\in\\{0,\dots,N-1\\}}$ for the function that
maps any $x\in[\alpha_{i}-1/4,\alpha_{i}+1/4]$ and $y\in[j-1/4,j+1/4]$ to
$V_{i,j}$, for all $i\in\\{1,\dots,n\\}$, $j\in\\{0,\dots,N-1\\}$.

There is some function in $\mathbb{RCD_{*}}+\xi$, that we write
$\operatorname{\mathcal{C}-send}(2^{m},{(\alpha_{i},j)\mapsto
V_{i,j}})_{i\in\\{1,\dots,n\\},j\in\\{0,\dots,N-1\\}}$, that maps any
$x\in[\alpha_{i}-1/4,\alpha_{i}+1/4]$ and $y\in[j-1/4,j+1/4]$ to a real at
distance at most $2^{-m}$ of $V_{i,j}$, for all $i\in\\{1,\dots,n\\}$,
$j\in\\{0,\dots,N-1\\}$.

#### Working with one step of a Turing machine

As the proof of [4, Lemmas 30] is done using all the functions provided by
these lemmas, we obtain:

######  Lemma 7.

We can construct some function $\overline{\textit{Next}}$ in
$\mathbb{RCD_{*}}+\xi$ that simulates one step of $M$, i.e. that computes the
Next function sending a configuration $\overline{C}$ of Turing machine $M$ to
$\overline{C}^{\prime}$, where $C^{\prime}$ is the next one:
$\|\textit{Next}(2^{m},2^{S},\overline{C})-\overline{C}^{\prime}\|\leq
2^{-m}$. Furthermore, it is robust to errors on its input, up to space $S$:
considering $\|\tilde{C}-\overline{C}\|\leq 4^{-(S+2)}$,
$\|\textit{Next}(2^{m},2^{S},\tilde{C})-\overline{C}^{\prime}\|\leq 2^{-m}$
remains true.

#### Converting integers an dyadics to words and conversely

The authors of [4] also defined some functions for converting integers and
dyadics to their encoding as words, and conversely. Namely, they consider the
following encoding: every digit in the binary expansion of dyadic $d$ is
encoded by a pair of symbols in the radix $4$ expansion of
$\overline{d}\in\mathcal{I}\cap[0,1]$: digit $0$ (respectively: $1$) is
encoded by $11$ (resp. $13$) if before the “decimal” point in $d$, and digit
$0$ (respectively: $1$) is encoded by $31$ (resp. $33$) if after. For example,
for $d=101.1$ in base $2$, $\overline{d}=0.13111333$ in base $4$. Conversely,
given $\overline{d}$, they provided a way to construct $d$. This corresponds
to [4, Lemmas 33 and 34]:

######  Lemma 8 (From $\mathbb{N}$ to $\mathcal{I}$).

We can construct some function
$\textit{Decode}\mathrel{\mathop{\ordinarycolon}}\mathbb{N}^{2}\to\mathbb{R}$
in $\mathbb{LDL}^{\circ}$ that maps $m$ and $n$ to some point at distance less
than $2^{-m}$ from $\gamma_{word}(\overline{n})$.

######  Lemma 9 (From $\mathcal{I}$ to $\mathbb{R}$, and multiplying in
parallel).

We can construct some function
$\textit{EncodeMul}\mathrel{\mathop{\ordinarycolon}}\mathbb{N}^{2}\times[0,1]\times\mathbb{R}\to\mathbb{R}$
in $\mathbb{LDL}^{\circ}$ that maps $m$, $2^{S}$,
$\gamma_{word}(\overline{d})$ and (bounded) $\lambda$ to some real at distance
at most $2^{-m}$ from $\lambda d$, whenever $\overline{d}$ is of length less
than $S$.

As for $\xi$, it is not clear that we have some real extensions of these
functions in $\mathbb{RCD_{*}}$: we write
$\mathbb{RCD_{*}}+\xi+\textit{Decode}+\textit{Encode}$ for the algebra where
some real extension of these functions is added as a base function.

###  3.6 Constructing the missing functions

We need a way to construct some substitute of “missing functions” ($\xi$,
Decode and EncodeMul). As all of them are defined using discrete ODEs, an idea
is to use a continuous ODE to simulate the respective discrete ODEs: we hence
revisit the construction of the ideal iteration trick of Section 3.3, dealing
with errors and not exact functions $\theta(z)$ and $r(x)$.

The key is to revisit Lemma 1, and do a more detailed analysis of possible
involved errors in dynamics of the form (9). This equation has been studied by
various authors in several articles, including [21, 22, 35, 13, 33]. We use
the following statement from [33, Lemma 4.5], [11, Lemma 5.2], obtained
basically by a case analysis of error propagations in Lemma 1.

######  Lemma 10 (Improved error analysis of ODE (9), [33, Lemma 4.5] [11,
5.2]).

Consider a point $b\in\mathbb{R}$, some $\gamma>0$ some reals $t_{0}<t_{1}$,
and a function
$\phi\mathrel{\mathop{\ordinarycolon}}\mathbb{R}\rightarrow\mathbb{R}$ with
the property that $\phi(t)\geq 0$ for all $t\geq t_{0}$ and
$\int_{t_{0}}^{t_{1}}\phi(t)dt>0$. Let $\rho,\delta\geq 0$ and let
$\overline{b},E\mathrel{\mathop{\ordinarycolon}}\mathbb{R}\rightarrow\mathbb{R}$
be functions such that that $|\overline{b}(t)-b|\leq\rho$ and
$|E(t)|\leq\delta$ for all $t\geq t_{0}$. Then the IVP defined by

$z^{\prime}=c(\overline{b}(t)-z)^{3}\phi(t)+E(t)$

with the initial condition $z\left(t_{0}\right)=\overline{z}_{0}$, where
$\gamma>0$ and $c\geq\frac{1}{2\gamma^{2}\int_{t_{0}}^{t_{1}}\phi(t)dt}$
satisfies

  1. 1.

$\left|z\left(t_{1}\right)-b\right|<\rho+\gamma+\delta\left(t_{1}-t_{0}\right)$,
independently of the initial condition $\overline{z}_{0}\in\mathbb{R}$

  2. 2.

$\min\left(\overline{z}_{0},b-\rho\right)-\delta\left(t_{1}-t_{0}\right)\leq
z(t)\leq\max\left(\overline{z}_{0},b+\rho\right)+\delta\left(t_{1}-t_{0}\right)$
for all $t\in\left[t_{0},t_{1}\right]$.

######  Proposition 2 (Simulating a discrete ODE by a continuous ODE).

Assumge $\mathbf{G}$ is almost constant around $\mathbb{N}\delta$ and $r$ is a
rounding function around $\mathbb{N}\delta$ for some $\delta>0$. Suppose that,
in (8), we replace function $\theta(z)$ and function $r(z)$ by some suitable
approximations: we take $\theta(x)=\operatorname{ReLU}(x)$,
$\theta_{\epsilon^{\prime}}(x)$, $r_{\epsilon^{\prime}(z)}$ such that
$\theta(z)=_{\epsilon^{\prime}}\theta(z)$, and
$r_{\epsilon^{\prime}}(x)=_{\epsilon^{\prime}}$, and take constant $c$ big
enough. Then the solution of the obtained ODE will continuously simulate the
discrete ODE (7), with same bounds as in the analysis in Section 3.3, i.e.
with error at most $\epsilon$ if $\epsilon^{\prime}$ is taken sufficiently
small. To guarantee $\epsilon=2^{-n}$, it is sufficient to take
$\epsilon^{\prime}=2^{-p(n)}$ and
$\theta_{\epsilon^{\prime}}(x)=\operatorname{ReLU\mathfrak{-s}}(2^{p(n)},x)$
for some polynomial $p$.

###### Proof.

The key is that involved errors propagate additively, from Lemma 10. Namely,
they are in $\mathcal{O}(\epsilon^{\prime})$, but they are then corrected from
the reasoning in Section 3.3: rounding function corrects errors or order
$\epsilon$ whenever its argument is at distance less than $1/4\delta$ of some
$n\delta$ exactly as in the reasoning in Section 3.3 (where $\delta=1$, even
if now it introduces some error $\epsilon^{\prime}$ at every step; but the
latter is corrected at the next step). Observe that the involved constant $c$,
is of order $2^{n}$.

We did up to that point the reasoning, assuming that parameter $\mathbf{x}$ is
fixed. But if we consider a function that either depends in a controlled way
on $t$ on that parameter, or that is such that a small perturbation on its
inputs does not change much its output (we mean in a controlled way, in the
way we consider the rounding function $r$), then the analysis remains
perfectly valid, even when this parameter may not be exact. ∎

We claim that for all $n\in\mathbb{N}$,
$\mathbf{y}_{1}(n,\mathbf{x})=_{\epsilon}\mathbf{y}_{2}(n,\mathbf{x})=_{\epsilon}\mathbf{f}(n,\mathbf{x})$,
and
$\mathbf{y}_{1}(t+\frac{1}{2},\mathbf{x})=_{\epsilon}\mathbf{y}_{2}(t,\mathbf{x})=_{\epsilon}\mathbf{f}(n,\mathbf{x})$
for all $t\in[n,n+\frac{1}{2}]$.

For $n=0$, initially
$\mathbf{f}(0,\mathbf{x})=\mathbf{y}_{1}(0,\mathbf{x})=\mathbf{y}_{2}(0,\mathbf{x})=\mathbf{g}(\mathbf{x})$.
For $t\in[n,n+1/2]$, we have $\theta(-\sin(2\pi t))=_{\epsilon^{\prime}}0$,
and hence $\mathbf{y}_{2}^{\prime}=_{\epsilon^{\prime}}0$, so $\mathbf{y}_{2}$
is kept close to value $\mathbf{g}(\mathbf{x})$ for $t\in[0,\frac{1}{2}]$,
with an error less than $\frac{1}{2}\epsilon^{\prime}$.

Consequently, for $t\in[0,1/2]$, $r(\mathbf{y}_{2})$ is kept close to a
constant value $\mathbf{g}(\mathbf{x})$, when an error less than
$\epsilon^{\prime}$, if we choose $\epsilon^{\prime}<\frac{1}{4}\delta$.
Meanwhile, $r(t)$ is also at a value close to $n$ with error less than
$\epsilon^{\prime}$.

Consequently, on this interval, if we write $C(t)=c\theta(\sin(2\pi t))$, then
the dynamics of $\mathbf{y}_{1}$ is given by a dynamic of the form of Lemma
10. This lemma states that $\mathbf{y}_{1}(t,\mathbf{x})$ will approach
$\mathbf{G}(\mathbf{g}(\mathbf{x}),0,\mathbf{x})=\mathbf{f}(1,\mathbf{x})$ on
this interval, with an error of order
$\epsilon^{\prime}+\epsilon^{\prime}+\frac{1}{2}\epsilon^{\prime}$.

Here the hypothesis that $\mathbf{G}$ is almost constant around
$\mathbb{N}\delta$ means that its value is guaranteed to be at
$\epsilon^{\prime}$ from $\mathbf{G}(\mathbf{g}(\mathbf{x}),0,\mathbf{x})$ on
the interval.

Thus,
$\mathbf{y}_{1}(\frac{1}{2},\mathbf{x})=_{\epsilon/2}\mathbf{f}(1,\mathbf{x})$,
if we choose $\frac{5}{2}\epsilon^{\prime}<\epsilon/2$. At $t=n+\frac{1}{2}$,
$\mathbf{y}_{1}$ will hence have simulated one step of discrete ODE (7), with
error less than $\epsilon/2$, and $\mathbf{y}_{2}$ will be close to
$\mathbf{g}(\mathbf{x})$ with error less than $\epsilon^{\prime}<\epsilon/2$.

Now, for $t\in[n+\frac{1}{2},n+1]$ the roles of $\mathbf{y}_{1}$ and
$\mathbf{y}_{2}$ are exchanged :
$\mathbf{y}_{1}^{\prime}(t,\mathbf{x})=_{\epsilon^{\prime}}0$, so
$\mathbf{y}_{1}$ is kept almost fixed, with a new error less than
$\frac{1}{2}\epsilon^{\prime}$. In the same time $\mathbf{y}_{2}$ approaches
$r(\mathbf{y}_{1})=\mathbf{f}(1,\mathbf{x})$ by Lemma 10, with some new error
of order less than $\frac{5}{2}\epsilon^{\prime}<\epsilon/2$.

Consequently, we get the property at rank $n+1$.

######  NB 6.

Observe that, somehow, the constructions always replace every function with a
function that does not change much locally (i.e., change in a controlled way).
This is the key that provides a robust ODE as in Definition 1, leading to
polynomial space complexity by Theorem 4.

In other words, whenever we have some discrete ODE as in (7) defining some
function $\mathbf{f}(t,\mathbf{x})$, we can construct some continuous ODE,
using only functions from $\mathbb{RCD_{*}}$, such that one of its projection
provides a function $\mathbf{f}(z,t,\mathbf{x})$, with the guarantee
$\mathbf{f}(2^{n},t,\mathbf{x})$ is $2^{-n}$ close to
$\mathbf{f}(n,\mathbf{x})$, whenever $t$ is close (at distance less than
$1/4$) to some integer $n$.

This works, as we can obtain such a $r_{\epsilon^{\prime}}(x)$ from the
functions from Corollary 1: Consider
$r(x,2^{m})=\sigma_{2}(2^{m},2^{n},x+\frac{1}{4})$ that works over $\lfloor
x\rfloor\in[-2^{n}+1,2^{n}]$, and observe that this is sufficient to apply the
trick for the required functions, from the form of the considered discrete ODE
in [4].

Except that we have a bootstrap problem: $\xi$ was defined using a discrete
ODE in [4], and as the functions from Corollary 1 are defined above using
$\xi$, we cannot apply this reasoning to get function $\xi$. But the point is
that for the special case of $\xi$, it is easy to construct a function in
$\mathbb{RCD}$ that corresponds to some real extension of $\xi$, as we have
functions such as $\sin(x)=\cos(\frac{\pi}{2}-x)$ and $\pi$.

######  Lemma 11.

Function $\xi$ has some real extension in $\mathbb{RCD_{*}}$.

###### Proof.

If we succeed to obtain a function $i(2^{m},2^{n},x)$ that values $\lfloor
x\rfloor$ whenever $x\in[\lfloor x\rfloor,\lfloor x\rfloor+\frac{3}{4}]$, we
are done, as we can then obtain $\xi(2^{m},2^{n},x)$ by considering
$\xi(2^{m},2^{n},x)=x+\frac{7}{8}-i(2^{m},2^{n},x+\frac{7}{8})$.

A possible solution is then the following: consider function
$R_{e}(x)\mathrel{\mathop{\ordinarycolon}}=\operatorname{\mathfrak{s}}(x,0,e/2)$,
and then $t_{e}(x)=(1-R_{e}(\sin(2\pi x)))((1-R_{e}(\sin(4\pi x)))$. If we put
aside some interval of width $e/2$ around $\frac{1}{2}$ and $\frac{7}{8}$
where it takes values in $[0,1]$, it values $0$ on $[\lfloor x\rfloor,\lfloor
x\rfloor+\frac{7}{8}]$, and then $1$ on $[\lfloor x\rfloor+\frac{7}{8},\lfloor
x\rfloor+1]$ (see following graphical illustration). We can then consider
$I_{e}(t)=8\int_{0}^{t}t_{e}(x)dx$ (i.e. the solution of ODE
$l^{\prime}_{e}=8t_{e}$), and then $i(t)=_{e.t}l_{e}(t)$. It is then
sufficient to replace $\operatorname{\mathfrak{s}}$ by
$\operatorname{\mathfrak{\mathcal{C}-s}}$, in the above expressions, in order
to control the error and make it smaller than $2^{-m}$.

Here is a graphical representation of $R_{\frac{1}{10}}(x)$:

Then of $R_{\frac{1}{10}}(\sin(2\pi x))$:

and $R_{\frac{1}{10}}(\sin(4\pi x))$:

We then get $t_{\frac{1}{10}}(x)$.

Its integral is then close to $\frac{1}{8}\lfloor x\rfloor$ on $[\lfloor
x\rfloor,\lfloor x\rfloor+\frac{3}{4}]$.

∎

Consequently, this is true that we can substitute a discrete ODE with a
continuous ODE for the required functions Decode and EncodeMul: just replace
$\xi$ in the involved schemas by the above function. Notice that we can also
easily get a real extension of the function that maps $n$ to $2^{n}$.

###  3.7 Working with all steps of a Turing machines

We can then go from one step of a Turing machine, to arbitrarily many steps.
We are following the idea of [4], but replacing discrete ODEs with continuous
ODEs.

######  Theorem 10.

Consider some Turing machine $M$ that computes some function
$f\mathrel{\mathop{\ordinarycolon}}\Sigma^{*}\to\Sigma^{*}$ in some polynomial
space $S(\ell(\omega))$ on input $\omega$. One can construct some function
$\tilde{f}\mathrel{\mathop{\ordinarycolon}}\mathbb{N}^{2}\times\mathbb{R}\to\mathbb{R}$
in $\mathbb{RCD}$ that does the same: we have
$\tilde{f}(2^{m},2^{S(\ell(\omega))},\gamma_{word}(\omega))$ that is at most
$2^{-m}$ far from $\gamma_{word}(f(\omega))$.

###### Proof.

We denote by $\mathcal{M}$ the Turing machine computing $f$. Similarly to the
arguments in [4], we can state that there exists a function $Exec$ solution of
a robust linear discrete ODE (E) that "computes" the execution of
$\mathcal{M}$, with $C_{init}$ the initial configuration :

$(E)\mathrel{\mathop{\ordinarycolon}}\left\\{\begin{aligned}
&Exec(2^{m},0,2^{S},C_{init})=C_{init}\\\ &\frac{\delta
Exec(2^{m},t,2^{S},C_{init})}{\delta
t}=Next(2^{m},2^{S},Exec(2^{m},t,2^{S},C_{init}))-Exec(2^{m},t,2^{S},C_{init})\end{aligned}\right.$

For any configuration $\overline{C}$ of $\mathcal{M}$, let write
$F(\overline{C}))=F(2^{m},2^{S},\overline{C})=Next(2^{m},2^{S},\overline{C})+\overline{C}$,
associated to the righthand side of the above discrete ODE. Denoting by
$\tilde{C}$ the errorless encoding of the configuration $C$, from the
constructions of [4] (Lemma 7), it is true that if
$\left|\overline{C}-\tilde{C}\right|\leq 4^{-(S+2)}$, then
$\left|F(\overline{C})-F(\tilde{C})\right|\leq 4^{-(S+2)}$. $F$ does not
change much locally on the space of configuration. Denoting by $S$ the space
of $\mathcal{M}$, and replacing $m$ by $m+2S+4$ as in [4], we have
$\left|Next(2^{m},2^{S},\overline{C})-\overline{C}\right|\leq 4^{-(S+2)}$. So
at each step of the TM, the error is fixed (and bounded). We can then apply
the above arguments (Proposition 2) to simulate continuously (E), with some
controlled error: all involved quantities have encoding polynomials in the
size of the inputs. ∎

##  4 Proof of Theorem 1

###### Proof.

$\subseteq$: In this direction, we just need to prove that $\mathbb{RCD}$
contains only functions over the reals that are computable in polynomial
space. Indeed, then for a function
$\mathbf{f}\mathrel{\mathop{\ordinarycolon}}\mathbb{R}^{d}\to\mathbb{R}^{d^{\prime}}$
sending every integer ${\mathbf{n}}\in\mathbb{N}^{d}$ to the vicinity of some
integer of $\mathbb{N}^{d}$, at a distance less than $1/4$, by approximating
its value with precision $1/4$ on its input arguments, and taking the closest
integer, we will get a function from the integers to the integers, that
corrresponds to $\operatorname{DP}(f)$, and that will be in
$\mathbf{FPSPACE}\cap\mathbb{N}^{\mathbb{N}}$.

This is indeed the case, since i) all the base functions of
$\overline{\mathbb{RCD}}$ are in $\mathbf{FPSPACE}$: they are even in
$\mathbf{FPTIME}$, see [46] ii) $\mathbb{R}^{\mathbb{R}}\cap\mathbf{FPSPACE}$
is stable under composition. iii) stability under robust ODE follows from
Theorem 4.

$\supseteq$: In the other direction, we use an argument similar to [4]:
namely, as the function is polynomial space computable, this means that there
is a polynomial space computable function
$g\mathrel{\mathop{\ordinarycolon}}\mathbb{N}^{d^{\prime\prime}+1}\to\\{1,3\\}^{*}$
so that on $\mathbf{m},2^{n}$, it provides the encoding
$\overline{\phi(\mathbf{m},n)}$ of some dyadic $\phi(\mathbf{m},n)$ with
$\|\phi(\mathbf{m},n)-\mathbf{f}(\mathbf{m})\|\leq 2^{-n}$ for all
$\mathbf{m}$. The problem is then to decode, compute and encode the result to
produce this dyadic. More precisely, from Theorem 10, we get $\tilde{g}$ with

$|\tilde{g}(2^{e},2^{p(max(\mathbf{m},n))},\textit{Decode}(2^{e},\mathbf{m},n))-\gamma_{word}(g(\mathbf{m},n))|\leq
2^{-e}$

for some polynomial $p$ corresponding to the time required to compute $g$, and
$e=\max(p(max(\mathbf{m},n)),n)$. Then we need to transform the value to the
correct dyadic: we mean
$\tilde{\mathbf{f}}(\mathbf{m},n)=\textit{EncodeMul}(2^{e},2^{t},\tilde{g}(2^{e},2^{t},\textit{Decode}(2^{e},\mathbf{m},n)),1)$,
where
$t=p(max(\mathbf{m},n))$, $e=\max(p(max(\mathbf{m},n)),n)$ provides a solution
with $\|\tilde{\mathbf{f}}(\mathbf{m},2^{n})-\mathbf{f}(\mathbf{m})\|\leq
2^{-n}.$ ∎

##  5 Proof of Theorem 2

###### Proof.

$\subseteq$: To prove that
$\overline{\mathbb{RCD}}\subseteq\mathbb{R}^{\mathbb{R}}\cap\mathbf{FPSPACE}$,
we only need to add to the previous arguments that
$\mathbb{R}^{\mathbb{R}}\cap\mathbf{FPSPACE}$ is also stable under $ELim$.

$\supseteq$: In this direction, we have the same issue as in [4]: the strategy
of decoding, working with the Turing machine, and encoding is not guaranteed
to work for all inputs. But, we can solve it by using an adaptative barycenter
technique as in [4].

We recall the principle here for a function whose domain is $\mathbb{R}$, but
it can be generalised to $\mathbb{R}^{d}$. The idea is to construct some
function
$\lambda\mathrel{\mathop{\ordinarycolon}}\mathbb{N}^{2}\times\mathbb{R}\rightarrow[0,1]$
definable in $\mathbb{RCD_{*}}$ as in Corollary 1, but with a continuous ODE :
Adapting the proof from [4] and using the simulation of $\xi$ in our
continuous framework, we can consider
$\lambda(2^{m},N,x)=\varPsi(\Xi(2^{m+1},N,x-9/8))$ where
$\varPsi(x)=\operatorname{\mathfrak{\mathcal{C}-s}}(2^{m+1},1/4,1/2,x)$. In
particular, by definition, $\lambda\in\mathbb{RCD_{*}}$. Thus, by Lemma 9, if
$\lambda(2^{m},N,x)=_{2^{-m}}0$, then $\sigma_{2}(2^{m},N,x)=_{2^{-m}}\lfloor
x\rfloor$. If $\lambda(2^{m},N,x)=_{2^{-m}}1$, then
$\sigma_{1}(2^{m},N,x)=_{2^{-m}}\lfloor x\rfloor$ and if
$\lambda(2^{m},N,x)\in(0,1)$, then $\sigma_{1}(2^{m},N,x)=_{2^{-m}}\lfloor
x\rfloor+1$ and furthermore $\sigma_{2}(2^{m},N,x)=_{2^{-m}}\lfloor x\rfloor$.
So,

$\lambda(\cdot,2^{n},x)\textit{Formula}_{1}(x,u,M,n)+(1-\lambda(\cdot,2^{n},n))\textit{Formula}_{2}(x,u,M,n)$

and we are sure to be close (up to some bounded error) to some $2^{-m}$
approximation of a function $f$. ∎

##  6 Conclusion

We characterised polynomial space using an algebraically defined class of
functions, using a finite set of basic functions, closure under composition,
and a schema for defining functions from robust ODEs. We proposed a concept of
robust ODEs solvable in polynomial space. As far as we know, this is an
original method for solving ODEs optimising space. It is based on classical
constructions such as Savitch’s theorem. We extended existing
characterisations to a characterisation of functions over the reals and not
only over the integers.

The interesting message from our statements is that we provide a clear and
simple concept associated with continuous ODEs for space: space corresponds to
the precision for numerically stable systems. Hence, compiled with [16], we
now know the length of solutions corresponds to time and precision to memory.

Considering future work: We have an algebraically defined class of functions.
It remains to know whether this could be transferred at the level of
polynomial ODE. We know that soluting of polynomial ODEs defined a very robust
class of functions, stable by many operations: sum, products, division, ODE
solving, etc: see [34, 12]. Hence, all the base functions we consider in our
algebraic class can be turned into polynomial ODEs, by adding some variables.
It would be interesting to understand if we could define space complexity
directly at the level of polynomial ODEs, using precision.

Recently, a characterisation of $\mathbf{PSPACE}$ was obtained for polynomial
ODEs using rather ad-hoc definitions in [33, 11] and working over a non-
compact space. Could our characterisation be put at this simplest class of
ODEs, but working with precision? The point is that the characterisation there
uses unbounded domains, hence, precision is harder to interpret in their
constructions, where the schemas are somehow done to control errors.

Of course, from our statements, adding any $\mathbf{FPSPACE}$-computable
function over the reals among the base functions would not change the class.
However, we did not intend to minimise the number of base functions. For
example, $\tanh(t)$ is solution of ODE $f^{\prime}=1+f^{2}$ and $\cos(t)$ can
be obtained by the two dimensional ODE $y_{1}^{\prime}=-y_{2}$,
$y_{2}^{\prime}=y_{1}$. Minimising the number of base functions is also left
for future work. We believe that even in this settings, proving space
complexity corresponds to precision is already significant, independently of
this question of a minimal set of base functions.

## References

  * [1] Oliver Aberth.  The failure in computable analysis of a classical existence theorem for differential equations.  Proceedings of the American Mathematical Society, 30:151–156, 1971\.
  * [2] Daniel Bertschinger, Christoph Hertrich, Paul Jungeblut, Tillmann Miltzow, and Simon Weber.  Training Fully Connected Neural Networks is $\exists\mathbb{R}$-Complete.  Preprint, https://doi. org/10.48550/arXiv, 2204, 2022.
  * [3] Manon Blanc and Olivier Bournez.  A characterization of polynomial time computable functions from the integers to the reals using discrete ordinary differential equations.  In Jérôme Durand-Lose and György Vaszil, editors, Machines, Computations, and Universality - 9th International Conference, MCU 2022, Debrecen, Hungary, August 31 - September 2, 2022, Proceedings, volume 13419 of Lecture Notes in Computer Science, pages 58–74. Springer, 2022.  MCU’22 Best Student Paper Award.  doi:10.1007/978-3-031-13502-6\\_4.
  * [4] Manon Blanc and Olivier Bournez.  A characterisation of functions computable in polynomial time and space over the reals with discrete ordinary differential equations: Simulation of turing machines with analytic discrete odes (MFCS’2023 best paper award).  In 48th International Symposium on Mathematical Foundations of Computer Science (MFCS 2023), volume 272. Schloss Dagstuhl-Leibniz-Zentrum für Informatik, 2023.
  * [5] Manon Blanc and Olivier Bournez.  Simulation of turing machines with analytic discrete ODEs: FPTIME and FPSPACE over the reals characterised with discrete ordinary differential equations.  arXiv preprint arXiv:2307.11747, 2023.
  * [6] Manon Blanc and Olivier Bournez.  Quantifiying the robustness of dynamical systems. relating time and space to length and precision.  In Computer Science Logic CSL’24, Naples, Italy, February 2024.
  * [7] Lenore Blum, Felipe Cucker, Mike Shub, and Steve Smale.  Complexity and Real Computation.  Springer, 1998.
  * [8] Lenore Blum, Mike Shub, and Steve Smale.  On a theory of computation and complexity over the real numbers; NP completeness, recursive functions and universal machines.  Bulletin of the American Mathematical Society, 21(1):1–46, July 1989.
  * [9] Olivier Bournez and Arnaud Durand.  Recursion schemes, discrete differential equations and characterization of polynomial time computation.  In Peter Rossmanith, Pinar Heggernes, and Joost-Pieter Katoen, editors, 44th Int Symposium on Mathematical Foundations of Computer Science, MFCS, volume 138 of LIPIcs, pages 23:1–23:14. Schloss Dagstuhl - Leibniz-Zentrum für Informatik, 2019.
  * [10] Olivier Bournez and Arnaud Durand.  A characterization of functions over the integers computable in polynomial time using discrete ordinary differential equations.  Computational Complexity, 32(2):7, 2023.
  * [11] Olivier Bournez, Riccardo Gozzi, Daniel S Graça, and Amaury Pouly.  A continuous characterization of PSPACE using polynomial ordinary differential equations.  Journal of Complexity, 77:101755, august 2023.  URL: https://www.sciencedirect.com/science/article/pii/S0885064X23000249?dgcid)=author.
  * [12] Olivier Bournez, Daniel Graça, and Amaury Pouly.  On the Functions Generated by the General Purpose Analog Computer.  Information and Computation, 257:34–57, 2017.  arXiv:1602.00546, doi:10.1016/j.ic.2017.09.015.
  * [13] Olivier Bournez, Daniel Graça, and Amaury Pouly.  Computing with polynomial ordinary differential equations.  Journal of Complexity, 36:106 – 140, 2016.  URL: http://www.sciencedirect.com/science/article/pii/S0885064X16300255, doi:http://dx.doi.org/10.1016/j.jco.2016.05.002.
  * [14] Olivier Bournez, Daniel S Graça, and Amaury Pouly.  On the complexity of solving initial value problems.  In Proceedings of the 37th International Symposium on Symbolic and Algebraic Computation, pages 115–121, 2012.
  * [15] Olivier Bournez, Daniel S. Graça, and Amaury Pouly.  Polynomial Time corresponds to Solutions of Polynomial Ordinary Differential Equations of Polynomial Length. The General Purpose Analog Computer and Computable Analysis are two efficiently equivalent models of computations (ICALP’2017 Track B best paper award).  In 43rd International Colloquium on Automata, Languages, and Programming, ICALP 2016, July 11-15, 2016, Rome, Italy, volume 55 of LIPIcs, pages 109:1–109:15. Schloss Dagstuhl - Leibniz-Zentrum fuer Informatik, 2016.
  * [16] Olivier Bournez, Daniel S. Graça, and Amaury Pouly.  Polynomial Time corresponds to Solutions of Polynomial Ordinary Differential Equations of Polynomial Length.  Journal of the ACM, 64(6):38:1–38:76, 2017.  doi:10.1145/3127496.
  * [17] Olivier Bournez and Amaury Pouly.  A universal ordinary differential equation.  In International Colloquium on Automata Language Programming, ICALP’2017, 2017.
  * [18] Olivier Bournez and Amaury Pouly.  A survey on analog models of computation.  In Handbook of Computability and Complexity in Analysis, pages 173–226. Springer, 2021.
  * [19] M. S. Branicky.  Universal computation and other capabilities of hybrid and continuous dynamical systems.  Theoretical Computer Science, 138(1):67–100, 6 February 1995.
  * [20] Vasco Brattka, Peter Hertling, and Klaus Weihrauch.  A tutorial on computable analysis.  In New computational paradigms, pages 425–491. Springer, 2008.
  * [21] Mark Braverman.  Hyperbolic Julia sets are poly-time computable.  Electronic Notes in Theoretical Computer Science, 120:17–30, 2005\.
  * [22] Manuel L. Campagnolo, Cristopher Moore, and José Félix Costa.  Iteration, inequalities, and differentiability in analog computers.  Journal of Complexity, 16(4):642–660, 2000.
  * [23] Tian Qi Chen, Yulia Rubanova, Jesse Bettencourt, and David K Duvenaud.  Neural ordinary differential equations.  In Advances in Neural Information Processing Systems, pages 6571–6583, 2018.
  * [24] P. Clote.  Computational models and function algebras.  In Edward R. Griffor, editor, Handbook of Computability Theory, pages 589–681. North-Holland, Amsterdam, 1998.
  * [25] Peter Clote and Evangelos Kranakis.  Boolean functions and computation models.  Springer Science & Business Media, 2013.
  * [26] Alan Cobham.  The intrinsic computational difficulty of functions.  In Y. Bar-Hillel, editor, Proceedings of the International Conference on Logic, Methodology, and Philosophy of Science, pages 24–30. North-Holland, Amsterdam, 1962.
  * [27] Pieter Collins and Daniel S Graça.  Effective computability of solutions of ordinary differential equations the thousand monkeys approach.  Electronic Notes in Theoretical Computer Science, 221:103–114, 2008\.
  * [28] Pieter Collins and Daniel S. Graça.  Effective Computability of Solutions of Differential Inclusions The Ten Thousand Monkeys Approach.  Journal of Universal Computer Science, 15(6):1162–1185, 2009.
  * [29] J.-P. Demailly.  Analyse Numérique et Équations Différentielles.  Presses Universitaires de Grenoble, 1996.
  * [30] R. L. Devaney.  An Introduction to Chaotic Dynamical Systems.  Addison-Wesley, 2nd edition, 1989.
  * [31] Kousha Etessami and Mihalis Yannakakis.  On the complexity of Nash equilibria and other fixed points.  SIAM Journal on Computing, 39(6):2531–2597, 2010.
  * [32] Francois Fages, Guillaume Le Guludec, Olivier Bournez, and Amaury Pouly.  Strong turing completeness of continuous chemical reaction networks and compilation of mixed analog-digital programs.  In Computational Methods in Systems Biology-CMSB 2017, 2017.  CMSB’2017 Best Paper Award.
  * [33] Riccardo Gozzi.  Analog Characterization of Complexity Classes.  PhD thesis, Instituto Superior Técnico, Lisbon, Portugal and University of Algarve, Faro, Portugal, 2022.
  * [34] Daniel S. Graça.  Computability with Polynomial Differential Equations.  PhD thesis, Instituto Superior Técnico, 2007.
  * [35] Daniel S. Graça, Manuel L. Campagnolo, and Jorge Buescu.  Robust simulations of Turing machines with analytic maps and flows.  In B. Cooper, B. Loewe, and L. Torenvliet, editors, Proceedings of CiE’05, New Computational Paradigms, volume 3526 of Lecture Notes in Computer Science, pages 169–179. Springer-Verlag, 2005.
  * [36] Daniel S. Graça and José Félix Costa.  Analog computers and recursive functions over the reals.  19(5):644–664, 2003.
  * [37] Daniel S. Graça, N. Zhong, and J. Buescu.  Computability, noncomputability and undecidability of maximal intervals of IVPs.  Transactions of the American Mathematical Society, 2006.  To appear.
  * [38] Daniel S. Graça and Ning Zhong.  Handbook of Computability and Complexity in Analysis, chapter Computability of Differential Equations.  Springer., 2018.
  * [39] Daniel S Graça and Ning Zhong.  Robust non-computability and stability of dynamical systems.  arXiv preprint arXiv:2305.14448, 2023.
  * [40] Morris W. Hirsch, Stephen Smale, and Robert Devaney.  Differential Equations, Dynamical Systems, and an Introduction to Chaos.  Elsevier Academic Press, 2003.
  * [41] Akitoshi Kawamura.  Lipschitz continuous ordinary differential equations are polynomial-space complete.  In 2009 24th Annual IEEE Conference on Computational Complexity, pages 149–160. IEEE, 2009.
  * [42] Akitoshi Kawamura, Hiroyuki Ota, Carsten Rösnick, and Martin Ziegler.  Computational complexity of smooth differential equations.  Logical Methods in Computer Science, 10, 2014.
  * [43] Akitoshi Kawamura, Florian Steinberg, and Holger Thies.  Parameterized complexity for uniform operators on multidimensional analytic functions and ode solving.  In International Workshop on Logic, Language, Information, and Computation, pages 223–236. Springer, 2018.
  * [44] Patrick Kidger.  On neural differential equations.  arXiv preprint arXiv:2202.02435, 2022.
  * [45] Ker-I Ko.  On the computational complexity of ordinary differential equations.  Information and Control, 58(1-3):157–194, July/August/September 1983.
  * [46] Ker-I Ko.  Complexity Theory of Real Functions.  Progress in Theoretical Computer Science. Birkhaüser, Boston, 1991\.
  * [47] E. N. Lorenz.  Deterministic nonperiodic flow.  Journal of the Atmospheric Science, 20:130–141, 1963.
  * [48] Webb Miller.  Recursive function theory and numerical analysis.  Journal of Computer and System Sciences, 4(5):465–472, October 1970\.
  * [49] John Milnor.  On the concept of attractor.  Communications in Mathematical Physics, 99:177–195, 1985.
  * [50] Amaury Pouly and Daniel Graça.  Computational complexity of solving polynomial differential equations over unbounded domains.  Theoretical Computer Science, 2016.
  * [51] Marian Boykan Pour-El and J. Ian Richards.  A computable ordinary differential equation which possesses no computable solution.  Annals of Mathematical Logic, 17:61–90, 1979.
  * [52] Cristobal Rojas and Mathieu Sablik.  On the algorithmic descriptive complexity of attractors in topological dynamics.  arXiv preprint arXiv:2311.15234, 2023.
  * [53] Keijo Ruohonen.  An effective Cauchy-Peano existence theorem for unique solutions.  International Journal of Foundations of Computer Science, 7(2):151–160, 1996.
  * [54] Claude E. Shannon.  Mathematical theory of the differential analyser.  Journal of Mathematics and Physics MIT, 20:337–354, 1941.
  * [55] Michael Sipser.  Introduction to the Theory of Computation.  PWS Publishing Company, 1997.
  * [56] Holger Thies.  Uniform computational complexity of ordinary differential equations with applications to dynamical systems and exact real arithmetic.  PhD thesis, University of Tokyo, Graduate School of Arts and Sciences, 2018.
  * [57] Warwick Tucker.  A rigorous ODE solver and Smale’s 14th problem.  Foundations of Computational Mathematics, 2:53–117, 2002.
  * [58] Alan Turing.  On computable numbers, with an application to the Entscheidungsproblem.  Proceedings of the London Mathematical Society, 42(2):230–265, 1936\.  Reprinted in Martin Davis. The Undecidable: Basic Papers on Undecidable Propositions, Unsolvable Problems and Computable Functions. Raven Press, 1965.
  * [59] Bernd Ulmann.  Analog computing.  Walter de Gruyter, 2013.
  * [60] Bernd Ulmann.  Analog and hybrid computer programming.  De Gruyter Oldenbourg, 2020.
  * [61] Veritasum.  Future computers will be radically different (analog computing).  Youtube video, 2022.  URL: https://www.youtube.com/watch?v=GVsUOuSjvcg.
  * [62] Klaus Weihrauch.  Computable Analysis: an Introduction.  Springer, 2000.

##  7 Appendix

###  7.1 Very basics of computable analysis

In very short, repeating [6], the idea behind classical computability and
complexity is to fix some representations of objects (such as graphs,
integers, etc.) using finite words over some finite alphabet, say
$\Sigma=\\{0,1\\}$ and to say that such an object is computable when such a
representation can be produced using a Turing machine. The computable analysis
is designed to be able to also talk about objects such as real numbers,
functions over the reals, closed subsets, compacts subsets, …, which cannot be
represented by finite words over $\Sigma$ (a clear reason for it is that such
words are countable while the set $\mathbb{R}$, for example, is not). However,
they can be represented by some infinite words over $\Sigma$ and the idea is
to fix such representations for these various objects, called _names_ , with
suitable computable properties. In particular, in all the following proposed
representations, it was proved that an object is computable iff it has some
computable representation.

######  NB 7.

Here the notion of computability involved is one of Type 2 Turing machines,
that is to say, computability over possibly infinite words: the idea is that
such a machine has some read-only input tape(s), that contains the input(s),
which can correspond to either a finite or infinite word(s), a read-write
working tape and one (or several) write-only output tape(s). It evolves as a
classical Turing machine, the only difference being that we consider it
outputs an infinite word when it writes forever the symbols of that word on
its (or one of its) write-only infinite output tape(s): see [62] for details.

A name for a point ${\mathbf{x}}\in\mathbb{R}^{d}$ is a sequence $(I_{n})$ of
nested open rational balls with $I_{n+1}\subseteq I_{n}$ for all
$n\in\mathbb{N}$ and $\\{x\\}=\bigcap_{n\in\mathbb{N}}I_{n}$. Such a name can
be encoded as an infinite sequence of symbols.

We call a real function
$f\mathrel{\mathop{\ordinarycolon}}\subseteq\mathbb{R}\rightarrow\mathbb{R}$
computable, iff some (Type 2 Turing) machine maps any name of any
$x\in\operatorname{dom}(f)$ to a name of $f(x)$. For real functions
$\mathbf{f}\mathrel{\mathop{\ordinarycolon}}\subseteq\mathbb{R}^{n}\rightarrow\mathbb{R}$
we consider machines reading $n$ names in parallel. A computable function is
necessarily continuous: see [62] for all details.

We also need the concept of polynomial time computable function in computable
analysis: see [46]. In short, a quickly converging name of
${\mathbf{x}}\in\mathbb{R}^{d}$ is a name of ${\mathbf{x}}$, with $I_{n}$ of
radius $<2^{-n}$. A function
$\mathbf{f}\mathrel{\mathop{\ordinarycolon}}\mathbb{R}^{d}\to\mathbb{R}^{d^{\prime}}$
is said to be computable in polynomial time, if there is some oracle TM $M$,
such that, for all ${\mathbf{x}}$, given any fast converging name of
${\mathbf{x}}$ as an oracle, given $n$, $M$ produces some open rational ball
of radius $<2^{-n}$ containing $\mathbf{f}({\mathbf{x}})$, in a time
polynomial in $n$.

A function
$\mathbf{f}\mathrel{\mathop{\ordinarycolon}}\mathbb{R}^{d}\to\mathbb{R}^{d^{\prime}}$
is computable in polynomial space if there exists an oracla TM $M$, such that,
for all ${\mathbf{x}}$, given any fast converging name of ${\mathbf{x}}$ as an
oracle, given $n$, $M$ produces some open rational ball of radius $<2^{-n}$
containing $\mathbf{f}({\mathbf{x}})$, in a space polynomial in $n$.

###  7.2 Some classical statements from numerical analysis

######  Lemma 12 (Discrete Grönwall’s lemma, e.g [29, page 213]).

Consider sequences $h_{n},\theta_{n}\geq 0$ and $\varepsilon_{n}\in\mathbb{R}$
such that

$\theta_{n+1}\leq(1+\Lambda h_{n})\theta_{n}+|\epsilon_{n}|$

Then

$\theta_{n}\leq e^{\Lambda(t_{n}-t_{0})}\theta_{0}+\sum_{0\leq i\leq
n-1}e^{\Lambda(t_{n}-t_{i+1})}|\epsilon_{i}|$

###### Proof.

By recurrence over $n$. For $n=0$, the inequality is
$\theta_{0}\leq\theta_{0}$.

Suppose now the inequality at order $n$. Observe that

$(1+\Lambda h_{n})\leq e^{\Lambda(t_{n+1}-t_{n})}$

By hypothesis, we have

$\displaystyle\theta_{n+1}$ $\displaystyle\leq$ $\displaystyle
e^{\Lambda(t_{n+1}-t_{n})}\theta_{n}+|\epsilon_{n}|$ $\displaystyle\leq$
$\displaystyle e^{\Lambda(t_{n+1}-t_{0})}\theta_{0}+\sum_{0\leq i\leq
n-1}e^{\Lambda(t_{n+1}-t_{i+1})}|\epsilon_{i}|+|\epsilon_{n}|$

The inequality at order $n+1$ follows.

∎