The computational complexity of the initial value problem for the three body problem

TL;DR

This paper proves that solving the initial value problem for the three body problem cannot be done in polynomial time, using analysis of complex oscillatory solutions in the Sitnikov problem to demonstrate computational intractability.

Contribution

It establishes the non-polynomial complexity of the three body problem's initial value problem through rigorous analysis of oscillatory solutions.

Findings

The IVP for the three body problem is not solvable in polynomial time.

Oscillatory solutions in the Sitnikov problem exhibit complex behavior.

Polynomial-time algorithms for the three body problem are unlikely to exist.

Abstract

The paper is concerned with the computational complexity of the initial value problem (IVP) for a system of ordinary dynamical equations. Formal problem statement is given, containing a Turing machine with an oracle for getting the initial values as real numbers. It is proven that the computational complexity of the IVP for the three body problem is not bounded by a polynomial. The proof is based on the analysis of oscillatory solutions of the Sitnikov problem that have complex dynamical behavior. These solutions contradict the existence of an algorithm that solves the IVP in polynomial time.

Equations

$$\begin{array}[]{lcl}\dot{\mathbf{x}}&=&\mathbf{f}(\mathbf{x})\\ \mathbf{x}(0)&=&\mathbf{x}_{0}\end{array}$$

$$\left.\begin{array}[]{rcl}\dot{\mathbf{p}}_{i}&=&\mathbf{v}_{i},\quad i=1..N\\ \dot{\mathbf{v}}_{i}&=&\sum\limits_{\begin{subarray}{c}j=1\\ j\neq i\end{subarray}}^{N}\mu_{j}\frac{\mathbf{p}_{j}-\mathbf{p}_{i}}{|\mathbf{p}_{j}-\mathbf{p}_{i}|^{3}},\quad i=1..N\end{array}\quad\right\}$$

$$\ddot{z}=-\frac{2\mu z}{{\sqrt{z^2+r(t{)}^2}}^3},$$

$$\begin{array}[]{rcl}r(t)&=&a(1-e\cos E(t))\\ E(t)-e\sin E(t)&=&\sqrt{\frac{\mu}{a^{3}}}(t-t_{0})\end{array}$$

$$\begin{array}[]{rcl}\dot{\mathbf{x}}&=&\mathbf{f}(\mathbf{x})=(0,0,0,v,\ddot{z},\dot{E})\\ \ddot{z}&=&-\frac{2\mu z}{\sqrt{z^{2}+a^{2}(1-e\cos E)^{2}}^{3}}\\ \dot{E}&=&\frac{\sqrt{\mu a}}{1-e\cos E}\end{array}$$

$$\partial v/\partial v$$

$$\partial \ddot{z}/\partial z$$

$$\partial \ddot{z}/\partial E$$

$$\partial \dot{E}/\partial E$$

$$\lfloor \frac{{\tau }_{k+1}-{\tau }_k}{P}\rfloor =s_k,\forall k\in \mathbb{Z}.$$

$$\ddot{x}=-q(t)x$$

$$\ddot{x}=-Q(t)x,$$

$$h=H(\tau )=\sqrt{{\left(\frac{2\mu {\tau }^2}{{\pi }^2}\right)}^{\frac{2}{3}}-a^2}$$

$$\ddot{z}=-\frac{2\mu z}{{\sqrt{{z^{\ast }(t)}^2+r(t{)}^2}}^3},$$

$$\ddot{z}=-Q(t)z.$$

$$Q(t)>\frac{2\mu }{{\sqrt{h^2+a^2}}^3}$$

$$q=2\mu /{\sqrt{h^2+a^2}}^3,$$

$$\ddot{z}=-qz.$$

$$z^{\ast \ast }(t)=v_0\sin (\sqrt{q}t)$$

Full text

The final publication is available at Springer via

http://doi.org/10.1007/s10958-017-3407-3

1 Introduction

The problem of numerical integration of ODE systems is undoubtedly one of the most popular problems in applied mathematics. There exists a huge number of algorithms and program packages for obtaining numerical solutions of systems of differential equations originating from math, physics, celestial mechanics and engineering. However, there is little available research in the area of computational complexity of the initial value problem itself (some results are obtained in [1] and [2]). In most other works, complexity of particular algorithms is analyzed, in terms of either the number of basic arithmetical operations performed on each step, or the number of calls to first or higher-order derivatives.

