AverageCase PolynomialTime Computability of Hamiltonian Dynamics
M F C S
AverageCase PolynomialTime Computability of Hamiltonian Dynamics
Martin Ziegler KAIST Daejeon 0 1
0 Holger Thies University of Tokyo Tokyo , Japan
1 Akitoshi Kawamura Kyushu University Fukuoka , Japan
We apply averagecase complexity theory to physical problems modeled by continuoustime dynamical systems. The computational complexity when simulating such systems for a bounded timeframe mainly stems from trajectories coming close to complex singularities of the system. We show that if for most initial values the trajectories do not come close to singularities the simulation can be done in polynomial time on average. For Hamiltonian systems we relate this to the volume of “almost singularities” in phase space and give some general criteria to show that a Hamiltonian system can be simulated efficiently on average. As an application we show that the planar circularrestricted threebody problem is averagecase polynomialtime computable. 2012 ACM Subject Classification Theory of computation → Computational complexity and cryptography Acknowledgements This work was supported by the Japan Society for the Promotion of Science (JSPS), CoretoCore Program (A. Advanced Research Networks), JSPS KAKENHI Grant Numbers JP18H03203 and JP18J10407, the Korean Ministry of Science and ICT grant NRF2016K1A3A7A03950702. We thank Florian Steinberg for discussions on averagecase complexity in analysis and the anonymous reviewers for many helpful comments.
and phrases Computable Analysis; Real computation; Dynamical systems; Averagecase complexity; Computation in physics

