Average-Case Polynomial-Time Computability of Hamiltonian Dynamics

LIPICS - Leibniz International Proceedings in Informatics, Aug 2018

We apply average-case complexity theory to physical problems modeled by continuous-time dynamical systems. The computational complexity when simulating such systems for a bounded time-frame 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 circular-restricted three-body problem is average-case polynomial-time computable.

A PDF file should load here. If you do not see its contents the file may be temporarily unavailable at the journal website or you do not have a PDF plug-in installed and enabled in your browser.

Alternatively, you can download the file locally and open with any standalone PDF reader:


Average-Case Polynomial-Time Computability of Hamiltonian Dynamics

M F C S Average-Case Polynomial-Time 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 average-case complexity theory to physical problems modeled by continuous-time dynamical systems. The computational complexity when simulating such systems for a bounded time-frame 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 circular-restricted three-body problem is average-case polynomial-time 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), Core-to-Core 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 average-case 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 continuous-time 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 time-frame. The extended Church-Turing 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 (non-quantum) 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 time-frame 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 Church-Turing 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 Church-Turing 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 continuous-time 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 right-hand side function f of the equation y˙ = f (y) is polynomial-time computable and Lipschitz-continuous, the unique solution y can be computed in PSPACE and can be hard for this class [6, 4]. On the other hand, for analytic right-hand side function the solution is also a polynomial-time 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 worst-case 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 average-case complexity theory, which often provides a more significant measure of the performance of an algorithm than worst-case complexity when the hard instances are rare. However, finding the right notion of average-case complexity poses some subtle difficulties. A structural theory of average-case complexity for discrete problems was introduced by Levin [13]. Schröder, Steinberg and Ziegler recently extended Levin’s definition of average-case 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 right-hand 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 right-hand 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 phase-space. Finally, we apply our theorem to show that a special case of the three-body problem, the planar circular restricted three-body 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 so-called oracle tape, is replaced by the value of ϕ(w) in a single time-step. 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 (worst-case) 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 (worst-case) time-complexity T : N → N if TM (x, n) ≤ T (n) for all x ∈ dom f and that f is polynomial-time computable if T is a polynomial. Thus, the time-complexity 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 time-bound 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 Average-case complexity for real functions Worst-case 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. Average-case 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 worst-case complexity we would like to have a class of average-case polynomial-time problems that we consider as “efficient on average”. This notion should fulfill some basic robustness criteria. For example, composition of an average-case polynomial-time computable function with a worst-case polynomial-time 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 polynomial-time on average if there is some ε > 0 such that the function If an algorithm runs in polynomial-time on average, there is a high probability that it runs in polynomial-time 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 non-empty 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 continuous-time 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 first-order 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 multi-index notation for tuples of non-negative 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) non-open 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 right-hand side function f is analytic. By the Cauchy-Kowalevski 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 average-case 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 C-polynomial-time 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 short-hand notation “f is α−1-polynomial-time computable on Dα”. Formally, this means that f is C-polynomial-time 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 right-hand 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 polynomial-time 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 C-polynomial-time 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 C-polynomial-time computable, and C(y0) is a derivative bound for f on K(y0) as in (3). Then there exists a C-polynomial-time 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 Average-case 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 average-case complexity. Such an average-case 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 phase-space 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 degree-of-freedom (time-independent) Hamiltonian system is a 2d-dimensional 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, real-valued 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 time-independent systems where the Hamiltonian is constant over time and therefore is sometimes also called the energy of the system. 4.1 Average-case 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 α−1-polynomial-time computable. Then 1. Y is α−1-polynomial-time 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 d-dimensional 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 phase-space 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 phase-space. Saari uses a similar idea to show that for the three-body 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 phase-space volume is in fact the only attribute we use about Hamiltonian systems, we can formulate the theorem in terms of the bigger class of volume-preserving (sometimes also called conservative) dynamical systems. I Theorem 7. Let y˙ = f (y) be a volume-preserving 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 α−1-polynomial-time 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 polynomial-time 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 volume-preserving 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 non-trivial 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 volume-preserving it holds λ(Btk (2α)) ≤ λ(NA(2α)) and thus the statement follows. J 4.2 Average-case complexity for the restricted three-body problem As an application of theorem 7 we show that the solution of the restricted three-body problem can be computed in polynomial-time on average. The classical three-body 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 two-body 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 three-body 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 three-body 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 non-modified 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 phase-space 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 phase-space 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 + |q|2 + 2 α on the velocity. Now assume kqk ≥ 2 then both d1 ≥ 1 and d2 ≥ 2 and thus kvk2 ≤ 2h + |q|2 + 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 right-hand 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 polynomial-time computable on average: I Theorem 11. Simulating the restricted three-body 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 average-case complexity in analysis to problems in classical physics. We gave some general conditions which show that a continuous-time 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 phase-space usually suffices. We applied our theory to the planar circular restricted three-body 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 three-body 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 Chi-Chih Yao. Classical physics and the Church–Turing Thesis. Journal of the ACM (JACM), 50(1):100–105, 2003. Lei Zhao. Quasi-Periodic Almost-Collision Orbits in the Spatial Three-Body 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 C-polynomial-time 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 time-interval. The polynomial-time 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 k-th 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 right-hand 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 derivative-bound 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 Lipschitz-continuous on K with Lipschitz-constant 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 while-loop 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 2012-Proceedings 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 . Ker-I Ko . On the computational complexity of ordinary differential equations . Information and Control , 58 ( 1-3 ): 157 - 194 , 1983 . Birkhäuser Boston Inc., Boston, MA, 1991 . doi: 10 .1007/978-1- 4684 -6802-1. Ker-I 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 . Ker-I 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 three-body 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 N-body dynamical systems , pages 11 - 33 . Springer, 1991 . R. H. Miller . Numerical difficulties with the gravitational n-body 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.uni-trier.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 four-body 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 . Average-Case Bit-Complexity Theory of Real Functions . In Mathematical Aspects of Computer and Information Sciences , pages 505 - 519 . Springer, Cham, 2015 . doi: 10 .1007/978-3- 319 -32859-1_ 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 .

This is a preview of a remote PDF: http://drops.dagstuhl.de/opus/volltexte/2018/9612/pdf/LIPIcs-MFCS-2018-30.pdf

Akitoshi Kawamura, Holger Thies, Martin Ziegler. Average-Case Polynomial-Time Computability of Hamiltonian Dynamics, LIPICS - Leibniz International Proceedings in Informatics, 2018, 30:1-30:17, DOI: 10.4230/LIPIcs.MFCS.2018.30