In this work, a formal statement is presented of the IVP for a system of ODEs. In that statement, the input data for a problem will be: the initial conditions, the point $t$ in time, and a precision $\varepsilon$. An algorithm is supposed to consume the input and produce the output (an approximate state of the system at time $t$) that matches the actual state of the system up to $\varepsilon$. It must be noted that since the said statement includes real numbers, we can not work with just data of finite length. While it is sufficient to treat $t$ and $\varepsilon$ as rationals, initial conditions are a different story: there is no prior knowledge of how many digits in them will be sufficient to ensure that the solution at $t>0$ will be obtained with precision $\varepsilon$.

There are several approaches to work around that difficulty. The first is to consider an infinite input tape (or several infinite tapes) whose cells contain the digits of the initial conditions. A Turing machine for the given IVP can read the digits on demand. The second approach is to have a Turing machine use an oracle that gives the needed digits on demand. The third approach is to use a secondary Turing machine that prints out the digits into the tape by request from the main Turing machine.

The third approach, as opposed to the first two, is that it would limit us to just the constructive real numbers. In this work, the second approach (with an oracle) is used. It differs from the first one in the conventions of complexity analysis: the calls to an oracle account for the time complexity of algorithm as a function from the (finite) input length, while in the first approach makes the input infinite, rendering the complexity analysis difficult.

Another obstacle in the formal statement of the problem is the following: even if the derivatives in the system of ODEs are known Lipschitz-continuous functions, the problem of existence of the ODE solution at point $t$ can be undecidable.

We will show that even in the provably decidable case, the complexity of the IVP can be non-polynomial. The proof is based on the investigation of systems with complex dynamical behavior. As a basic example, we will use the classical Sitnikov problem for a three-body gravitational system, where two bodies follow elliptic orbits on a plane, and the third body stays on the line perpendicular to that plane. In the general case, the third body does unending oscillations with arbitrary amplitudes.

Instead of the oscillating solution of Sitnikov problem, we could use other dynamical systems exhibiting complex behavior, like, for instance, a neighborhood of some homoclinic solution. Their computational complexity would have turned non-polynomial, too.

In this work, we do not use a natural representation of such solution in the terms of symbolic dynamics [6]. Rather, to prove the absence of a polynomial algorithm for our formal IVP statement, it is sufficient to show that for a certain neighborhood of initial conditions in phase space, the number of algorithmically distinguishable trajectories is exponential in $t$.

2 Turing machine for the initial value problem

We estimate the computational complexity of the initial value problem for the dynamical system

[TABLE]

where $\mathbf{x}\in D$, $\mathbf{x}_{0}\in D$ is a real vector, and $\mathbf{f}:D\to\mathbb{R}^{n}$ is a computable real vector-valued function (open set $D\subseteq\mathbb{R}^{n}$ is the phase space of the system).

This work deals with the case when the solution $\mathbf{x}^{*}(t):\mathbb{R}\rightarrow D$:

exists on the whole $\mathbb{R}$; 2. 2.

is unique; 3. 3.

is a computable real vector-valued function.

Solutions that do not extend to $\mathbb{R}$ are called singular. The problem of determining the singularity of a solution is undecidable (see section 3.2). Uniqueness of a solution, it it exists, is guaranteed given that the function $\mathbf{f}$ is locally Lipschitz-continuous in every point in $D$. (The proof of that fact can be found e.g. in [7, p. 15].) However, the local Lipschitz-continuity does not imply the existence of the solution on $\mathbb{R}$.

If $\mathbf{f}$ is defined on $D$ when $D=\mathbb{R}^{n}$ and is (globally) Lipschitz-continuous, then the solution on $\mathbb{R}$ does exist for all $\mathbf{x}_{0}$ and is unique due to the Cauchy-Lipschitz theorem.

If $\mathbf{f}$ is continuous at every point in $D$, then every unique solution is computable by a (non-practical) combinatorial algorithm [8]. In particular, that holds for any computable $\mathbf{f}$, since every computable function is continuous.