Many phenomena in nature can be modeled by continuoustime dynamical systems. Analyzing
such phenomena is usually done by simulating the evolution of a system with digital computers.
It is therefore crucial to better understand the computational properties of dynamical systems.
One of the most basic questions one can ask about such a system is given the state of the
system at some time t0 what will be the state at some time t > t0, i.e., to simulate the
evolution for a finite timeframe.
The extended ChurchTuring thesis is the statement that any physical computation device
can be simulated efficiently with a Turing machine. While the thesis might be incompatible
with quantum computation, at least for classical (nonquantum) physical computation the
assumption seems reasonable. Thus, it should be possible to simulate trajectories of dynamical
systems for problems in (classical) physics for a bounded timeframe efficiently as nature
already provides an efficient “computation device”. However, accurate numerical simulation
of dynamical systems can be very hard. For example, Miller pointed out the difficulties
when integrating the famous gravitational N body problem [16]. The N body problem is
the problem of predicting the motion of N point masses under their mutual gravitational
attraction. For a moderate number N of point masses, the main difficulty arises because
close encounters of particles cause instabilities. Indeed, two particles colliding leads to a
singularity in the analytic function describing the dynamics. On the other hand, Saari
could show that at least for N ≤ 4 singularities are rare in the sense that the set of initial
values leading to singularities has Lebesgue measure zero [
21, 22
] (for N > 4 this is an open
problem). Thus, a possible resolution to the inconsistency with the extended ChurchTuring
thesis might be that hard instances are extremely rare in nature and therefore it is possible
to do the simulation efficiently on typical inputs. For further discussion on the extended
ChurchTuring thesis in classical physics, see e.g. [26].
The theory of real number computation based on a realistic model of approximation
is known as computable analysis. Using this model, computational complexity theory
can be applied to real functions [7, 8]. Simulating a continuoustime dynamical system
corresponds to solving an initial value problem (IVP) for systems of ordinary differential
equations. There are already several results on the computational complexity of solving
IVPs in computable analysis. In particular, if the righthand side function f of the equation
y˙ = f (y) is polynomialtime computable and Lipschitzcontinuous, the unique solution y can
be computed in PSPACE and can be hard for this class [6, 4]. On the other hand, for analytic
righthand side function the solution is also a polynomialtime computable function [18, 9].
However, this formulation does not really capture the notion of what is usually understood
by simulating a dynamical system, as it is assumed that the solution exists on the whole
time interval and only takes values in a known compact set, and there are several hidden
factors depending on the function and the initial value that heavily influence the efficiency
in practice.
Instead of fixing the initial value, it is therefore more natural to study the complexity of
the function mapping initial values and time to the corresponding solution. This, however,
poses the problem that the worstcase complexity for most interesting systems is unbounded.
Indeed, it is quite obvious that the simulation should take longer the closer a trajectory
approaches a singularity of the system as then higher precision is required. Nonetheless, the
system might behave well for most initial values in the sense that the trajectory does stay
far away from any singularities. We would then expect efficient simulation to be possible on
typical inputs.
In this paper we want to formalize this intuition. In classical complexity theory, the
notion of being efficiently computable on typical inputs is usually captured by averagecase
complexity theory, which often provides a more significant measure of the performance of an
algorithm than worstcase complexity when the hard instances are rare. However, finding the
right notion of averagecase complexity poses some subtle difficulties. A structural theory
of averagecase complexity for discrete problems was introduced by Levin [13]. Schröder,
Steinberg and Ziegler recently extended Levin’s definition of averagecase complexity to
problems on real numbers [23].
In this paper we use their definition to show that many physical problems are indeed
efficiently solvable on average. We first give a short introduction to the model of computation
in Section 2 and define the most important notions that we need in the rest of the paper. In
Section 3 we formalize a parameterized result for IVPs with analytic righthand side. We
show that restricted to a compact domain where the dynamics take place the solution can be
computed in time polynomial in the output precision and a natural parameter depending on
the function and the domain. In Section 4 we use this to show that if the “probability of
trajectories to get close to complex singularities of the system” is small and if the righthand
side function can be evaluated efficiently on points not close to singularities, the simulation
can be done in polynomial time on average. We then focus on a special case of dynamical
systems that play an important role in classical physics, the Hamiltonian systems, and
show that there is a simple way to bound the above probability in terms of the volume of
singularities in phasespace. Finally, we apply our theorem to show that a special case of the
threebody problem, the planar circular restricted threebody problem, can be simulated in
polynomial time on average.
2
Computability and complexity in analysis
Computable Analysis extends classical computability theory to deal with uncountable
quantities such as real numbers. The theory already dates back to Turing [24] with later important
contributions for example by Grzegorczyk [3] and Lacombe [12]. The rigorous study of
complexity in this model was initiated by Ko and Friedman [8]. In this section we briefly
summarize the most important definitions. For a more detailed overview the reader is referred,
e.g., to Weihrauch’s monograph [25].
2.1
Computing real numbers and functions
Classically, computability is typically defined using Turing machines or an equivalent model
of computation. Turing machines compute functions from finite strings to finite strings, that
is, functions F : Σ∗ → Σ∗ for some fixed finite alphabet Σ. While computations over discrete
objects like natural numbers, rational numbers or graphs can be defined by choosing an
appropriate encoding for these objects as finite strings, the set of reals is uncountable and it
is therefore impossible to find such an encoding. On the other hand, any real number can
be approximated with arbitrarily small error by rational numbers. A real number x can
therefore be encoded by a function that gives arbitrary exact approximations of x. More
formally: A function ϕ : Σ∗ → Σ∗ is called a name for a real number x ∈ R if it maps strings
of length n to the binary encoding of some integer z such that x − 2zn ≤ 2−n. Any real
number has infinitely many different names. We say x is computable if it has a computable
name.
Similarly, a function f : R → R is computable if there is a computable function mapping
names of x to names of f (x), i.e., an algorithm that computes approximations of the output
f (x) with arbitrary precision while having access to arbitrary exact approximations of the
input x. Such computation on names can be formally defined using oracle machines. Oracle
machines can make queries to an oracle, a function ϕ : Σ∗ → Σ∗, during the computation.
When making such a query the string w on a special tape, the socalled oracle tape, is
replaced by the value of ϕ(w) in a single timestep. We denote the oracle machine M with
oracle ϕ by M ϕ. M ϕ again computes a function Σ∗ → Σ∗. A computable real function can
then be defined as follows: A partial function f :⊆ R → R is computable if there is an oracle
machine M such that whenever ϕ is a name for x ∈ dom f the machine M ϕ computes a name
for f (x). The machine can behave arbitrarily on elements outside the domain of the function.
The definition can be easily generalized to multidimensional functions f :⊆ Rn → Rm by
allowing multiple oracle and output tapes.
The running time TM (ϕ, w) for an oracle machine M with oracle ϕ : Σ∗ → Σ∗ and
input w ∈ Σ∗ is defined as the number of steps the machine makes. If the machine M
computes a real function f : R → R we can thus define the (worstcase) running time
TM (x, n) ∈ N ∪{∞} of M for an argument x ∈ R and output precision n ∈ N as an upper
bound of TM (ϕ, 0n) over all names ϕ for x. We further say that a computable function
f : R → R has (worstcase) timecomplexity T : N → N if TM (x, n) ≤ T (n) for all x ∈ dom f
and that f is polynomialtime computable if T is a polynomial. Thus, the timecomplexity
of a real function is given in terms of the number of steps necessary to approximate the
solution up to precision 2−n depending on the output precision n ∈ N.
Most work in real complexity theory restricts dom f to be a fixed compact set as this
makes sure that a (finite) complexity bound exists. A more general theory taking into account
the “size” of an oracle can be found in [5].
In this paper we sometimes use the notation that f is computable on some maximal
domain dom f and uniformly allows some timebound on a restricted subset of the domain.
By this we mean that there is a single algorithm to compute f on its domain and the running
time of this algorithm can be bounded on this subset.
2.2
Averagecase complexity for real functions
Worstcase complexity as defined above has the disadvantage that the cost can be dominated
by a small number of instances. A problem that can be solved efficiently on most inputs
might thus still be hard from the perspective of computational complexity. Averagecase
complexity studies the average running time of an algorithm over all inputs and often gives a
more realistic measure for the efficiency of an algorithm. Similarly to the class P in worstcase
complexity we would like to have a class of averagecase polynomialtime problems that we
consider as “efficient on average”. This notion should fulfill some basic robustness criteria.
For example, composition of an averagecase polynomialtime computable function with a
worstcase polynomialtime computable function should still be computable in polynomial
time on average. However, already in the discrete case this robustness property is not fulfilled
by the most intuitive definitions. A definition that resolves this problem was given by Levin
[13]. Schröder, Steinberg and Ziegler [23] recently extended Levin’s notion to real number
computations:
I Definition 1. Let (X, Σ, μ) denote a probability space and T : X ×N → [0, ∞] a measurable
function. We say that T is polynomialtime on average if there is some ε > 0 such that the
function
If an algorithm runs in polynomialtime on average, there is a high probability that it runs in
polynomialtime for a randomly selected input as by Markov’s inequality for any k > 0 it is
c
Pr[ T (x, n) ≥ kn 1ε ] ≤ kε .
In this paper we only consider the case where X ⊆ Rd, Σ is the Borel algebra over Rd and
μ(A) = λλ((XA)) where λ denotes the Lebesgue measure on Rd.
1 Z
n 7→ n X (T (x, n))ε dμ
is bounded by some constant c > 0.
(1)
Complexity of simulating dynamical systems
Dynamical systems are used to describe the evolution of a system over time. A dynamical
system on the reals can be formally defined as follows:
I Definition 2. A dynamical system is a triple (X, T, Φ) of a nonempty set X ⊂ Rd called
the phase space, a time set T ⊆ R and a (partial) function Φ :⊆ X × T → X called evolution
operator satisfying
1. Φ(x, 0) = x, and
2. Φ(Φ(x, t1), t2) = Φ(x, t1 + t2)
for x ∈ X and all t1, t2 ∈ T such that (x, t1), (Φ(x, t1), t2) ∈ dom Φ.
A point x ∈ X is also called a state of the system and Φ(x0, t) the state at time t (w.r.t. the
initial value x0). For a fixed initial value x0 the function Φ(x0, ·) is called trajectory through
x0. In this paper we only consider continuoustime systems where the time set is some real
interval and the evolution operator is continuously differentiable with respect to time. The
dynamics of such a system can be described by the solution of a system of autonomous
firstorder ODEs
y˙ = f (y), y(t0) = y0
(2)
for some function f : D ⊆ Rd → Rd and y0 ∈ Rd.
We call Equation (2) an initial value problem (IVP) with initial value y0. The
corresponding solution function y : R → Rd of the IVP gives a trajectory of the system. We
sometimes also write y(t; t0, y0) for the solution function with initial value y(t0) = y0 to make
the dependency on the initial value explicit.
In the following, we use multiindex notation for tuples of nonnegative integers β =
(β1, . . . , βd), i.e., for β ∈ Nd and x ∈ Rd it is β = β1 + · · · + βd, β! = β1! . . . βd!, xβ =
x1β1 . . . xdβd and Dβf = D1β1 . . . Ddβd f .
A function f : U → R for U ⊆ Rd open is called analytic if for each x0 ∈ U there is an
open neighborhood V ⊆ U of x0 such that for all x ∈ V the power series expansion
X
β∈Nd aβ(x − x0)β
converges to f (x). A function f : D → R on a (possibly) nonopen domain D is called
analytic if it can be extended to an analytic function on some open domain U ⊇ D.
We further say a function f : D → Rm, f (x) = (f1(x), . . . , fm(x)) is analytic if each of
the component functions f1, . . . , fd : D → R are analytic. Any function analytic on some
subset D ⊆ Rd can be extended to a complex analytic function on an open set G ⊆ Cd. In
this paper we only consider IVPs where the righthand side function f is analytic. By the
CauchyKowalevski theorem the solution function y of such an IVP is again analytic.
3.1
Parameterized complexity for analytic initial value problems
Most results on IVPs in computable analysis assume the initial value to be fixed. In this work,
however, we want to consider the average complexity over all initial values. We therefore first
need to formalize how exactly the complexity of the solution depends on the chosen initial
value. More formally, assume that we allow initial values from some set A ⊆ Rd. We want to
give a complexity bound for the solution function Y :⊆ A × [
0, 1
] → Rd, (y0, t) 7→ y(t; 0, y0)
that maps initial values at time 0 and time t ∈ [
0, 1
] to the solution of the initial value
problem at time t. Note that the restriction that the time is between 0 and 1 is somehow
arbitrary. However, for the averagecase analysis in the following chapter it is reasonable to
assume that the time is bounded as there is no natural choice for a distribution on time.
The domain of Y are tuples (y0, t) ∈ A × [
0, 1
] where the solution with initial value y0
exists up to time t. As this domain is not necessarily compact, we can not expect to find a
simple complexity bound on the whole domain. However, we can subdivide the domain into
subsets depending on some natural parameter of the system and give a complexity not only
in terms of the output precision n but also in terms of an additional integer parameter. That
is, we want to say that there is a uniform algorithm that computes the function on its whole
domain and if we restrict the domain to some subset this algorithm runs in time polynomial
in the output precision n and some parameter depending on the subset. We formalize this in
the following definition.
I Definition 3. Let f : D ⊆ Rd → Re be a computable real function and let D = Si∈I Di
be some (fixed) partition of its domain with index set I. We say f is Cpolynomialtime
computable for a function C : I → N if there is an oracle machine M that
1. computes f on all inputs x ∈ D, and
2. there is a polynomial p : N → N such that on inputs x ∈ Di the machine outputs a 2−n
approximation of f (x) after at most p(C(i) + n) steps.
We will usually partition the domain in terms of some positive real parameter α and use the
shorthand notation “f is α−1polynomialtime computable on Dα”. Formally, this means
that f is Cpolynomialtime computable w.r.t. the partition of the domain D = Sα>0 Dα
and the function C : R+ → N, C(α) = dα−1e.
Assume we are given an initial value problem (2) and the righthand side function f is
analytic on some compact K ⊆ dom f . Then there is some integer CK such that the bound
holds for all β ∈ Nd and x ∈ K [
10
]. On the other hand, if for some y0 ∈ dom f the solution
exists up to time T then y([0, T ]) is a compact subset of the domain. Thus, restricting f to
this set suffices to compute the solution.
We usually can not assume that f is polynomialtime computable on its whole domain
since its values can be unbounded. On the other hand on any compact subset K ⊆ dom f
of the domain, CK is an upper bound for the values of f . It is therefore a reasonable
assumption that f is Cpolynomialtime computable where C denotes a function that assigns
each compact subset K ⊆ dom f an integer CK as in Equation (3).
For simplicity let us from now assume that T = 1, i.e., we only consider initial values for
which the solution exists up to time 1. The following theorem characterizes the complexity
of the solution in terms of the parameterization given by C.
I Theorem 4. Assume f : D ⊆ Rd → Rd is analytic and computable. Let D0 ⊆ D be
some subset of initial values y0 ∈ D such that the solution y to the IVP (2) with initial
value y0 exists for all t ∈ [
0, 1
]. For each y0 ∈ D0 let further K(y0) := y([
0, 1
]; 0, y0) and
D0 := Sy0∈D0 K(y0). Assume there is a function C : D0 → N such that
f restricted to D0 is Cpolynomialtime computable, and
C(y0) is a derivative bound for f on K(y0) as in (3).
Then there exists a Cpolynomialtime computable function Z : Rd × [
0, 1
] × N → Rd such
that Z(y0, t, C(y0)) = Y (y0, t) for all y0 ∈ D0 and t ∈ [
0, 1
]. That is, Z computes Y when
given C(y0) as additional information.
Dβ f (x) ≤ CKβ+1β!
(3)
Note that in general it is not possible to effectively get the parameter C(y0) from the function
[2, 20]. Thus, in the above theorem we have to provide it as an additional input. For most
natural systems, however, it is usually easy to get an upper bound for the parameter. We
therefore assume that we can effectively get the parameter for the rest of the paper.
The main idea of the proof is that a simple power series based approach suffices to compute
a local solution on some small time interval [t0, t0 + δ] in time polynomial in n + C(y0)
(see e.g. [17]). To get a solution on a bigger interval this algorithm can be iterated several
times. To show theorem 4 it therefore suffices to show that polynomially in C(y0) many
iterations suffice and that it suffices to compute the intermediate values in each iteration
with polynomial precision. The proof of the theorem is very similar to the proof of a recent
result by Bournez, Graça and Pouly [
1
] for polynomial ODEs over unbounded time. We
therefore chose to omit the details here. However, as both the statement and notation of our
theorem are quite different from theirs, we included a proof in the appendix for completeness.
4
Averagecase complexity for dynamical systems
While theorem 4 gives a very general characterization for the complexity of simulating initial
value problems, it is usually not very useful as the complexity bound can differ for each
initial value y0. In general there is no way around this as the trajectories for distinct initial
values can be quite different. On the other hand, it might be the case that most trajectories
behave nicely in the sense that they stay far away from singularities of the system. In that
case, we would like to say that the simulation can “usually” be done efficiently. This notion
can be made formal using averagecase complexity.
Such an averagecase analysis, however, requires that we can get a bound on the probability
for a trajectory to come close to a singularity. Usually, there is no easy way to assign such a
probability. On the other hand, if the system has the special property that the volume of
subsets of phasespace is preserved over time then a small volume of singularities in phase
space indicates that the set of initial values coming close to singularities is also small. A
large class of dynamical systems that have this property are Hamiltonian systems.
A d degreeoffreedom (timeindependent) Hamiltonian system is a 2ddimensional
dynamical system y˙ = f (y) where the state can be split into two variables y(t) = (q(t), p(t)) for
smooth functions q, p : R → Rd that satisfy the system of 2d first order ordinary differential
equations
q˙(t) = ∂H(p, q)
∂p
,
p˙(t) = −
∂H(p, q)
∂q
(4)
for some smooth, realvalued function H : R2d → R. The functions q and p are called the
position and momentum of the system and H is called the Hamiltonian. We only consider
timeindependent systems where the Hamiltonian is constant over time and therefore is
sometimes also called the energy of the system.
4.1
Averagecase complexity for Hamiltonian system
In this chapter we define some criteria under which a given Hamiltonian system can be
simulated in polynomial time on average. Here, it is more natural to formalize the results
in terms of distance to complex singularities of a complex analytic extension instead of a
derivative bound. We first (independently of the system being Hamiltonian or not) formalize
the notion of a point being close to a (complex) singularity of an analytic function.
α
S
N (α)
N (2α)
Figure 1 An αsingularity is a state of the system where the distance to a singularity is at most
α. If the distance at time t = 0 is at least 2α and there is an αsingularity at time t > 0 there has
to be some time t0 when y(t0) is 2α close to the singularity. While the distance is between 2α and α
the velocity is bounded which can be used to bound the minimum time the particle has to spend in
the green region.
I Definition 5. Let f : D ⊆ Rd → Rd analytic. For any α > 0 and A ⊆ D we define the
following sets:
1. G(α) is the set of points x ∈ D such that there is a complex analytic extension fe of f
that is analytic on all z ∈ Cd with x − z ≤ α,
2. N (α) := D \ G(α),
3. RA ⊆ D is the set of points reachable from A in time t ≤ 1, i.e., x ∈ RA if there is y0 ∈ A
and t ∈ [
0, 1
] such that y(t; 0, y0) = x,
4. GA(α) := G(α) ∩ RA and NA(α) := N (α) ∩ RA,
5. A(α) := {y0 ∈ A : y(t; 0, y0) ∈ G(α) for all t ∈ [
0, 1
]},
6. B(α) := A \ A(α).
We call a point x ∈ N (α) an αsingularity, x ∈ A(α) an αgood initial value (αbad for
x ∈ B(α)) and points in RA reachable points.
If we restrict the initial values to be contained in a set A and the volume of B(α) is small
in comparison to A and if the growth of f behaves “nicely” on G(α) then the simulation can
be done in polynomial time on average:
I Theorem 6. Let f : D ⊆ Rd → Rd analytic and computable. Consider the solution
function Y to the IVP (2) restricted to a bounded subset of initial values A ⊆ D. Assume
that there is some polynomial m : N → N such that f˜ ≤ m(α−1) holds for any point in
GA( α2 ) and that f on GA(α) is α−1polynomialtime computable. Then
1. Y is α−1polynomialtime computable on initial values in A(α), and
2. If there is additionally some constant γ > 0 such that λ(λB(A(α))) ≤ αγ then Y can be
computed in polynomial time on average (w.r.t. the restriction of ddimensional Lebesgue
measure to A).
Proof. We show how to apply theorem 4 to prove the first part of the theorem. Let us first
show that the polynomial m can be used to get a derivative bound on G(α).
Assume f˜(z) ≤ M (α) for all z ∈ G( α2 ) for a function M with M (α) ≥ α2 . For z0 ∈ G(α),
consider the polydisc D = Qid=1 Di with Di := {z ∈ C : z0,i − z ≤ α2 }. Now, f˜ is analytic
on D by the definition of G(α). Furthermore f˜ is bounded on D by M (α). By Cauchy’s
integral formula it follows
Dβf (z0) =
Therefore, Dβf (z0) ≤ β!M (α)β+1 for all z0 ∈ G(α). As we assume that the set of initial
values is bounded and t ∈ [
0, 1
], it further follows that the absolute value of points in GA(α)
is bounded by a linear function in M (α). Thus we can apply theorem 4 using the partition of
the domain into sets Kα = GA(α) which concludes the proof of the first part of the theorem.
For the second part let μ be the restriction of the Lebesgue measure to A, i.e., μ(B) = λλ((BA))
for all measurable subsets B ⊆ A. Let further Ai = A(2−(i+1)) \ A(2−i), i.e., Ai contains
the initial values where the minimum distance to a complex singularity for any t ∈ [
0, 1
] is
between 2−i and 2−(i+1). Since Ai ⊆ B(2−i) it holds λ(Ai) ≤ λ(B(2−i)) and by assumption
μ(Ai) ≤ 2−iγ. On the other hand Ai is contained in A(2−(i+1)) thus by the first part of the
theorem and the assumption on the time bound of f it follows that Y is computable on
initial values in Ai in time u(n + 2i) for a polynomial u. Since A = Si∞=0 Ai it holds
Let m be the highest coefficient in the polynomial u then for ε ≤ 2γm the sum converges and
thus RA n1 T (x, n)εdμ is bounded. J
Theorem 6 holds even if the system is not Hamiltonian, but usually there is no simple way to
get a bound for the measure for the second part of the theorem. For Hamiltonian systems, on
the other hand, we can exploit the fact that they preserve phasespace volume over time (this
is known as Liouville’s theorem). Thus, there is a relation between the volume λ(A(α)) of
initial values leading to almost singularities and the volume λ(NA(α)) of almost singularities
in phasespace.
Saari uses a similar idea to show that for the threebody problem the set of initial
values leading to collisions has measure 0. However, for our application we need a stronger
quantification of how the measure correlates to the distance of the particles. The main
difficulty is that almost singularities can occur at any time t ∈ [
0, 1
]. Thus, theoretically
we would have to consider infinitely many “copies” of NA(α), one for each possible time.
Saari manages to replace this by only countably infinite many copies which suffices to show
that the set of initial values leading to collisions has measure zero. However, this approach
does does not suffice to give a bound for αsingularities with α > 0. We therefore need a
slightly more complicated idea to show that finitely many copies suffice in our case. Note,
however, that Saari’s result holds for unbounded time while our idea is based on the time
being bounded. In fact, a generalization of our result to unbounded time turns out to be
false [27].
Let us now state our main theorem. As the property of preserving phasespace volume is
in fact the only attribute we use about Hamiltonian systems, we can formulate the theorem in
terms of the bigger class of volumepreserving (sometimes also called conservative) dynamical
systems.
I Theorem 7. Let y˙ = f (y) be a volumepreserving dynamical system with f : D ⊆ Rd → Rd
computable and analytic and let A ⊆ D be a bounded subset of initial values. Assume that
f is α−1polynomialtime computable on αgood initial values y0 ∈ GA(α),
for all α > 0 the measure of αbad subsets of phase space is bounded by λ(NA(α)) ≤ αγ
for some γ > 0,
fe is bounded by Mα := M (dα−1e) for a polynomial M : N → N on N (2α) \ N (α).
Then
1. λ(B(α)) ≤ 2Mαα + 1 λ(NA(2α))
2. The solution function Y :⊆ A × [
0, 1
] → Rd is polynomialtime computable on average.
Proof. We only prove the first part. The second part then follows by applying theorem 6.
For any t ∈ [
0, 1
] let Bt(α) ⊆ A be the subset of initial values that lead to an αsingularity
at time t, i.e., y0 ∈ Bt(α) iff y(t; 0, y0) ∈ N (α). As the system is volumepreserving it holds
λ(Bt(α)) ≤ λ(NA(α)) for all t and α. Obviously, B(α) = St∈[
0,1
] Bt(α).
We now show that we can replace the infinite union by a union of finitely many slightly
bigger sets. Let y(0) ∈ Bt(α) for some t ∈ [
0, 1
], i.e., y(0) is an initial value such that
y(t) ∈ N (α). If y(0) ∈ N (2α) then it holds that y(0) ∈ B0(2α) by definition. Otherwise there
has to be some time t0 > 0 where y(t0) ∈ N (2α) \ N (α) for the first time (see Fig. 1). As for
any z ∈ N (2α) \ N (α) by assumption f˜(z) ≤ Mα. It follows that for any s ∈ [t0, t0 + Mαα ] the
difference of the state at time t0 and the state at time s can be at most α, i.e., y(s) ∈ G(α). In
particular, if y(0) ∈ G(2α) and for some t ∈ [
0, 1
] it holds that y(t) ∈ N (α) then there is some
nontrivial time interval [t1, t2] of length at least Mαα such that y(s) ∈ N (2α) \ N (α) for all
s ∈ [t1, t2]. Therefore, [t1, t2] contains some multiple t∗ of 2Mαα . In particular, y(t∗) ∈ N (2α)
and thus by definition y(0) ∈ Bt∗ (2α). This shows that for the finite subset {tk}k ⊂ [
0, 1
] of
multiples of 2Mα α it holds that B(α) ⊆ Stk Btk (2α). By applying the fact that the system is
volumepreserving it holds λ(Btk (2α)) ≤ λ(NA(2α)) and thus the statement follows. J
4.2
Averagecase complexity for the restricted threebody problem
As an application of theorem 7 we show that the solution of the restricted threebody problem
can be computed in polynomialtime on average. The classical threebody problem is the
problem of predicting the motion of three point masses under their mutual gravitational
attraction. An important special case is that we assume that one particle P3 has much
smaller mass than the other two P1 and P2. In that case P3 does not influence the motion
of P1 and P2 significantly. The problem can therefore be simplified by assuming that P3
is massless. The motion of the heavy particles can then be seen as a twobody problem.
A further simplification can be achieved by assuming that the heavy particles move on a
circular orbit around their common center of mass and P3 only moves in the plane defined
by P1 and P2. This problem is known as the planar circular restricted threebody problem. In
spite of being much simpler than the general problem, the restricted problem shares many of
the properties of the N body problem that makes it interesting for our analysis. In particular,
the restricted problem is a Hamiltonian dynamical system where the equation of motion is
given by an analytic initial value problem. There is no general closed form solution and the
motion can be chaotic [15]. The restricted threebody problem has been studied extensively
by mathematicians and engineers.
By choosing appropriate units of measurement and a rotating coordinate system, the
system can be brought in a simpler, normalized form only depending on a single parameter
μ ∈ (0, 0.5]. The masses of the particles P1 and P2 in the new units are given by μ and 1 − μ,
respectively. The position of the heavy particles in the rotating coordinate system remains
1 − μ
fixed at (−μ, 0) and (1 − μ, 0). We use q = (q1, q2) to describe the coordinates of P3 relative
to that coordinate system (see Figure 2). A full derivation for the transformations as well as
formulas to translate a solution for the normalized system to a solution for the nonmodified
system can e.g. be found in [11].
In Hamiltonian form the IVP can be written in terms of position q ∈ R2 and moment
p ∈ R2 as
H(p, q) = 21 kpk2 + q2p1 − q1p2 − dμ1 −
where d1 := p(q1 + μ)2 + q22 and d2 := p(q1 + μ − 1)2 + q22. This defines a dynamical
system with phase space Γ ⊆ R4. We sometimes also consider the velocity of the particle
given by v(t) = (p1(t) + q2(t), p2(t) − q1(t)) instead of the moment.
We assume that we start the simulation with initial values in the set A of points q0, p0
with q0 ≤ 1 and p0 ≤ 1. We first make the additional assumption that the energy is
bounded by some constant h > 0, i.e., it holds H(q0, p0) ≤ h. We denote the subset of
initial values satisfying this bound by Ah.
In this problem, an αsingularity corresponds to the situation where P3 gets close to either
P1 or P2. Note that the absolute value of the moments tend to infinity when approaching a
singularity. A bound on the volume of αsingularities in phasespace is given by the following
lemma.
I Lemma 8. For any α ∈ [0, 0.5], the Lebesgue measure of NAh (α) is given by λ(NAh (α)) ≤
8π2hα2.
Proof. We show that for the set N1(α) of points coming close to P1 it holds that λ(N1(α)) ≤
π2hα2. Then the claim follows as NAh (α) = N1(α) ∪ N2(α). We first change to a new
coordinate system (x, v) in terms of position and velocity of the particle such that P1 is at
the origin (note that this change preserves volume). The Hamiltonian of the problem in
these coordinates can be written as
1 1
E(x, v) = 2 kvk2 − 2 kxk2 + x1μ −
μ
kxk −
1 − μ − 0.5μ2
d2
Let Γ ⊆ R4 the phasespace of the problem. Γ can be parameterized in terms of the position
and the energy E by
Φ : (E, r, ϕ, ψ) 7→
with
r cos (ϕ), r sin (ϕ), R(E, r, ϕ) cos (ψ), R(E, r, ϕ) sin (ψ)
r
R(E, r, ϕ) :=
2E + μ2 − 2μr cos(ϕ) +
2μ + 2 − 2μ + r2
r d2
It is N1(α) ⊆ Φ(G) with G := [−h, h] × [0, α] × [0, 2π) × [0, 2π). Thus the volume can be
bounded by
Z
Φ(G)
Z
G
λ(N1(h, α)) ≤
dq1 dq2 dv1 dv2 =
 det(D Φ) dr dh dϕ dψ = 4π2 · α2 · h.
