Higher-Order.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Higher-Order Equation of Lagrange in Terms of Differential Forms #\n",
    "\n",
    "We show that the \"two equations of Lagrange\" given by Oziewicz and Ramirez are not independent.\n",
    "The second equation is exactly $v$ times the first equation.\n",
    "\n",
    "This worksheet: <font size=2>[https://cloud.sagemath.com/projects/b04b5777-e269-4c8f-a4b8-b21dbe1c93c6/files/Two%20Equations%20of%20Lagrange.sagews](https://cloud.sagemath.com/projects/b04b5777-e269-4c8f-a4b8-b21dbe1c93c6/files/Two%20Equations%20of%20Lagrange.sagews)</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%display typeset"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider vectors and forms on a 6-dimensional manifold"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "M = Manifold(5,'M')\n",
    "coord.<t, x, v, a, s> = M.chart()\n",
    "f = M.scalar_field(function('f')(*list(coord)))\n",
    "basis = coord.frame()\n",
    "[Dt,Dx,Dv,Da,Ds] = [basis[i] for i in range(M.dim())]\n",
    "cobasis=coord.coframe()\n",
    "[dt,dx,dv,da,ds] = [cobasis[i] for i in range(M.dim())]\n",
    "d=xder\n",
    "def ev(N): return (lambda x: N.contract(x))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## For the most general Lagrangian ##"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\begin{array}{llcl} & M & \\longrightarrow & \\mathbb{R} \\\\ & \\left(t, x, v, a, s\\right) & \\longmapsto & L\\left(t, x, v, a, s\\right) \\end{array}</script></html>"
      ],
      "text/plain": [
       "M --> R\n",
       "(t, x, v, a, s) |--> L(t, x, v, a, s)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "L = M.scalar_field(function('L')(*list(coord))); L.display()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$f$ is assumed to be a dependent variable"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "f = M.scalar_field(var('f'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Kinematics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\frac{\\partial}{\\partial t } + v \\frac{\\partial}{\\partial x } + a \\frac{\\partial}{\\partial v } + s \\frac{\\partial}{\\partial a } + f \\frac{\\partial}{\\partial s }</script></html>"
      ],
      "text/plain": [
       "d/dt + v d/dx + a d/dv + s d/da + f d/ds"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "N = Dt + v*Dx + a*Dv + s*Da + f*Ds; N.display()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\mathrm{True}</script></html>"
      ],
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "N(M.scalar_field(t))==1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\mathrm{True}</script></html>"
      ],
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "N(M.scalar_field(x))==v"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\mathrm{True}</script></html>"
      ],
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "N(M.scalar_field(v))==a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\mathrm{True}</script></html>"
      ],
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "N(M.scalar_field(a))==s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\mathrm{True}</script></html>"
      ],
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "N(M.scalar_field(s))==f"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define auxillary fields"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\begin{array}{llcl} & M & \\longrightarrow & \\mathbb{R} \\\\ & \\left(t, x, v, a, s\\right) & \\longmapsto & \\frac{\\partial\\,L}{\\partial s} \\end{array}</script></html>"
      ],
      "text/plain": [
       "M --> R\n",
       "(t, x, v, a, s) |--> d(L)/ds"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "r=Ds(L); r.display()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\begin{array}{llcl} & M & \\longrightarrow & \\mathbb{R} \\\\ & \\left(t, x, v, a, s\\right) & \\longmapsto & -v \\frac{\\partial^2\\,L}{\\partial x\\partial s} - a \\frac{\\partial^2\\,L}{\\partial v\\partial s} - s \\frac{\\partial^2\\,L}{\\partial a\\partial s} - f \\frac{\\partial^2\\,L}{\\partial s^2} - \\frac{\\partial^2\\,L}{\\partial t\\partial s} + \\frac{\\partial\\,L}{\\partial a} \\end{array}</script></html>"
      ],
      "text/plain": [
       "M --> R\n",
       "(t, x, v, a, s) |--> -v*d^2(L)/dxds - a*d^2(L)/dvds - s*d^2(L)/dads - f*d^2(L)/ds^2 - d^2(L)/dtds + d(L)/da"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "q=Da(L)-N(r); q.display()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\begin{array}{llcl} & M & \\longrightarrow & \\mathbb{R} \\\\ & \\left(t, x, v, a, s\\right) & \\longmapsto & v^{2} \\frac{\\partial^3\\,L}{\\partial x^2\\partial s} + a^{2} \\frac{\\partial^3\\,L}{\\partial v^2\\partial s} + 2 \\, a s \\frac{\\partial^3\\,L}{\\partial v\\partial a\\partial s} + 2 \\, a f \\frac{\\partial^3\\,L}{\\partial v\\partial s^2} + s^{2} \\frac{\\partial^3\\,L}{\\partial a^2\\partial s} + 2 \\, f s \\frac{\\partial^3\\,L}{\\partial a\\partial s^2} + f^{2} \\frac{\\partial^3\\,L}{\\partial s^3} + {\\left(2 \\, a \\frac{\\partial^3\\,L}{\\partial x\\partial v\\partial s} + 2 \\, s \\frac{\\partial^3\\,L}{\\partial x\\partial a\\partial s} + 2 \\, f \\frac{\\partial^3\\,L}{\\partial x\\partial s^2} + 2 \\, \\frac{\\partial^3\\,L}{\\partial t\\partial x\\partial s} - \\frac{\\partial^2\\,L}{\\partial x\\partial a}\\right)} v + 2 \\, a \\frac{\\partial^3\\,L}{\\partial t\\partial v\\partial s} + 2 \\, s \\frac{\\partial^3\\,L}{\\partial t\\partial a\\partial s} + 2 \\, f \\frac{\\partial^3\\,L}{\\partial t\\partial s^2} + a \\frac{\\partial^2\\,L}{\\partial x\\partial s} - a \\frac{\\partial^2\\,L}{\\partial v\\partial a} + s \\frac{\\partial^2\\,L}{\\partial v\\partial s} - s \\frac{\\partial^2\\,L}{\\partial a^2} + \\frac{\\partial^3\\,L}{\\partial t^2\\partial s} - \\frac{\\partial^2\\,L}{\\partial t\\partial a} + \\frac{\\partial\\,L}{\\partial v} \\end{array}</script></html>"
      ],
      "text/plain": [
       "M --> R\n",
       "(t, x, v, a, s) |--> v^2*d^3(L)/dx^2ds + a^2*d^3(L)/dv^2ds + 2*a*s*d^3(L)/dvdads + 2*a*f*d^3(L)/dvds^2 + s^2*d^3(L)/da^2ds + 2*f*s*d^3(L)/dads^2 + f^2*d^3(L)/ds^3 + (2*a*d^3(L)/dxdvds + 2*s*d^3(L)/dxdads + 2*f*d^3(L)/dxds^2 + 2*d^3(L)/dtdxds - d^2(L)/dxda)*v + 2*a*d^3(L)/dtdvds + 2*s*d^3(L)/dtdads + 2*f*d^3(L)/dtds^2 + a*d^2(L)/dxds - a*d^2(L)/dvda + s*d^2(L)/dvds - s*d^2(L)/da^2 + d^3(L)/dt^2ds - d^2(L)/dtda + d(L)/dv"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p=Dv(L)-N(q); p.display()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First Equation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\begin{array}{llcl} & M & \\longrightarrow & \\mathbb{R} \\\\ & \\left(t, x, v, a, s\\right) & \\longmapsto & v^{3} \\frac{\\partial^4\\,L}{\\partial x^3\\partial s} + a^{3} \\frac{\\partial^4\\,L}{\\partial v^3\\partial s} + 3 \\, a^{2} s \\frac{\\partial^4\\,L}{\\partial v^2\\partial a\\partial s} + 3 \\, a^{2} f \\frac{\\partial^4\\,L}{\\partial v^2\\partial s^2} + 3 \\, a s^{2} \\frac{\\partial^4\\,L}{\\partial v\\partial a^2\\partial s} + 6 \\, a f s \\frac{\\partial^4\\,L}{\\partial v\\partial a\\partial s^2} + 3 \\, a f^{2} \\frac{\\partial^4\\,L}{\\partial v\\partial s^3} + s^{3} \\frac{\\partial^4\\,L}{\\partial a^3\\partial s} + 3 \\, f s^{2} \\frac{\\partial^4\\,L}{\\partial a^2\\partial s^2} + 3 \\, f^{2} s \\frac{\\partial^4\\,L}{\\partial a\\partial s^3} + f^{3} \\frac{\\partial^4\\,L}{\\partial s^4} + {\\left(3 \\, a \\frac{\\partial^4\\,L}{\\partial x^2\\partial v\\partial s} + 3 \\, s \\frac{\\partial^4\\,L}{\\partial x^2\\partial a\\partial s} + 3 \\, f \\frac{\\partial^4\\,L}{\\partial x^2\\partial s^2} + 3 \\, \\frac{\\partial^4\\,L}{\\partial t\\partial x^2\\partial s} - \\frac{\\partial^3\\,L}{\\partial x^2\\partial a}\\right)} v^{2} + 3 \\, a^{2} \\frac{\\partial^4\\,L}{\\partial t\\partial v^2\\partial s} + 6 \\, a s \\frac{\\partial^4\\,L}{\\partial t\\partial v\\partial a\\partial s} + 6 \\, a f \\frac{\\partial^4\\,L}{\\partial t\\partial v\\partial s^2} + 3 \\, s^{2} \\frac{\\partial^4\\,L}{\\partial t\\partial a^2\\partial s} + 6 \\, f s \\frac{\\partial^4\\,L}{\\partial t\\partial a\\partial s^2} + 3 \\, f^{2} \\frac{\\partial^4\\,L}{\\partial t\\partial s^3} + 3 \\, a^{2} \\frac{\\partial^3\\,L}{\\partial x\\partial v\\partial s} + 3 \\, a s \\frac{\\partial^3\\,L}{\\partial x\\partial a\\partial s} + 3 \\, a f \\frac{\\partial^3\\,L}{\\partial x\\partial s^2} - a^{2} \\frac{\\partial^3\\,L}{\\partial v^2\\partial a} + 3 \\, a s \\frac{\\partial^3\\,L}{\\partial v^2\\partial s} - 2 \\, a s \\frac{\\partial^3\\,L}{\\partial v\\partial a^2} + 3 \\, f s \\frac{\\partial^3\\,L}{\\partial v\\partial s^2} - s^{2} \\frac{\\partial^3\\,L}{\\partial a^3} + f s \\frac{\\partial^3\\,L}{\\partial a^2\\partial s} + 2 \\, f^{2} \\frac{\\partial^3\\,L}{\\partial a\\partial s^2} + {\\left(3 \\, a^{2} \\frac{\\partial^4\\,L}{\\partial x\\partial v^2\\partial s} + 6 \\, a s \\frac{\\partial^4\\,L}{\\partial x\\partial v\\partial a\\partial s} + 6 \\, a f \\frac{\\partial^4\\,L}{\\partial x\\partial v\\partial s^2} + 3 \\, s^{2} \\frac{\\partial^4\\,L}{\\partial x\\partial a^2\\partial s} + 6 \\, f s \\frac{\\partial^4\\,L}{\\partial x\\partial a\\partial s^2} + 3 \\, f^{2} \\frac{\\partial^4\\,L}{\\partial x\\partial s^3} + 6 \\, a \\frac{\\partial^4\\,L}{\\partial t\\partial x\\partial v\\partial s} + 6 \\, s \\frac{\\partial^4\\,L}{\\partial t\\partial x\\partial a\\partial s} + 6 \\, f \\frac{\\partial^4\\,L}{\\partial t\\partial x\\partial s^2} + 3 \\, a \\frac{\\partial^3\\,L}{\\partial x^2\\partial s} - 2 \\, a \\frac{\\partial^3\\,L}{\\partial x\\partial v\\partial a} + 3 \\, s \\frac{\\partial^3\\,L}{\\partial x\\partial v\\partial s} - 2 \\, s \\frac{\\partial^3\\,L}{\\partial x\\partial a^2} + f \\frac{\\partial^3\\,L}{\\partial x\\partial a\\partial s} + 3 \\, \\frac{\\partial^4\\,L}{\\partial t^2\\partial x\\partial s} - 2 \\, \\frac{\\partial^3\\,L}{\\partial t\\partial x\\partial a} + \\frac{\\partial^2\\,L}{\\partial x\\partial v}\\right)} v + 3 \\, a \\frac{\\partial^4\\,L}{\\partial t^2\\partial v\\partial s} + 3 \\, s \\frac{\\partial^4\\,L}{\\partial t^2\\partial a\\partial s} + 3 \\, f \\frac{\\partial^4\\,L}{\\partial t^2\\partial s^2} + 3 \\, a \\frac{\\partial^3\\,L}{\\partial t\\partial x\\partial s} - 2 \\, a \\frac{\\partial^3\\,L}{\\partial t\\partial v\\partial a} + 3 \\, s \\frac{\\partial^3\\,L}{\\partial t\\partial v\\partial s} - 2 \\, s \\frac{\\partial^3\\,L}{\\partial t\\partial a^2} + f \\frac{\\partial^3\\,L}{\\partial t\\partial a\\partial s} - a \\frac{\\partial^2\\,L}{\\partial x\\partial a} + s \\frac{\\partial^2\\,L}{\\partial x\\partial s} + a \\frac{\\partial^2\\,L}{\\partial v^2} + {\\left(a f + 3 \\, s^{2}\\right)} \\frac{\\partial^3\\,L}{\\partial v\\partial a\\partial s} + 2 \\, f \\frac{\\partial^2\\,L}{\\partial v\\partial s} - f \\frac{\\partial^2\\,L}{\\partial a^2} + \\frac{\\partial^4\\,L}{\\partial t^3\\partial s} - \\frac{\\partial^3\\,L}{\\partial t^2\\partial a} + \\frac{\\partial^2\\,L}{\\partial t\\partial v} - \\frac{\\partial\\,L}{\\partial x} \\end{array}</script></html>"
      ],
      "text/plain": [
       "M --> R\n",
       "(t, x, v, a, s) |--> v^3*d^4(L)/dx^3ds + a^3*d^4(L)/dv^3ds + 3*a^2*s*d^4(L)/dv^2dads + 3*a^2*f*d^4(L)/dv^2ds^2 + 3*a*s^2*d^4(L)/dvda^2ds + 6*a*f*s*d^4(L)/dvdads^2 + 3*a*f^2*d^4(L)/dvds^3 + s^3*d^4(L)/da^3ds + 3*f*s^2*d^4(L)/da^2ds^2 + 3*f^2*s*d^4(L)/dads^3 + f^3*d^4(L)/ds^4 + (3*a*d^4(L)/dx^2dvds + 3*s*d^4(L)/dx^2dads + 3*f*d^4(L)/dx^2ds^2 + 3*d^4(L)/dtdx^2ds - d^3(L)/dx^2da)*v^2 + 3*a^2*d^4(L)/dtdv^2ds + 6*a*s*d^4(L)/dtdvdads + 6*a*f*d^4(L)/dtdvds^2 + 3*s^2*d^4(L)/dtda^2ds + 6*f*s*d^4(L)/dtdads^2 + 3*f^2*d^4(L)/dtds^3 + 3*a^2*d^3(L)/dxdvds + 3*a*s*d^3(L)/dxdads + 3*a*f*d^3(L)/dxds^2 - a^2*d^3(L)/dv^2da + 3*a*s*d^3(L)/dv^2ds - 2*a*s*d^3(L)/dvda^2 + 3*f*s*d^3(L)/dvds^2 - s^2*d^3(L)/da^3 + f*s*d^3(L)/da^2ds + 2*f^2*d^3(L)/dads^2 + (3*a^2*d^4(L)/dxdv^2ds + 6*a*s*d^4(L)/dxdvdads + 6*a*f*d^4(L)/dxdvds^2 + 3*s^2*d^4(L)/dxda^2ds + 6*f*s*d^4(L)/dxdads^2 + 3*f^2*d^4(L)/dxds^3 + 6*a*d^4(L)/dtdxdvds + 6*s*d^4(L)/dtdxdads + 6*f*d^4(L)/dtdxds^2 + 3*a*d^3(L)/dx^2ds - 2*a*d^3(L)/dxdvda + 3*s*d^3(L)/dxdvds - 2*s*d^3(L)/dxda^2 + f*d^3(L)/dxdads + 3*d^4(L)/dt^2dxds - 2*d^3(L)/dtdxda + d^2(L)/dxdv)*v + 3*a*d^4(L)/dt^2dvds + 3*s*d^4(L)/dt^2dads + 3*f*d^4(L)/dt^2ds^2 + 3*a*d^3(L)/dtdxds - 2*a*d^3(L)/dtdvda + 3*s*d^3(L)/dtdvds - 2*s*d^3(L)/dtda^2 + f*d^3(L)/dtdads - a*d^2(L)/dxda + s*d^2(L)/dxds + a*d^2(L)/dv^2 + (a*f + 3*s^2)*d^3(L)/dvdads + 2*f*d^2(L)/dvds - f*d^2(L)/da^2 + d^4(L)/dt^3ds - d^3(L)/dt^2da + d^2(L)/dtdv - d(L)/dx"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(N(p)-Dx(L)).display()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Second Equation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\begin{array}{llcl} & M & \\longrightarrow & \\mathbb{R} \\\\ & \\left(t, x, v, a, s\\right) & \\longmapsto & v^{4} \\frac{\\partial^4\\,L}{\\partial x^3\\partial s} + {\\left(3 \\, a \\frac{\\partial^4\\,L}{\\partial x^2\\partial v\\partial s} + 3 \\, s \\frac{\\partial^4\\,L}{\\partial x^2\\partial a\\partial s} + 3 \\, f \\frac{\\partial^4\\,L}{\\partial x^2\\partial s^2} + 3 \\, \\frac{\\partial^4\\,L}{\\partial t\\partial x^2\\partial s} - \\frac{\\partial^3\\,L}{\\partial x^2\\partial a}\\right)} v^{3} + {\\left(3 \\, a^{2} \\frac{\\partial^4\\,L}{\\partial x\\partial v^2\\partial s} + 6 \\, a s \\frac{\\partial^4\\,L}{\\partial x\\partial v\\partial a\\partial s} + 6 \\, a f \\frac{\\partial^4\\,L}{\\partial x\\partial v\\partial s^2} + 3 \\, s^{2} \\frac{\\partial^4\\,L}{\\partial x\\partial a^2\\partial s} + 6 \\, f s \\frac{\\partial^4\\,L}{\\partial x\\partial a\\partial s^2} + 3 \\, f^{2} \\frac{\\partial^4\\,L}{\\partial x\\partial s^3} + 6 \\, a \\frac{\\partial^4\\,L}{\\partial t\\partial x\\partial v\\partial s} + 6 \\, s \\frac{\\partial^4\\,L}{\\partial t\\partial x\\partial a\\partial s} + 6 \\, f \\frac{\\partial^4\\,L}{\\partial t\\partial x\\partial s^2} + 3 \\, a \\frac{\\partial^3\\,L}{\\partial x^2\\partial s} - 2 \\, a \\frac{\\partial^3\\,L}{\\partial x\\partial v\\partial a} + 3 \\, s \\frac{\\partial^3\\,L}{\\partial x\\partial v\\partial s} - 2 \\, s \\frac{\\partial^3\\,L}{\\partial x\\partial a^2} + f \\frac{\\partial^3\\,L}{\\partial x\\partial a\\partial s} + 3 \\, \\frac{\\partial^4\\,L}{\\partial t^2\\partial x\\partial s} - 2 \\, \\frac{\\partial^3\\,L}{\\partial t\\partial x\\partial a} + \\frac{\\partial^2\\,L}{\\partial x\\partial v}\\right)} v^{2} + {\\left(a^{3} \\frac{\\partial^4\\,L}{\\partial v^3\\partial s} + 3 \\, a^{2} s \\frac{\\partial^4\\,L}{\\partial v^2\\partial a\\partial s} + 3 \\, a^{2} f \\frac{\\partial^4\\,L}{\\partial v^2\\partial s^2} + 3 \\, a s^{2} \\frac{\\partial^4\\,L}{\\partial v\\partial a^2\\partial s} + 6 \\, a f s \\frac{\\partial^4\\,L}{\\partial v\\partial a\\partial s^2} + 3 \\, a f^{2} \\frac{\\partial^4\\,L}{\\partial v\\partial s^3} + s^{3} \\frac{\\partial^4\\,L}{\\partial a^3\\partial s} + 3 \\, f s^{2} \\frac{\\partial^4\\,L}{\\partial a^2\\partial s^2} + 3 \\, f^{2} s \\frac{\\partial^4\\,L}{\\partial a\\partial s^3} + f^{3} \\frac{\\partial^4\\,L}{\\partial s^4} + 3 \\, a^{2} \\frac{\\partial^4\\,L}{\\partial t\\partial v^2\\partial s} + 6 \\, a s \\frac{\\partial^4\\,L}{\\partial t\\partial v\\partial a\\partial s} + 6 \\, a f \\frac{\\partial^4\\,L}{\\partial t\\partial v\\partial s^2} + 3 \\, s^{2} \\frac{\\partial^4\\,L}{\\partial t\\partial a^2\\partial s} + 6 \\, f s \\frac{\\partial^4\\,L}{\\partial t\\partial a\\partial s^2} + 3 \\, f^{2} \\frac{\\partial^4\\,L}{\\partial t\\partial s^3} + 3 \\, a^{2} \\frac{\\partial^3\\,L}{\\partial x\\partial v\\partial s} + 3 \\, a s \\frac{\\partial^3\\,L}{\\partial x\\partial a\\partial s} + 3 \\, a f \\frac{\\partial^3\\,L}{\\partial x\\partial s^2} - a^{2} \\frac{\\partial^3\\,L}{\\partial v^2\\partial a} + 3 \\, a s \\frac{\\partial^3\\,L}{\\partial v^2\\partial s} - 2 \\, a s \\frac{\\partial^3\\,L}{\\partial v\\partial a^2} + 3 \\, f s \\frac{\\partial^3\\,L}{\\partial v\\partial s^2} - s^{2} \\frac{\\partial^3\\,L}{\\partial a^3} + f s \\frac{\\partial^3\\,L}{\\partial a^2\\partial s} + 2 \\, f^{2} \\frac{\\partial^3\\,L}{\\partial a\\partial s^2} + 3 \\, a \\frac{\\partial^4\\,L}{\\partial t^2\\partial v\\partial s} + 3 \\, s \\frac{\\partial^4\\,L}{\\partial t^2\\partial a\\partial s} + 3 \\, f \\frac{\\partial^4\\,L}{\\partial t^2\\partial s^2} + 3 \\, a \\frac{\\partial^3\\,L}{\\partial t\\partial x\\partial s} - 2 \\, a \\frac{\\partial^3\\,L}{\\partial t\\partial v\\partial a} + 3 \\, s \\frac{\\partial^3\\,L}{\\partial t\\partial v\\partial s} - 2 \\, s \\frac{\\partial^3\\,L}{\\partial t\\partial a^2} + f \\frac{\\partial^3\\,L}{\\partial t\\partial a\\partial s} - a \\frac{\\partial^2\\,L}{\\partial x\\partial a} + s \\frac{\\partial^2\\,L}{\\partial x\\partial s} + a \\frac{\\partial^2\\,L}{\\partial v^2} + {\\left(a f + 3 \\, s^{2}\\right)} \\frac{\\partial^3\\,L}{\\partial v\\partial a\\partial s} + 2 \\, f \\frac{\\partial^2\\,L}{\\partial v\\partial s} - f \\frac{\\partial^2\\,L}{\\partial a^2} + \\frac{\\partial^4\\,L}{\\partial t^3\\partial s} - \\frac{\\partial^3\\,L}{\\partial t^2\\partial a} + \\frac{\\partial^2\\,L}{\\partial t\\partial v} - \\frac{\\partial\\,L}{\\partial x}\\right)} v \\end{array}</script></html>"
      ],
      "text/plain": [
       "M --> R\n",
       "(t, x, v, a, s) |--> v^4*d^4(L)/dx^3ds + (3*a*d^4(L)/dx^2dvds + 3*s*d^4(L)/dx^2dads + 3*f*d^4(L)/dx^2ds^2 + 3*d^4(L)/dtdx^2ds - d^3(L)/dx^2da)*v^3 + (3*a^2*d^4(L)/dxdv^2ds + 6*a*s*d^4(L)/dxdvdads + 6*a*f*d^4(L)/dxdvds^2 + 3*s^2*d^4(L)/dxda^2ds + 6*f*s*d^4(L)/dxdads^2 + 3*f^2*d^4(L)/dxds^3 + 6*a*d^4(L)/dtdxdvds + 6*s*d^4(L)/dtdxdads + 6*f*d^4(L)/dtdxds^2 + 3*a*d^3(L)/dx^2ds - 2*a*d^3(L)/dxdvda + 3*s*d^3(L)/dxdvds - 2*s*d^3(L)/dxda^2 + f*d^3(L)/dxdads + 3*d^4(L)/dt^2dxds - 2*d^3(L)/dtdxda + d^2(L)/dxdv)*v^2 + (a^3*d^4(L)/dv^3ds + 3*a^2*s*d^4(L)/dv^2dads + 3*a^2*f*d^4(L)/dv^2ds^2 + 3*a*s^2*d^4(L)/dvda^2ds + 6*a*f*s*d^4(L)/dvdads^2 + 3*a*f^2*d^4(L)/dvds^3 + s^3*d^4(L)/da^3ds + 3*f*s^2*d^4(L)/da^2ds^2 + 3*f^2*s*d^4(L)/dads^3 + f^3*d^4(L)/ds^4 + 3*a^2*d^4(L)/dtdv^2ds + 6*a*s*d^4(L)/dtdvdads + 6*a*f*d^4(L)/dtdvds^2 + 3*s^2*d^4(L)/dtda^2ds + 6*f*s*d^4(L)/dtdads^2 + 3*f^2*d^4(L)/dtds^3 + 3*a^2*d^3(L)/dxdvds + 3*a*s*d^3(L)/dxdads + 3*a*f*d^3(L)/dxds^2 - a^2*d^3(L)/dv^2da + 3*a*s*d^3(L)/dv^2ds - 2*a*s*d^3(L)/dvda^2 + 3*f*s*d^3(L)/dvds^2 - s^2*d^3(L)/da^3 + f*s*d^3(L)/da^2ds + 2*f^2*d^3(L)/dads^2 + 3*a*d^4(L)/dt^2dvds + 3*s*d^4(L)/dt^2dads + 3*f*d^4(L)/dt^2ds^2 + 3*a*d^3(L)/dtdxds - 2*a*d^3(L)/dtdvda + 3*s*d^3(L)/dtdvds - 2*s*d^3(L)/dtda^2 + f*d^3(L)/dtdads - a*d^2(L)/dxda + s*d^2(L)/dxds + a*d^2(L)/dv^2 + (a*f + 3*s^2)*d^3(L)/dvdads + 2*f*d^2(L)/dvds - f*d^2(L)/da^2 + d^4(L)/dt^3ds - d^3(L)/dt^2da + d^2(L)/dtdv - d(L)/dx)*v"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(Dt(L)-(N(L)-(p*a+q*s+r*f)-(v*N(p)+a*N(q)+s*N(r)))).display()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Multiplying by $v$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\mathrm{True}</script></html>"
      ],
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "v*(N(p)-Dx(L)) == (Dt(L)-(N(L)-(p*a+q*s+r*f)-(v*N(p)+a*N(q)+s*N(r))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Checking the calculations from the paper ##"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "L = M.scalar_field(function('L')(*list(coord)))\n",
    "p = M.scalar_field(function('p')(*list(coord)))\n",
    "q = M.scalar_field(function('q')(*list(coord)))\n",
    "r = M.scalar_field(function('r')(*list(coord)))\n",
    "t=M.scalar_field(t)\n",
    "x=M.scalar_field(x)\n",
    "v=M.scalar_field(v)\n",
    "a=M.scalar_field(a)\n",
    "s=M.scalar_field(s)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Action differential Form"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left( -v p\\left(t, x, v, a, s\\right) - a q\\left(t, x, v, a, s\\right) - s r\\left(t, x, v, a, s\\right) + L\\left(t, x, v, a, s\\right) \\right) \\mathrm{d} t + p\\left(t, x, v, a, s\\right) \\mathrm{d} x + q\\left(t, x, v, a, s\\right) \\mathrm{d} v + r\\left(t, x, v, a, s\\right) \\mathrm{d} a</script></html>"
      ],
      "text/plain": [
       "(-v*p(t, x, v, a, s) - a*q(t, x, v, a, s) - s*r(t, x, v, a, s) + L(t, x, v, a, s)) dt + p(t, x, v, a, s) dx + q(t, x, v, a, s) dv + r(t, x, v, a, s) da"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "alpha = L*dt + p*(dx-v*dt) + q*(dv-a*dt) + r*(da-s*dt)\n",
    "alpha.display()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\mathrm{True}</script></html>"
      ],
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "alpha == L*dt+p*dx+q*dv+r*da-(p*v+q*a+r*s)*dt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left( v \\frac{\\partial\\,p}{\\partial x} + a \\frac{\\partial\\,q}{\\partial x} + s \\frac{\\partial\\,r}{\\partial x} - \\frac{\\partial\\,L}{\\partial x} + \\frac{\\partial\\,p}{\\partial t} \\right) \\mathrm{d} t\\wedge \\mathrm{d} x + \\left( v \\frac{\\partial\\,p}{\\partial v} + a \\frac{\\partial\\,q}{\\partial v} + s \\frac{\\partial\\,r}{\\partial v} + p\\left(t, x, v, a, s\\right) - \\frac{\\partial\\,L}{\\partial v} + \\frac{\\partial\\,q}{\\partial t} \\right) \\mathrm{d} t\\wedge \\mathrm{d} v + \\left( v \\frac{\\partial\\,p}{\\partial a} + a \\frac{\\partial\\,q}{\\partial a} + s \\frac{\\partial\\,r}{\\partial a} + q\\left(t, x, v, a, s\\right) - \\frac{\\partial\\,L}{\\partial a} + \\frac{\\partial\\,r}{\\partial t} \\right) \\mathrm{d} t\\wedge \\mathrm{d} a + \\left( v \\frac{\\partial\\,p}{\\partial s} + a \\frac{\\partial\\,q}{\\partial s} + s \\frac{\\partial\\,r}{\\partial s} + r\\left(t, x, v, a, s\\right) - \\frac{\\partial\\,L}{\\partial s} \\right) \\mathrm{d} t\\wedge \\mathrm{d} s + \\left( -\\frac{\\partial\\,p}{\\partial v} + \\frac{\\partial\\,q}{\\partial x} \\right) \\mathrm{d} x\\wedge \\mathrm{d} v + \\left( -\\frac{\\partial\\,p}{\\partial a} + \\frac{\\partial\\,r}{\\partial x} \\right) \\mathrm{d} x\\wedge \\mathrm{d} a -\\frac{\\partial\\,p}{\\partial s} \\mathrm{d} x\\wedge \\mathrm{d} s + \\left( -\\frac{\\partial\\,q}{\\partial a} + \\frac{\\partial\\,r}{\\partial v} \\right) \\mathrm{d} v\\wedge \\mathrm{d} a -\\frac{\\partial\\,q}{\\partial s} \\mathrm{d} v\\wedge \\mathrm{d} s -\\frac{\\partial\\,r}{\\partial s} \\mathrm{d} a\\wedge \\mathrm{d} s</script></html>"
      ],
      "text/plain": [
       "(v*d(p)/dx + a*d(q)/dx + s*d(r)/dx - d(L)/dx + d(p)/dt) dt/\\dx + (v*d(p)/dv + a*d(q)/dv + s*d(r)/dv + p(t, x, v, a, s) - d(L)/dv + d(q)/dt) dt/\\dv + (v*d(p)/da + a*d(q)/da + s*d(r)/da + q(t, x, v, a, s) - d(L)/da + d(r)/dt) dt/\\da + (v*d(p)/ds + a*d(q)/ds + s*d(r)/ds + r(t, x, v, a, s) - d(L)/ds) dt/\\ds + (-d(p)/dv + d(q)/dx) dx/\\dv + (-d(p)/da + d(r)/dx) dx/\\da - d(p)/ds dx/\\ds + (-d(q)/da + d(r)/dv) dv/\\da - d(q)/ds dv/\\ds - d(r)/ds da/\\ds"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "d(alpha).display()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\mathrm{True}</script></html>"
      ],
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "d(alpha) == d(L).wedge(dt) + d(p).wedge(dx) + d(q).wedge(dv) + d(r).wedge(da) - d(p*v + q*a + r*s).wedge(dt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\mathrm{True}</script></html>"
      ],
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ev(alpha)(N)==L"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}p\\left(t, x, v, a, s\\right) \\mathrm{d} t\\wedge \\mathrm{d} x + q\\left(t, x, v, a, s\\right) \\mathrm{d} t\\wedge \\mathrm{d} v + r\\left(t, x, v, a, s\\right) \\mathrm{d} t\\wedge \\mathrm{d} a</script></html>"
      ],
      "text/plain": [
       "p(t, x, v, a, s) dt/\\dx + q(t, x, v, a, s) dt/\\dv + r(t, x, v, a, s) dt/\\da"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Omega = -(p*dx+q*dv+r*da).wedge(d(t)); Omega.display()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\mathrm{True}</script></html>"
      ],
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "alpha == L*dt + ev(N)(Omega)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Equation of Motion ($E=0$)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left( -v^{2} \\frac{\\partial\\,p}{\\partial x} - a^{2} \\frac{\\partial\\,q}{\\partial v} - a s \\frac{\\partial\\,q}{\\partial a} - a f \\frac{\\partial\\,q}{\\partial s} - a s \\frac{\\partial\\,r}{\\partial v} - s^{2} \\frac{\\partial\\,r}{\\partial a} - f s \\frac{\\partial\\,r}{\\partial s} - {\\left(a \\frac{\\partial\\,p}{\\partial v} + s \\frac{\\partial\\,p}{\\partial a} + f \\frac{\\partial\\,p}{\\partial s} + a \\frac{\\partial\\,q}{\\partial x} + s \\frac{\\partial\\,r}{\\partial x} - \\frac{\\partial\\,L}{\\partial x} + \\frac{\\partial\\,p}{\\partial t}\\right)} v - a p\\left(t, x, v, a, s\\right) - s q\\left(t, x, v, a, s\\right) - f r\\left(t, x, v, a, s\\right) + a \\frac{\\partial\\,L}{\\partial v} + s \\frac{\\partial\\,L}{\\partial a} + f \\frac{\\partial\\,L}{\\partial s} - a \\frac{\\partial\\,q}{\\partial t} - s \\frac{\\partial\\,r}{\\partial t} \\right) \\mathrm{d} t + \\left( v \\frac{\\partial\\,p}{\\partial x} + a \\frac{\\partial\\,p}{\\partial v} + s \\frac{\\partial\\,p}{\\partial a} + f \\frac{\\partial\\,p}{\\partial s} - \\frac{\\partial\\,L}{\\partial x} + \\frac{\\partial\\,p}{\\partial t} \\right) \\mathrm{d} x + \\left( v \\frac{\\partial\\,q}{\\partial x} + a \\frac{\\partial\\,q}{\\partial v} + s \\frac{\\partial\\,q}{\\partial a} + f \\frac{\\partial\\,q}{\\partial s} + p\\left(t, x, v, a, s\\right) - \\frac{\\partial\\,L}{\\partial v} + \\frac{\\partial\\,q}{\\partial t} \\right) \\mathrm{d} v + \\left( v \\frac{\\partial\\,r}{\\partial x} + a \\frac{\\partial\\,r}{\\partial v} + s \\frac{\\partial\\,r}{\\partial a} + f \\frac{\\partial\\,r}{\\partial s} + q\\left(t, x, v, a, s\\right) - \\frac{\\partial\\,L}{\\partial a} + \\frac{\\partial\\,r}{\\partial t} \\right) \\mathrm{d} a + \\left( r\\left(t, x, v, a, s\\right) - \\frac{\\partial\\,L}{\\partial s} \\right) \\mathrm{d} s</script></html>"
      ],
      "text/plain": [
       "(-v^2*d(p)/dx - a^2*d(q)/dv - a*s*d(q)/da - a*f*d(q)/ds - a*s*d(r)/dv - s^2*d(r)/da - f*s*d(r)/ds - (a*d(p)/dv + s*d(p)/da + f*d(p)/ds + a*d(q)/dx + s*d(r)/dx - d(L)/dx + d(p)/dt)*v - a*p(t, x, v, a, s) - s*q(t, x, v, a, s) - f*r(t, x, v, a, s) + a*d(L)/dv + s*d(L)/da + f*d(L)/ds - a*d(q)/dt - s*d(r)/dt) dt + (v*d(p)/dx + a*d(p)/dv + s*d(p)/da + f*d(p)/ds - d(L)/dx + d(p)/dt) dx + (v*d(q)/dx + a*d(q)/dv + s*d(q)/da + f*d(q)/ds + p(t, x, v, a, s) - d(L)/dv + d(q)/dt) dv + (v*d(r)/dx + a*d(r)/dv + s*d(r)/da + f*d(r)/ds + q(t, x, v, a, s) - d(L)/da + d(r)/dt) da + (r(t, x, v, a, s) - d(L)/ds) ds"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "E = ev(N)(d(alpha))\n",
    "E.display()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Rewriting it in various ways."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\mathrm{True}</script></html>"
      ],
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "E == N(L)*dt - N(t) * d(L) + N(p)*dx - N(x)*d(p) + N(q)*dv - N(v)*d(q) + N(r)*da - N(a)*d(r) - N(p*v+q*a+r*s)*dt + N(t)*d(p*v+q*a+r*s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\mathrm{True}</script></html>"
      ],
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "E == N(L)*dt - d(L) + N(p)*dx - v*d(p) + N(q)*dv - a*d(q) + N(r)*da - s*d(r) - N(p*v+q*a+r*s)*dt + d(p*v+q*a+r*s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\mathrm{True}</script></html>"
      ],
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "E == N(L-p*v-q*a-r*s)*dt - d(L - p*v - q*a- r*s) - v*d(p) - a*d(q) - s*d(r) + N(p)*dx + N(q)*dv + N(r)*da"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\mathrm{True}</script></html>"
      ],
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "d(p*v) == v*d(p) + p*d(v)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\mathrm{True}</script></html>"
      ],
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "d(q*a) == a*d(q) + q*d(a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\mathrm{True}</script></html>"
      ],
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "d(r*s) == s*d(r) + r*d(s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\mathrm{True}</script></html>"
      ],
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "E == N(L - p*v - q*a - r*s)*dt - d(L) + p*dv + q*da + r*ds + N(p)*dx + N(q)*dv + N(r)*da"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\mathrm{True}</script></html>"
      ],
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "N(p*v) == p*N(v) + v*N(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\mathrm{True}</script></html>"
      ],
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "E == N(L)*dt - (p*a+q*s+r*f)*dt -(v*N(p) + a*N(q) + s*N(r) )*dt - d(L) + p*dv + q*da + r*ds + N(p)*dx + N(q)*dv + N(r)*da"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\mathrm{True}</script></html>"
      ],
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "E == ( N(L) - (p*a + q*s + r*f) - (v*N(p) + a*N(q) + s*N(r)))*dt + N(p)*dx + (N(q)+p)*dv + (N(r)+q)*da + r*ds - d(L)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\mathrm{True}</script></html>"
      ],
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "d(L) == Dt(L)*dt + Dx(L)*dx + Dv(L)*dv + Da(L)*da + Ds(L)*ds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\begin{array}{llcl} & M & \\longrightarrow & \\mathbb{R} \\\\ & \\left(t, x, v, a, s\\right) & \\longmapsto & \\frac{\\partial\\,L}{\\partial s} \\end{array}</script></html>"
      ],
      "text/plain": [
       "M --> R\n",
       "(t, x, v, a, s) |--> d(L)/ds"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "r=Ds(L); r.display()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\begin{array}{llcl} & M & \\longrightarrow & \\mathbb{R} \\\\ & \\left(t, x, v, a, s\\right) & \\longmapsto & -v \\frac{\\partial^2\\,L}{\\partial x\\partial s} - a \\frac{\\partial^2\\,L}{\\partial v\\partial s} - s \\frac{\\partial^2\\,L}{\\partial a\\partial s} - f \\frac{\\partial^2\\,L}{\\partial s^2} - \\frac{\\partial^2\\,L}{\\partial t\\partial s} + \\frac{\\partial\\,L}{\\partial a} \\end{array}</script></html>"
      ],
      "text/plain": [
       "M --> R\n",
       "(t, x, v, a, s) |--> -v*d^2(L)/dxds - a*d^2(L)/dvds - s*d^2(L)/dads - f*d^2(L)/ds^2 - d^2(L)/dtds + d(L)/da"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "q=Da(L)-N(r); q.display()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\begin{array}{llcl} & M & \\longrightarrow & \\mathbb{R} \\\\ & \\left(t, x, v, a, s\\right) & \\longmapsto & v^{2} \\frac{\\partial^3\\,L}{\\partial x^2\\partial s} + a^{2} \\frac{\\partial^3\\,L}{\\partial v^2\\partial s} + 2 \\, a s \\frac{\\partial^3\\,L}{\\partial v\\partial a\\partial s} + 2 \\, a f \\frac{\\partial^3\\,L}{\\partial v\\partial s^2} + s^{2} \\frac{\\partial^3\\,L}{\\partial a^2\\partial s} + 2 \\, f s \\frac{\\partial^3\\,L}{\\partial a\\partial s^2} + f^{2} \\frac{\\partial^3\\,L}{\\partial s^3} + {\\left(2 \\, a \\frac{\\partial^3\\,L}{\\partial x\\partial v\\partial s} + 2 \\, s \\frac{\\partial^3\\,L}{\\partial x\\partial a\\partial s} + 2 \\, f \\frac{\\partial^3\\,L}{\\partial x\\partial s^2} + 2 \\, \\frac{\\partial^3\\,L}{\\partial t\\partial x\\partial s} - \\frac{\\partial^2\\,L}{\\partial x\\partial a}\\right)} v + 2 \\, a \\frac{\\partial^3\\,L}{\\partial t\\partial v\\partial s} + 2 \\, s \\frac{\\partial^3\\,L}{\\partial t\\partial a\\partial s} + 2 \\, f \\frac{\\partial^3\\,L}{\\partial t\\partial s^2} + a \\frac{\\partial^2\\,L}{\\partial x\\partial s} - a \\frac{\\partial^2\\,L}{\\partial v\\partial a} + s \\frac{\\partial^2\\,L}{\\partial v\\partial s} - s \\frac{\\partial^2\\,L}{\\partial a^2} + \\frac{\\partial^3\\,L}{\\partial t^2\\partial s} - \\frac{\\partial^2\\,L}{\\partial t\\partial a} + \\frac{\\partial\\,L}{\\partial v} \\end{array}</script></html>"
      ],
      "text/plain": [
       "M --> R\n",
       "(t, x, v, a, s) |--> v^2*d^3(L)/dx^2ds + a^2*d^3(L)/dv^2ds + 2*a*s*d^3(L)/dvdads + 2*a*f*d^3(L)/dvds^2 + s^2*d^3(L)/da^2ds + 2*f*s*d^3(L)/dads^2 + f^2*d^3(L)/ds^3 + (2*a*d^3(L)/dxdvds + 2*s*d^3(L)/dxdads + 2*f*d^3(L)/dxds^2 + 2*d^3(L)/dtdxds - d^2(L)/dxda)*v + 2*a*d^3(L)/dtdvds + 2*s*d^3(L)/dtdads + 2*f*d^3(L)/dtds^2 + a*d^2(L)/dxds - a*d^2(L)/dvda + s*d^2(L)/dvds - s*d^2(L)/da^2 + d^3(L)/dt^2ds - d^2(L)/dtda + d(L)/dv"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p=Dv(L)-N(q); p.display()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## For example: The Schiff and Poirier Lagrangian ##"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\begin{array}{llcl} & M & \\longrightarrow & \\mathbb{R} \\\\ & \\left(t, x, v, a, s\\right) & \\longmapsto & \\frac{4 \\, m^{2} v^{6} - 8 \\, m v^{4} V\\left(x\\right) + 5 \\, a^{2} {\\hbar}^{2} - 2 \\, {\\hbar}^{2} s v}{8 \\, m v^{4}} \\end{array}</script></html>"
      ],
      "text/plain": [
       "M --> R\n",
       "(t, x, v, a, s) |--> 1/8*(4*m^2*v^6 - 8*m*v^4*V(x) + 5*a^2*hbar^2 - 2*hbar^2*s*v)/(m*v^4)"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "hbar = var('hbar',latex_name='\\hbar')\n",
    "m = var('m')\n",
    "V = M.scalar_field(function('V')(var('x')))\n",
    "L = 1/2*m*v^2  - V - hbar^2/4/m*(s/v^3-5/2*a^2/v^4); L.display()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\begin{array}{llcl} & M & \\longrightarrow & \\mathbb{R} \\\\ & \\left(t, x, v, a, s\\right) & \\longmapsto & -\\frac{{\\hbar}^{2}}{4 \\, m v^{3}} \\end{array}</script></html>"
      ],
      "text/plain": [
       "M --> R\n",
       "(t, x, v, a, s) |--> -1/4*hbar^2/(m*v^3)"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "r=Ds(L); r.display()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\begin{array}{llcl} & M & \\longrightarrow & \\mathbb{R} \\\\ & \\left(t, x, v, a, s\\right) & \\longmapsto & \\frac{a {\\hbar}^{2}}{2 \\, m v^{4}} \\end{array}</script></html>"
      ],
      "text/plain": [
       "M --> R\n",
       "(t, x, v, a, s) |--> 1/2*a*hbar^2/(m*v^4)"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "q=Da(L)-N(r); q.display()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\begin{array}{llcl} & M & \\longrightarrow & \\mathbb{R} \\\\ & \\left(t, x, v, a, s\\right) & \\longmapsto & \\frac{4 \\, m^{2} v^{6} - 2 \\, a^{2} {\\hbar}^{2} + {\\hbar}^{2} s v}{4 \\, m v^{5}} \\end{array}</script></html>"
      ],
      "text/plain": [
       "M --> R\n",
       "(t, x, v, a, s) |--> 1/4*(4*m^2*v^6 - 2*a^2*hbar^2 + hbar^2*s*v)/(m*v^5)"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p=Dv(L)-N(q); p.display()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First Equation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}a m + \\frac{5 \\, a^{3} {\\hbar}^{2}}{2 \\, m v^{6}} - \\frac{2 \\, a {\\hbar}^{2} s}{m v^{5}} + \\frac{f {\\hbar}^{2}}{4 \\, m v^{4}} + D[0]\\left(V\\right)\\left(x\\right)</script></html>"
      ],
      "text/plain": [
       "a*m + 5/2*a^3*hbar^2/(m*v^6) - 2*a*hbar^2*s/(m*v^5) + 1/4*f*hbar^2/(m*v^4) + D[0](V)(x)"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(N(p)-Dx(L)).expr().expand()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Second Equation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}a m v + v D[0]\\left(V\\right)\\left(x\\right) + \\frac{5 \\, a^{3} {\\hbar}^{2}}{2 \\, m v^{5}} - \\frac{2 \\, a {\\hbar}^{2} s}{m v^{4}} + \\frac{f {\\hbar}^{2}}{4 \\, m v^{3}}</script></html>"
      ],
      "text/plain": [
       "a*m*v + v*D[0](V)(x) + 5/2*a^3*hbar^2/(m*v^5) - 2*a*hbar^2*s/(m*v^4) + 1/4*f*hbar^2/(m*v^3)"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(Dt(L)-(N(L)-(p*a+q*s+r*f)-(v*N(p)+a*N(q)+s*N(r)))).expr().expand()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "SageMath 7.2",
   "language": "",
   "name": "sagemath"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}