In [9], it is proven that the solution of an IVP is computable with a modification of Picard–Lindelöf method, if $\mathbf{f}$ is Lipschitz-continuous on $D$. This important fact is quite non-trivial, despite the existence of hundreds numerical integrators for ODE. The vast majority of these integrators suffer from saturation: the step size being small enough, the error grows upon further decrease of the step size. Therefore, these integrators can not in principle obtain a solution up to an arbitrary precision [10].

To summarize: with $D=\mathbb{R}^{n}$ and Lipschitz-continuous $\mathbf{f}$, the solution of (1) with any $\mathbf{x}_{0}$ exists on $\mathbb{R}$, is unique and computable. It follows independently from [8] and [9]. In both sources, the computability is proven for the solution being the function of $\mathbf{x}_{0}$ and $t$, rather than just $t$.

In this work, we limit ourselves with the study of a particular instance of the three-body problem (see Section 3.3). The subject for study is the asymptotic dependence of the computational complexity of the solution $\mathbf{x}^{*}(t)$ on the value of $t$; the dependence on the precision of $t$ is not considered. In the text that follows, $t$ in the IVP is treated as rational, while $\mathbf{x}_{0}$ is a real vector. The complexity analysis of another special case of IVP, where $t\in\mathbb{R}$, is given in [1].

Definition 1**.**

The solution function of an initial value problem (1) is the function $S(\mathbf{x}_{0},t):D\times\mathbb{Q}\to D$, where $S|_{\mathbf{x}=\mathbf{x}_{0}}:\mathbb{Q}\to D$ is a computable real vector-valued function, whose closure on the real axis is the solution of (1).

Definition 2**.**

Turing machine that computes the solution function of an IVP is a Turing machine that accepts rational $t$ and $\varepsilon$ as input; has an oracle $\varphi$ that instruments $\mathbf{x}_{0}$ as a computable real vector; and produces the value of the solution $\mathbf{x}(t)$ corresponding to given $\mathbf{x}_{0}$ and $t$, with the precision $\varepsilon$.

It should be noted that in terms of complexity theory, the IVP belongs to the class of function problems, as opposed to more studied decision problems. The job of the oracle in the Turing machine is to write into its tape the representation of $\mathbf{x}_{0}$ up to an arbitrary precision, specified by the machine itself. It is obvious that the time required by the Turing machine includes the time to read the oracle tape.

Definition 3**.**

The IVP (1) has polynomial complexity if there exists a Turing machine from the definition 2 that computes its solution function in time bounded by $\mathcal{P}(\textrm{LENGTH}(t),\textrm{LENGTH}(\varepsilon))$, where $\mathcal{P}$ is an arbitrary polynomial.

Remark. Without loss of generality, it can be assumed that $\varepsilon=2^{-l}$, hence $\textrm{LENGTH}(\varepsilon)=l$.

Definition 4**.**

Suppose A and B are two IVPs. A is called polynomially reducible to B if there exist the following functions, computable in polynomial time: $G:D^{(A)}\to D^{(B)}$ and $H:D^{(B)}\to D^{(A)}$, so that for any initial state $\mathbf{x}_{0}^{(A)}\in D^{(A)}$ and a corresponding solution $\mathbf{x}^{*(A)}(t)$ the following holds: $\mathbf{x}^{*(A)}(t)=H(\mathbf{x}^{*(B)}(t))$, where $\mathbf{x}^{*(B)}(t)$ is a solution of B with initial state $\mathbf{x}_{0}^{(B)}=G(\mathbf{x}_{0}^{(A)})$.

Statement. If IVP A is polynomially reducible to IVP B, and B has polynomial complexity, then A has polynomial complexity as well.

3 Analysis of the computational complexity of the IVP for the three-body problem

3.1 $N$-body problem

Gravitational $N$-body problem is concerned with the Newtonian motion of $N$ point-masses in three dimensions. The system of ODEs for this problem is the following::

[TABLE]

where $\mu_{i}\in\mathbb{R}$, $\mu_{i}\geq 0$, $\mathbf{p}_{i}\in\mathbb{R}^{3}$, $\mathbf{v}_{i}\in\mathbb{R}^{3}$.

