Wolfram Mathematica package for deriving Lagrangian and Hamiltonian (classical mechanics)

demo.nb

(* Content-type: application/vnd.wolfram.mathematica *)

(*** Wolfram Notebook File ***)
(* http://www.wolfram.com/nb *)

(* CreatedBy='Mathematica 9.0' *)

(*CacheID: 234*)
(* Internal cache information:
NotebookFileLineBreakTest
NotebookFileLineBreakTest
NotebookDataPosition[       157,          7]
NotebookDataLength[     10994,        332]
NotebookOptionsPosition[     10018,        296]
NotebookOutlinePosition[     10530,        316]
CellTagsIndexPosition[     10487,        313]
WindowFrame->Normal*)

(* Beginning of Notebook Content *)
Notebook[{
Cell[BoxData[
 RowBox[{"Import", "[", 
  RowBox[{
   RowBox[{"NotebookDirectory", "[", "]"}], "<>", "\"\<Hamiltonian.m\>\""}], 
  "]"}]], "Input",
 CellChangeTimes->{{3.748247581702917*^9, 3.7482476035718985`*^9}}],

Cell[CellGroupData[{

Cell["Examples", "Section",
 CellChangeTimes->{{3.694028635610221*^9, 3.6940286381088347`*^9}}],

Cell[CellGroupData[{

Cell["Pendulum", "Subsection",
 CellChangeTimes->{{3.694025777955056*^9, 3.694025780543663*^9}, {
  3.694028616484562*^9, 3.694028617731206*^9}, {3.694028682178193*^9, 
  3.694028683121877*^9}}],

Cell["with control", "Text",
 CellChangeTimes->{{3.6940286841239758`*^9, 3.694028688955155*^9}}],

Cell[BoxData[{
 RowBox[{
  RowBox[{"x1", "[", "\[Phi]_", "]"}], ":=", 
  RowBox[{"{", 
   RowBox[{
    RowBox[{"l", " ", 
     RowBox[{"Sin", "[", "\[Phi]", "]"}]}], ",", " ", 
    RowBox[{"l", " ", 
     RowBox[{"Cos", "[", "\[Phi]", "]"}]}]}], "}"}]}], "\[IndentingNewLine]", 
 RowBox[{
  RowBox[{"T", "[", 
   RowBox[{"\[Phi]_", ",", " ", "d\[Phi]_"}], "]"}], ":=", 
  RowBox[{
   FractionBox["1", "2"], "m", 
   RowBox[{"(", 
    RowBox[{
     SuperscriptBox[
      RowBox[{"(", 
       RowBox[{
        RowBox[{
         RowBox[{
          RowBox[{"x1", "'"}], "[", "\[Phi]", "]"}], "\[LeftDoubleBracket]", 
         "1", "\[RightDoubleBracket]"}], "d\[Phi]"}], ")"}], "2"], "+", 
     SuperscriptBox[
      RowBox[{"(", 
       RowBox[{
        RowBox[{
         RowBox[{
          RowBox[{"x1", "'"}], "[", "\[Phi]", "]"}], "\[LeftDoubleBracket]", 
         "2", "\[RightDoubleBracket]"}], "d\[Phi]"}], ")"}], "2"]}], 
    ")"}]}]}], "\[IndentingNewLine]", 
 RowBox[{
  RowBox[{"V", "[", 
   RowBox[{"\[Phi]_", ",", "d\[Phi]_"}], "]"}], ":=", 
  RowBox[{
   RowBox[{"m", " ", "g", " ", 
    RowBox[{
     RowBox[{"x1", "[", "\[Phi]", "]"}], "\[LeftDoubleBracket]", "2", 
     "\[RightDoubleBracket]"}]}], "-", 
   RowBox[{"\[Phi]", " ", 
    RowBox[{"u", "[", "t", "]"}]}]}]}], "\[IndentingNewLine]", 
 RowBox[{
  RowBox[{"L", "[", 
   RowBox[{"q_", ",", "dq_", ",", "t_"}], "]"}], ":=", 
  RowBox[{
   RowBox[{"T", "[", 
    RowBox[{"q", ",", " ", "dq"}], "]"}], "-", 
   RowBox[{"V", "[", 
    RowBox[{"q", ",", " ", "dq"}], "]"}]}]}]}], "Input",
 CellChangeTimes->{{3.694025129997781*^9, 3.6940251532591066`*^9}, {
   3.6940253703549137`*^9, 3.694025370522032*^9}, 3.6940256610311203`*^9, {
   3.694025772390231*^9, 3.694025798225137*^9}, {3.6940277113560925`*^9, 
   3.6940277151687593`*^9}, {3.6940319940450907`*^9, 
   3.6940319941631527`*^9}, {3.6940877619419765`*^9, 3.6940877988728275`*^9}}],

Cell[BoxData[{
 RowBox[{
  RowBox[{"L", "[", 
   RowBox[{"\[Phi]", ",", "d\[Phi]", ",", "t"}], "]"}], "//", 
  "Simplify"}], "\[IndentingNewLine]", 
 RowBox[{
  RowBox[{"MakeLagrangeEquation", "[", 
   RowBox[{"L", ",", "1"}], "]"}], "//", "Simplify"}], "\[IndentingNewLine]", 
 RowBox[{"GeneralizedMomentumInverze", "[", 
  RowBox[{"L", ",", "1"}], "]"}], "\[IndentingNewLine]", 
 RowBox[{
  RowBox[{
   RowBox[{"MakeHamiltonian", "[", 
    RowBox[{"L", ",", "1"}], "]"}], "[", 
   RowBox[{"q", ",", "p", ",", "t"}], "]"}], "//", 
  "Simplify"}], "\[IndentingNewLine]", 
 RowBox[{
  RowBox[{"MakeHamiltonEquation", "[", 
   RowBox[{
    RowBox[{"MakeHamiltonian", "[", 
     RowBox[{"L", ",", "1"}], "]"}], ",", "1"}], "]"}], "//", 
  "Simplify"}]}], "Input",
 CellChangeTimes->{{3.694027759220313*^9, 3.6940277596073446`*^9}, {
  3.6940280940254135`*^9, 3.694028100148511*^9}, {3.6940281984350615`*^9, 
  3.69402820022972*^9}, {3.6940737980361958`*^9, 3.6940737985993853`*^9}}]
}, Open  ]],

Cell[CellGroupData[{

Cell["Pendulum on rails", "Subsection",
 CellChangeTimes->{{3.694028624441499*^9, 3.6940286337421274`*^9}, {
  3.694029987905763*^9, 3.6940299883666463`*^9}}],

Cell["with control", "Text",
 CellChangeTimes->{{3.6940286841239758`*^9, 3.694028688955155*^9}}],

Cell[BoxData[{
 RowBox[{
  RowBox[{"x1", "[", 
   RowBox[{"s_", ",", " ", "\[Phi]_"}], "]"}], ":=", 
  RowBox[{"{", 
   RowBox[{"s", ",", " ", "0"}], "}"}]}], "\[IndentingNewLine]", 
 RowBox[{
  RowBox[{"x2", "[", 
   RowBox[{"s_", ",", " ", "\[Phi]_"}], "]"}], ":=", 
  RowBox[{"{", " ", 
   RowBox[{
    RowBox[{"s", "+", 
     RowBox[{"l", " ", 
      RowBox[{"Sin", "[", "\[Phi]", "]"}]}]}], ",", " ", 
    RowBox[{"l", " ", 
     RowBox[{"Cos", "[", "\[Phi]", "]"}]}]}], "}"}]}]}], "Input"],

Cell[BoxData[
 RowBox[{
  RowBox[{"(", 
   RowBox[{
    RowBox[{
     RowBox[{
      RowBox[{"Table", "[", 
       RowBox[{
        RowBox[{
         FractionBox["1", "2"], 
         RowBox[{"(", 
          RowBox[{
           RowBox[{
            RowBox[{"#", ".", "#"}], "&"}], "@", 
           RowBox[{"(", 
            RowBox[{"D", "[", 
             RowBox[{
              RowBox[{"y", "[", 
               RowBox[{
                RowBox[{"s", "[", "t", "]"}], ",", 
                RowBox[{"\[Phi]", "[", "t", "]"}]}], "]"}], ",", "t"}], "]"}],
             ")"}]}], ")"}]}], ",", 
        RowBox[{"{", 
         RowBox[{"y", ",", 
          RowBox[{"{", 
           RowBox[{"x1", ",", "x2"}], "}"}]}], "}"}]}], "]"}], ".", 
      RowBox[{"{", 
       RowBox[{"M", ",", "m"}], "}"}]}], "//", "Simplify"}], "//", "Expand"}],
    ")"}], "/.", 
  RowBox[{"{", 
   RowBox[{
    RowBox[{
     RowBox[{"s", "[", "t", "]"}], "\[Rule]", "s"}], ",", 
    RowBox[{
     RowBox[{"\[Phi]", "[", "t", "]"}], "\[Rule]", "\[Phi]"}], ",", " ", 
    RowBox[{
     RowBox[{
      SuperscriptBox["s", "\[Prime]",
       MultilineFunction->None], "[", "t", "]"}], "\[Rule]", "ds"}], ",", " ", 
    RowBox[{
     RowBox[{
      SuperscriptBox["\[Phi]", "\[Prime]",
       MultilineFunction->None], "[", "t", "]"}], "\[Rule]", "d\[Phi]"}]}], 
   "}"}]}]], "Input",
 CellChangeTimes->{{3.6940322289095736`*^9, 3.6940322793996468`*^9}, {
  3.694032314781993*^9, 3.694032316135379*^9}, {3.6940337341917825`*^9, 
  3.694033802541934*^9}}],

Cell[BoxData[{
 RowBox[{
  RowBox[{"T", "[", 
   RowBox[{
   "s_", ",", " ", "\[Phi]_", ",", " ", "ds_", ",", " ", "d\[Phi]_", ",", " ",
     "t_"}], "]"}], ":=", 
  RowBox[{
   FractionBox[
    RowBox[{
     SuperscriptBox["ds", "2"], " ", "m"}], "2"], "+", 
   RowBox[{
    FractionBox["1", "2"], " ", 
    SuperscriptBox["d\[Phi]", "2"], " ", 
    SuperscriptBox["l", "2"], " ", "m"}], "+", 
   FractionBox[
    RowBox[{
     SuperscriptBox["ds", "2"], " ", "M"}], "2"], "+", 
   RowBox[{"ds", " ", "d\[Phi]", " ", "l", " ", "m", " ", 
    RowBox[{"Cos", "[", "\[Phi]", "]"}]}]}]}], "\[IndentingNewLine]", 
 RowBox[{
  RowBox[{
   RowBox[{"V", "[", 
    RowBox[{"q1_", ",", " ", "q2_", ",", " ", "t_"}], "]"}], ":=", 
   RowBox[{
    RowBox[{"-", "g"}], " ", 
    RowBox[{
     RowBox[{"{", 
      RowBox[{"M", ",", "m"}], "}"}], ".", 
     RowBox[{"Table", "[", 
      RowBox[{
       RowBox[{
        RowBox[{"y", "[", 
         RowBox[{"q1", ",", "q2"}], "]"}], "\[LeftDoubleBracket]", "2", 
        "\[RightDoubleBracket]"}], ",", 
       RowBox[{"{", 
        RowBox[{"y", ",", 
         RowBox[{"{", 
          RowBox[{"x1", ",", "x2"}], "}"}]}], "}"}]}], "]"}]}]}]}], 
  RowBox[{"(*", 
   RowBox[{
    RowBox[{"-", "q1"}], " ", 
    RowBox[{"u", "[", "t", "]"}]}], "*)"}]}], "\[IndentingNewLine]", 
 RowBox[{
  RowBox[{"L", "[", 
   RowBox[{
   "q1_", ",", " ", "q2_", ",", " ", "dq1_", ",", " ", "dq2_", ",", " ", 
    "t_"}], "]"}], ":=", 
  RowBox[{
   RowBox[{"T", "[", 
    RowBox[{
    "q1", ",", " ", "q2", ",", " ", "dq1", ",", " ", "dq2", ",", " ", "t"}], 
    "]"}], "-", 
   RowBox[{"V", "[", 
    RowBox[{"q1", ",", " ", "q2", ",", " ", "t"}], "]"}]}]}]}], "Input",
 CellChangeTimes->{{3.694028712564621*^9, 3.6940287128394337`*^9}, {
   3.6940287563897595`*^9, 3.6940287567088594`*^9}, {3.694031867813525*^9, 
   3.6940318679419327`*^9}, {3.69403230417896*^9, 3.694032306721613*^9}, {
   3.694033691187278*^9, 3.6940337208295918`*^9}, {3.6940337871644015`*^9, 
   3.694033817716612*^9}, {3.6940340444165525`*^9, 3.6940340490490417`*^9}, {
   3.694087830331758*^9, 3.694087833594532*^9}, 3.694088051803752*^9}],

Cell[BoxData[
 RowBox[{"ClearAll", "[", 
  RowBox[{"m", ",", "M", ",", "g", ",", "l"}], "]"}]], "Input",
 CellChangeTimes->{{3.6940880636240325`*^9, 3.6940880687932277`*^9}}],

Cell[BoxData[{
 RowBox[{
  RowBox[{"L", "[", 
   RowBox[{"s", ",", "\[Phi]", ",", 
    OverscriptBox["s", "."], ",", 
    OverscriptBox["\[Phi]", "."], ",", "t"}], "]"}], "//", 
  "Simplify"}], "\[IndentingNewLine]", 
 RowBox[{
  RowBox[{"MakeLagrangeEquation", "[", 
   RowBox[{"L", ",", "2"}], "]"}], "//", "Simplify"}], "\[IndentingNewLine]", 
 RowBox[{"GeneralizedMomentumInverze", "[", 
  RowBox[{"L", ",", "2"}], "]"}], "\[IndentingNewLine]", 
 RowBox[{
  RowBox[{
   RowBox[{"MakeHamiltonian", "[", 
    RowBox[{"L", ",", "2"}], "]"}], "[", 
   RowBox[{"s", ",", "\[Phi]", ",", "ps", ",", "p\[Phi]", ",", "t"}], "]"}], "//",
   "Simplify"}], "\[IndentingNewLine]", 
 RowBox[{
  RowBox[{"MakeHamiltonEquation", "[", 
   RowBox[{
    RowBox[{"MakeHamiltonian", "[", 
     RowBox[{"L", ",", "2"}], "]"}], ",", "2", ",", " ", "q", ",", "p"}], 
   "]"}], "//", "Simplify"}]}], "Input",
 CellChangeTimes->{{3.694028773220317*^9, 3.6940288068994613`*^9}, {
  3.6940341097825117`*^9, 3.694034137866268*^9}, {3.6940341966446342`*^9, 
  3.6940342259321556`*^9}, {3.69407023086001*^9, 3.6940702842028427`*^9}, {
  3.6940711637869034`*^9, 3.6940711643181934`*^9}, {3.6940713368086443`*^9, 
  3.6940713368617744`*^9}, {3.6940737442957954`*^9, 3.694073745318322*^9}}]
}, Open  ]]
}, Open  ]]
},
WindowSize->{936, 816},
WindowMargins->{{Automatic, -105}, {Automatic, -6}},
PrintingCopies->1,
PrintingPageRange->{32000, 32000},
PrintingOptions->{"Magnification"->1.,
"PaperOrientation"->"Portrait",
"PaperSize"->{595.3199999999999, 841.92}},
FrontEndVersion->"11.0 for Microsoft Windows (64-bit) (July 28, 2016)",
StyleDefinitions->"Default.nb"
]
(* End of Notebook Content *)

(* Internal cache information *)
(*CellTagsOutline
CellTagsIndex->{}
*)
(*CellTagsIndex
CellTagsIndex->{}
*)
(*NotebookFileOutline
Notebook[{
Cell[557, 20, 214, 5, 30, "Input"],
Cell[CellGroupData[{
Cell[796, 29, 95, 1, 63, "Section"],
Cell[CellGroupData[{
Cell[916, 34, 194, 3, 43, "Subsection"],
Cell[1113, 39, 96, 1, 30, "Text"],
Cell[1212, 42, 1907, 53, 106, "Input"],
Cell[3122, 97, 979, 24, 107, "Input"]
}, Open  ]],
Cell[CellGroupData[{
Cell[4138, 126, 158, 2, 43, "Subsection"],
Cell[4299, 130, 96, 1, 30, "Text"],
Cell[4398, 133, 495, 15, 50, "Input"],
Cell[4896, 150, 1519, 46, 68, "Input"],
Cell[6418, 198, 2132, 59, 90, "Input"],
Cell[8553, 259, 174, 3, 30, "Input"],
Cell[8730, 264, 1260, 28, 109, "Input"]
}, Open  ]]
}, Open  ]]
}
]
*)

Hamiltonian.m

(* ::Package:: *)

BeginPackage["Hamiltonian`"]

arg::usage = "Generates symbolic argument list of given length, optionally given symbol name.";
arg[n_, symbol_: q] := 
 Table[Symbol[ToString[symbol] <> ToString[i]], {i, 1, n}]

KineticEnergy::usage = "Generates kinetic energy formula from parameterized coordinates.";
KineticEnergy[X_, M_, n_, qsym_: q, tsym_: t] := 
 Module[{Args = arg[n, qsym]},
  (1/2 (#.# & /@ D[X @@ (#[tsym] & /@ Args), tsym]).M) /. 
   Flatten@Table[{a[tsym] -> a, 
      Derivative[1][a][tsym] -> 
       ToExpression["d" <> ToString[a]]}, {a, Args}]
  ]

GravitationalPotential::usage = "Generates gravitational potential formula from parameterized coordinates. The gravitational constant is g."
GravitationalPotential[X_, M_, n_, qsym_: q] := g M.(Last /@ (X @@ arg[n, qsym]))

MakeLagrangeEquation::usage = "Generates Langrange equations of the second kind from a given Lagrangian.
Explicit time dependence should be included, even if zero!";
MakeLagrangeEquation[L_, n_, qsym_: q, tsym_: t] := 
 Module[{Args = Join[arg[2 n, qsym], {tsym}, 1], Subs},
  Subs = Join @@ Table[{Args[[j]] -> Args[[j]][tsym], Args[[n + j]] -> Derivative[1][Args[[j]]][tsym]}, {j, 1, n, 1}];
  Table[
   D[((Derivative @@ UnitVector[2 n + 1, n + j])[L]) @@ Args /. Subs, tsym] -
     (((Derivative @@ UnitVector[2 n + 1, j])[L]) @@ Args /. Subs) == 0,
    {j, 1, n}]
  ]

GeneralizedMomentumInverze::usage = "Calculates dq in terms of generalized momentums from a given Lagrangian.";
GeneralizedMomentumInverze[L_, n_, qsym_: q, psym_: p, tsym_: t] := 
 Module[{Args = Join[arg[2 n, qsym], {tsym}, 1]},
  Simplify[
   #[[2]] & /@ Solve[
      Table[
       arg[n, psym][[j]] == ((Derivative @@ UnitVector[2 n + 1, n + j])[L]) @@ Args,
       {j, 1, n, 1}],
       Args[[n + 1 ;; 2 n]]
       ][[1]]
   ]
  ]
MakeHamiltonian::usage = "Calculates Hamiltoninan from Lagrangian (Legendre transform).";
MakeHamiltonian[L_, n_, qsym_: q, psym_: p, tsym_: t] := Module[
  {
     LArgs = Join[arg[2 n, qsym], {tsym}, 1],
     HArgs = Join[arg[n, qsym], GeneralizedMomentumInverze[L, n, qsym, psym, tsym], {tsym}, 1]
   },
  Function @@ {Join[arg[n, qsym], arg[n, psym], {tsym}, 1], 
    GeneralizedMomentumInverze[L, n, qsym, psym, tsym].arg[n, psym] - 
     L @@ HArgs}
  ]

MakeHamiltonEquation::usage = "Generates Hamiltonian equations from a given Hamiltonian.
Explicit time dependence should be included, even if zero!";
MakeHamiltonEquation[H_, n_, qsym_: q, psym_: p, tsym_: t] := Module[
  {args = Join[arg[n, qsym], arg[n, psym], {tsym}, 1], rules},
  rules = (Rule @@ # & /@ (Partition[Riffle[args[[;; -2]], #[tsym] & /@ args[[;; -2]]], 2]));
  Flatten[Table[
    {args[[n + k]]'[t] == -D[H @@ args, args[[k]]] /. rules, 
     args[[k]]'[t] == (D[H @@ args, args[[n + k]]] /. rules)},
    {k, 1, n, 1}], 1]
  ]

EndPackage[]