Where will the cat go? (orbital mechanics)

Question

A nearly massless cat is dropped in space (don't worry, with a space suit and everything) at the point (x, y, z) with velocity (vx, vy, vz). There is an fixed, infinitely dense planet (with volume of 0) at the point (0, 0, 0) and it attracts objects at distance r with acceleration 1/r^2. According to Newtonian gravity, where does the object go after time t?

Nearly massless in this case means you are outputting the value of lim (mass --> 0) <position of cat>. The mass is affected by the planet's gravity, but the planet is not affected by the cat's gravity. In other words, the central body is fixed.

This is somewhat similar to Code Golf: What is the spaceship's fate? [floating point version], but this is different because it is measuring accuracy.

You may implement a solution based on a simulation, which must run in less than 3 seconds, OR you may implement a program that gives an exact value (must also run in less than 3 seconds). See scoring details below. If you implement a simulation, it does not have to be exact, but your score will be lower because of the inaccuracy.

Input: x y z vx vy vz t, not necessarily integers representing x, y, z coordinates, velocity in the x, y, and z directions and time respectively. It is guaranteed that the speed of the cat is strictly less than the escape velocity at that altitude. Input may be taken from anywhere, including parameters to a function. The program must run in less than three seconds on my laptop for t < 2^30, which means, if you are running a simulation, you must adjust your timestep accordingly. If you are planning on hitting the 3 second limit for every test case, make sure there is a tuneable parameter that can make it more accurate/less accurate for speed gains, so that I can make it run in three seconds on my computer.

Output: x y z, the position after time t.

Since the two-body problem can be solved perfectly, it is in theory possible to get a perfect, correct answer.

Scoring: For any test case, the error is defined to be the distance between your output and the "true" output. The true output is defined to be the one that the test case snippet generates. If the error is less than 10^(-8), the error is rounded down to zero. Your score is the average error on 100 (or more) random test cases. If you write a perfectly accurate answer, you should get a score of 0; lowest score wins, and ties will be broken by code length.

Test cases:

1 0 0 0 -1 0 1000000000 --> 0.83789 -0.54584 0

In this case , the orbit is perfectly circular with period 2*pi, so after circling 159154943 times, the cat ends up at approximately (0.83789, -0.54584). This isn't a test case your code will be tested on; if you submit a perfectly accurate answer, however, you may want to test it on this.

The snippet below generates random additional test cases and will be used to judge submissions; let me know if there is a bug with this:

var generateTestCase = function() {
  var periapsis = Math.random() * 10 + 0.5;
  var circularVelocity = Math.sqrt(1 / periapsis);
  var escapeVelocity = Math.sqrt(2) * circularVelocity;
  var v = Math.random() * (escapeVelocity - circularVelocity) + circularVelocity;
  var a = 1 / (2 / periapsis - v * v);
  var ecc = (a - periapsis) * (2 / periapsis - v * v);
  var semiLatus = a * (1 - ecc * ecc);
  var angM = Math.sqrt(semiLatus);
  var period = 2 * Math.PI * a * Math.sqrt(a);
  var getCoordAfterTime = function(t) {
    var meanAnomaly = t / (a * Math.sqrt(a));

    function E(e, M, n) {
      var result = M;
      for (var k = 1; k <= n; ++k) result = M + e * Math.sin(result);
      return result;
    }
    var eAnomaly = E(ecc, meanAnomaly, 10000);
    var radius = a * (1 - ecc * Math.cos(eAnomaly));
    var theta = Math.acos((semiLatus / radius - 1) / ecc);
    if (t > period / 2) theta = -theta;
    return {
      r: radius,
      t: theta,
      x: radius * Math.cos(theta),
      y: radius * Math.sin(theta)
    };
  }
  var t1 = Math.random() * period;
  var p1 = getCoordAfterTime(t1);
  var v1 = Math.sqrt(2 / p1.r - 1 / a);
  var velocityAngle1 = Math.acos(angM / p1.r / v1);
  var horizontal = Math.atan(-p1.x / p1.y);
  if (t1 < period / 2) {
    velocityAngle1 *= -1;
    horizontal += Math.PI;
  }
  var vx = v1 * Math.cos(horizontal + velocityAngle1);
  var vy = v1 * Math.sin(horizontal + velocityAngle1);
  var t2 = Math.random() * period;
  var p2 = getCoordAfterTime(t2);
  var v2 = Math.sqrt(2 / p2.r - 1 / a);
  var dummyPeriods = Math.floor(Math.random() * 100) * period + period;
  var rotate = function(x, y, z, th1, th2, th3) {
    var x1 = x * Math.cos(th1) - y * Math.sin(th1);
    var y1 = x * Math.sin(th1) + y * Math.cos(th1);
    var z1 = z;
    var x2 = x1 * Math.cos(th2) + z1 * Math.sin(th2);
    var y2 = y1;
    var z2 = -x1 * Math.sin(th2) + z1 * Math.cos(th2);
    var x3 = x2;
    var y3 = y2 * Math.cos(th3) - z2 * Math.sin(th3);
    var z3 = y2 * Math.sin(th3) + z2 * Math.cos(th3);
    return [x3, y3, z3];
  }
  var th1 = Math.random() * Math.PI * 2,
    th2 = Math.random() * 2 * Math.PI,
    th3 = Math.random() * 2 * Math.PI;
  var input = rotate(p1.x, p1.y, 0, th1, th2, th3)
    .concat(rotate(vx, vy, 0, th1, th2, th3))
    .concat([t2 - t1 + dummyPeriods]);
  var output = rotate(p2.x, p2.y, 0, th1, th2, th3);
  return input.join(" ") + " --> " + output.join(" ");
}
var $ = function(x) {
  return document.querySelector(x);
}
$("#onetestcase").onclick = function() {
  $("#result1").innerHTML = generateTestCase();
}
$("#huntestcases").onclick = function() {
  var result = [];
  for (var i = 0; i < 100; ++i) result.push(generateTestCase());
  $("#result100").innerHTML = result.join("\n");
}

<div>
  <button id="onetestcase">Generate</button>
  <pre id="result1"></pre>
</div>
<div>
  <button id="huntestcases">Generate 100 test cases</button>
  <pre id="result100"></pre>
</div>

Answer 1

Python 3.5 + NumPy, exact, 186 bytes

from math import*
def o(r,v,t):
 d=(r@r)**.5;W=2/d-v@v;U=W**1.5;b=[0,t*U+9]
 while 1:
  a=sum(b)/2;x=1-cos(a);y=sin(a)/U;k=r@v*x/W+d*y*W
  if a in b:return k*v-r*x/W/d+r
  b[k+a/U-y>t]=a

This is an exact solution, using a formula I engineered based on Jesper Göranssonhis, “Symmetries of the Kepler problem”, 2015. It uses a binary search to solve the transcendental equation Ax + B cos x + C sin x = D, which has no closed-form solution.

The function expects the position and velocity to be passed in as NumPy arrays:

>>> from numpy import array
>>> o(array([1,0,0]),array([0,-1,0]),1000000000)
array([ 0.83788718, -0.54584345,  0.        ])
>>> o(array([-1.1740058273269156,8.413493259550673,0.41996042044140003]),array([0.150014367067652,-0.09438816345868332,0.37294941703455975]),7999.348650387233)
array([-4.45269544,  6.93224929, -9.27292488])

Answer 2

Javascript

This is just to get the ball rolling, since no one seems to be posting answers. Here is a very naive, simple way that can be improved a lot:

function simulate(x, y, z, vx, vy, vz, t) {
  var loops = 1884955; // tune this parameter
  var timestep = t / loops;
  for (var i = 0; i < t; i += timestep) {
    var distanceSq = x*x + y*y + z*z; // distance squared from origin
    var distance = Math.sqrt(distanceSq);
    var forceMag = 1/distanceSq; // get the force of gravity
    var forceX = -x / distance * forceMag;
    var forceY = -y / distance * forceMag;
    var forceZ = -z / distance * forceMag;
    vx += forceX * timestep;
    vy += forceY * timestep;
    vz += forceZ * timestep;
    x += vx * timestep;
    y += vy * timestep;
    z += vz * timestep;
  }
  return [x, y, z];
}

Testing:

simulate(1, 0, 0, 0, -1, 0, Math.PI*2) --> [0.9999999999889703, -0.0000033332840909716455, 0]

Hey, that's pretty good. It has an error of about 3.333*10^(-6) which isn't enough for it to be rounded down... it's close.

Just for fun:

console.log(simulate(1, 0, 0, 0, -1, 0, 1000000000))
--> [-530516643639.4616, -1000000000.0066016, 0]

Oh well; so it's not the best.

And on a random test case from the generator:

simulate(-1.1740058273269156,8.413493259550673,0.41996042044140003,0.150014367067652,-0.09438816345868332,0.37294941703455975,7999.348650387233)
-->    [-4.528366392498373, 6.780385554803544, -9.547824236472668]
Actual:[-4.452695438880813, 6.932249293597744, -9.272924876103785]

With an error of only approximately 0.32305!

This can be improved a lot by using something like Verlet integration or some fancy algorithm. In fact, those algorithms may even get perfect scores, despite being simulations.