With $N=3$, the initial state of the system is given by a 21-vector $\mathbf{x}_{0}=(\mu_{1},\mu_{2},\mu_{3},p_{1,1},\ldots,p_{3,3},v_{1,1},\ldots,v_{3,3})$, while the system (2) defines a computable real vector-valued function $\dot{\mathbf{x}}={\mathbf{f}}(\mathbf{x})$. (The first three variables do not depend on $\mathbf{x}$ or $t$.)

3.2 Known results

The classical two-body problem ($N=2$) has a solution in algebraic functions of initial state and $t$. Depending on the configuration of the system, the two bodies follow either a Keplerian orbit (a parabola, hyperbola, or ellipse) or move along a line. The detailed description of the solutions can be found in multiple sources. Given those algebraic solutions, it is not difficult to show that the IVP for a nonsingular two-body problem has polynomial complexity.

With $N=3$ the problem does not have a generic algebraic solution, as proven by Poincaré. However, Sundman in 1912 derived a solution in the form of converging series. Unfortunately, the estimate of the number of terms required to calculate the series at point $t$ with a sensible precision is exponential in $t$ [11]. Merman improved Sundman’s result and found other series [12], though still exponential in $t$.

In practical tasks related to the $N$-body problem (in particular, in ephemeris astronomy) algorithms of numerical integration are used to obtain approximate solutions. The time complexity of such algorithms has a fundamental lower bound of $O(t)$, hence it can not be upper-bounded by a polynomial of $\textrm{LENGTH}(t)$.

The bottom line is that the known algorithms for the IVP for the three-body problem are non-polynomial. However, that does not disprove the polynomial complexity of the problem.

On a different note, let us show that there is a singular solution of the $N$-body problem that has a nonsingular one in any neighborhood. Let $N=2$. Two bodies collide if they are thrown upon each other along a straight line, while a smallest deviation from the straight line will prevent the collision (if the velocity is big enough). This implies the undecidability of the problem of determination of singularity with computable real $\mathbf{x}_{0}$: it requires the solution of equality relation of real numbers which does not exist.

3.3 Sitnikov problem

From now on, we will focus on a special case of the three-body problem, where two of the bodies are of equal positive mass, while the third body is massless and lies on a line, perpendicular to the plane of the motion of the first two bodies and passing through their center of mass (Fig. 1). Hence, the two bodies follow the unperturbed (Keplerian) orbit; in this problem, the elliptic orbit is the case.

Let us place the center of mass at the origin, and the $Z$ axis along the line where the third body is. Let us denote $r(t)$ the distance from the first body (and the second, as their trajectories are symmetric) to the origin.

Following Newtonian laws (2), the coordinate of the third body, denoted as $z$, obeys the following differential equation:

[TABLE]

where $\mu$ is the gravitational constant of the first and second bodies. Periodic function $r(t)$ comes from the solution of the two-body problem:

[TABLE]

$a$ (semimajor axis), $e$ (eccentricity) and $t_{0}$ (epoch) are constants that can be calculated from the initial state of the two bodies. $E(t)$ is the eccentric anomaly angle. The period of $r(t)$ is $P=2\pi\sqrt{\frac{a^{3}}{\mu}}$.

The initial values in the Sitnikov problem are:

• •

$a>0$, $e\in(0..1)$, $\mu>0$ — parameters of the orbit of the two bodies;

• •

$z_{0}=z(0)$ — initial position of the third body in the $Z$ axis.

• •

$v_{0}=\dot{z}(0)$ – initial velocity of the third body in the $Z$ axis.

• •

$\phi=E(0)$, $0\leq\phi<2\pi$ — initial value of the eccentric anomaly of the orbit of the two bodies.

The state vector of the system is accordingly $\mathbf{x}=(a,e,\mu,z,v,E)$. $a$, $e$ and $\mu$ do not depend on time; $\dot{z}=v$; $\dot{v}=\ddot{z}$ from (3); $\dot{E}$ follows from (4):

[TABLE]

Statement. IVP for the Sitnikov problem (3) is polynomially reducible to the IVP for the three-body problem (2).

