wolfram Mathematica package to calculate geodetics from Riemannian metric

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  ]]
}
]
*)

Geodetics.m

BeginPackage["Geodetics`"]

ChristoffelSymbols::usage = "Calculates Christoffel symbols of the second kind.
Provide the Riemannian metric as a matrix vaued function and the number of parameters.";
ChristoffelSymbols[G_, n_] := 
 Function[{Args(*=Table[Subscript[x, q], {q, n}]*)},
  Module[{ginv=Inverse[G@Args]},
   Table[ 
    1/2 Sum[
            ginv[[m, k]] (D[G[Args][[k, i]], Args[[j]]] + 
            D[G[Args][[k, j]], Args[[i]]] - 
            D[G[Args][[i, j]], Args[[k]]])
        , {k, n}]
    , {i, n}, {j, n}, {m, n}]
  ]
 ]

CurvatureTensor[G_, n_] := Function[{Args},
  Module[{Ch = ChristoffelSymbols[G, n][Args]},
   Table[D[Ch[[i, k, l]], Args[[j]]] - D[Ch[[i, j, l]], Args[[k]]] + 
     Sum[Ch[[j, s, l]] Ch[[i, k, s]] - 
       Ch[[k, s, l]] Ch[[i, j, s]], {s, n}], {i, n}, {j, n}, {k, 
     n}, {l, n}]]
  ]
  
RicciTensor[G_, n_] := 
 Function[{Args}, 
  Module[{CT = CurvatureTensor[G, n][Args]}, 
   Table[Sum[CT[[i, k, j, k]], {k, n}], {i, n}, {j, n}]
  ]
 ]
 
RicciScalar[G_, n_] := Function[{Args}, 
  Module[{RT = RicciTensor[G, n][Args], ginv=Inverse[G@Args]},
    Sum[ginv[[i,j]]RT[[i,j]],{i,n},{j,n}]]
  ]

GeodeticEquation[G_, vars_, p_] := Module[
  {n = Length[vars], 
   X = Table[Subscript[x, i], {i, Length[vars]}]},
  Table[vars[[i]]''[p] +
  Sum[
    ChristoffelSymbols[G, Length[vars]][#[p]& /@ vars][[j, k, i]] vars[[j]]'[p] vars[[k]]'[p]
    , {k, n}, {j, n}] == 0
  , {i, n}]
 ]
 
EndPackage[]