The last equality holds since det(D Φ) = r.
It further holds that a bound on the energy implies that both the position and the velocity of
the particle are bounded as long as the particle does not get close to one of the singularities.
I Lemma 9. For any q, v ∈ G(α) ∩ Ah it holds kqk ≤ 2h + 10 and kvk ≤ q(2h + 11)2 + α1 .
Proof. Assume d1 ≥ α and d2 ≥ α. The Hamiltonian in terms of the velocity is given by
H(v, q) = 21 kvk2 − dμ1 −
1
− 2 kqk2.
Thus the bound on the Hamiltonian implies the bound
kvk2 ≤ 2h + q2 +
2
α
on the velocity.
Now assume kqk ≥ 2 then both d1 ≥ 1 and d2 ≥ 2 and thus kvk2 ≤ 2h + q2 + 2 holds.
Let h0 = max(h, 2) then kvk ≤ h0 + k k
q . Now consider for n ≥ 1 the subsets Un where
n − 1 ≤ kqk ≤ n. For n ≥ 3, (q, v) ∈ Un implies that kvk ≤ h0 + n. In particular the minimum
time to get from Un to Un+1 is tn := h01+n . As for the initial values kqk ≤ 2 the minimum
time needed to reach a state where kqk ≥ N is bounded by
TN := XN 1
i=3 h0 + n
= HN+h0+3 − Hh0+3
where HN := PkN=1 k1 denotes the N th harmonic number. By the well known bound
γ + log(N ) ≤ HN ≤ γ + log(N + 1)
it holds
TN ≥ log
N + h0 + 3
h0 + 4
.
kvk ≤
r
2h + (2h + 10)2 +
2
α
from which the claim follows.
As the time is bounded by 1 it follows that kqk ≤ 2h0 + 6 ≤ 2h + 10. Inserting this back
into (8) yields
(6)
(7)
J
(8)
From this it can be easily followed that the righthand side function of the ODE system is
computable in time polynomial in dα−1e and h and that a polynomial bound in those terms
holds for complex analytic extensions of the function on N (2α) \ N (α).
As for initial values it holds that the magnitude of both position and moment are bounded
by 1, it follows from the Equation for the Hamiltonian (5) that high energy can only be due
to particles being already close at time 0. In particular for αgood initial values the following
bound on the Hamiltonian holds.
I Lemma 10. Assume kqk ≤ 1, kpk ≤ 1 and q, p ∈/ N (α) then H(p, q) ≤ 3 + α2 .
Thus the above polynomial bounds can all be stated only in terms of α. In particular all
conditions in theorem 7 are satisfied and thus the problem is polynomialtime computable
on average:
I Theorem 11. Simulating the restricted threebody problem for time t ≤ 1 for initial values
p0, q0 with p0 ≤ 1, q0 ≤ 1 is possible in polynomial time on average.
5
Conclusion
We applied a recent definition of averagecase complexity in analysis to problems in classical
physics. We gave some general conditions which show that a continuoustime system can be
computed efficiently on average. For the important special case of Hamiltonian systems, we
showed that a simpler approach based on the volume of almost singularities in phasespace
usually suffices. We applied our theory to the planar circular restricted threebody problem
and showed that it indeed can be computed in polynomial time on average. The same can
easily be done for some other simple dynamical systems.
However, our theory does not easily generalize to some other more complicated systems
like the classical N body problem as bounding the volume of singularities in phase space is
more complicated for these systems. While we think that most systems in nature can indeed
be simulated efficiently and that a similar result holds for, e.g., the general N body problem
it is unlikely that our approach can be easily adapted to this case as for N > 4 it is not even
known if the set of initial conditions leading to singularities has measure zero.
Nonetheless, we hope that we can at least extend our theory to the general threebody
problem and some other interesting problems in future work.
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
Andrew ChiChih Yao. Classical physics and the Church–Turing Thesis. Journal of the
ACM (JACM), 50(1):100–105, 2003.
Lei Zhao. QuasiPeriodic AlmostCollision Orbits in the Spatial ThreeBody Problem.
Communications on Pure and Applied Mathematics, 68(12):2144–2176, 2015.
Proof of Theorem 4
The proof consists of two parts. We first show that the local solution on some small interval
with length polynomial in C(y0)−2 can be computed in C(y0)polynomial time. We then
show that iterating the algorithm for the local solution until reaching time t = 1 can still be
done within this complexity bound.
A.1
Computing a local solution
We formulate the complexity result for the local solution as the following Lemma.
I Lemma 12. Let D ⊆ Rd analytic and computable. Then there is an algorithm that given
an initial y0 ∈ Rd, time t ∈ R and an integer C ∈ N such that C is a derivative bound for f
on the closed ball BC (y0) of radius C1 around y0 returns the solution y(t) to the IVP (2) for
1
any t with t − t0 ≤ 2(d+1)C2 . Assume f is Cpolynomialtime computable on BC (y0) then
so is y for all t that satisfy the above bound.
Proof. Cauchy’s existence theorem guarantees that the solution exists and is analytic on
this timeinterval. The polynomialtime result has in principle already been shown, e.g., by
Moiske and Müller [17]. The main idea is to first compute the coefficients of the power series
of f around t0. Müller [19] showed that with the above bounds for any β ∈ Nd the coefficient
aβ is computable in time polynomial in n + C + β. Further, if the power series of f around
y0 can be computed in time Tf (n + β) then the kth coefficient of the power series (bk)k∈N
of any of the d components y1, . . . , yd of solution function y around t0 can be computed from
the power series of f in time O(kdTf (n + kd)).
For any i = 1, . . . , d, the solution yi(t) for small enough t can then be approximated by
summing up the truncated power series PN
i=0 bi(t − t0)i for some N ∈ N. It can be shown
1
that bi ≤ (d + 1)iC2i. Therefore for any t with t − t0 ≤ 2(d+1)C2 it suffices to sum up
O(n) coefficients to make the truncation error less than 2−(n+1). To additionally make the
approximation error less than 2−(n+1) and thereby the total error less than 2−n it suffices to
approximate each bi with error less than 2−(2n+1). J
A.2
Extending to a global solution
We now prove the theorem by iterating the local solution algorithm. Let us first summarize
the assumptions:
1. The righthand side function f : D ⊆ Rd → Rd is analytic and computable,
2. The solution y(t) := y(t; 0, y0) exists for all time t ∈ [
0, 1
],
3. Restricted to the set K := y([
0, 1
]; 0, y0) the integer C is a derivativebound for f , i.e.,
Dβf (x) ≤ Cβ+1β! for all x ∈ K,
4. The algorithm computing f gives a 2−n approximation of f (x) in time poly(n + C) on
any x ∈ K.
Note that each x ∈ K is bounded by x ≤ C as it is in the image of y.
We want to show that iteratively using the local solution algorithm yields a
polynomialtime algorithm mapping initial values y0 ∈ K and time t ∈ [
0, 1
] to the solution y(t) at time t.
Algorithm 1 Solving Initial Value Problems.
function solve_ivp(f, y0, t, n, C)
tcurr ← 0
ycurr ← Approx(y0, m)
1
h ← 2(d+1)C2
while tcurr + h ≤ t do
ycurr ← LocalSolution(ycurr, h, m, C)
tcurr ← tcurr + h
return LocalSolution(f, ycurr, t − tcurr, m, C)
For this it suffices to show that we can always reach time 1 in polynomially many iterations
and that it suffices to approximate all intermediate values with polynomial precision. For
simplicity we fix the index i of the solution function yi and simply denote it by y. By
lemma 12 we can assume that there is an algorithm LocalSolution that computes for inputs
1
y0 ∈ K, t ∈ [
0, 1
] and C, n ∈ N with t ≤ 2(d+1)C2 a rational approximation q ∈ Q such
that y(t; t0, y0) − q ≤ 2−n. For a suitable choice of the intermediate precision parameter m
Algorithm 1 computes an approximation q to the solution at any time t ∈ [
0, 1
] such that
y(t; t0, y0) − q ≤ 2−n. It remains to show that m can be chosen to be polynomial in n + C.
The algorithm makes at most 2(d + 1)C2 steps to reach any time t ∈ [
0, 1
].
Let us first consider the error when computing the local solution yi+1 from yi. There
are two types of errors. The first one is due to the local solution algorithm only returning
an approximation to the real solution with precision 2−m. The second kind of error arises
because of the accumulated error in yi: Instead of solving the IVP problem with initial value
yi we use an approximation zi of yi as initial value. To bound the second kind of error we
use the following fact:
I Lemma 13. Assume f : Rd → Rd fulfills the conditions above. Then each component
function fi : Rd → R is Lipschitzcontinuous on K with Lipschitzconstant L = √dC2.
Proof. This simply follows from the fact that all partial derivatives of f are bounded by
C2. J
This can be used to get a bound on the local error.
I Lemma 14. Assume f fulfills the conditions above. Let y0, z0 ∈ K be initial condition
1
with ky0 − z0k∞ ≤ ε. Then for all t < 2(d+1)C2 it is y(t; y0) − y(t; z0) ≤ 2ε.
Proof. It holds that y(t; y0) = y0 + R0t f (y(τ )). Thus
y(t; y0) − y(t; z0)
Z t
≤ y0 − z0∞ +
0
Z t
t0
≤ ε +
L y(τ ; y0) − y(τ ; z0) dτ
y(t; y0) − y(t; z0) ≤ εeLt.
f (y(τ ; y0)) − f (y(τ ; z0)) dτ
By Grönwall’s Lemma [14, Chapter 1, Theorem 8.1] it follows that
1 √d
For t ≤ 2(d+1)C2 ≤ 2(d+1)L it then holds
J
It is now easy to compute the total error for the algorithm depending on the parameter m
for the intermediate precision.
I Lemma 15. Assume algorithm 1 needs N iterations to reach the final time t. The total
error E is bounded by E ≤ 2N+1−m.
Proof. Let Ek be the error of ycurr after k iterations of the whileloop in Algorithm 1. It is
E0 ≤ 2−m and Ek+1 ≤ 2Ek + 2−m. It follows by induction that EN ≤ 2N+1−m − 2−m. J
Since the maximum number of iterations is bounded by N = 2(d + 1)C2 to achieve precision
2−n it suffices to choose m ≥ n + 2(d + 1)C2 + 1 which proves the theorem.
1 2 3 4 Olivier Bournez , Daniel S. Graça, and Amaury Pouly . On the complexity of solving initial value problems . In ISSAC 2012Proceedings of the 37th International Symposium on Symbolic and Algebraic Computation , pages 115  121 . ACM, New York, 2012 .
doi:10.1145/2442829 .2442849.
Daniel S Graça , Ning Zhong, and Jorge Buescu . Computability, noncomputability and undecidability of maximal intervals of IVPs . Transactions of the American Mathematical Society , 361 ( 6 ): 2913  2927 , 2009 .
A. Grzegorczyk . On the definitions of computable real continuous functions . Fund . Math., 44 : 61  71 , 1957 .
Akitoshi Kawamura. Lipschitz Continuous Ordinary Differential Equations are PolynomialSpace Complete . Computational Complexity , 19 ( 2 ): 305  332 , 2010 .
Akitoshi Kawamura and Stephen Cook . Complexity theory for operators in analysis . ACM Transactions on Computation Theory , 4 ( 2 ):Article 5, 2012 .
KerI Ko . On the computational complexity of ordinary differential equations . Information and Control , 58 ( 13 ): 157  194 , 1983 .
Birkhäuser Boston Inc., Boston, MA, 1991 . doi: 10 .1007/9781 4684 68021.
KerI Ko and Harvey Friedman . Computational complexity of real functions . Theoretical Computer Science , 20 ( 3 ): 323  352 , 1982 . doi: 10 .1016/S0304 3975 ( 82 ) 80003  0 .
KerI Ko and Harvey Friedman . Computing power series in polynomial time . Advances in Applied Mathematics , 9 ( 1 ): 40  50 , 1988 .
Hikosaburo Komatsu . A characterization of real analytic functions . Proceedings of the Japan Academy , 36 ( 3 ): 90  93 , 1960 .
Wang Sang Koon , Martin W. Lo, Jerrold E. Marsden , and Shane D. Ross . Dynamical systems, the threebody problem and space mission design . World Scientific , 2000 .
Daniel Lacombe . Sur les possibilités d'extension de la notion de fonction récursive aux fonctions d'une ou plusieurs variables réelles . In Le raisonnement en mathématiques et en sciences expérimentales, Colloques Internationaux du Centre National de la Recherche Scientifique , LXX, pages 67  75 . Editions du Centre National de la Recherche Scientifique, Paris, 1958 .
Leonid A. Levin . Average case complete problems . SIAM Journal on Computing , 15 ( 1 ): 285  286 , 1986 .
X. Mao. Stochastic Differential Equations and Their Applications. Ellis Horwood series in mathematics and its applications . Horwood Pub., 1997 .
Andrea Milani . Chaos in the three body problem . In Predictability, stability, and chaos in Nbody dynamical systems , pages 11  33 . Springer, 1991 .
R. H. Miller . Numerical difficulties with the gravitational nbody problem . In Dale G. Bettis, editor, Proceedings of the Conference on the Numerical Solution of Ordinary Differential Equations , pages 260  275 , Berlin, Heidelberg, 1974 . Springer Berlin Heidelberg.
In Proc. 22 JAIIO  PANEL '93, Part 2 , pages 283  293 , Buenos Aires , 1993 . URL: http: //www.unitrier.de/mueller.
Norbert Th. Müller . Uniform Computational Complexity of Taylor Series. In Proc. 14th International Colloquium on Automata, Languages, and Programming , volume 267 of LNCS , pages 435  444 . Springer, 1987 .
Norbert Th. Müller . Polynomial time computation of Taylor series . Proc. 22 JAIIOPANEL'93 , Part 2 , Buenos Aires , 259 : 281 , 1993 .
Norbert Th. Müller . Constructive Aspects of Analytic Functions . In Proc. Workshop on Computability and Complexity in Analysis , volume 190 of InformatikBerichte , pages 105  114 . FernUniversität Hagen, 1995 .
Donald G. Saari . Improbability of Collisions in Newtonian Gravitational Systems. II. Transactions of the American Mathematical Society , 181 : 351  368 , 1973 . doi: 10 .2307/1996638.
Donald G. Saari . A global existence theorem for the fourbody problem of Newtonian mechanics . Journal of Differential Equations , 26 ( 1 ): 80  111 , 1977 . doi: 10 .1016/ 0022  0396 ( 77 ) 90100  0 .
Matthias Schröder , Florian Steinberg, and Martin Ziegler . AverageCase BitComplexity Theory of Real Functions . In Mathematical Aspects of Computer and Information Sciences , pages 505  519 . Springer, Cham, 2015 . doi: 10 .1007/9783 319 328591_ 43 .
Proceedings of the London Mathematical Society , 2 ( 1 ): 230  265 , 1936 . doi: 10 .1112/plms/ s2 42 .1.230.
Klaus Weihrauch . Computable Analysis . Springer, Berlin/Heidelberg, 2000 .