The study of the trajectories of $z(t)$ in this system was started by Kolmogorov, while Sitnikov was the first to prove the existence of the oscillatory motions in this system [13]. His proof was also the first proof of this kind for three-body systems in general.

Theorem 1**.**

In the Sitnikov problem, there are no singularities, and the function $\mathbf{f}$ is Lipschitz-continuous on the whole domain.

Proof.

From Eqs. (3) and (4), along with the fact that $r(t)>0$, instantly follows that $\mathbf{f}$ is defined and continuous with any $z,v,E\in\mathbb{R}$.

Let us prove the Lipschitz-continuity of $\mathbf{f}$ by showing that all its partial derivatives w.r.t. $\mathbf{x}$ are bounded. We write down those derivatives, skipping the zero ones:

[TABLE]

(Notion $w=\sqrt{z^{2}+a^{2}(1-e\cos E)^{2}}$ is used for brevity.)

It is evident that all those functions are defined and continuous for any $z,v,E\in\mathbb{R}$ (for (9) it is important that $0<e<1$). The boundedness of (6) and (9) is trivial. The boundedness of (7) follows from the fact that it approaches zero as ${z\to\pm\infty}$: $\frac{1}{w^{3}}\to 0$ and $\frac{z^{2}}{w^{5}}\to 0$. Similarly, (8) is bounded because $\frac{z}{w^{3}}\to 0$ at ${z\to\pm\infty}$. ∎

Existence, uniqueness, and computability of the solution of the IVP for the Sitnikov problem follow from Theorem 1 and the references given in Section 2.

For the rest of the article, we consider the Sitnikov problem with $z_{0}=0$, omitting the solutions where the third body never crosses the plane.

3.4 Combinatorial properties of the solutions of the Sitnikov problem

Sitnikov’s result about the oscillatory motion was significantly extended by Alexeyev, who not only discovered the existence of all the classes of final motions in this problem, but also proved the following [3, 4, 5]:

Theorem 2**.**

For any sufficiently small eccentricity $e>0$ there exists an $m(e)$ such that for any double-infinite sequence $\{s_{n}\}_{n\in\mathbb{Z}},s_{n}\geq m$ there exists a solution $z(t)$ of the equation (3) whose roots satisfy the equation

[TABLE]

The shortened version of the original theorem is given, excluding the finite and semi-infinite sequences. Alexeyev also proved a generalization of his theorem to the case when the third body has a nonzero mass. A simpler proof was later obtained by Moser [14].

In what follows, we restrict our analysis to $t\geq 0,k\geq 0$ ($\tau_{0}=0$).

Lemma 1**.**

Let $C(T)$ be the set of (finite) sequences of the form $(s_{1},\ldots,s_{k})$, $s_{i}\geq m>1,\ s_{i}\mod 2=0,\ m\mod 2=0$, for each of which any sequence $(\tau_{0},\ldots,\tau_{k+1})$ satisfying (10) lies in the interval $[0,T]$ (i.e. $\tau_{k+1}\leq T$). $|C(T)|$ has an asymptotic lower bound exponential in $T$.

Proof.

Obviously, $C((m+1)P)=1$. For some $T\geq(m+1)P$, let us consider the interval $[T,T+(m+1)P]$. Any sequence $(s_{1},\ldots s_{k})\in C(T)$ can be extended to a sequence from $C(T+(m+1)P)$ by the following ways:

• •

$(s_{1},\ldots,s_{k},m)\in C(T+(m+1)P)$

• •

$(s_{1},\ldots,s_{k}+2i)\in C(T+(m+1)P),\ \forall 0<i\leq m/2$

Consequently, $|C(T+(m+1)P)|\geq(m/2+1)|C(T)|$, and that implies $|C(T)|\geq(m/2+1)^{\frac{T}{(m+1)P}}$ for sufficiently large $T$. If $m>0$, this bound is exponential in $T$. ∎

3.5 Computational complexity of the IVP for the Sitnikov problem

We give two lemmas that describe important properties of $z(t)$. The first lemma gives a lower bound of $|z(t)|$ between two roots separated by a certain distance. In the proof of the lemma, the Sturm’s comparison theorem is used:

Theorem 3** (Sturm’s comparison theorem).**

Consider two equations:

[TABLE]

and

