notebook/models/05-kepler_problem.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "source": [
        "The Hamiltonian $\\mathcal {H}$ and the angular momentum $L$ for the Kepler problem are\n\n$$\\mathcal {H} = \\frac{1}{2}(\\dot{q}^2_1+\\dot{q}^2_2)-\\frac{1}{\\sqrt{q^2_1+q^2_2}},\\quad\nL = q_1\\dot{q_2} - \\dot{q_1}q_2$$\n\nAlso, we know that\n\n$${\\displaystyle {\\frac {\\mathrm {d} {\\boldsymbol {p}}}{\\mathrm {d} t}}=-{\\frac {\\partial {\\mathcal {H}}}{\\partial {\\boldsymbol {q}}}}\\quad ,\\quad {\\frac {\\mathrm {d} {\\boldsymbol {q}}}{\\mathrm {d} t}}=+{\\frac {\\partial {\\mathcal {H}}}{\\partial {\\boldsymbol {p}}}}}$$"
      ],
      "metadata": {}
    },
    {
      "outputs": [],
      "cell_type": "code",
      "source": [
        "using OrdinaryDiffEq, LinearAlgebra, ForwardDiff, Plots; gr()\nH(q,p) = norm(p)^2/2 - inv(norm(q))\nL(q,p) = q[1]*p[2] - p[1]*q[2]\n\npdot(dp,p,q,params,t) = ForwardDiff.gradient!(dp, q->-H(q, p), q)\nqdot(dq,p,q,params,t) = ForwardDiff.gradient!(dq, p-> H(q, p), p)\n\ninitial_position = [.4, 0]\ninitial_velocity = [0., 2.]\ninitial_cond = (initial_position, initial_velocity)\ninitial_first_integrals = (H(initial_cond...), L(initial_cond...))\ntspan = (0,20.)\nprob = DynamicalODEProblem(pdot, qdot, initial_velocity, initial_position, tspan)\nsol = solve(prob, KahanLi6(), dt=1//10);"
      ],
      "metadata": {},
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "source": [
        "Let's plot the orbit and check the energy and angular momentum variation. We know that energy and angular momentum should be constant, and they are also called first integrals."
      ],
      "metadata": {}
    },
    {
      "outputs": [],
      "cell_type": "code",
      "source": [
        "plot_orbit(sol) = plot(sol,vars=(3,4), lab=\"Orbit\", title=\"Kepler Problem Solution\")\n\nfunction plot_first_integrals(sol, H, L)\n    plot(initial_first_integrals[1].-map(u->H(u[2,:], u[1,:]), sol.u), lab=\"Energy variation\", title=\"First Integrals\")\n    plot!(initial_first_integrals[2].-map(u->L(u[2,:], u[1,:]), sol.u), lab=\"Angular momentum variation\")\nend\nanalysis_plot(sol, H, L) = plot(plot_orbit(sol), plot_first_integrals(sol, H, L))"
      ],
      "metadata": {},
      "execution_count": null
    },
    {
      "outputs": [],
      "cell_type": "code",
      "source": [
        "analysis_plot(sol, H, L)"
      ],
      "metadata": {},
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "source": [
        "Let's try to use a Runge-Kutta-Nyström solver to solve this problem and check the first integrals' variation."
      ],
      "metadata": {}
    },
    {
      "outputs": [],
      "cell_type": "code",
      "source": [
        "sol2 = solve(prob, DPRKN6())  # dt is not necessary, because unlike symplectic\n                              # integrators DPRKN6 is adaptive\n@show sol2.u |> length\nanalysis_plot(sol2, H, L)"
      ],
      "metadata": {},
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "source": [
        "Let's then try to solve the same problem by the `ERKN4` solver, which is specialized for sinusoid-like periodic function"
      ],
      "metadata": {}
    },
    {
      "outputs": [],
      "cell_type": "code",
      "source": [
        "sol3 = solve(prob, ERKN4()) # dt is not necessary, because unlike symplectic\n                            # integrators ERKN4 is adaptive\n@show sol3.u |> length\nanalysis_plot(sol3, H, L)"
      ],
      "metadata": {},
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "source": [
        "We can see that `ERKN4` does a bad job for this problem, because this problem is not sinusoid-like.\n\nOne advantage of using `DynamicalODEProblem` is that it can implicitly convert the second order ODE problem to a *normal* system of first order ODEs, which is solvable for other ODE solvers. Let's use the `Tsit5` solver for the next example."
      ],
      "metadata": {}
    },
    {
      "outputs": [],
      "cell_type": "code",
      "source": [
        "sol4 = solve(prob, Tsit5())\n@show sol4.u |> length\nanalysis_plot(sol4, H, L)"
      ],
      "metadata": {},
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "source": [
        "#### Note\n\nThere is drifting for all the solutions, and high order methods are drifting less because they are more accurate.\n\n### Conclusion\n\n---\n\nSymplectic integrator does not conserve the energy completely at all time, but the energy can come back. In order to make sure that the energy fluctuation comes back eventually, symplectic integrator has to have a fixed time step. Despite the energy variation, symplectic integrator conserves the angular momentum perfectly.\n\nBoth Runge-Kutta-Nyström and Runge-Kutta integrator do not conserve energy nor the angular momentum, and the first integrals do not tend to come back. An advantage Runge-Kutta-Nyström integrator over symplectic integrator is that RKN integrator can have adaptivity. An advantage Runge-Kutta-Nyström integrator over Runge-Kutta integrator is that RKN integrator has less function evaluation per step. The `ERKN4` solver works best for sinusoid-like solutions.\n\n## Manifold Projection\n\nIn this example, we know that energy and angular momentum should be conserved. We can achieve this through mainfold projection. As the name implies, it is a procedure to project the ODE solution to a manifold. Let's start with a base case, where mainfold projection isn't being used."
      ],
      "metadata": {}
    },
    {
      "outputs": [],
      "cell_type": "code",
      "source": [
        "using DiffEqCallbacks\n\nplot_orbit2(sol) = plot(sol,vars=(1,2), lab=\"Orbit\", title=\"Kepler Problem Solution\")\n\nfunction plot_first_integrals2(sol, H, L)\n    plot(initial_first_integrals[1].-map(u->H(u[1:2],u[3:4]), sol.u), lab=\"Energy variation\", title=\"First Integrals\")\n    plot!(initial_first_integrals[2].-map(u->L(u[1:2],u[3:4]), sol.u), lab=\"Angular momentum variation\")\nend\n\nanalysis_plot2(sol, H, L) = plot(plot_orbit2(sol), plot_first_integrals2(sol, H, L))\n\nfunction hamiltonian(du,u,params,t)\n    q, p = u[1:2], u[3:4]\n    qdot(@view(du[1:2]), p, q, params, t)\n    pdot(@view(du[3:4]), p, q, params, t)\nend\n\nprob2 = ODEProblem(hamiltonian, [initial_position; initial_velocity], tspan)\nsol_ = solve(prob2, RK4(), dt=1//5, adaptive=false)\nanalysis_plot2(sol_, H, L)"
      ],
      "metadata": {},
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "source": [
        "There is a significant fluctuation in the first integrals, when there is no mainfold projection."
      ],
      "metadata": {}
    },
    {
      "outputs": [],
      "cell_type": "code",
      "source": [
        "function first_integrals_manifold(residual,u)\n    residual[1:2] .= initial_first_integrals[1] - H(u[1:2], u[3:4])\n    residual[3:4] .= initial_first_integrals[2] - L(u[1:2], u[3:4])\nend\n\ncb = ManifoldProjection(first_integrals_manifold)\nsol5 = solve(prob2, RK4(), dt=1//5, adaptive=false, callback=cb)\nanalysis_plot2(sol5, H, L)"
      ],
      "metadata": {},
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "source": [
        "We can see that thanks to the manifold projection, the first integrals' variation is very small, although we are using `RK4` which is not symplectic. But wait, what if we only project to the energy conservation manifold?"
      ],
      "metadata": {}
    },
    {
      "outputs": [],
      "cell_type": "code",
      "source": [
        "function energy_manifold(residual,u)\n    residual[1:2] .= initial_first_integrals[1] - H(u[1:2], u[3:4])\n    residual[3:4] .= 0\nend\nenergy_cb = ManifoldProjection(energy_manifold)\nsol6 = solve(prob2, RK4(), dt=1//5, adaptive=false, callback=energy_cb)\nanalysis_plot2(sol6, H, L)"
      ],
      "metadata": {},
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "source": [
        "There is almost no energy variation but angular momentum varies quite bit. How about only project to the angular momentum conservation manifold?"
      ],
      "metadata": {}
    },
    {
      "outputs": [],
      "cell_type": "code",
      "source": [
        "function angular_manifold(residual,u)\n    residual[1:2] .= initial_first_integrals[2] - L(u[1:2], u[3:4])\n    residual[3:4] .= 0\nend\nangular_cb = ManifoldProjection(angular_manifold)\nsol7 = solve(prob2, RK4(), dt=1//5, adaptive=false, callback=angular_cb)\nanalysis_plot2(sol7, H, L)"
      ],
      "metadata": {},
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "source": [
        "Again, we see what we expect."
      ],
      "metadata": {}
    },
    {
      "outputs": [],
      "cell_type": "code",
      "source": [
        "using SciMLTutorials\nSciMLTutorials.tutorial_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])"
      ],
      "metadata": {},
      "execution_count": null
    }
  ],
  "nbformat_minor": 2,
  "metadata": {
    "language_info": {
      "file_extension": ".jl",
      "mimetype": "application/julia",
      "name": "julia",
      "version": "1.4.2"
    },
    "kernelspec": {
      "name": "julia-1.4",
      "display_name": "Julia 1.4.2",
      "language": "julia"
    }
  },
  "nbformat": 4
}