{ "cells": [ { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "\n", "# `TaylorIntegration.jl`\n", "\n", "### Taylor's method for ODE initial value problems in Julia\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Jorge A. Pérez and Luis Benet\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "-" } }, "source": [ "### Instituto de Ciencias Físicas\n", "### Universidad Nacional Autónoma de México\n" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "-" } }, "source": [ "**Thanks:** DGAPA-UNAM: IG-100616, and CONACYT." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# Motivation" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "-" } }, "source": [ "Taylor's method provides a *high accuracy method* for the integration of initial value problems. \n", "\n", "The error in each step of the integration can be made close to the round-off error." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "### 99942 Apophis \n", "\n", "![Apophis](2004MN4_Sormano.gif)\n", "\n", "Credit: Osservatorio Astronomico Sormano via [Wikipedia](https://en.wikipedia.org/wiki/99942_Apophis)\n", "\n", "**Semi-major axis**: 0.92228 AU; **eccentricity**: 0.19108\n", "\n", "**Aphelion**: 1.09851 AU; **perihelion**: 0.74605 AU\n", "\n", "**Important close approach** April 13th, 2029." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# Taylor's method" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "-" } }, "source": [ "Consider the initial value problem\n", "\n", "\\begin{equation*}\n", "\\dot{x} = f(x, t), \\qquad x(0) = x_0.\n", "\\end{equation*}" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "-" } }, "source": [ "Approximate locally the solution by a (Taylor) polynomial in $t$ \n", "\n", "\\begin{equation*}\n", "x(t) = \\sum_{n=0}\\, x_{[n]}\\, t^n.\n", "\\end{equation*}" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "Using the Taylor expansion of $f(x(t),t)$ around $t=0$\n", "\n", "\\begin{equation*}\n", "f(x(t), t) = \\sum_{n=0}\\, f_{[n]}\\, t^n,\n", "\\end{equation*}\n", "\n", "we arrive at the **recursion relation**\n", "\n", "\\begin{equation*}\n", "x_{[n+1]} = \\frac{f_{[n]}}{n+1}, \\qquad x_{[0]} = x_0.\n", "\\end{equation*}\n", "\n", "We use [TaylorSeries.jl](https://github.com/JuliaDiff/TaylorSeries.jl) to obtain the expansion of $f(x(t),t)$ around a given value of $t$." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "To assure convergence of the series, the coefficients $x_{[n]}$ **must** be computed to a *large enough* order.\n", "\n", "The time step is determined through the last computed coefficients, $ | x_{[n]} h^n | \\leq \\epsilon$, then\n", "\n", "\\begin{equation*}\n", "h \\leq \\min\\Big\\{ \\Big(\\frac{\\epsilon}{x_{[n]}}\\Big)^{1/n}, \\,\\Big(\\frac{\\epsilon}{x_{[n-1]}}\\Big)^{1/(n-1)}\\Big\\}.\n", "\\end{equation*}\n", "\n", "Note that higher order of the Taylor expansion provides a larger time step." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "## A blackboard example\n", "\n", "Let's consider the equation \n", "\n", "\\begin{equation*}\n", "\\dot{x} = -x, \\qquad x(0) = 1\n", "\\end{equation*}\n", "\n", "whose solution is $x(t) = \\exp(-t)$.\n", "\n", "Clearly, $x_{[0]}=1$ and $f_{[n]} = -x_{[n]}$. Then,\n", "\n", "\\begin{eqnarray*}\n", "x_{[k+1]} & = & \\frac{f_{[k]}}{k+1} = -\\frac{x_{[k]}}{k+1} = (-1)^2\\frac{x_{[k-1]}}{k (k+1)}\\\\\n", "& = & \\dots = (-1)^{k+1} \\frac{x_{[0]}}{(k+1)!} = \\frac{(-1)^{k+1}}{(k+1)!},\n", "\\end{eqnarray*}\n", "\n", "which is the (k+1)-Taylor coefficient of $x(t) = \\exp(-t)\\,$.\n" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# Jet transport" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "*Jet transport* is a tool that allows the propagation under the flow of a small neighborhood in phase space around a given initial condition, instead of propagating single initial conditions only.\n", "\n", "To compute the propagation of $\\mathbf{x}_0+\\delta\\mathbf{x}$, where $\\delta\\mathbf{x}$ are independent small displacements in phase space, one has to solve high-order variational equations. The idea is to treat $\\mathbf{x}_0+\\delta\\mathbf{x}$ as a (truncated) polynomial in the $\\delta\\mathbf{x}$ variables.\n" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "-" } }, "source": [ "Jet transport works with *any* ordinary differential equations integrator, provided the computations can be done using *polynomial algebra*." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "## An example\n", "\n", "\n", "From D. Pérez-Palau et al, *Celest. Mech, Dyn. Astron.* **123**, 239-262 (2015)." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "-" } }, "source": [ "Let us consider the differential equations for the harmonic oscillator:\n", "\n", "\\begin{eqnarray*}\n", "\\dot{x} & = & y, \\\\ \n", "\\dot{y} & = & -x, \n", "\\end{eqnarray*}\n", "\n", "with the initial condition $\\mathbf{x}_0=[x_0, y_0]^T$." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "Euler's method corresponds to:\n", "\n", "\\begin{equation*}\n", "\\mathbf{x}_{n+1} = \\mathbf{x}_n + h \\,\\mathbf{f}(\\mathbf{x}_n).\n", "\\end{equation*}\n", "\n", "Instead of the initial condition, we consider the polynomial\n", "\n", "\\begin{equation*}\n", "P_{0,\\mathbf{x}_0}\\,(\\delta\\mathbf{x}) = [x_0+\\delta x, y_0 + \\delta y]^T.\n", "\\end{equation*}" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "Then,\n", "\n", "\\begin{eqnarray*}\n", "\\mathbf{x}_1 &=& P_{h, \\mathbf{x}_0}(\\delta\\mathbf{x}) = \n", "P_{0,\\mathbf{x}_0}(\\delta\\mathbf{x}) + h\\, \\mathbf{f}(P_{0,\\mathbf{x}_0}(\\delta\\mathbf{x}))\\\\ & = & \n", "\\left(\n", "\\begin{array}{c}\n", "x_0 + h y_0 \\\\ y_0 - h x_0\n", "\\end{array}\n", "\\right) + \\left(\n", "\\begin{array}{cc}\n", "1 & h \\\\ -h & 1\n", "\\end{array}\n", "\\right) \\left(\n", "\\begin{array}{c}\n", "\\delta x\\\\ \\delta y\n", "\\end{array}\n", "\\right).\n", "\\end{eqnarray*}" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "At each step of Taylor's method, we can write\n", "\n", "\\begin{equation*}\n", "\\phi(t_n+h; \\mathbf{x}_n) = \\sum_{i=0}^p \\mathbf{x}^{(i)}(\\mathbf{x}_n, t_n) \\,h^i.\n", "\\end{equation*}\n" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# The Lorenz system\n", "\n", "\\begin{eqnarray*}\n", "\\dot{x_1} & = & σ (x_2-x_1), \\\\\n", "\\dot{x_2} & = & x_1 (ρ-x_3)-x_2, \\\\\n", "\\dot{x_3} & = & x_1 x_2-β x_3.\n", "\\end{eqnarray*}" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "\u001b[1m\u001b[36mINFO: \u001b[39m\u001b[22m\u001b[36mRecompiling stale cache file /Users/Jorge/.julia/lib/v0.6/TaylorSeries.ji for module TaylorSeries.\n", "\u001b[39m" ] } ], "source": [ "#Pkg.add(\"TaylorIntegration.jl\")\n", "\n", "# Some experimental things used are here\n", "# Pkg.checkout(\"TaylorIntegration\", \"jp/poincare\")\n", "using TaylorIntegration, BenchmarkTools\n", "\n", "using Plots" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [ { "data": { "text/plain": [ "lorenz! (generic function with 1 method)" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "#Lorenz system parameters and ODE:\n", "const σ, β, ρ = 16.0, 4.0, 45.92\n", "\n", "function lorenz!(t, x, dx)\n", " dx[1] = σ*(x[2]-x[1])\n", " dx[2] = x[1]*(ρ-x[3])-x[2]\n", " dx[3] = x[1]*x[2]-β*x[3]\n", " nothing\n", "end" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [ { "data": { "text/plain": [ "100.0" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "const x0 = [19.0, 20.0, 50.0] #the initial condition\n", "\n", "const t0 = 0.0 #the initial time\n", "\n", "const tmax = 100.0 #final time of integration" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " 0.679648 seconds (3.37 M allocations: 341.860 MiB, 6.05% gc time)\n" ] } ], "source": [ "@time tv, xv = taylorinteg(lorenz!, x0, t0, tmax, 28, 1e-20; \n", " maxsteps=10000);" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [ { "data": { "text/html": [ "\n", "\n", " \n", "

Plotly javascript loaded.

\n", "

To load again call

init_notebook(true)

\n", " " ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
\n", "\n", "\n" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "plotlyjs()\n", "plot(xv[:,1], xv[:,2], xv[:,3])\n", "title!(\"Lorenz attractor\")" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# Lyapunov spectrum of Lorenz system" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " 14.530125 seconds (227.63 M allocations: 15.303 GiB, 23.01% gc time)\n" ] } ], "source": [ "@time tv, xv, λv = \n", " liap_taylorinteg(lorenz!, x0, t0, tmax, 28, 1e-20; maxsteps=2000000);" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "#Plot the Lyapunov spectrum vs time\n", "plot(tv, λv[:,1], label=\"λ₁\"); plot!(tv, λv[:,2], label=\"λ₂\")\n", "plot!(tv, λv[:,3], label=\"λ₃\"); xlabel!(\"time\"); ylabel!(\"λᵢ, i=1,2,3\")\n", "title!(\"Lyapunov spectrum vs time\")" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "In the Lorenz system, it is known that\n", "\n", "\\begin{equation*}\n", "\\sum_{i=1}^3 \\lambda_i = \\textrm{trace}(DF) = -\\sigma -1 -\\beta.\n", "\\end{equation*}" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [ { "data": { "text/plain": [ "3-element Array{Float64,1}:\n", " 1.53183 \n", " -0.00994121\n", " -22.5219 " ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "#the Lyapunov spectrum at the end of the integration\n", "λv[end,:]" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [ { "data": { "text/plain": [ "-21.0" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "lorenztr = -σ - 1 - β" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [ { "data": { "text/plain": [ "-20.99999999999998" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sum(λv[end,:])" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "slideshow": { "slide_type": "skip" } }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "#Σᵢλᵢ should be equal to lorenztr...\n", "plotlyjs()\n", "plot(tv[2:end], \n", " log10.(abs.(λv[2:end,1]+λv[2:end,2]+λv[2:end,3]-lorenztr)),\n", " label=\"log₁₀|Σᵢλᵢ-Tr(J)|\", xscale=:log10)\n", "ylims!(-15,-12)\n", "xlabel!(\"time\")\n", "ylabel!(\"log₁₀|Σᵢλᵢ-Tr(J)|\")\n", "title!(\"Σᵢλᵢ: computed minus expected (log-log)\")" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# Hénon-Heiles Hamiltonian\n", "\n", "\\begin{equation*}\n", "H = \\frac{1}{2}(p_x^2 + p_y^2) + \\frac{1}{2}(x^2+y^2) + (\\frac{x^3}{3} - x y^2).\n", "\\end{equation*}" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [ { "data": { "text/plain": [ "henonheiles! (generic function with 1 method)" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "#the Hénon-Heiles system\n", "function henonheiles!(t, x, dx)\n", " dx[1] = x[3]\n", " dx[2] = x[4]\n", " dx[3] = -x[1]-(x[2]^2-x[1]^2)\n", " dx[4] = -x[2]-2x[2]*x[1]\n", " nothing\n", "end" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "collapsed": true, "slideshow": { "slide_type": "skip" } }, "outputs": [], "source": [ "#select initial energy and initial condition\n", "const E0, x02 = 0.1025, [0.45335, 0.0, 0.0, 0.0]\n", "\n", "#potential and Hamiltonian\n", "VV(x,y) = 0.5*( x^2 + y^2 )+( x*y^2 - x^3/3)\n", "H(x,y,p,q) = 0.5*(p^2+q^2) + VV(x, y)\n", "H(x) = H(x...)\n", "\n", "#select py0 such that E=E0\n", "py(x,E) = sqrt(2(E-VV(x[1],x[2]))-x[3]^2)\n", "function py!(x,E)\n", " mypy = py(x,E)\n", " x[4] = mypy\n", " x\n", "end\n", "py!(x02,E0)\n", "nothing" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " 1.948519 seconds (24.29 M allocations: 2.321 GiB, 16.60% gc time)\n" ] } ], "source": [ "#integrate\n", "tv2 = 0.0:0.1:3000.0\n", "@time xv2 = taylorinteg(henonheiles!, x02, tv2, 25, 1e-20, maxsteps=30000 );" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "plotlyjs()\n", "path3d(xv2[:,1],xv2[:,2],xv2[:,3], ratio=:equal, label=\"(x(t), y(t), pₓ(t))\", width=3, zlab=\"pₓ\")\n", "xlabel!(\"x\")\n", "ylabel!(\"y\")\n", "title!(\"Hénon-Heiles orbit in x-y-pₓ space\")" ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "collapsed": true, "slideshow": { "slide_type": "subslide" } }, "outputs": [], "source": [ "#Is energy preserved during the integration?\n", "\n", "Ev = Array{Float64}(length(tv2))\n", "\n", "for i in eachindex(tv2)\n", " Ev[i] = H(xv2[i,:])\n", "end\n", "\n", "δEv = (Ev-E0)/E0/eps();" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "plot(tv2,δEv, label=\"δE(t)\")\n", "title!(\"Energy variations as a function of time\")\n", "xlabel!(\"time\")\n", "ylabel!(\"δE(t)=(E(t)-E₀)/E₀ (machine-epsilons)\")" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "Using `TaylorIntegration.jl`, we can numerically show that at each time-step this property is preserved when the system is Hamiltonian, essentially to machine accuracy!\n", "\n", "$$\n", "D( g^t(\\boldsymbol{x}_0) )^\\mathrm{T} \\cdot J \\cdot D( g^t(\\boldsymbol{x}_0) ) = J\n", "$$" ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "slideshow": { "slide_type": "skip" } }, "outputs": [ { "data": { "text/plain": [ "vareqs_taylorinteg (generic function with 1 method)" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "using TaylorSeries\n", "\n", "function vareqs_taylorinteg{T<:Number}(f, q0::Array{T,1}, t0::T, tmax::T,\n", " order::Int, abstol::T; maxsteps::Int=500)\n", " # Allocation\n", " const tv = Array{T}(maxsteps+1)\n", " dof = length(q0)\n", " const xv = Array{T}(dof, maxsteps+1)\n", " const jt = eye(T, dof)\n", " const vT = zeros(T, order+1)\n", " vT[2] = one(T)\n", "\n", " # NOTE: This changes GLOBALLY internal parameters of TaylorN\n", " global _δv = set_variables(\"δ\", order=1, numvars=dof)\n", "\n", " # Initial conditions\n", " @inbounds tv[1] = t0\n", " for ind in 1:dof\n", " @inbounds xv[ind,1] = q0[ind]\n", " end\n", " const x0 = vcat(q0, reshape(jt, dof*dof))\n", " nx0 = dof*(dof+1)\n", " t00 = t0\n", "\n", " # Initialize the vector of Taylor1 expansions\n", " const x = Array{Taylor1{T}}(nx0)\n", " for i in eachindex(x0)\n", " @inbounds x[i] = Taylor1( x0[i], order )\n", " end\n", "\n", " #Allocate auxiliary arrays\n", " const dx = Array{Taylor1{T}}(nx0)\n", " const xaux = Array{Taylor1{T}}(nx0)\n", " const δx = Array{TaylorN{Taylor1{T}}}(dof)\n", " const dδx = Array{TaylorN{Taylor1{T}}}(dof)\n", " const jac = Array{Taylor1{T}}(dof,dof)\n", " for i in eachindex(jac)\n", " @inbounds jac[i] = zero(x[1])\n", " end\n", "\n", " #auxiliary arrays for symplectic structure tests\n", " const δSv = Array{T}(maxsteps+1); δSv[1] = zero(T)\n", " auxJn = Int(dof/2)\n", " const J_n = vcat( hcat(zeros(auxJn,auxJn), eye(auxJn,auxJn)), hcat(-eye(auxJn,auxJn), zeros(auxJn,auxJn)) )\n", "\n", " # Integration\n", " nsteps = 1\n", " while t0 < tmax\n", " δt = TaylorIntegration.liap_taylorstep!(f, x, dx, xaux, δx, dδx, jac, t0, tmax, x0, order, abstol, vT)\n", " for ind in eachindex(jt)\n", " @inbounds jt[ind] = x0[dof+ind]\n", " end\n", " t0 += δt\n", " tspan = t0-t00\n", " nsteps += 1\n", " @inbounds tv[nsteps] = t0\n", " @inbounds for ind in 1:dof\n", " xv[ind,nsteps] = x0[ind]\n", " end\n", " δSv[nsteps] = norm( jt'*J_n*jt-J_n, Inf)\n", " if nsteps > maxsteps\n", " warn(\"\"\"\n", " Maximum number of integration steps reached; exiting.\n", " \"\"\")\n", " break\n", " end\n", " end\n", "\n", " return view(tv,1:nsteps), view(transpose(xv),1:nsteps,:), view(δSv,1:nsteps)\n", "end" ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " 1.893656 seconds (6.18 M allocations: 426.140 MiB, 13.35% gc time)\n" ] } ], "source": [ "#vareqs_taylorinteg integrates the original system and its 1st-order variational equations\n", "@time tvS, xvS, δSv = vareqs_taylorinteg(henonheiles!, x02, 0.0, 30.0, 25, 1e-20, maxsteps=30000 ); " ] }, { "cell_type": "code", "execution_count": 20, "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [ { "data": { "text/plain": [ "1.6272260836797863e-16" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "#δS = Dgᵗ'*J_4*Dgᵗ-J_4\n", "#the error in the preservation of the symplectic structure \n", "# during each integration time step is at most:\n", "norm(δSv,Inf)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# Poincaré map for the Hénon-Heiles system" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this example, we will construct a Poincaré map associated to $y=0$, $\\dot y>0$ and $E=0.1025$ for the Hénon-Heiles system:" ] }, { "cell_type": "code", "execution_count": 21, "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [ { "data": { "text/plain": [ "g (generic function with 1 method)" ] }, "execution_count": 21, "metadata": {}, "output_type": "execute_result" } ], "source": [ "g(t, x, dx) = x[2] # y=0 section" ] }, { "cell_type": "code", "execution_count": 22, "metadata": { "slideshow": { "slide_type": "skip" } }, "outputs": [ { "data": { "text/plain": [ "g (generic function with 3 methods)" ] }, "execution_count": 22, "metadata": {}, "output_type": "execute_result" } ], "source": [ "#specialized method of g for Taylor1's\n", "function g(t,x::Array{Taylor1{T},1},dx::Array{Taylor1{T},1}) where T<:Real\n", " if x[4][1]>zero( T )\n", " return x[2]\n", " else\n", " return zero(x[1])\n", " end\n", "end\n", "\n", "#specialized method of g for Taylor1{TaylorN}'s\n", "function g(t,x::Array{Taylor1{TaylorN{T}},1},dx::Array{Taylor1{TaylorN{T}},1}) where T<:Real\n", " if x[4][1][1][1]>zero( T )\n", " return x[2]\n", " else\n", " return zero(x[1])\n", " end\n", "end" ] }, { "cell_type": "code", "execution_count": 23, "metadata": { "collapsed": true, "slideshow": { "slide_type": "skip" } }, "outputs": [], "source": [ "tvSv2 = Array{Float64,1}[]\n", "xvSv2 = Array{Float64,2}[]\n", "gvSv2 = Array{Float64,1}[];" ] }, { "cell_type": "code", "execution_count": 24, "metadata": { "scrolled": false, "slideshow": { "slide_type": "-" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " 34.062830 seconds (406.59 M allocations: 43.390 GiB, 25.35% gc time)\n" ] } ], "source": [ "nconds2 = 1000 #number of initial conditions\n", "@time for i in 1:nconds2\n", " rand1 = rand(); rand2 = rand()\n", " x_ini = py!(x02+0.005*[sqrt(rand1)*cos(2pi*rand2),0.0,sqrt(rand1)*sin(2pi*rand2),0.0],E0)\n", " tv_i, xv_i, tvS_i, xvS_i, gvS_i = TaylorIntegration.poincare(henonheiles!, g, x_ini, 0.0, 135.0, 25, 1e-25, maxsteps=30000);\n", " push!(tvSv2, vcat(0.0, tvS_i))\n", " push!(xvSv2, vcat(x_ini', xvS_i))\n", " push!(gvSv2, vcat(0.0, gvS_i))\n", "end" ] }, { "cell_type": "code", "execution_count": 25, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "\u001b[1m\u001b[36mINFO: \u001b[39m\u001b[22m\u001b[36mSaved animation to /Users/Jorge/TaylorIntegration.jl/examples/JuliaCon2017/henonheilespoincaremap5.gif\n", "\u001b[39m" ] }, { "data": { "text/html": [ "\" />" ], "text/plain": [ "Plots.AnimatedGif(\"/Users/Jorge/TaylorIntegration.jl/examples/JuliaCon2017/henonheilespoincaremap5.gif\")" ] }, "execution_count": 25, "metadata": {}, "output_type": "execute_result" } ], "source": [ "pyplot()\n", "poincareani5 = @animate for i=1:21\n", " scatter(map(x->x[i,1], xvSv2), map(x->x[i,3], xvSv2), label=\"$(i-1)-th iterate\", m=(1,stroke(0)), ratio=:equal)\n", " xlims!(0.08,0.48)\n", " ylims!(-0.13,0.13)\n", " xlabel!(\"x\")\n", " ylabel!(\"pₓ\")\n", " title!(\"Hénon-Heiles Poincaré map near a period 5 orbit\")\n", "end\n", "gif(poincareani5, \"./henonheilespoincaremap5.gif\", fps = 2)" ] }, { "cell_type": "code", "execution_count": 26, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "\u001b[1m\u001b[36mINFO: \u001b[39m\u001b[22m\u001b[36mSaved animation to /Users/Jorge/TaylorIntegration.jl/examples/JuliaCon2017/henonheilespoincaremap5z.gif\n", "\u001b[39m" ] }, { "data": { "text/html": [ "\" />" ], "text/plain": [ "Plots.AnimatedGif(\"/Users/Jorge/TaylorIntegration.jl/examples/JuliaCon2017/henonheilespoincaremap5z.gif\")" ] }, "execution_count": 26, "metadata": {}, "output_type": "execute_result" } ], "source": [ "pyplot()\n", "poincareani5z = @animate for i=0:4\n", " scatter(map(x->x[5i+1,1], xvSv2), map(x->x[5i+1,3], xvSv2), label=\"$(5i)-th iterate\", m=(5,stroke(0)))\n", " xlims!(0.4475,0.45875)\n", " ylims!(-0.02,0.02)\n", " xlabel!(\"x\")\n", " ylabel!(\"pₓ\")\n", " title!(\"Poincaré map near a period 5 orbit (zoom)\")\n", "end\n", "gif(poincareani5z, \"./henonheilespoincaremap5z.gif\", fps = 2)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "## Hénon-Heiles: Poincaré map using jet transport techniques" ] }, { "cell_type": "code", "execution_count": 27, "metadata": { "collapsed": true, "slideshow": { "slide_type": "-" } }, "outputs": [], "source": [ "ξ = set_variables(\"ξ\", numvars=length(x02), order=3)\n", "x0TN = x02+ξ\n", "py!(x0TN,E0);" ] }, { "cell_type": "code", "execution_count": 28, "metadata": { "scrolled": true, "slideshow": { "slide_type": "-" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " 29.846923 seconds (411.18 M allocations: 27.331 GiB, 27.19% gc time)\n" ] } ], "source": [ "@time tvTN, xvTN, tvSTN, xvSTN, gvSTN = \n", "TaylorIntegration.poincare(henonheiles!, g, x0TN, 0.0, 135.0, 25, 1e-20, maxsteps=30000);" ] }, { "cell_type": "code", "execution_count": 29, "metadata": { "slideshow": { "slide_type": "skip" } }, "outputs": [], "source": [ "#some auxiliaries:\n", "xvSTNaa = Array{Array{TaylorN{Float64},1}}( length(tvSTN)+1 );\n", "xvSTNaa[1] = x0TN\n", "for ind in 2:length(tvSTN)+1\n", " whatever = xvSTN[ind-1,:]\n", " xvSTNaa[ind] = whatever\n", "end\n", "tvSTNaa = union([zero(tvSTN[1])],tvSTN);\n", "#length(tvSTNaa), length(xvSTNaa)" ] }, { "cell_type": "code", "execution_count": 30, "metadata": { "slideshow": { "slide_type": "skip" } }, "outputs": [], "source": [ "myrnd = 0:0.01:1 #rand(100)\n", "myrnd2 = 0:0.01:1 #rand(100)\n", "npoints = length(myrnd)\n", "ncrosses = length(tvSTN)\n", "xS = Array{Float64}(ncrosses+1,npoints)\n", "pS = Array{Float64}(ncrosses+1,npoints)\n", "HS = Array{Float64}(ncrosses+1,npoints)\n", "gS = Array{Float64}(ncrosses+1,npoints)\n", "\n", "myrad=0.005\n", "\n", "#myrnd = rand(npoints)\n", "\n", "ξx = myrad*cos.(2pi*myrnd)\n", "ξp = myrad*sin.(2pi*myrnd)\n", "#ξx = myrad.*(sqrt.(myrnd2)).*cos.(2pi*myrnd)\n", "#ξp = myrad.*(sqrt.(myrnd2)).*sin.(2pi*myrnd)\n", "#ξx = -myrad+2myrad*rand(npoints)\n", "#ξp = zeros(npoints) #-myrad+2myrad*rand(npoints)\n", "\n", "for indpoint in 1:npoints\n", " xS[1,indpoint] = x02[1]+ξx[indpoint]\n", " pS[1,indpoint] = x02[3]+ξp[indpoint]\n", " mycond = [ξx[indpoint], 0.0, ξp[indpoint], 0.0]\n", " HS[1,indpoint] = evaluate(H(xvSTNaa[1]),mycond)\n", " gS[1,indpoint] = evaluate(g(tvSTNaa[1],xvSTNaa[1],xvSTNaa[1]),mycond)\n", " for indS in 2:ncrosses+1 \n", " temp = evaluate(xvSTNaa[indS],mycond)\n", " xS[indS,indpoint] = temp[1]\n", " pS[indS,indpoint] = temp[3]\n", " HS[indS,indpoint] = evaluate(H(xvSTNaa[indS]),mycond)\n", " gS[indS,indpoint] = evaluate(g(tvSTNaa[indS],xvSTNaa[indS],xvSTNaa[indS]),mycond)\n", " end\n", "end" ] }, { "cell_type": "code", "execution_count": 31, "metadata": { "scrolled": false, "slideshow": { "slide_type": "subslide" } }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "\u001b[1m\u001b[36mINFO: \u001b[39m\u001b[22m\u001b[36mSaved animation to /Users/Jorge/TaylorIntegration.jl/examples/JuliaCon2017/henonheilespoincaremap5jt.gif\n", "\u001b[39m" ] }, { "data": { "text/html": [ "\" />" ], "text/plain": [ "Plots.AnimatedGif(\"/Users/Jorge/TaylorIntegration.jl/examples/JuliaCon2017/henonheilespoincaremap5jt.gif\")" ] }, "execution_count": 31, "metadata": {}, "output_type": "execute_result" } ], "source": [ "pyplot()\n", "poincareani5jt = @animate for i=1:21\n", " plot(xS[i,:],pS[i,:], ratio=:equal, label=\"$(i-1)-th iterate\")\n", " xlims!(0.08,0.48)\n", " ylims!(-0.13,0.13)\n", " xlabel!(\"x\")\n", " ylabel!(\"pₓ\")\n", " title!(\"Poincaré map with jet transport near a period 5 orbit\")\n", "end\n", "gif(poincareani5jt, \"./henonheilespoincaremap5jt.gif\", fps = 2)" ] }, { "cell_type": "code", "execution_count": 32, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "\u001b[1m\u001b[36mINFO: \u001b[39m\u001b[22m\u001b[36mSaved animation to /Users/Jorge/TaylorIntegration.jl/examples/JuliaCon2017/poincareani5zJTvMC.gif\n", "\u001b[39m" ] }, { "data": { "text/html": [ "\" />" ], "text/plain": [ "Plots.AnimatedGif(\"/Users/Jorge/TaylorIntegration.jl/examples/JuliaCon2017/poincareani5zJTvMC.gif\")" ] }, "execution_count": 32, "metadata": {}, "output_type": "execute_result" } ], "source": [ "pyplot()\n", "poincareani5zJTvMC = @animate for i=0:4\n", " scatter(map(x->x[5i+1,1], xvSv2), map(x->x[5i+1,3], xvSv2), m=(4,stroke(0)),label=\"$(5i)-th it. MC\")\n", " plot!(xS[5i+1,:],pS[5i+1,:], width=:3,label=\"$(5i)-th it. 3rd-order Jet transport\")\n", " xlims!(0.4475,0.45875)\n", " ylims!(-0.02,0.02)\n", " xlabel!(\"x\")\n", " ylabel!(\"pₓ\")\n", " title!(\"Poincaré map: 3rd-order jet transport vs Monte Carlo\")\n", "end\n", "gif(poincareani5zJTvMC, \"./poincareani5zJTvMC.gif\", fps = 2)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# Apophis close-approach to Earth on 2029\n", "\n", "In this example, we will show a numerical integration of asteroid Apophis orbit using `TaylorIntegration.jl`. This asteroid will have a close-approach to on 2029. Our model includes the following interactions:\n", "\n", "1. Massive bodies: Sun, eight planets, Moon and dwarf planet Ceres\n", "2. Newtonian accelerations, plus relativistic corrections\n", "3. Effects due to polar flattening of Earth and Sun\n", "4. Apophis non-gravitational accelerations\n", "5. Earth's axial precession" ] }, { "cell_type": "code", "execution_count": 33, "metadata": { "collapsed": true, "slideshow": { "slide_type": "skip" } }, "outputs": [], "source": [ "const celobj_symbols=(\"s\",\"☿\", \"♀\", \"e\", \"m\", \"♂\", \"♃\", \"♄\",\"♅\", \"♆\", \"Ce\", \"Ap\")\n", "const N = length(celobj_symbols)\n", "const au = 149597870.700\n", "const ld = 384400.0/au\n", "const RE=6378.1363/au\n", "const su = 1 #Sun's index\n", "const ea = 4 #Earth's index\n", "const mo = 5 #Moon's index\n", ";" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "skip" } }, "source": [ "We are interested in calculating the time and distance of Apophis at each of its close approaches to Earth. That is, let $d_{EA}(t)$ be the Earth-Apophis distance at time $t$:\n", "\n", "$$\n", "d_{EA}(t) = \\sqrt{(x_{E}(t)-x_{A}(t))^2+(y_{E}(t)-y_{A}(t))^2+(z_{E}(t)-z_{A}(t))^2}\n", "$$\n", "\n", "Then, we are looking to find the values $t^\\star$ such that\n", "\n", "$$\n", "\\frac{d}{dt} d_{EA}(t^\\star)=0\n", "$$" ] }, { "cell_type": "code", "execution_count": 34, "metadata": { "slideshow": { "slide_type": "skip" } }, "outputs": [ { "data": { "text/plain": [ "d_EA2 (generic function with 1 method)" ] }, "execution_count": 34, "metadata": {}, "output_type": "execute_result" } ], "source": [ "#the square of the Earth-Apophis distance\n", "d_EA(t,x,dx) = (x[3N-2]-x[3ea-2])^2+(x[3N-1]-x[3ea-1])^2+(x[3N]-x[3ea])^2\n", "\n", "#the derivative of the squared Earth-Apophis distance\n", "d_EA2(t,x,dx) = (x[3N-2]-x[3ea-2])*(x[6N-2]-x[3(N+ea)-2])+(x[3N-1]-x[3ea-1])*(x[6N-1]-x[3(N+ea)-1])+(x[3N]-x[3ea])*(x[6N]-x[3(N+ea)])" ] }, { "cell_type": "code", "execution_count": 35, "metadata": { "slideshow": { "slide_type": "skip" } }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "WARNING: redefining constant t0\n", "WARNING: redefining constant tmax\n" ] }, { "data": { "text/plain": [ "(30.78713210130048, 2007-12-07T00:00:00, 2038-09-20T00:00:00)" ] }, "execution_count": 35, "metadata": {}, "output_type": "execute_result" } ], "source": [ "const order = 30\n", "const abstol = 1.0E-30\n", "const T = 365.25\n", "const t0 = Dates.datetime2julian(DateTime(2007,12,7,0,0,0))\n", "const tmax = t0+30.7871321013T\n", "(tmax-t0)/T, Dates.julian2datetime(t0), Dates.julian2datetime(tmax)" ] }, { "cell_type": "code", "execution_count": 36, "metadata": { "collapsed": true, "slideshow": { "slide_type": "skip" } }, "outputs": [], "source": [ "#Un-comment lines below if instead of running the integrations, the solution is being read from the saved files:\n", "\n", "qN = readdlm(\"data/Apophis.dat\", ' ')\n", "\n", "qN_S = readdlm(\"data/Apophis_S.dat\", ' ')\n", "\n", "#times of evaluation of solution:\n", "tv = qN[:,1]\n", "t = (tv -t0)/T\n", "#dynamical state at each solution evaluation:\n", "xv = qN[:,2:end]\n", "\n", "#times of passage through surface of section:\n", "tvS = qN_S[:,1]\n", "tS = (tvS -t0)/T\n", "\n", "#dynamical state at surface of section crossings:\n", "xvS = qN_S[:,2:end-1]\n", "\n", "#residuals of Newton-Raphson iteration for each crossing:\n", "gvS = qN_S[:,end];" ] }, { "cell_type": "code", "execution_count": 37, "metadata": { "slideshow": { "slide_type": "skip" } }, "outputs": [ { "data": { "text/plain": [ "(true, true)" ] }, "execution_count": 37, "metadata": {}, "output_type": "execute_result" } ], "source": [ "#=\n", "Make solution vector more user-friendly: for a given suffix, e.g. \"σσ\", xσσ, yσσ, zσσ are arrays \n", "which represent the position of the body σσ at each time step, with the corresponding suffix defined \n", "by celobj_symbols. Similarly, uσσ , vσσ , wσσ represent the velocities. E.g., for the Earth: xe, ye, ze, \n", "etc. For Jupiter: x♃, y♃, z♃, etc...\n", "=#\n", "\n", "i=1\n", "for op = celobj_symbols\n", " for var = (:x, :y, :z)\n", " ex = Symbol(var,op)\n", " @eval $ex = xv[:,$i]\n", " # println(i,\" \",ex)\n", " i+=1\n", " end\n", "end\n", "\n", "j=1 #this index runs over number-of-body\n", "\n", "for op = celobj_symbols\n", " for var = (:u, :v, :w)\n", " #`xv[:,$i]` represents a canonical momentum.\n", " #In order to get the velocity, we must divide by the mass\n", " #when integrating backwards in time, a minus sign\n", " #must be added for the velocities: $ex = xv[:,$i] ---> $ex = -xv[:,$i], etc.\n", " ex = Symbol(var,op)\n", " @eval $ex = xv[:,$i]\n", " i+=1\n", " end\n", " j+=1\n", "end\n", "#if everything was done right index-wise, both expressions below should evaluate to zero:\n", "i == 6N+1, j == N+1" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "# The orbit of Apophis in space" ] }, { "cell_type": "code", "execution_count": 38, "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n" ] }, "execution_count": 38, "metadata": {}, "output_type": "execute_result" } ], "source": [ "plotlyjs()\n", "plot(xAp-xs,yAp-xs,zAp-xs, label=\"Apophis\")\n", "plot!(xe-xs,ye-xs,ze-xs, label=\"Earth\")\n", "plot!(x♀-xs,y♀-xs,z♀-xs, label=\"Venus\")\n", "plot!(x☿-xs,y☿-xs,z☿-xs, label=\"Mercury\")\n", "scatter!([0], [0], [0], markercolor=:yellow, markerstrokewidth = 0.1, markersize=2.0, label=\"Sun\")\n", "title!(\"Heliocentric orbits: Inner S.S. + Apophis\")\n", "xlabel!(\"X axis\")\n", "ylabel!(\"Y axis\")" ] }, { "cell_type": "code", "execution_count": 39, "metadata": { "scrolled": false, "slideshow": { "slide_type": "subslide" } }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n" ] }, "execution_count": 39, "metadata": {}, "output_type": "execute_result" } ], "source": [ "plot(xAp-xe,yAp-ye,zAp-ze, label=\"Apophis\")\n", "plot!(xm-xe,ym-ye,zm-ze, label=\"Moon\")\n", "scatter!([0], [0], [0], markercolor=:blue, markerstrokewidth = 0, markersize=0.5, label=\"Earth\")\n", "title!(\"Moon and Apophis geocentric orbits\")\n", "xlabel!(\"X axis\")\n", "ylabel!(\"Y axis\")" ] }, { "cell_type": "code", "execution_count": 40, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n" ] }, "execution_count": 40, "metadata": {}, "output_type": "execute_result" } ], "source": [ "plot(xAp-xe,yAp-ye,zAp-ze, label=\"Apophis\")\n", "plot!(xm-xe,ym-ye,zm-ze, label=\"Moon\")\n", "scatter!([0], [0], [0], markercolor=:blue, markerstrokewidth = 0, markersize=2, label=\"Earth\")\n", "xlims!(-10ld,10ld)\n", "ylims!(-10ld,10ld)\n", "zlims!(-10ld,10ld)\n", "title!(\"Moon and Apophis geocentric orbits (zoom)\")\n", "xlabel!(\"X axis\")\n", "ylabel!(\"Y axis\")" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "## Apophis close-encounter with Earth in 2013: comparison vs radar observations" ] }, { "cell_type": "code", "execution_count": 41, "metadata": { "collapsed": true, "slideshow": { "slide_type": "skip" } }, "outputs": [], "source": [ "# range data from NEODyS\n", "neodys_range_data = [15391792.43632,\n", "14605784.93200,\n", "14526476.23426,\n", "14500639.40073,\n", "14457708.50168,\n", "14460886.85185,\n", "14582360.35599,\n", "14702834.53945,\n", "23669285.49162,\n", "24520157.18491,\n", "35258719.15715]/au\n", "\n", "# range statistical uncertainties data\n", "neodys_range_error = [0.44969,\n", "0.44969,\n", "0.44969,\n", "0.44969,\n", "0.14990,\n", "0.22484,\n", "0.14990,\n", "0.14990,\n", "0.14990,\n", "0.14990,\n", "0.29979]/au\n", "\n", "# the times reported for the ranging data\n", "neodys_times_table = [\"2012-12-22T11:00:00\",\n", "\"2013-01-03T10:00:00\",\n", "\"2013-01-05T10:50:00\",\n", "\"2013-01-06T08:20:00\",\n", "\"2013-01-09T08:00:00\",\n", "\"2013-01-10T08:00:00\",\n", "\"2013-01-14T09:50:00\",\n", "\"2013-01-16T06:30:00\",\n", "\"2013-02-18T01:36:00\",\n", "\"2013-02-20T01:27:00\",\n", "\"2013-03-15T23:59:00\"]\n", "\n", "#some little processing\n", "neodys_datetime = map(x->Dates.DateTime(x,\"yyyy-mm-ddTHH:MM:SS\"), neodys_times_table);\n", "tv_neodys_2013 = Dates.datetime2julian.(neodys_datetime);\n", "\n", "indrange2013 = 1313:1375\n", "ranges2013 = Array{Float64}(length(indrange2013))\n", "for i in eachindex(indrange2013)\n", " ranges2013[i] = sqrt(d_EA(tv[indrange2013[i]], xv[indrange2013[i],:], [0]))\n", "end\n", "# ranges2013\n", ";" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "skip" } }, "source": [ "Then, we plot the data from the integration, and compare them to radar astrometry data performed at Arecibo, reported in NEODyS. Reported error bars are not discernible at the plotted scale:" ] }, { "cell_type": "code", "execution_count": 42, "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n" ] }, "execution_count": 42, "metadata": {}, "output_type": "execute_result" } ], "source": [ "plot((tv[indrange2013]-t0)/T, ranges2013/ld, label=\"range (integrated)\")\n", "scatter!((tv_neodys_2013-t0)/T, neodys_range_data/ld, label=\"range (radar astrometry)\", markershape=:cross, markerstrokewidth=0)\n", "xlabel!(\"Julian years since Dec 7th, 2007 00:00:00.0 TDB (yr)\")\n", "title!(\"Arecibo 2013 observations vs integration\")\n", "ylabel!(\"Geocentric range (lunar distances, ld)\")" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "## Apophis close-encounter with Earth in 2029: comparison vs JPL prediction" ] }, { "cell_type": "code", "execution_count": 43, "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [ { "data": { "text/plain": [ "(2029-04-13T21:46:12.187, 37964.016066159034)" ] }, "execution_count": 43, "metadata": {}, "output_type": "execute_result" } ], "source": [ "#our prediction\n", "moidApr2029 = 49\n", "\n", "g_star = sqrt(d_EA(tvS[moidApr2029], xvS[moidApr2029,:], xvS[moidApr2029,:]))\n", "\n", "Dates.julian2datetime( tvS[moidApr2029] ), au*g_star" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "-" } }, "source": [ "since Earth's equatorial radius is roughly equal to 6378 km, according to our integration, Apophis will pass at a distance of about 6-1=**5 Earth radii from Earth's surface in 2029!!!**" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "-" } }, "source": [ "||JPL HORIZONS|`TaylorIntegration.jl`|\n", "|-|-|-|\n", "|**time**|April 13, 2029, 21:46:12.012 TDB|April 13, 2029, 21:46:12.187 (TDB)|\n", "|**min. distance**|37,952±255 km|37,964 km|" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "skip" } }, "source": [ "Now, we plot the Apophis-Earth distance versus time:" ] }, { "cell_type": "code", "execution_count": 44, "metadata": { "collapsed": true, "slideshow": { "slide_type": "skip" } }, "outputs": [], "source": [ "d_Ap_Ea = sqrt.((xAp-xe).^2+(yAp-ye).^2+(zAp-ze).^2);\n", "d_Ap_Ea_S = sqrt.( (xvS[:,3N-2].-xvS[:,3ea-2]).^(2).+(xvS[:,3N-1].-xvS[:,3ea-1]).^(2).+(xvS[:,3N].-xvS[:,3ea]).^2 );" ] }, { "cell_type": "code", "execution_count": 45, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n" ] }, "execution_count": 45, "metadata": {}, "output_type": "execute_result" } ], "source": [ "plot(t, d_Ap_Ea/ld, label=\"d(Ea,Ap)\")\n", "scatter!(tS[1:2:end] , d_Ap_Ea_S[1:2:end]/ld, label=\"MOID\", markershape=:cross, markerstrokewidth=0)\n", "plot!([t[1],tS[moidApr2029],t[end]],[1,1,1], label=\"Moon's orbit\")\n", "plot!([t[1],tS[moidApr2029],t[end]],[RE/ld,RE/ld,RE/ld], label=\"Earth's surface\")\n", "xlabel!(\"time (Julian years since Dec 7, 2007)\")\n", "ylabel!(\"Apophis geocentric range (lunar distances)\")\n", "title!(\"Apophis 2029 close-approach to Earth\")\n", "ylims!(0,1000)\n", "xaxis!((0,30),0:2:30)" ] }, { "cell_type": "code", "execution_count": 46, "metadata": { "slideshow": { "slide_type": "skip" } }, "outputs": [], "source": [ "midnight2029Apr13 = Dates.datetime2julian(Dates.DateTime(2029,4,13,0,0,0))\n", "tv_since_midnight2029Apr13 = tv-midnight2029Apr13\n", "t_since_midnight2029Apr13 = tv_since_midnight2029Apr13/T\n", "tvS_since_midnight2029Apr13 = tvS-midnight2029Apr13\n", "tS_since_midnight2029Apr13 = tvS_since_midnight2029Apr13/T\n", "\n", "tplotfirst = 0.0\n", "tplotlast = 1.8\n", "tplotstep = 0.25\n", "timesxaxis = string.( Dates.julian2datetime.( midnight2029Apr13+(tplotfirst:tplotstep:tplotlast) ) );" ] }, { "cell_type": "code", "execution_count": 47, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n" ] }, "execution_count": 47, "metadata": {}, "output_type": "execute_result" } ], "source": [ "plot(tv_since_midnight2029Apr13, d_Ap_Ea/ld, label=\"d(Ea,Ap)\")\n", "scatter!(tvS_since_midnight2029Apr13[1:2:end] , d_Ap_Ea_S[1:2:end]/ld, label=\"MOID\", markershape=:cross, markerstrokewidth=0)\n", "plot!([tv_since_midnight2029Apr13[1],tvS_since_midnight2029Apr13[moidApr2029],tv_since_midnight2029Apr13[end]],[1,1,1], label=\"Moon's orbit\")\n", "plot!([tv_since_midnight2029Apr13[1],tvS_since_midnight2029Apr13[moidApr2029],tv_since_midnight2029Apr13[end]],[RE/ld,RE/ld,RE/ld], label=\"Earth's surface\")\n", "xlabel!(\"time (TDB)\")\n", "ylabel!(\"Apophis geocentric range (lunar distances)\")\n", "title!(\"Apophis 2029 close-approach to Earth\")\n", "xaxis!((tplotfirst,tplotlast),(tplotfirst:tplotstep:tplotlast,timesxaxis), xrotation=-30)\n", "yaxis!((0,1.2),0:0.2:1.5)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# Conclusions\n", "\n", "`TaylorIntegration.jl` provides a highly accurate Taylor method to integrate ODEs,\n", "including whole Lyapunov spectrum calculations and jet transport. \n", "\n", "It interacts nicely with other numeric types (`BigFloat`, `ArbFloats`, `Interval`).\n", "\n", "Slow, but not too slow for what it achieves. There is room for performance improvement!" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "## To Do's\n", "\n", "- Add `JuliaDiffEq` common interface bindings ([#19](https://github.com/PerezHz/TaylorIntegration.jl/issues/19))\n", "\n", "- Add nice API for Monte-Carlo calculations\n", "\n", "- Improve performance \n", "\n", "- Poincaré surfaces of section (WIP)\n", "\n", "- Validated numerical integration of ODEs (WIP)\n" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "\n", "# Thank you !" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# Towards validated Taylor integrations" ] }, { "cell_type": "code", "execution_count": 48, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "WARNING: Method definition string(ArbFloats.ArbFloat{P}) in module ArbFloats at /Users/Jorge/.julia/v0.6/ArbFloats/src/basics/string.jl:69 overwritten at /Users/Jorge/.julia/v0.6/ArbFloats/src/basics/string.jl:220.\n" ] } ], "source": [ "using IntervalArithmetic, ArbFloats\n", "@format full;" ] }, { "cell_type": "code", "execution_count": 49, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " 0.981991 seconds (432.90 k allocations: 23.762 MiB, 1.40% gc time)\n" ] } ], "source": [ "@time tvi, xvi = taylorinteg( (t,x)->t^2, 0..0, 0..0, 1..1, 20, 1.e-20); " ] }, { "cell_type": "code", "execution_count": 50, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "true" ] }, "execution_count": 50, "metadata": {}, "output_type": "execute_result" } ], "source": [ "@interval(1/3) ∈ xvi[end]" ] }, { "cell_type": "code", "execution_count": 51, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " 0.496704 seconds (213.06 k allocations: 11.361 MiB)\n" ] } ], "source": [ "@time tva, xva = taylorinteg( (t,x)->t^2, ArbFloat(0), ArbFloat(0), ArbFloat(1), 20, ArbFloat(1.e-20)); " ] }, { "cell_type": "code", "execution_count": 52, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "2-element SubArray{ArbFloats.ArbFloat{116},1,Array{ArbFloats.ArbFloat{116},1},Tuple{UnitRange{Int64}},true}:\n", " 0 \n", " 0.3333333333333333333333333.." ] }, "execution_count": 52, "metadata": {}, "output_type": "execute_result" } ], "source": [ "xva" ] }, { "cell_type": "code", "execution_count": 53, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " 0.171007 seconds (72.40 k allocations: 3.875 MiB)\n" ] } ], "source": [ "@time tvi, xvi = taylorinteg( (t,x)->sin(t), 0..0, 0..0, @interval(pi), 20, 1.e-20); " ] }, { "cell_type": "code", "execution_count": 54, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Interval(1.999999999999998, 2.0000000000000027)" ] }, "execution_count": 54, "metadata": {}, "output_type": "execute_result" } ], "source": [ "xvi[end]" ] }, { "cell_type": "code", "execution_count": 55, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "6" ] }, "execution_count": 55, "metadata": {}, "output_type": "execute_result" } ], "source": [ "length(xvi)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] } ], "metadata": { "celltoolbar": "Slideshow", "kernelspec": { "display_name": "Julia 0.6.1-pre", "language": "julia", "name": "julia-0.6" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "0.6.1" }, "livereveal": { "height": 768, "scroll": true, "slideNumber": false, "transition": "fade", "width": 1024 } }, "nbformat": 4, "nbformat_minor": 2 }