{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Taylor integration of the Kepler problem\n", "\n", "Here, we try to reproduce __exactly__ the [Kepler problem integration example](http://nbviewer.jupyter.org/github/JuliaDiff/TaylorSeries.jl/blob/master/examples/1-KeplerProblem.ipynb) made by Luis Benet in [JuliaDiff/TaylorSeries.jl](https://github.com/JuliaDiff/TaylorSeries.jl).\n", "\n", "The Kepler problem is the basis of planetary motion; it describes the motion of a secondary body (e.g., a planet, asteroid, comet, etc.), around a primary body (e.g., the Sun). In cartesian coordinates over the orbital plane, the Hamiltonian for the Kepler problem reads:\n", "\n", "$$\n", "\\begin{align}\n", "H_{\\mathrm{Kepler}} &= \\frac{1}{2\\mu}(p_x^2+p_y^2)-\\frac{\\mu}{\\sqrt{x^2+y^2}}\n", "\\end{align}\n", "$$\n", "\n", "where $\\mu=G(m_1+m_2)$, $G$ is the gravitational constant, $m_1$ is the mass of the primary body and $m_2$ is the mass of the secondary body. If we write $\\vec r_1 = (x_1,y_1)$ and $\\vec r_2 = (x_2,y_2)$, respectively, for the position of the primary and secondary body, then the vector $\\vec r = \\vec r_2-\\vec r_1$, with coordinates $\\vec r = (x, y)=(x_2-x_1,y_2-y_1)$, represents the position of the secondary body relative to the primary body; i.e., $\\vec r$ represents the so-called relative coordinates. In terms of the vector $\\vec{r}$, the position $\\vec{r}=0$ corresponds to the position of the primary body.\n", "\n", "Using Hamilton equations, we can obtain the equations of motion for the Kepler problem:\n", "\n", "$$\n", "\\begin{align}\n", "\\dot x &= u \\\\\n", "\\dot y &= v \\\\\n", "\\dot u &= -\\frac{\\mu x}{(x^2+y^2)^{3/2}}\\\\\n", "\\dot v &= -\\frac{\\mu y}{(x^2+y^2)^{3/2}}\n", "\\end{align}\n", "$$\n", "\n", "Note that the canonical momenta are $p_x$ and $p_y$, while $u$ and $v$ are, respectively, the $x$ and $y$ components of the velocity; i.e., $p_x = \\mu u$ and $p_y = \\mu v$. \n", "\n", "\n", "First of all, we shall include all relevant packages:" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "Fontconfig warning: ignoring UTF-8: not a valid region tag\n" ] }, { "data": { "text/plain": [ "Plots.PyPlotBackend()" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "using TaylorIntegration, Plots, LaTeXStrings\n", "pyplot()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Some parameters necessary for the integration:\n", "\n", "+ $\\mu$: the gravitational parameter\n", "+ `q0`: the initial condition (we will select an initial condition which corresponds to elliptical motion)\n", "+ `order`: the order of the Taylor expansion\n", "+ `t_max`: the final time of the integration\n", "+ `abs_tol`: the absolute tolerance\n", "+ `n_iter`: the number of time-steps" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "500000" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "const μ = 1.0\n", "const q0 = [0.19999999999999996, 0.0, 0.0, 3.0] # a initial condition for elliptical motion\n", "const order = 28\n", "const t0 = 0.0\n", "const t_max = 10000*(2π) # we are just taking a wild guess about the period ;)\n", "const abs_tol = 1.0E-20\n", "const steps = 500000" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As usual, we write down the equations of motion into a `function`, which here we will name `kepler!`. `q` represents the system state; i.e., the set of values of the dynamical variables at a given instant; `dq` represents the time derivatives of the components of `q`." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "kepler! (generic function with 1 method)" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "#an auxiliary array which helps with optimization:\n", "const r_p3d2 = Array{TaylorSeries.Taylor1{Float64}}(1)\n", "\n", "#the equations of motion for the Kepler problem:\n", "function kepler!(t, q, dq)\n", " r_p3d2[1] = (q[1]^2+q[2]^2)^(3/2)\n", " \n", " dq[1] = q[3]\n", " dq[2] = q[4]\n", " dq[3] = -μ*q[1]/r_p3d2[1]\n", " dq[4] = -μ*q[2]/r_p3d2[1]\n", " \n", " nothing\n", "end" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The Taylor integration:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " 43.060857 seconds (464.94 M allocations: 39.523 GB, 12.87% gc time)\n" ] } ], "source": [ "t, q = taylorinteg(kepler!, q0, t0, 0.01, order, abs_tol, maxsteps=2); #warm-up lap\n", "@time t, q = taylorinteg(kepler!, q0, t0, t_max, order, abs_tol, maxsteps=steps);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The final state:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(10000.0,[0.2,2.72854e-9,-2.27378e-8,3.0])" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "t[end]/(2pi), q[end,:]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's extract the values of $x$, $y$, $u$ and $v$ for each time-step:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "x, y, u, v = view(q,:,1), view(q,:,2), view(q,:,3), view(q,:,4);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Orbital motion" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The initial conditions we selected correspond to a elliptical orbit with a relatively high eccentricity: $e=0.8$ for the initial condition `q0=[0.19999999999999996, 0.0, 0.0, 3.0]`. How does the motion of the planet/asteroid/comet looks like? Well, let's plot its orbit over the $x-y$ plane (the orbit is shown in blue; the position of the primary body is shown as a yellow dot):" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "scatter(\n", "\n", "[0.0], [0.0],\n", "label = L\"\\mathrm{Center\\; of\\; attraction}\",\n", "ms=10\n", "\n", ")\n", "scatter!(\n", "\n", "x[1:75:end], y[1:75:end],\n", "title = L\"\\mathrm{Orbital\\;\\;motion}\", \n", "xaxis = (L\"x\", (-2.0, 0.5), -2.0:0.5:2.0), \n", "yaxis = (L\"y\", (-0.8, 0.8), -2.0:0.5:2.0),\n", "label = L\"\\mathrm{Keplerian\\; ellipse}\",\n", "ms=0.5\n", "\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Energy conservation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Below, we write the energy function of the Kepler problem; i.e., the Hamiltonian $H_{\\mathrm{Kepler}}$ in terms of $x$, $y$, $u$ and $v$:" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "H (generic function with 1 method)" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "H(x_, y_, u_, v_) = 0.5*(u_^2+v_^2)-μ/sqrt(x_^2+y_^2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, using the function `H` defined above, we calculate the energy during each time-step:" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "E = H.(x, y, u, v);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$E_0=H_{\\mathrm{Kepler}}(x_0,y_0,u_0,v_0)$, the initial value of the energy, is:" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "-0.5000000000000009" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "E0=E[1]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We define $\\delta E$ as the relative error in the energy; i.e.:\n", "\n", "$$\n", "\\begin{align}\n", "\\delta E(t) &= \\frac{E(t)-E_0}{E_0}\n", "\\end{align}\n", "$$" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "436153-element Array{Float64,1}:\n", " -0.0 \n", " 1.77636e-15\n", " 1.77636e-15\n", " 8.88178e-16\n", " 8.88178e-16\n", " 2.66454e-15\n", " 3.55271e-15\n", " 2.66454e-15\n", " 3.55271e-15\n", " 2.66454e-15\n", " 3.9968e-15 \n", " 2.66454e-15\n", " 3.10862e-15\n", " ⋮ \n", " -7.12763e-14\n", " -7.19425e-14\n", " -7.14984e-14\n", " -7.10543e-14\n", " -7.19425e-14\n", " -7.10543e-14\n", " -7.10543e-14\n", " -6.92779e-14\n", " -7.01661e-14\n", " -6.92779e-14\n", " -6.92779e-14\n", " -7.10543e-14" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "δE = (E-E0)/(E0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The analytical solution preserves the energy; i.e., $\\delta E(t)=0$ for the analytical solution. Thus, we expect our solution to be close to zero. So, even if the $\\delta E$ is not perfectly zero during the whole integration, it is comparable to Julia's `Float64` machine-epsilon.\n", "\n", "Below, we plot $\\delta E$ in units of Julia's `Float64` machine-epsilon as a function of time, $t$. In Julia, the machine-epsilon for the `Float64` type has a value" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "2.220446049250313e-16" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "eps(Float64)" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "plot(\n", "\n", "t/(2π), δE/eps(Float64),\n", "title = L\"\\mathrm{Energy\\;\\; relative\\;\\; error\\;\\; vs\\;\\; time}\", \n", "xaxis = (L\"t\\mathrm{\\,(orbital\\;\\;periods)}\"),\n", "yaxis = (L\"\\delta E \\mathrm{\\;\\;(machine\\;\\;epsilons)}\"),\n", "label = L\"\\mathrm{Energy\\;\\; relative\\;\\; error\\;\\; vs\\;\\; time}\"\n", "\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, how does the energy error distribute around zero?" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "\u001b[1m\u001b[34mINFO: binning = 20\n", "\u001b[0m" ] }, { "data": { "text/html": [ "" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "histogram(\n", "\n", "δE/eps(Float64),\n", "title = L\"\\mathrm{Distribution\\;\\;of\\;\\;energy\\;\\;relative\\;\\;error}\", \n", "xaxis = (L\"\\delta E\"),\n", "yaxis = (L\"\\mathrm{Number\\;\\;of\\;\\;}\\delta E\\mathrm{\\;\\;values\\;\\;within\\;\\;bin\\;\\;range}\"),\n", "leg = false,\n", "nbins=20\n", "\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The energy error behaves roughly as a random walk, which means that the numerical error in the energy is dominated by rounding errors." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Angular momentum conservation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Analogously to the energy, we now focus on the angular momentum, which is preserved by the analytical solution too. The value of the angular momentum is given by\n", "$$\n", "\\begin{align}\n", "L &= xv-yu\n", "\\end{align}\n", "$$\n", "We write the angular momentum function as:" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "ang_mom (generic function with 1 method)" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ang_mom(x_, y_, u_, v_) = x_.*v_-y.*u_" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So the angular momentum during each time-step is:" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [], "source": [ "L = ang_mom(x, y, u, v);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$L_0=x_0v_0-y_0u_0$, the initial value of the angular momentum, is:" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.5999999999999999" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "L0 = L[1]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We define $\\delta L$ as the relative error in the angular momentum; i.e.:\n", "\n", "$$\n", "\\begin{align}\n", "\\delta L(t) &= \\frac{L(t)-L_0}{L_0}\n", "\\end{align}\n", "$$" ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "collapsed": true }, "outputs": [], "source": [ "δL = (L-L0)/L0;" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Just as we did for the energy, we now plot $\\delta L$ in units of `eps(Float64)` vs $t$:" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "plot(\n", "\n", "t/(2π), δL/eps(Float64),\n", "title = L\"\\mathrm{Angular\\;\\;momentum\\;\\; relative\\;\\; error\\;\\; vs\\;\\; time}\", \n", "xaxis = (L\"t\\mathrm{\\,(orbital\\;\\;periods)}\"),\n", "yaxis = (L\"\\delta L \\mathrm{\\;\\;(machine\\;\\;epsilons)}\"),\n", "label = L\"\\mathrm{Angular\\;\\;momentum\\;\\; relative\\;\\; error\\;\\; vs\\;\\; time}\",\n", "leg=false,\n", "color=:orange\n", "\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Again, we see that the angular momentum relative error is comparable to `eps(Float64)`; the maximum variation is ~$150$ machine epsilons. How does the distribution of this error look like?" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "\u001b[1m\u001b[34mINFO: binning = 20\n", "\u001b[0m" ] }, { "data": { "text/html": [ "" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "histogram(\n", "\n", "δL/eps(Float64),\n", "title = L\"\\mathrm{Distribution\\;\\;of\\;\\;angular\\;\\;momentum\\;\\;relative\\;\\;error}\", \n", "xaxis = (L\"\\delta L\"),\n", "yaxis = (L\"\\mathrm{Number\\;\\;of\\;\\;}\\delta L\\mathrm{\\;\\;values\\;\\;within\\;\\;bin\\;\\;range}\"),\n", "leg = false,\n", "nbins=20,\n", "color=:orange\n", "\n", ")\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The distribution of $\\delta L$ shows a peak near zero and roughly symmetrical around this value. This means that the error in the angular momentum is dominated also by rounding errors of floating-point arithmetic." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Lastly, we reproduce __exactly__ the last plot shown in the original version of this example, authored by Luis Benet and included in `JuliaDiff/TaylorSeries.jl`'s [Kepler problem integration example](http://nbviewer.jupyter.org/github/JuliaDiff/TaylorSeries.jl/blob/master/examples/1-KeplerProblem.ipynb) jupyter notebook." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A $\\delta E$, $\\delta L$ plot vs $t$:" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ] }, "execution_count": 21, "metadata": {}, "output_type": "execute_result" } ], "source": [ "plot(\n", "\n", "[t/(2π), t/(2π)], [δE/eps(Float64), δL/eps(Float64)],\n", "title = L\"\\delta E\\mathrm{\\;\\;(blue),\\;\\;}\\delta L\\mathrm{\\;\\;(green)\\;\\;vs\\;\\;time}\", \n", "xaxis = (L\"t\\mathrm{\\,(orbital\\;\\;periods)}\"),\n", "yaxis = (L\"\\delta E\\mathrm{,\\;\\;}\\delta L\\;\\;\\mathrm{(machine\\;\\;epsilons)}\"),\n", "leg=false,\n", "\n", ")" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Julia 0.5.3-pre", "language": "julia", "name": "julia-0.5" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "0.5.3" } }, "nbformat": 4, "nbformat_minor": 1 }