[TABLE]

where $q$ and $Q$ are continuous functions. Let a nonzero solution of (11) $x(t)$ has roots $a$ and $b$, and $Q(t)>q(t)$ on $t\in[a,b]$. Then any solution of (12) has a root on $(a,b)$.

Lemma 2**.**

Let $z^{*}(t)$ be a solution of the Sitnikov problem (3) with initial values $a,e,\mu,\phi,v_{0}$. According to the previous assumptions, let $z^{*}(0)=0$. To be specific, we consider $v_{0}>0$ (the case of negative $v_{0}$ is a mirroring of that). Let $\tau$ be the smallest positive root of $z^{*}$. Then $\exists t\in(0,\tau):z^{*}(t)\geq h$, where

[TABLE]

Proof by contradiction.

Suppose $z^{*}(t)<h$, $0\leq t\leq\tau$. Since $z^{*}$ is the solution of (3), then it is also the solution of the following equation:

[TABLE]

where the factor of $z$ depends only on $t$, but not on $z$. Let us denote this factor $Q(t)$:

[TABLE]

Since $z^{*}(t)<h$ by the assumption, and $r(t)\leq a$, then

[TABLE]

Denoting

[TABLE]

we write a differential equation

[TABLE]

Since $q>0$ the equation (17) is the equation of a harmonic oscillator. We examine its solution $z^{**}$ for initial conditions $z(0)=0,\dot{z}(0)=v_{0}$:

[TABLE]

By the Sturm’s comparison theorem, between two roots of $z^{**}$—0 and $\pi/\sqrt{q}$—there exist roots of any solution of (15), including $z^{*}$. Since $\tau$ was chosen as the smallest positive root of $z^{*}$, it must be that $\tau<\pi/\sqrt{q}$. However, by construction of $q$ (16) and $h$ (13) it follows $\tau=\pi/\sqrt{q}$, hence the contradiction. ∎

Lemma 3**.**

Consider a nonnegative function $z(t)$, continuous and convex on $[t_{1},t_{2}]$; let $z(t_{1})=z(t_{2})=0$; let at some $t\in[t_{1},t_{2}]$ $z(t)>h>0$. Then $\exists t_{a},t_{b}\in[t_{1},t_{2}]:(t_{b}-t_{a})>\frac{3}{4}(t_{2}-t_{1}),\forall t\in(t_{a},t_{b})\ z(t)>h/4$.

Proof.

$z(t)$ has one (strict) maximum at $(t_{1},t_{2})$, let us say that $t_{3}$ is the point where the maximum is reached. Let us place points (Fig. 2): A$(t_{1},0)$, B$(t_{3},z(t_{3}))$, C$(t_{2},0)$. Let the line $z=h/4$ cross $AB$ at point $D$ and $BC$ at point $E$. Similarly, let the same line cross the $z(t)$ curve at $F$ and $G$.

Since $z(t)$ is convex, it lies above ABC, with the exception of A, B and C themselves (Fig. 2). Consequently, $|\textrm{DE}|<|\textrm{FG}|$. At the same time, from the similarity of triangles it follows that $\frac{|\textrm{DE}|}{|\textrm{AC}|}=1-\frac{h/4}{z(t_{3})}$. Since $z(t_{3})>h$ and $|\textrm{AC}|=(t_{2}-t_{1})$, we get $|\textrm{FG}|>\frac{3}{4}(t_{2}-t_{1})$. The horizontal coordinates of $F$ and $G$ are the desired $t_{a}$ and $t_{b}$. ∎

Theorem 4**.**

The time complexity of an initial value problem for the Sitnikov problem with any fixed value of eccentricity does not have a polynomial upper bound.

Proof by contradiction.

Suppose that there exists a Turing machine $M$ that calculates the solution function of the IVP for the Sitnikov problem in time $\mathcal{P}(\textrm{LENGTH}(t),\textrm{LENGTH}(\varepsilon))$, where $\mathcal{P}$ is arbitrary polynomial.

