{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# N body simulation" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "\n", "\n", "\n", "\n" ], "text/plain": [ "" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from IPython.core.display import HTML\n", "css_file = 'https://raw.githubusercontent.com/ngcm/training-public/master/ipython_notebook_styles/ngcmstyle.css'\n", "HTML(url=css_file)" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": true }, "outputs": [], "source": [ "import numpy\n", "from matplotlib import pyplot\n", "%matplotlib inline\n", "from matplotlib import rcParams\n", "rcParams['font.family'] = 'serif'\n", "rcParams['font.size'] = 16" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In [Greg's original post](http://software-carpentry.org/blog/2014/10/why-we-dont-teach-testing.html) he specifically references the [Python 3 version](http://benchmarksgame.alioth.debian.org/u32/program.php?test=nbody&lang=python3&id=1) of the [n-body benchmark](http://benchmarksgame.alioth.debian.org/u32/performance.php?test=nbody). In particular, he asks how to test the `advance` function. \n", "\n", "There's a number of useful and important comments underneath the original post, in particular [this thread on details of the algorithm and the problems with the code](http://software-carpentry.org/blog/2014/10/why-we-dont-teach-testing.html#comment-1662471640), but here we'll start from scratch." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's reproduce that code here (please note the [OSI BSD license](http://benchmarksgame.alioth.debian.org/license.html))." ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# The Computer Language Benchmarks Game\n", "# http://benchmarksgame.alioth.debian.org/\n", "#\n", "# originally by Kevin Carson\n", "# modified by Tupteq, Fredrik Johansson, and Daniel Nanz\n", "# modified by Maciej Fijalkowski\n", "# 2to3\n", "\n", "import sys \n", "\n", "def combinations(l):\n", " result = []\n", " for x in range(len(l) - 1):\n", " ls = l[x+1:]\n", " for y in ls:\n", " result.append((l[x],y))\n", " return result\n", "\n", "PI = 3.14159265358979323\n", "SOLAR_MASS = 4 * PI * PI\n", "DAYS_PER_YEAR = 365.24\n", "\n", "BODIES = {\n", " 'sun': ([0.0, 0.0, 0.0], [0.0, 0.0, 0.0], SOLAR_MASS),\n", "\n", " 'jupiter': ([4.84143144246472090e+00,\n", " -1.16032004402742839e+00,\n", " -1.03622044471123109e-01],\n", " [1.66007664274403694e-03 * DAYS_PER_YEAR,\n", " 7.69901118419740425e-03 * DAYS_PER_YEAR,\n", " -6.90460016972063023e-05 * DAYS_PER_YEAR],\n", " 9.54791938424326609e-04 * SOLAR_MASS),\n", "\n", " 'saturn': ([8.34336671824457987e+00,\n", " 4.12479856412430479e+00,\n", " -4.03523417114321381e-01],\n", " [-2.76742510726862411e-03 * DAYS_PER_YEAR,\n", " 4.99852801234917238e-03 * DAYS_PER_YEAR,\n", " 2.30417297573763929e-05 * DAYS_PER_YEAR],\n", " 2.85885980666130812e-04 * SOLAR_MASS),\n", "\n", " 'uranus': ([1.28943695621391310e+01,\n", " -1.51111514016986312e+01,\n", " -2.23307578892655734e-01],\n", " [2.96460137564761618e-03 * DAYS_PER_YEAR,\n", " 2.37847173959480950e-03 * DAYS_PER_YEAR,\n", " -2.96589568540237556e-05 * DAYS_PER_YEAR],\n", " 4.36624404335156298e-05 * SOLAR_MASS),\n", "\n", " 'neptune': ([1.53796971148509165e+01,\n", " -2.59193146099879641e+01,\n", " 1.79258772950371181e-01],\n", " [2.68067772490389322e-03 * DAYS_PER_YEAR,\n", " 1.62824170038242295e-03 * DAYS_PER_YEAR,\n", " -9.51592254519715870e-05 * DAYS_PER_YEAR],\n", " 5.15138902046611451e-05 * SOLAR_MASS) }\n", "\n", "\n", "SYSTEM = list(BODIES.values())\n", "PAIRS = combinations(SYSTEM)\n", "\n", "\n", "def advance(dt, n, bodies=SYSTEM, pairs=PAIRS):\n", "\n", " for i in range(n):\n", " for (([x1, y1, z1], v1, m1),\n", " ([x2, y2, z2], v2, m2)) in pairs:\n", " dx = x1 - x2\n", " dy = y1 - y2\n", " dz = z1 - z2\n", " mag = dt * ((dx * dx + dy * dy + dz * dz) ** (-1.5))\n", " b1m = m1 * mag\n", " b2m = m2 * mag\n", " v1[0] -= dx * b2m\n", " v1[1] -= dy * b2m\n", " v1[2] -= dz * b2m\n", " v2[0] += dx * b1m\n", " v2[1] += dy * b1m\n", " v2[2] += dz * b1m\n", " for (r, [vx, vy, vz], m) in bodies:\n", " r[0] += dt * vx\n", " r[1] += dt * vy\n", " r[2] += dt * vz" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In previous notebooks we compared the results of the algorithm with the expected behaviour within the context of a *model*. The model only considered the behaviour of the numerical algorithm: we took it for granted that the correct physics had been implemented. In this case, the `advance` function is doing two things, neither of which are immediately obvious:\n", "\n", "1. it evaluates the forces due to gravitational interaction between all of the particles, and\n", "2. it evolves the system of particles forwards in time.\n", "\n", "In testing whether the results are *close enough*, we need to test within the context of *two* models: we need to check that the correct physics is implemented, and we need to check that the numerical evolution behaves as expected." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## The models" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The physics that we're modelling is the [classical n body problem](http://en.wikipedia.org/wiki/N-body_problem) of massive particles interacting only through gravity. The [general formulation](http://en.wikipedia.org/wiki/N-body_problem#General_formulation) gives the equations of motion\n", "\n", "$$\n", "\\begin{equation}\n", " m_i \\frac{d^2 \\vec{r}_i}{d t^2} = \\sum_{j \\ne i} \\frac{ G m_i m_j (\\vec{r}_j - \\vec{r}_i) }{|\\vec{r}_j - \\vec{r}_i|^3}\n", "\\end{equation}\n", "$$\n", "\n", "where $m_i, \\vec{r}_i$ are the masses and positions of the $i^{\\text{th}}$ particle. This model is Newton's laws (which are equivalent to the conservation of energy and momentum) plus Newton's law of gravity." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The numerical method we need to model uses the [*semi-implicit Euler method*](http://en.wikipedia.org/wiki/Semi-implicit_Euler_method)1 to update in time, and a [direct particle-particle method](http://en.wikipedia.org/wiki/N-body_problem#Few_bodies) to compute the spatial interactions. The semi-implicit Euler method is first order accurate (with respect to the timestep $\\Delta t$), and the particle-particle interactions are exact." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "That is, we can think of the algorithm as solving the ODE\n", "\n", "$$\n", "\\begin{equation}\n", " \\frac{d y}{dt} = F(t, y)\n", "\\end{equation}\n", "$$\n", "\n", "where the vector $y$ contains the positions $\\vec{r}$ and velocities $\\vec{v}$ of each body. Due to the physics of the problem (ie, there's a total energy that we want to conserve) there's an advantage in considering the positions and velocities separately, but when it comes to measure the error we don't need to do that: we just treat the `advance` function as an algorithm that takes the full vector $y_n$ as input and returns $y_{n+1}$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "What this means in practical terms is that we expect the error in time to have similar convergence rate behaviour to the standard Euler method discussed before: the convergence rate should be $1$, which can be checked both by the self-convergence being within $0.585$ and $1.585$, and also that the error in the convergence rate goes down by at least $1/3$ as the base resolution of $\\Delta t$ is reduced by $2$. Checking the local truncation error would require more detailed analysis, but in principle should follow the same steps as before." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "What's more interesting here is the additional checks from the complexities of the particle-particle interaction. We can use our previous techniques to check that the time update algorithm is behaving as expected. However, the particle-particle interaction itself contains many steps, and those also need testing.\n", "\n", "Let's check the output on (a copy of) the default input, taking a *single* step." ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[([12.905190971986693, -15.102456648865521, -0.22341579268383283], [1.0821409847561863, 0.8694752833109706, -0.010821379117708814], 0.0017237240570597112), ([8.333218184594687, 4.1430349688596575, -0.40343728476075313], [-1.0148533649892928, 1.8236404735352374, 0.00861323535682711], 0.011286326131968767), ([4.847339930856543, -1.1321630550460413, -0.10387091638393478], [0.5908488391821772, 2.8156988981387063, -0.02488719128116714], 0.03769367487038949), ([15.389485801827046, -25.913363875177772, 0.17891118741867673], [0.9788686976130451, 0.5950734810190818, -0.034758553169444574], 0.0020336868699246304), ([1.5983797301314366e-07, -3.0187796864956455e-08, -3.730115986360403e-09], [1.5983797301314366e-05, -3.0187796864956456e-06, -3.7301159863604027e-07], 39.47841760435743)]\n", "[([12.894369562139131, -15.111151401698631, -0.22330757889265573], [1.0827910064415354, 0.8687130181696082, -0.010832637401363636], 0.0017237240570597112), ([8.34336671824458, 4.124798564124305, -0.4035234171143214], [-1.0107743461787924, 1.8256623712304119, 0.008415761376584154], 0.011286326131968767), ([4.841431442464721, -1.1603200440274284, -0.10362204447112311], [0.606326392995832, 2.81198684491626, -0.02521836165988763], 0.03769367487038949), ([15.379697114850917, -25.919314609987964, 0.17925877295037118], [0.979090732243898, 0.5946989986476762, -0.034755955504078104], 0.0020336868699246304), ([0.0, 0.0, 0.0], [0.0, 0.0, 0.0], 39.47841760435743)]\n" ] } ], "source": [ "import copy\n", "bodies = copy.deepcopy(SYSTEM)\n", "advance(dt=0.01, n=1, bodies=bodies, pairs=combinations(bodies))\n", "print(bodies)\n", "print(SYSTEM)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Comparing these two is going to be a real pain. So, let's create a utility function that will give the difference between two configurations of bodies. It will return the differences in the locations and the velocities as arrays, but will also include the \"total\" error as a single number - the sum of the absolute values of all errors. If this norm is zero then the two configurations are identical." ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def differences(bodies1, bodies2):\n", " \"\"\"\n", " Compare two configurations.\n", " \"\"\"\n", " assert len(bodies1) == len(bodies2), \"Configurations must have same number of bodies! {}, {}\".format(len(bodies1), len(bodies2))\n", " N = len(bodies1)\n", " d_positions = numpy.zeros((N, 3))\n", " d_velocities = numpy.zeros((N, 3))\n", " norm_difference = 0.0\n", " for n in range(N):\n", " d_positions[n, :] = numpy.array(bodies1[n][0]) - numpy.array(bodies2[n][0])\n", " d_velocities[n, :] = numpy.array(bodies1[n][1]) - numpy.array(bodies2[n][1])\n", " norm_difference += numpy.sum(numpy.abs(d_positions[n, :])) + numpy.sum(numpy.abs(d_velocities[n, :]))\n", " return norm_difference, d_positions, d_velocities" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false }, "outputs": [], "source": [ "norm1, d_x, d_v = differences(bodies, SYSTEM)" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "array([[ -6.50021685e-04, 7.62265141e-04, 1.12582837e-05],\n", " [ -4.07901881e-03, -2.02189770e-03, 1.97473980e-04],\n", " [ -1.54775538e-02, 3.71205322e-03, 3.31170379e-04],\n", " [ -2.22034631e-04, 3.74482371e-04, -2.59766537e-06],\n", " [ 1.59837973e-05, -3.01877969e-06, -3.73011599e-07]])" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "d_v" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we have to think about how to test the particle-particle interaction. The parameter we can change in order to do the test is the configuration `bodies` (and the associated list of `pairs`, which we'll always take to be `combinations(bodies)`)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So we can either change the positions, or the velocities, of the bodies, or their masses (or some combination of the three). What changes make meaningful tests?\n", "\n", "Well, the coordinates chosen should not matter, *if* we change them in a self-consistent way, *and* interpret the results appropriately. For example, if we *ran time backwards* by changing the direction of time (`dt` $\\to$ `-dt`) and changing the sign of the velocities, then the results should be identical (after re-flipping the direction of time)." ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def flip_time(bodies):\n", " \"\"\"\n", " Flip the time by flipping the velocity.\n", " \"\"\"\n", " \n", " for i in range(len(bodies)):\n", " for j in range(3):\n", " bodies[i][1][j] *= -1.0" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Norm of differences is 0.0.\n" ] } ], "source": [ "bodies_flip_time = copy.deepcopy(SYSTEM)\n", "flip_time(bodies_flip_time)\n", "advance(dt=-0.01, n=1, bodies=bodies_flip_time, pairs=combinations(bodies_flip_time))\n", "flip_time(bodies_flip_time)\n", "norm_flip_time, d_x_flip_time, d_v_flip_time = differences(bodies, bodies_flip_time)\n", "print(\"Norm of differences is {}.\".format(norm_flip_time))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As expected, the positions are identical. Note that the velocities are significantly different, because they still have the opposite sign, but if we flip the velocities back as well, then they will be identical." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Why did we do this? Because one of the first steps in debugging, and hence in testing, should be to [check the conservation laws](http://software-carpentry.org/v5/novice/python/05-defensive.html#know-what-its-supposed-to-do). This is just an example of a conservation law resulting from a symmetry: the consequence of [Emmy Noether's famous theorem](http://en.wikipedia.org/wiki/Noether's_theorem). There is a [detailed list][1] of symmetries that a physical theory could have, each of which has an associated conservation law. Not all physical models obey all symmetries, and not every numerical simulation obeys all the symmetries of the physical model (maybe due to boundary conditions, or numerical approximations). However, in the n-body case we expect a large number of these symmetries to hold, which we can check.\n", "\n", "These include:\n", "\n", "* discrete parity symmetry (sending $\\vec{r} \\to -\\vec{r}$)\n", "* discrete time symmetry (sending $t \\to -t$)\n", "* discrete charge symmetry (not relevant here - no charges)\n", "* continuous space translation ($\\vec{r} \\to \\vec{r} + \\vec{a}$)\n", "* continuous space rotation\n", "* continuous time translation (this is obeyed by the continuum model, but not by the numerical model)\n", "* continuous conformal (like) symmetry ($t \\to b t$, $\\vec{r} \\to b \\vec{r}$, and $m \\to b m$ all together).\n", "\n", "[1]: http://en.wikipedia.org/wiki/Symmetry_(physics)#Conservation_laws_and_symmetry" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the first test where we \"ran time in reverse\" we were checking discrete time symmetry. We can now check the other symmetries systematically." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Discrete space symmetry" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We should be able to flip individual coordinates: if $y \\to -y$ and $v_y \\to -v_y$, for example, then the evolution of the other coordinates and velocities should be unaffected (and if this flip is reversed after evolution, the position should be identical)." ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def flip_coordinate(bodies, coord):\n", " \"\"\"\n", " Flip a single coordinate direction\n", " \"\"\"\n", " \n", " for i in range(len(bodies)):\n", " bodies[i][0][coord] *= -1.0\n", " bodies[i][1][coord] *= -1.0" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Norm of differences is 0.0 (flipped coordinate 0).\n", "Norm of differences is 0.0 (flipped coordinate 1).\n", "Norm of differences is 0.0 (flipped coordinate 2).\n" ] } ], "source": [ "for coord in range(3):\n", " bodies_flip_coord = copy.deepcopy(SYSTEM)\n", " flip_coordinate(bodies_flip_coord, coord)\n", " advance(dt=0.01, n=1, bodies=bodies_flip_coord, pairs=combinations(bodies_flip_coord))\n", " flip_coordinate(bodies_flip_coord, coord)\n", " norm_flip_coord, d_x_flip_coord, d_v_flip_coord = differences(bodies, bodies_flip_coord)\n", " print(\"Norm of differences is {} (flipped coordinate {}).\".format(norm_flip_coord, coord))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Continuous space translation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next, note that only the *separation* between bodies is meant to matter. So if we arbitrarily translate all coordinates by some amount whilst leaving the velocities untouched, and reverse this after evolution, the result should be identical." ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def translate_coordinate(bodies, shift):\n", " \"\"\"\n", " Translate the coordinates of all bodies\n", " \"\"\"\n", " \n", " for i in range(len(bodies)):\n", " for n in range(3):\n", " bodies[i][0][n] += shift[n]" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Norm of differences is 3.090701803516006e-15.\n" ] } ], "source": [ "shift = 10.0*(-0.5+numpy.random.rand(3)) # Random coordinate shift in [-5, 5]\n", "bodies_shift = copy.deepcopy(SYSTEM)\n", "translate_coordinate(bodies_shift, shift)\n", "advance(dt=0.01, n=1, bodies=bodies_shift, pairs=combinations(bodies_shift))\n", "translate_coordinate(bodies_shift, -shift)\n", "norm_shift, d_x_shift, d_v_shift = differences(bodies, bodies_shift)\n", "print(\"Norm of differences is {}.\".format(norm_shift))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this case the repeated operations introduce some floating point round-off error. This should be related to the round off introduced by the shift, and bounded by the number of bodies ($5$) multiplied by the total round-off introduced:" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "5.5511151231257827e-15" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "numpy.sum(numpy.abs(numpy.spacing(shift)))*5" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The precise differences are:" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[ -1.77635684e-15 0.00000000e+00 -1.11022302e-16]\n", " [ 0.00000000e+00 0.00000000e+00 -1.11022302e-16]\n", " [ 8.88178420e-16 0.00000000e+00 8.32667268e-17]\n", " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", " [ -6.62471880e-17 -4.71076778e-17 7.50034688e-18]]\n", "[[ 0. 0. 0.]\n", " [ 0. 0. 0.]\n", " [ 0. 0. 0.]\n", " [ 0. 0. 0.]\n", " [ 0. 0. 0.]]\n" ] } ], "source": [ "print(d_x_shift)\n", "print(d_v_shift)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We see that the some of the bodies have errors around the expected floating point limit, and these dominate the resulting error (which is well within the expected bound)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Continuous space rotation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Yet another coordinate transformation can be tried without affecting the results: a rotation. We'll pick three random Euler angles and [set up the rotation matrix](http://en.wikipedia.org/wiki/Rotation_matrix#In_three_dimensions). In general, we should do a translation first to ensure no body is at the origin (here, the Sun will be), to avoid specializing the problem, but given the generality of the code that's a very minor worry." ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def rotate_bodies(bodies, angles, invert=False):\n", " \"\"\"\n", " Rotate the coordinates of all bodies.\n", " \"\"\"\n", " \n", " Rx = numpy.array([[1.0, 0.0, 0.0],\n", " [0.0, numpy.cos(angles[0]), -numpy.sin(angles[0])],\n", " [0.0, numpy.sin(angles[0]), numpy.cos(angles[0])]])\n", " Ry = numpy.array([[numpy.cos(angles[1]), 0.0, numpy.sin(angles[1])],\n", " [0.0, 1.0, 0.0],\n", " [-numpy.sin(angles[1]), 0.0, numpy.cos(angles[1])]])\n", " Rz = numpy.array([[numpy.cos(angles[2]), -numpy.sin(angles[2]), 0.0],\n", " [numpy.sin(angles[2]), numpy.cos(angles[2]), 0.0],\n", " [0.0, 0.0, 1.0]])\n", " if invert:\n", " R = numpy.dot(numpy.dot(Rx, Ry), Rz)\n", " else:\n", " R = numpy.dot(Rz, numpy.dot(Ry, Rx))\n", " for i in range(len(bodies)):\n", " x = numpy.array(bodies[i][0])\n", " v = numpy.array(bodies[i][1])\n", " xp = numpy.dot(R, x)\n", " vp = numpy.dot(R, v)\n", " for n in range(3):\n", " bodies[i][0][n] = xp[n]\n", " bodies[i][1][n] = vp[n]" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Norm of differences is 2.7932524155836443e-14.\n" ] } ], "source": [ "angles = numpy.pi/4.0*numpy.random.rand(3) # Random Euler angles in [0, pi/4]\n", "bodies_rotate = copy.deepcopy(SYSTEM)\n", "rotate_bodies(bodies_rotate, angles)\n", "advance(dt=0.01, n=1, bodies=bodies_rotate, pairs=combinations(bodies_rotate))\n", "rotate_bodies(bodies_rotate, -angles, invert=True)\n", "norm_rotate, d_x_rotate, d_v_rotate = differences(bodies, bodies_rotate)\n", "print(\"Norm of differences is {}.\".format(norm_rotate))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The size of the differences is considerably larger, thanks to the number of operations performed. With three rotation matrices combined, each containing nine entries, applied to five bodies, we would expect a total error of the order of:" ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "4.496403249731884e-14" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "3*9*5*numpy.sum(numpy.abs(numpy.spacing(angles)))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Continuous conformal (like) symmetry" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We could also try scaling the coordinates and the masses. As the [force goes as](http://en.wikipedia.org/wiki/N-body_problem#General_formulation) $M^2/L^2$, if we scale time, position and mass by the same amount then the results should not change. Note that because we scale *both* length and time, the velocity does not change." ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def scale_bodies(bodies, scale):\n", " \"\"\"\n", " Scale coordinates and masses.\n", " \"\"\"\n", " \n", " bodies_scale = []\n", " for (x, v, m) in bodies:\n", " new_x = copy.deepcopy(x)\n", " new_v = copy.deepcopy(v)\n", " new_m = m * scale\n", " for i in range(3):\n", " new_x[i] *= scale\n", " bodies_scale.append((new_x, new_v, new_m))\n", " \n", " return bodies_scale" ] }, { "cell_type": "code", "execution_count": 20, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Norm of differences is 0.0.\n" ] } ], "source": [ "scale = 2.0\n", "bodies_scale = scale_bodies(SYSTEM, scale)\n", "advance(dt=0.02, n=1, bodies=bodies_scale, pairs=combinations(bodies_scale))\n", "bodies_rescale = scale_bodies(bodies_scale, 1.0/scale)\n", "norm_scale, d_x_scale, d_v_scale = differences(bodies, bodies_rescale)\n", "print(\"Norm of differences is {}.\".format(norm_scale))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This scaling does not *prove* that the force goes as $M^2/L^2$; it only shows that the force contains $M$ and $L$ to the same power. To show that it's got the appropriate form we should either [compare to an oracle or test simpler cases](http://software-carpentry.org/v5/novice/python/05-defensive.html#know-what-its-supposed-to-do). The [Java n-body code](http://benchmarksgame.alioth.debian.org/u32/benchmark.php?test=nbody&lang=java) is (likely) the oracle it was tested against; testing against a simpler case will be done later." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Convergence - continuous time translation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Above, we noted that continuous time translation (which corresponds to conservation of energy) may not be perfectly retained by our *numerical method*, even though it is perfectly preserved by the theory. However, as we're solving an ODE in time, we can use the techniques applied to the phugoid model to ensure the numerical method is converging at the expected rate. As we're expecting the method to converge at first order, we need the error in the measured self-convergence rate to reduce by at least $1/3$ when the resolution (in time) is reduced by $2$." ] }, { "cell_type": "code", "execution_count": 21, "metadata": { "collapsed": true }, "outputs": [], "source": [ "T = 10.0 # The base resolution will take 1000 steps\n", "dt_values = numpy.array([0.01*2**(-i) for i in range(4)])\n", "bodies_list = []\n", "for i, dt in enumerate(dt_values):\n", " bodies_loop = copy.deepcopy(SYSTEM)\n", " advance(dt=dt, n=int(T/dt), bodies=bodies_loop, pairs=combinations(bodies_loop))\n", " bodies_list.append(bodies_loop)\n" ] }, { "cell_type": "code", "execution_count": 22, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Convergence rate (base dt=0.01) is 0.9964571407422125 (error 0.003542859257787523).\n", "Convergence rate (base dt=0.005) is 0.9982227199049917 (error 0.0017772800950083267).\n", "Is the convergence rate close enough to 1 for the answers to match? True\n", "Does the convergence of the convergence rate show it's close enough to 1? True\n" ] } ], "source": [ "convergence_rate = numpy.zeros((2,))\n", "for i in range(len(convergence_rate)):\n", " numerator, diff_x, diff_v= differences(bodies_list[i], bodies_list[i+1])\n", " denominator, diff_x, diff_v = differences(bodies_list[i+1], bodies_list[i+2])\n", " convergence_rate[i] = numpy.log(numerator/denominator)/numpy.log(2.0)\n", " print(\"Convergence rate (base dt={}) is {} (error {}).\".format(\n", " dt_values[i], convergence_rate[i], numpy.abs(convergence_rate[i]-1.0)))\n", "print(\"Is the convergence rate close enough to 1 for the answers to match? {}\".format(\n", " numpy.log(1.0+0.5)/numpy.log(2.0) < convergence_rate[1] < numpy.log(2.0**2-1.0)/numpy.log(2.0)))\n", "print(\"Does the convergence of the convergence rate show it's close enough to 1? {}\".format(\n", " numpy.abs(convergence_rate[1]-1.0) < 2.0/3.0*numpy.abs(convergence_rate[0]-1.0)))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This shows that the time evolution is close enought to first order as expected. We haven't explicitly shown that it's the semi-implicit Euler method, as explicitly calculating the local truncation error would be a lot of work, and in this case we're really not interested in the specific time integration scheme - just that it converges." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Checking the specific force law" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In order to check that we've really implemented Newton's law of gravity, we need to test a very simplified situation where we can predict the precise answer. Let's set up a two body problem starting from rest.\n", "\n", "Once we've done that, we can vary one parameter and check that the acceleration behaves as expected. The two cases we have are varying a mass, so that if $m_1 \\to b m_1$ then\n", "\n", "$$\n", "\\begin{equation}\n", " \\vec{a}_2(b) = b \\vec{a_2}(1)\n", "\\end{equation}\n", "$$\n", "\n", "and varying separation $|\\vec{r}_1 - \\vec{r}_2|$ so that $|\\vec{r}_1 - \\vec{r}_2| \\to c |\\vec{r}_1 - \\vec{r}_2|$\n", "\n", "$$\n", "\\begin{equation}\n", " \\vec{a}_2(c) = c^{-2} \\vec{a}(1).\n", "\\end{equation}\n", "$$\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Over a single timestep $\\Delta t$ we can approximate the acceleration as\n", "\n", "$$\n", "\\begin{equation}\n", " \\vec{a}_2(t=0) \\simeq \\frac{\\vec{v}_2(t=\\Delta t) - \\vec{v}_2(t=0)}{\\Delta t}.\n", "\\end{equation}\n", "$$\n", "\n", "For *this specific algorithm* we have that \n", "\n", "$$\n", "\\begin{equation}\n", " \\vec{v}_2(t=\\Delta t) = \\vec{v}_2(t=0) + \\Delta t \\, \\vec{a}_2(t=0),\n", "\\end{equation}\n", "$$\n", "\n", "which means that the approximation for the acceleration is, in fact, exact:\n", "\n", "$$\n", "\\begin{equation}\n", " \\vec{a}_2(t=0) = \\frac{\\vec{v}_2(t=\\Delta t) - \\vec{v}_2(t=0)}{\\Delta t}.\n", "\\end{equation}\n", "$$\n", "\n", "For a general algorithm we would have a first order in $\\Delta t$ error in our approximation to the acceleration, so would have to use Richardson extrapolation methods as before in order to calculate the acceleration, and hence check the force." ] }, { "cell_type": "code", "execution_count": 23, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The acceleration (hence velocity) should decrease as separation^2\n", "So these two numbers should match: 0.009869604401089358 0.009869604401089358 (difference: 0.0)\n", "The acceleration (hence velocity) should increase as mass^1\n", "So these two numbers should match: 0.009869604401089358 0.009869604401089358 (difference: 0.0)\n" ] } ], "source": [ "body_1 = ([1.0, 0.0, 0.0], [0.0, 0.0, 0.0], 0.1*SOLAR_MASS)\n", "body_2 = ([-1.0, 0.0, 0.0], [0.0, 0.0, 0.0], 0.1*SOLAR_MASS)\n", "body_1_2_separation = ([3.0, 0.0, 0.0], [0.0, 0.0, 0.0], 0.1*SOLAR_MASS)\n", "body_1_2_mass = ([1.0, 0.0, 0.0], [0.0, 0.0, 0.0], 0.2*SOLAR_MASS)\n", "two_bodies = [copy.deepcopy(body_1), copy.deepcopy(body_2)]\n", "two_bodies_2_separation = [copy.deepcopy(body_1_2_separation), copy.deepcopy(body_2)]\n", "two_bodies_2_mass = [copy.deepcopy(body_1_2_mass), copy.deepcopy(body_2)]\n", "advance(dt=0.01, n=1, bodies=two_bodies, pairs=combinations(two_bodies))\n", "advance(dt=0.01, n=1, bodies=two_bodies_2_separation, pairs=combinations(two_bodies_2_separation))\n", "advance(dt=0.01, n=1, bodies=two_bodies_2_mass, pairs=combinations(two_bodies_2_mass))\n", "print(\"The acceleration (hence velocity) should decrease as separation^2\")\n", "print(\"So these two numbers should match: {} {} (difference: {})\".format(\n", " two_bodies[1][1][0], \n", " two_bodies_2_separation[1][1][0]*4.0,\n", " numpy.abs(two_bodies[1][1][0]-two_bodies_2_separation[1][1][0]*4.0)))\n", "print(\"The acceleration (hence velocity) should increase as mass^1\")\n", "print(\"So these two numbers should match: {} {} (difference: {})\".format(\n", " two_bodies[1][1][0], \n", " two_bodies_2_mass[1][1][0]/2.0,\n", " numpy.abs(two_bodies[1][1][0]-two_bodies_2_mass[1][1][0]/2.0)))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Of course, we could have just got (un)lucky in comparing single data points when we vary either $b$ or $c$. Instead we can choose a set of points at random, and fit to a more general model; for example, fit to\n", "\n", "$$\n", "\\begin{equation}\n", " \\vec{a}_2 (b) = \\sum_{i=0}^4 p_i b^i\n", "\\end{equation}\n", "$$\n", "\n", "or\n", "\n", "$$\n", "\\begin{equation}\n", " \\vec{a}_2 (c) = \\sum_{i=0}^4 p_i c^{-i}:\n", "\\end{equation}\n", "$$\n" ] }, { "cell_type": "code", "execution_count": 24, "metadata": { "collapsed": false }, "outputs": [], "source": [ "bodies_list_separation = []\n", "separation_scale = numpy.random.rand(5)\n", "separations = separation_scale + 1.0\n", "separation_v = numpy.zeros_like(separations)\n", "for i, scale in enumerate(separation_scale):\n", " body_1_separation = ([scale, 0.0, 0.0], [0.0, 0.0, 0.0], 0.1*SOLAR_MASS)\n", " two_bodies_separation = [copy.deepcopy(body_1_separation), copy.deepcopy(body_2)]\n", " advance(dt=0.01, n=1, bodies=two_bodies_separation, pairs=combinations(two_bodies_separation))\n", " bodies_list_separation.append(two_bodies_separation)\n", " separation_v[i] = two_bodies_separation[1][1][0]\n", " \n", "bodies_list_mass = []\n", "mass_scale = numpy.random.rand(5)\n", "masses = 0.1 * mass_scale\n", "mass_v = numpy.zeros_like(masses)\n", "for i, scale in enumerate(mass_scale):\n", " body_1_mass = ([1.0, 0.0, 0.0], [0.0, 0.0, 0.0], 0.1*scale*SOLAR_MASS)\n", " two_bodies_mass = [copy.deepcopy(body_1_mass), copy.deepcopy(body_2)]\n", " advance(dt=0.01, n=1, bodies=two_bodies_mass, pairs=combinations(two_bodies_mass))\n", " bodies_list_mass.append(two_bodies_mass)\n", " mass_v[i] = two_bodies_mass[1][1][0]" ] }, { "cell_type": "code", "execution_count": 25, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "We expect the third-to-last (separation^{-2}) coefficient to dominate:\n", "Coefficients from separation: [ 2.80087953e-13 -8.47221417e-13 3.94784176e-02 -4.76829993e-13\n", " 8.87847308e-14]\n", "Coefficient of separation^{-2}: 0.039478\n", "Largest other coefficient: 8.4722e-13\n", "We expect the second-to-last (mass^{1}) coefficient to dominate:\n", "Coefficients from mass: [ 1.30669862e-11 -2.85143237e-12 2.14919187e-13 9.86960440e-02\n", " 6.96273251e-17]\n", "Coefficient of mass^{1}: 0.098696\n", "Largest other coefficient: 1.3067e-11\n" ] } ], "source": [ "p_separation = numpy.polyfit(1./separations, separation_v, len(separations)-1)\n", "p_mass = numpy.polyfit(masses, mass_v, len(masses)-1)\n", "print(\"We expect the third-to-last (separation^{-2}) coefficient to dominate:\")\n", "print(\"Coefficients from separation: {}\".format(p_separation))\n", "print(\"Coefficient of separation^{{-2}}: {:.5g}\".format(p_separation[-3]))\n", "print(\"Largest other coefficient: {:.5g}\".format(numpy.max(numpy.abs(numpy.delete(p_separation,-3)))))\n", "print(\"We expect the second-to-last (mass^{1}) coefficient to dominate:\")\n", "print(\"Coefficients from mass: {}\".format(p_mass))\n", "print(\"Coefficient of mass^{{1}}: {:.5g}\".format(p_mass[-2]))\n", "print(\"Largest other coefficient: {:.5g}\".format(numpy.max(numpy.abs(numpy.delete(p_mass,-2)))))" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false }, "source": [ "We we see that the expected coefficient dominates, but we're not getting a fit that's absolutely perfect. This is likely due to the limitations of the fitting algorithm." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So, we have an algorithm that\n", "\n", "* obeys all the expected symmetries;\n", "* depends on the size of the separation and the mass only;\n", "* scales as expected with mass and separation.\n", "\n", "This strongly suggests that it implements Newton's laws, and the force term that's implemented has the form\n", "\n", "$$\n", " \\vec{F} = G \\frac{m_1 m_2 (\\vec{r}_2 - \\vec{r}_1)}{|\\vec{r}_2 - \\vec{r}_1|^3}.\n", "$$\n", "\n", "The only remaining question: is the value of $G$ correct?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Our previous calculations give us the value of $G$ internal to the algorithm; we just take our value for the acceleration $\\vec{a}_2 = \\vec{F}/m_2$ and take out the masses and separations:" ] }, { "cell_type": "code", "execution_count": 26, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Algorithm G, in units of AU^3 SolarMass^{-1} years^{-2}: 39.47841760435743\n" ] } ], "source": [ "dt = 0.01\n", "separation = 2.0\n", "mass = 0.1\n", "G = two_bodies[1][1][0]/mass*separation**2/dt \n", "print(\"Algorithm G, in units of AU^3 SolarMass^{{-1}} years^{{-2}}: {}\".format(G))" ] }, { "cell_type": "code", "execution_count": 27, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Accepted value of G, in units of m^3 kg^-1 s^-2: 6.67384e-11\n" ] } ], "source": [ "from scipy import constants\n", "print(\"Accepted value of G, in units of {}: {}\".format(\n", " constants.unit('Newtonian constant of gravitation'), \n", " constants.G))" ] }, { "cell_type": "code", "execution_count": 28, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Algorithm G, in units of m^3 kg^-1 s^-2: 6.67253235346714e-11\n" ] } ], "source": [ "solar_mass_in_kg = 1.9891e30\n", "year_in_seconds = 3.15569e7\n", "AU_in_metres = 149597870700.0\n", "print(\"Algorithm G, in units of m^3 kg^-1 s^-2: {}\".format(\n", " G*(AU_in_metres**3/(solar_mass_in_kg)/year_in_seconds**2)))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So the constants agree to four significant figures; without more in-depth knowledge of the original implementation, this might just reflect experimental uncertainties in the knowledge of $G$ at the time. Given the [experimental uncertainty](http://en.wikipedia.org/wiki/Gravitational_constant#Laws_and_constants) in the value of the gravitational constant, I would say this check is close enough." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### Footnotes" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
    \n", "
  1. \n", "

    The method that got me to this conclusion is a classic *just Google it* style, as follows. Looking at the original n-body example there's no useful documentation, but the [about section](http://benchmarksgame.alioth.debian.org/u32/performance.php?test=nbody#about) says that is uses a simple [symplectic integrator](http://en.wikipedia.org/wiki/Symplectic_integrator). The [first order example](http://en.wikipedia.org/wiki/Symplectic_integrator#A_first-order_example) is the [semi-implicit Euler method](http://en.wikipedia.org/wiki/Semi-implicit_Euler_method), which appears to match the code. ↩

    \n", "

    2. \n", "
" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.4.2" } }, "nbformat": 4, "nbformat_minor": 0 }