Where the body will be after time T
A body moves around an attracting centre (a planet around the Sun, a station around the Earth), and nothing but that centre pulls it. The acceleration points to the centre and equals μ/r², where r is the distance and μ the centre’s gravitational parameter:
$$\mathbf{r}'' = -\frac{\mu}{|\mathbf{r}|^3}\,\mathbf{r}.$$Given the position and velocity at time 0, find the position at time T. Any units will do as long as they are consistent: kilometres and seconds, astronomical units and days.
The second example is the Earth on 17 October 1926 (relative to the Sun, from JPL Horizons data), a year later. The hidden tests include Halley’s Comet from the same day, the ISS on 17 October 2026 and the asteroid Apophis; all in the two-body problem, that is, without the pull of the other planets.
Input
The first line holds the number of queries q (1 ≤ q ≤ 20). Each of the following lines holds eight numbers: μ, x, y, z, vx, vy, vz, T. The orbit is an ellipse (the total energy is negative), the eccentricity is at most 0.996, and T is at most twenty orbital periods.
Output
For each query, three numbers: x, y, z at time T. Each is accepted with an error of up to 10⁻⁶ of the absolute value of the correct answer (for answers below one, with an absolute error of 10⁻⁶).
Examples
Example 1
2 1 1 0 0 0 1 0 3.1415926535897931 1 1 0 0 0 1 0 6.2831853071795862
-1 -2.953085909e-12 0 1 -1.129566406e-11 0
Example 2
1 0.00029591310798513925 0.91110800447119999 0.40365312087548 7.96530645313e-05 -0.0072492298972944994 0.015666005457084169 2.515955426257e-06 365.25
0.9110958776 0.4036793268 7.965727319e-05
Limits
- Time · Python
- 8000 ms
- Time · Go
- 2000 ms
- Memory
- 256 MB
- Judging
- 9 tests
Walkthrough
Open it step by step: each step gives away a little more. You can jump straight to the solution, but a problem solved with one hint stays with you longer.
1Hint 1
Rewrite the second-order equation as a first-order system: the state is six numbers, position and velocity, and the derivative of the state is velocity and acceleration. After that any method from the article will do. Which step should you choose?
2Hint 2
Look at the third test: Halley’s Comet crawls far from the Sun for decades and turns round at the Sun within a few weeks. A single step for the whole orbit is either too fine far out or too coarse at the Sun.
3The idea
A step that adjusts itself is needed: the Dormand–Prince pair estimates the error of each step from the difference between the fifth- and fourth-order answers and shortens the step where the solution changes fast. The tolerance must have a large margin: step errors accumulate, and the answer is checked to 10⁻⁶ after dozens of revolutions.
4Walkthrough
A Dormand–Prince step: seven slopes from the tableau, a fifth-order answer, and an error estimate, the difference of the two orders’ weights multiplied by the same slopes. Divide the error of each coordinate by a scale: the initial distance for positions and the initial speed for velocities. If the largest ratio is at most a tolerance of 10⁻¹², accept the step, and the seventh slope, taken at the end of the step, becomes the first slope of the next. Either way the new step is the old one times 0.9·(1/error)1/5, but by at most a factor of 5 up and 0.2 down at a time. Shorten the last step to land exactly on T.
Why RK4 with a constant step is not enough. Two thousand steps per query are generous for the Earth over a year, but for Halley’s Comet over 60 years they are steps of 11 days, and at the perihelion the comet turns by tens of degrees in 11 days, so the method skips the turn: the answer’s coordinates drift by tens of per cent. The ISS makes 16 revolutions a day, and two thousand steps is 125 per revolution: an error of about 2·10⁻⁴, hundreds of times more than allowed.
There is also a way with no steps at all: in the two-body problem the orbit is an ellipse, and the position at any moment comes from Kepler’s equation, solved by Newton’s method or bisection (see the article on bisection and Newton’s method). Such a solution is more accurate and faster, but works only for two bodies; the numerical one works for any forces.
4+Solution
5Complexity
The number of steps depends on the orbit: on an almost circular one, a few hundred per revolution at a tolerance of 10⁻¹²; on an elongated one, most of them go on the neighbourhood of the pericentre. Each step is six acceleration evaluations, O(1) each.
6Common mistakes
- A constant step. Passes the first two tests and fails on the comet, the station and random elongated orbits.
- Euler. Even a million steps do not give six digits: its error falls only as the first power of the step, and the orbit’s energy grows.
- A loose tolerance. A tolerance of 10⁻⁶ per step does not mean 10⁻⁶ in the answer: the errors of hundreds of steps add up and inflate. Take one several orders of magnitude stricter than the answer requires.
- The error scale. An absolute tolerance in kilometres and in astronomical units are different things. Normalise the error by the size of the orbit.