We examine the solutions at the interval $t\in[0,T],\ T\in\mathbb{N}$. From Lemma 1 and Alexeyev’s theorem, the number $C(T)$ of different solutions $z(t)$, forming different sequences $(s_{1},\ldots,s_{k})$ with $s_{k}\mod 2=0,s_{k}\geq m\ (m\mod 2=0)$, has a lower bound of $(m/2+1)^{\frac{T}{(m+1)P}}$, where $m$ depends only on $e$. (The Alexeyev’s theorem allows zero and odd $m$, but we can round the $m$ up to be a nonzero even number, without trouble to the theorem.

We build an algorithm for recovery of the sequence $(s_{1},\ldots,s_{k})$ that corresponds to a solution $z(t)$ for some initial values, using our supposedly existing Turing machine $M$. We choose the parameters $\delta\in\mathbb{Q},\delta<mP/2$ and $\varepsilon=2^{-l}(l\in\mathbb{N}),\varepsilon<h/4$, where $h=H(mP)$. (Note that $P$ is a computable real number.)

Let us build on $[0,T]$ a uniform grid with a step $\delta$; on each node $\{t_{i}=i\delta,0<i\leq\lfloor T/\delta\rfloor\}$ we can compute the state of the system up to the precision $\varepsilon$. The grid has the following important properties:

• •

If $|z(t_{i})|>h/4$, then from Lemma 3 follows that the closest root to $t_{i}$ lies no farther than $mP/4$.

• •

From above it follows that two neighbor nodes can not both have $|z|<h/4$

• •

Calculated $z(t_{i})$ can be divided into three classes: positive ($z>0$ for sure), negative ($z<0$ for sure) and undefined (the sign of $z$ is not determined within the given precision).

• •

Positive and negative nodes can go any number in a row, while there can be only one undefined node in a row.

• •

From the estimate of the distance between roots, it is evident that if there are no nodes between a positive node and a negative node, or if there is (one) undefined node, then $z$ has exactly one root in between.

Given that the $s_{k}$ are even, it is easily seen that $p$ nodes in a row of the same sign correspond to $s_{k}=\lceil(p+1)/2\rceil$; undefined nodes do not correspond to any $s_{k}$.

It is not important how long it took to recover the sequence of $s_{k}$. What matters is that all the ‘‘calls’’ to out Turing machine $M$ have used the same oracle for the computation of the (same) initial state. But, as we supposed, $M$ did not have a chance to read more than $\mathcal{P}(\textrm{LENGTH}(t_{i}),\textrm{LENGTH}(\varepsilon))$ digits from the oracle tape for any $t_{i}$, which is no more than $P(\log_{2}T,l)$; hence, basing on what it had read, it can possibly generate no more than $2^{\mathcal{P}(\log_{2}T,l)}$ different outcomes. At the same time, we proved that our algorithm recovers any of at least $(m/2+1)^{T/((m+1)P)}$ sequences, which (as $m>0$) is not bounded by the said polynomial. ∎

4 Conclusion and future work

In this work we examined the theoretical complexity of the initial value problem. We have shown that the lower time bound of that complexity can not be polynomial for the three-body problem (instantly meaning the absence of such a bound for the $N$-body problem). The choice of the three-body problem and oscillatory trajectories is not principal. We believe that similar results can be obtained in other systems, where, with the help of methods of symbolic dynamics, complex dynamical behavior can be shown and analyzed. We already mentioned homoclinic trajectories, discovered by Poincaré for the three body problem. It seems appropriate to quote his work ‘‘New methods of celestial mechanics’’ [15] here:

‘‘One is struck by the complexity of this figure I am not even attempting to draw. Nothing can give us a better idea of the complexity of the three-body problem and of all problems of dynamics where there is no holomorphic integral and Bohlin’s series diverge.’’

On a different note, for the integrable dynamical systems—those who have computable integrals of motion with good complexity bounds in $t$ and $\varepsilon$— it is possible to derive complexity bounds for the initial value problem in our formal statement. Those bounds will be polynomial by $\log(t)$ and $\log(1/\varepsilon)$. That can point to a link between computational complexity of the IVP and integrability.

On another different note, in this work the computational complexity of the IVP is examined at the ‘‘macro level’’ (rational $t\to\infty$), but what is left aside is the ‘‘micro level’’ (real $t$), where the precision of $t$ plays an important role [1]. Another work is planned devoted to that case.