1999 Mathematica Developer Conference

Applications of Mathematica 
October 23, 1999

Simulation of Forces : 
Solution of the Inverse Dynamics 
in Classical Mechanics

Robert Kragler 
FH Ravensburg-Weingarten /University of Applied Sciences
kragler@fh-weingarten.de
A.N. Prokopenya and N.I. Chopchits
Brest Polytechnic Institute, Belarus (box@aprokop.belpak.brest.by)

Initializations

Working Directory "\\Math\\WnMath30\\Vortrag" and Settings

Uses StyleSheet Classroom.nb and working directory "\\Math\\WnMath30\\Vortrag"

actDirectory = "\\Math\\WnMath30"<>"\\Vortrag";

SetOptions["stdout", PageWidth -> 80]; \(
\)Off[SetDirectory::cdir]; \(
\)actDir  =  SetDirectory[actDirectory];

If[ StringPosition[actDir, ":", 1] == {{2, 2}},  
      actDir  =  StringDrop[actDir, 2] ];

subDir  =  actDir<>"\\Vortrag1";  \(
\){ actDir, subDir, $TopDirectory} // ColumnForm

($Path = $Path ⋃ {actDir, subDir}) // ColumnForm

PackageInfo

\(PackageInfo::usage\  = 

        "PackageInfo["context`"] provides basic information for all

names (symbols) in a given context. In lexicographic order all

all names used will be listed with their supplementary information

in order to give a survey of all procedures available in a package."\);

Attributes[ Info ]       =  HoldFirst;
Info[ symb_Symbol ]  := 
        ToExpression["? " <> ToString[symb] ];
Info[ name_String ]  :=  ToExpression["? " <> name ];

PackageInfo[ contxt_String] :=Scan[ Info[#]&,  Names[contxt <> "*"]];

openNotebookButton & filesInDirectory

Off[ General::spell]; \(
\)Off[General::spell1]; \(
\)Notebooks[];

PathToFileName[path_]:=StringDrop[path, Last[Last[StringPosition[path, $PathnameSeparator]]]]
        
PathToOpenBox[path_]:=ButtonBox[PathToFileName[path], ButtonFunction :> NotebookOpen[path]]
        
openNotebookButton[pattern_, depth_: 3, directory_: Automatic]:=
		\(Module[{temp, dir}, 
				\ \ \ \ \ \ \ \ \ \ \ \ 
        dir\  = If[\ directory\  === \ Automatic, \ 

            					\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Nest[\ ParentDirectory, \ Directory[], depth - 1], 

            				\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ directory]\ ; 
		\ 			\ \ \ \ \ \ \ \ 
        NotebookWrite[\ SelectedNotebook[], 
						
          Cell[\ BoxData[\ 
					\ \ \ \ \ \ \ \ \ \ 
              GridBox[
						\ \ \ \ \ \ \ \ \ \ \ \ 
                Partition[\ 
								\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
                  Map[\ PathToOpenBox, 
                    temp = FileNames[pattern, dir, depth]\ ], 

                  						\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 	\ \ \ \ \ Min[2, Length[temp]]

                  						\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ], 

                					\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ RowSpacings\  -> 1, 

                					\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 	ColumnSpacings -> 2

                					\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ]
		\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ], 
					\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
            "Text", Active -> True], 

          			\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ After]

        	\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ]\)

filesInDirectory[path_] := 
	\(Module[\ {}, \ 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
        fileList\  = \ FileNames["*.nb", path, 1]\ ; 

        \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
        extractFile\  := 
          StringDrop[\ #, 1\  + \ \(StringPosition[#, path]\)[\([1, 2]\)]\ ]&
            \ \ \ ; \ 
\ \ \ \ \ \ \ 
        Table[\ extractFile\ @\ fileList[\([i]\)], 
          \ {i, 1, Length[fileList]}\ ]\ \ \ ]\)

openFileBrowser

Calling openFileBrowser[ ]corresponds to Open in the  File menu.

openFileBrowser[] :=
 \(If[TrueQ[Head[$FrontEnd] === FrontEndObject], 

          \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
          \(Function[link, \ LinkWrite[link, \ FrontEnd`FileBrowse[False]]; 

              \ \ \ \ \ \ \ \ \ \ \ \ \ \ LinkRead[link]\ ]\)[
            First[$FrontEnd]]\ , \ $Failed
\ \ \ \ ]\);

Load Mathematica Packages

Needs["Graphics`Colors`"]; \(
\)Needs["Graphics`ArgColors`"]; \(
\)Needs["Graphics`Arrow`"]; \(
\)PackageInfo["Graphics`ArgColors`"];

Indexed Variables  (Subscript[a, 0]) , (OverVector[V]), (Subscript[(OverVector[V]), n])

Needs["Utilities`Notation`"] ;

In order to use symbols with indices for names of functions or as arguments of functions it is required to declare these quantities using the procedure Symbolize (siehe TMJ 7(1) : 35) from the Package Utilities`Notation  .

Off[General::spell1]; \(
\)Off[Symbolize::boxSymbolExists]

Symbolize[\(_\__\)];         
Symbolize[\(_\&⇀\)];           
Symbolize[\(_\&⇀\__\)];

On[Symbolize::boxSymbolExists]

The correct notation for indexed vectors  e.g.  (Subscript[(Overscript[V, ⇀]), n])  must be of the following form :
			(Subscript[(OverVector[V]), n])   ⟺   V⎵Overscript⎵RightVector⎵Subscript⎵n
Examples of the correct internal representation of combined symbols according to the definitions above are :

(Subscript[a, 0])  // FullForm\(
\)(Overscript[V, ⇀])    // FullForm\(
\)(Subscript[(Overscript[V, ⇀]), n])  // FullForm

Cell-Group with Icon ⊳ resp. ▽

This cell when unformatted with Ctrl+⇧+E  shows the mechanism of the underlying Triangle icon which appears when this cell is the head of a cell group

Cell["Turning triangle icon                          ", "Text",
  ShowCellBracket->True,
  ShowGroupOpenCloseIcon->True,
  CellMargins->{{Inherited, 44.5625}, {Inherited, Inherited}},
  PageBreakBelow->False,
  TextAlignment->Left,
  TextJustification->0]

Set notebook options

Off[General::Spell1];  Off[General::Spell]; \(
\)$DefaultFont  =  {"Helvetica-Bold", 10}; \(
\)th  =  AbsoluteThickness[2]; \(
\)SetOptions[Plot, PlotStyle -> {th}]; \(
\)SetOptions[ListPlot, PlotStyle -> {Red, AbsolutePointSize[2]}];

Information on important System Variables

\(StringForm["Mathematica Version : `` Release : ``", \ \ $Version, 

      	\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ $ReleaseNumber\ ]\)

$ContextPath  // ColumnForm

\(StringForm["Notebooks in `` : 
            `` ", 

      	\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Directory[], 
      FileNames["*.nb"] // ColumnForm]\)

datum = Date[ ];
tag   = StringForm["on  ``.``.``   ",datum[[1]],datum[[2]],datum[[3]] ];
If[ datum[[5]] < 10,
    zeit = StringForm["at ``:0`` h", datum[[4]], datum[[5]] ] ,
    zeit = StringForm["at ``:`` h" , datum[[4]], datum[[5]] ] ];

Print["Initializations finished ",tag,zeit ]

Motivation

In general, studying the equations of motion of any physical system it is necessary to solve differential equations. However, very seldomly the corresponding solutions can be found in analytical form. Thus, usually only systems which are simple enough will be analysed in a physics course at universities. Situation changed considerably, however, with the availability of computer algebra systems. Using Mathematica in many cases we can easily find either analytical or numerical solutions of the equations of motion and visualize them even in the case of an intricate physical system. The analysis of such systems therefore promotes better understanding of the underlying physical principles and may help to develop physical intuition of the students.

Generation of Figures

Figure.1.

\(HB = "Helvetica-Bold"; 

    Module[{t, y, y1, δ}, \ \ 
	\ \ \ \ \ \ \ \ \ \ \ \ 
      y[t_] = 5  Exp[\(-0.1\) t\^2] + 1; \ y1[t_] = \(y'\)[t]; 
      δ = ArcTan[y1[2]]; 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ 
      Plot[y[t], {t, 0, 6}, 
	\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
        PlotRange -> {{\(-2.5\), 7}, {\(-2\), 8}}, Ticks -> None, 

        	\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
        AxesOrigin -> {\(-1\), 0}, \ AspectRatio -> Automatic, 

        	\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
        PlotStyle -> {{Red, th}}, 

        		\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
        Prolog -> {
            	{Dashing[{0.05, 0.03}], 

              	\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 				
              Line[{{\(-1\), y[2.5]}, {4.5, y[2.5]}}]}, 

            				\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {Green, Disk[{2.5, y[2.5]}, 0.3]}\ }, 

        	\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
        Epilog -> {\ Arrow[{\(-1\), \(-0.5\)}, {\(-1\), 8}], \ 
            Arrow[{\(-0.5\), 0}, {7, 0}], 

            				\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 			\ \ \ \ \ \ \ \ \ Arrow[{2.5, y[2.5]}, {4, y[2.5] - 1.5/y1[2.5]}], 	\ \ \ \ 

            				\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Arrow[{2.5, y[2.5]}, {2.5, y[2.5] - 3}], 

            		\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 	{Dashing[{0.05, 0.03}], Line[{{0, 0}, {0, y[0]}}]}, 

            		\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 	{Dashing[{0.05, 0.03}], Line[{{2.5, 0}, {2.5, y[2.5]}}]\ }, 

            			\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {Dashing[{0.05, 0.03}], Line[{{\(-1\), y[0]}, {0, y[0]}}]\ }, 

            			\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Line[{{0.5, y[2.5] - 2\ y1[2.5]}, {4.5, y[2.5] + 2\ y1[2.5]}}], 

            				\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Circle[{2.5, y[2.5]}, 0.8, {0, π + δ}], 

            		\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 	Text[FontForm["x", {HB, 12}], {6.6, \(-0.4\)}], 

            			\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Text[FontForm["y", {HB, 12}], {\(-1.4\), 7.5}], 

            			\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Text[FontForm["O", {HB, 14}], {\(-1.4\), \(-0.4\)}], 

            			\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
            Text[FontForm[\*""(Subscript[y, 0])"", {HB, 12}], {\(-1.4\), 
                y[0]}], 

            	\ 		\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Text[FontForm["y", {HB, 12}], {\(-1.4\), y[2.5]}], 

            		\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
            Text[FontForm[\*""(Subscript[x, 0])"", {HB, 12}], {0, \(-0.4\)}], 
            
	\ 		\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Text[FontForm["x", {HB, 12}], {2.5, \(-0.4\)}], 

            		\ 	\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
            Text[FontForm[\*""(Overscript[N,  ->]) "", {HB, 12}], {
                4.4, 4.5}], 

            	\ 		\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
            Text[FontForm[\*""m(Overscript[g,  ->]) "", {HB, 12}], {
                3.2, 1}], 

            	\ 		\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Text[FontForm["α", {HB, 14}], {3, 4.8}], 

            \ 			\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Text[FontForm["y= f(x)", {HB, 12}], {6.2, y[6] + 0.7}], 

            \ 			\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Text[FontForm["Fig. 1.", {HB, 12}], {4, \(-1.7\)}]

            \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ }\ ]; 

      	\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ]; 
    \)

Figure 2.

\(\(\(Module[\ {x, y, φ}, 
			\ \ \ \ \ \ \ \ \ 
      x[φ_] = 
        3\/\@\(3 + Cos[2  φ + \(3  π\)\/4]\)\ Cos[φ]; 
      y[φ_] = 
        \(-\(3\/\@\(3 + Cos[2  φ + \(3  π\)\/4]\)\)\)\ 
          Sin[φ]; 
\ \ \ \ \ \ \ \ \ \ \ \ \ 
      ParametricPlot[{x[φ], y[φ]}, 
        \ {φ, 0, 2  π}, 

        			\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ PlotStyle -> {{Red, th}}, 

        	\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ AspectRatio -> Automatic, Ticks -> None, 

        \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 	\ \ \ \ \ \ PlotRange -> {{\(-2.3\), 4}, {\(-3.5\), 3}}, 

        	\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Prolog -> 	{{Green, Disk[{x[0], 0}, 0.2]}}, 

        \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Epilog 
          -> {\ Arrow[{0, 0}, {4, 0}], \ Arrow[{0, 0}, {0, 3}], 

            		\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {Dashing[{0.02, 0.015}], 

              		\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 		\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 	
              Line[{{0, 0}, {x[\(-\(π\/4\)\)], y[\(-\(π\/4\)\)]}}]
                \ }, 
	\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 	
            Circle[{0, 0}, 0.6, {0, π\/4}], 

            		\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 			\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
            Arrow[{0.58, 0.1}, {0.6, 0}, HeadLength ->  .035], 

            			\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {Blue, Disk[{0, 0}, 0.3]}, 

            	\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 				{Magenta, th, Circle[{0, 0}, 0.3]}, 

            				\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 	{Magenta, 

              				\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 	
              Rectangle[{\(-0.45\), \(-2.5\)}, {\(-0.15\), \(-2.1\)}]}, 

            				\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 	{Magenta, th, 
              Line[{{\(-0.3\), \(-2.1\)}, {\(-0.3\), 0}}]}, 

            \ 		\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 			
            Arrow[{x[0], \(-0.4\)}, {x[0] + 0.6, \(-0.4\)}], \ 

            			\ 		\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
            Arrow[{x[0], \(-0.4\)}, {x[0] - 0.6, \(-0.4\)}], 

            				\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {th, 
              Line[{{x[0] - 0.2, 0}, {x[0] + 0.2, 0}}]}, 

            					\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
            Circle[{0, 0}, 2.3, {\(5  π\)\/6, \(19  π\)\/20}], 

            	\ 				\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
            Arrow[{2.3  Cos[\(19  π\)\/20], 2.3  Sin[\(19  π\)\/20]}, 
              
			\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 		{2.3  Cos[\(29  π\)\/30], 2.3  Sin[\(29  π\)\/30]}], 

            \ 				\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 	
            Text[FontForm["φ", {HB, 14}], {x[0]/2 - 0.2, 0.4}], 
            
		\ 		\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 	
            Text[FontForm["m", {HB, 12}], {\(-0.8\), \(-2.3\)}], 

            				\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 	
            Text[FontForm["R", {HB, 12}], {\(-0.4\), 0.4}], 

            \ 				\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 	
            Text[FontForm["C", {HB, 16}], {\(-1.5\), 0.8}], 

            		\ 		\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 	
            Text[FontForm["M", {HB, 12}], {x[0] + 0.3, 0.4}], 

            	\ 				\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
            Text[FontForm["x", {HB, 12}], {3.8, 0.3}], 

            \ 				\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 	
            Text[FontForm["y", {HB, 12}], {\(-0.3\), 2.8}], 

            		\ 			\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
            Text[FontForm["r", {HB, 12}], {x[0]/2, \(-0.2\)}], 

            \ 				\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 	
            Text[FontForm["Fig. 2.", {HB, 12}], {2, \(-3.2\)}]

            			\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ }\ \ ]; 

      		\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ]; \)		\)\)

The Problem of Simulation of Forces

In classical mechanics the motion of a particle is determined by Newton's second law

m ((((d )^2)(Overscript[r, ->]))/(d (t^( 2))))=(Overscript[F,  ->])

where the position of the particle (Overscript[r,  ->])  is given in Cartesian coordinates. 
Here it is not assumed that this expression is a definition of the force (Overscript[F,  ->])  as it is e.g. in the case of (Overscript[p,  ->])  = m (Overscript[v,  ->]) , defining the momentum (Overscript[p,  ->])  of the particle, because forces can be determined independently on the lhs of Newton's second law.

Traditionally all forces considered in classical mechanics depend on coordinates and velocity of the particle, (Overscript[F,  ->])  = 
            (Overscript[F,  ->]) ((Overscript[r,  ->]) ,  (Overscript[v,  ->]) ) , and henceforth are well known. If, however, the particle interacts with another body whose law of motion is given then the force (Overscript[F,  ->])  may even explicitly depend on time. Thus, in the general case the force is a function of three variables : (Overscript[F,  ->])  = 
            (Overscript[F,  ->]) ((Overscript[r,  ->]) ,  (Overscript[v,  ->]) ,  t) and Newton's second law combined with initial conditions completely determines the motion of the particle.

On the other hand, if the law of motion of the particle (Overscript[r,  ->])  = (Overscript[r,  ->]) (t) is given then the resulting force (Overscript[F,  ->])  acting on the particle can be found as

(Overscript[(F(t)),  ->]) =m ((((d )^2)(Overscript[r,  ->]) )/(d (t^( 2))))

and may be presented as a function of time (Overscript[F,  ->]) (t) or as a function of velocity \(TraditionalForm
      \`\(F\& -> \)(t(\(v\& -> \))) = 
        \(F\& -> \)(\(v\& -> \))\), of coordinates or some combinations of these variables.

The very first example will demonstrate that simulation of forces which occur in the equations of motion is not a trivial problem.

Example 1 : Constant Force for particle moving along a straight line

For example, a constant force in the case of motion of a particle along a straight line can be written in the form:

(Subscript[F, x])=((m ((v^2)  -  (Subscript[B, 1])))/(2 ( x  -  (Subscript[B, 2]))))

where (Subscript[B, 1]), (Subscript[B, 2]) are some constants; (Subscript[B, 1])=((Subscript[v, 0])^2) is the initial velocity squared and (Subscript[B, 2])=(Subscript[x, 0]) is the initial coordinate. Indeed, the solution of the equation m (((d^2) x)/(d (t^2))) = m a = const with initial conditions \(TraditionalForm\`\(x\& . \)(0) = v\_0, \ \ x(0) = x\_0\) is well known

\(Clear["Global`*"]; 
Needs["Graphics`Colors`"]; \)

(Subscript[eq, 1]) = { x''[t] == a,  x'[0] == (Subscript[v, 0]),  x[0] == (Subscript[x, 0])}
            ;
(Subscript[sol, 1]) = 
        DSolve[(Subscript[eq, 1]), x[t], t][[1
          ]]

Eliminating time t from the equations x(t) = x,  \(\(x\& . \)(t) = v\) we obtain the following expression for the acceleration  a :

(Subscript[eq, 2])={( x[t]/.(Subscript[sol, 1]))==x,
		          D[ x[t]/.(Subscript[sol, 1]),t]==v}//Eliminate[#, t]&// Solve[#, a][[1]]&

It is now easy to test that in the case of a force  (Subscript[F, x]) =  ((m( (v^2) - (Subscript[v, 0]^2) ))/(2 (x - (Subscript[x, 0]))))  the solutions of the equation of motion coincide with the results obtained for constant force  (Subscript[F, x])=m a=const .

(Subscript[eq, 3])=x''[t] == 
            ((x'[t]^2 - ((Subscript[v, 0])^2))/(2 (x[t] - (Subscript[x, 0]))));
(Subscript[sol, 3]) = 
        DSolve[ (Subscript[eq, 3]), x[t], t][[1
            ]]   /. 
           {C[1] -> (Subscript[c, 1]), C[2] -> (Subscript[c, 2]) }

Equating coefficients of same power of time t in the solutions (Subscript[sol, 1]) and (Subscript[sol, 3]) and solving the set of equations obtained we find the constants (Subscript[c, 1]) and (Subscript[c, 2]).

\(eq\_4 = {
        Coefficient[sol\_1[[1, 2]], 
            t\^2] == \ 
          Coefficient[sol\_3[[1, 2]], 
            t\^2], \ 
		\ \ \ \ \ \ \ \ 
        Coefficient[sol\_1[[1, 2]], t
            \ \ ] == \ 
          Coefficient[sol\_3[[1, 2]], t
            \ \ ]}; 

    sol\_4 = \(Solve[eq\_4, {c\_1, c\_2}]\)[[1
        ]]\)

Substituting the values of (Subscript[c, 1]) and (Subscript[c, 2]) found into (Subscript[sol, 3])  we reduce it to (Subscript[sol, 1]) where we started from.

x[t] /. (Subscript[sol, 3]) /. (Subscript[sol, 4])  // Expand

Resume

Thus, it is evident that for given law of motion (Overscript[r,  ->])  = (Overscript[r,  ->]) (t) we can find the force (Overscript[F, ⇀]) which acts on the particle as (Overscript[F,  ->])  = 
            m ((((d )^2) (Overscript[r,  ->]) )/(d (t^2))). However, this force will not be defined as a function of coordinates and velocity in a unique way; it is not obvious what is the source for that particular force. 
Therefore the question arises whether it is possible to simulate a force (Overscript[F, ⇀])(t) with arbitrary time dependence using some standard forces and ideal constraints. Essentially this turns out to be a problem of simulation of forces. Subsequently, we will consider several examples which  solve this problem.

Example 2 : Particle on a Track

First we consider a well-known example of simulation of forces in the case of a one-dimensional motion of the particle on a given track y = f(x) under the influence of gravitational force.

Let a particle of mass m slide without friction on a given track y = f(x) in the xy-plane (Fig.1). There are two possibilities to describe the motion of the particle :

Representation in Cartesian coordinates x(t), y(t)

We may use Cartesian coordinates x(t) and y(t) with x being the independent coordinate.

The forces acting on this particle are the gravitational force m (Overscript[g,  ->])  and the normal force (Overscript[N, ->]). In this case Newton's second law  assumes the form

m (Overscript[x, ¨]) = N sinα ,     
      m (Overscript[y, ¨]) = -m g - N cosα

Taking into account the equation \(tan\ α = \(\(d\ y\)\/\(d\ x\) = y\& . \/x\& . \)\) and eliminating N/m we obtain

(OverDot[y]) (Overscript[y, ¨])+(OverDot[x]) (Overscript[x, ¨])+g (OverDot[y])=0

Integrating eqn. (1) we get the following equation which is the conservation law for the total energy of the system

(m/2) \((x\& . \^2 + y\& . \^2)\)+m g y=E=(m/2) \((x\& . \_0\^2 + y\& . \_0\^2)\)+m g (Subscript[y, 0]),

Hence, the differential equation for the function  y = y(t) can be rewritten in the form

((d y)/(d t))=±(Sqrt[ ((Subscript[(OverDot[x]), 0])^2)+((Subscript[(OverDot[y]), 0])^2)+2  g ((Subscript[y, 0]) - y)-((OverDot[x])^2)])

The solution of this equation  y = y(t)  together with x = x(t)  determines the particle's trajectory  y = y(x) in parametric form.

Thus, if the particle slides along the resulting curve under the influence of a gravitational force its coordinate x will change according to the law x = x(t). Of course, a solution of equations (1) and (2) only exists in some special cases of the function x(t).

Parametrization using arc length s(t)

Here we use as parametrization the arc length s(t) along the track.

If, instead of  x(t) the function  s = s(t) is given (where s is the arc length of the curve from some initial point) the problem of the simulation of forces is solved in a similar way. From the conservation of total energy of the system we obtain

\(\(v\^2 = 
      \(s\& . \^2 = 
        v\& . \_0\^2 + 
          2\ g\ \((y\_0 - y)
              \)\)\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \)\)

Hence, the coordinate y(t) of the particle changes according to the following equation

y (t)=(Subscript[y, 0])+((((Subscript[v, 0])^2)-((OverDot[s])^2))/(2  g))

The function x(t) which determines the curve in parametric form can be found as solution of the following differential equation

(OverDot[x])=±(Sqrt[((OverDot[s])^2)-((OverDot[y])^2)])=±(OverDot[s]) (Sqrt[1-(((((Overscript[s, ¨])/g)))^2)])

where \(TraditionalForm\`s\& . \) (Overscript[s, ¨]) = -g \(TraditionalForm\`y\& . \) is used which follows from the conservation law of total energy above.

Solution with Mathematica

a) Uniformly accelerated motion

At first we consider the simplest case of a uniformly accelerated motion of the particle.

(1) Representation in Cartesian coordinates x(t),y(t)

This means the Cartesian coordinate x of the particle changes according to the well-known equation : x=((a (t^( 2)))/2).

We define the equation for function y(t) and try to find its solution.

\(Clear["Global`*"]; \ 
Needs["Graphics`Colors`"]; \)

(Subscript[x, 1])[t_] = ((a (t^2))/2);

eq=x'[t]^2 + y'[t]^2 + 2 g  y[t]==(Subscript[E, 0]);

It is evident that without any restrictions we may choose the initial conditions such as (Subscript[y, 0]) = 0 ,  \(TraditionalForm\`y\& . \_0 = v\_0\). So we rewrite the equation above as

initial = {y[0] -> 0, y'[0] -> (Subscript[v, 0])};  
(Subscript[eq, 1])=eq  /.x->((Subscript[x, 1])[#]&);  
(Subscript[eq, 1]) = 
        (Subscript[eq, 1]) /. 
          Solve[ (Subscript[eq, 1]) /. t -> 0 /. initial, (Subscript[E, 0]) ]
              [[1]]

Obviously DSolve can not find the solution of (Subscript[eq, 1]). Thus, according to the standard procedure for solving differential equations which can not be solved directly we transform the differential equation using the following replacement rule y'(t)-> p(t).

(Subscript[eq, 2])=(Subscript[eq, 1]) /.y'[t] -> p[t];
\(eq\_3 = 
        \(p'\)[t] == 
          \((\(p'\)[t] /. 
		\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
              \(Solve[\ 
                  p[t] == \ 
                    D[\ y[t] /. \ 
                        \(Solve[eq\_2, y[t]]\)[[1
                          ]], t], 

                  				\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
                  \(p'\)[t]]\)[[1]])
            \)\)

Equation (Subscript[eq, 3]) cannot be solved with DSolve too. But (Subscript[eq, 3]) is a homogeneous equation which will be rewritten with the replacement rule p(t) -> t · z(t) into

(Subscript[eq, 4]) = 
      Solve[ (Subscript[eq, 3])  /. p -> (# z[#]&),  z'[t] ]
          [[1]]

The resulting equation (Subscript[eq, 4]) has separable variables and may have several solutions. The solution corresponds to the case where the numerator of rhs (Subscript[eq, 4]) equals zero.

(Subscript[sol, 4]) = 
      Solve[ (Subscript[eq, 4])[[1, 2, 4]] == 0, 
        z[t] ]

The corresponding function y(t) is

{(Subscript[y, 11])[t_] = 
        (∫(
                t z[t] /. (Subscript[sol, 4])[[1]])
                 \[DifferentialD]t) ,  
      (Subscript[y, 12])[t_] = 
        (∫(
                t z[t] /. (Subscript[sol, 4])[[2]])
                 \[DifferentialD]t)}

Then, the equation for the trajectory y = y(x) we are looking for is found by eliminating time t

\({sol\_51 = 
          \(Solve[\ Eliminate[{x\_1[t] == x, y\_11[t] == y}, t], y]
              \)[[1]]\ , 
\ \ 
        sol\_52 = 
          \(Solve[\ Eliminate[{x\_1[t] == x, y\_12[t] == y}, t], y]
              \)[[1]]\ } // Flatten\)

But there are additional solutions (in (Subscript[eq, 4])). In order to find them we regard t as a function of z,        i.e. t = t(z).

(Subscript[eq, 6])=t'[z] == 
          (1/(Subscript[eq, 4])[[1, 2]]  /. 
                z[t] -> z /. t -> t[z])

and obtain the solution

(Subscript[sol, 6]) = 
      (DSolve[(Subscript[eq, 6]), t[z], z][[1
              ]] // Simplify ) /. 
        C[1] -> (Subscript[c, 1    ])

The constant (Subscript[c, 1]) can be defined with the condition: \(\(y\& . \)(t\_0) = v\_0\). We assume that for (Subscript[t, 0]) the parameter z((Subscript[t, 0]))=(Subscript[z, 0]). Afterwards we determine the value of (Subscript[z, 0]) corresponding to the starting point of the particle's motion (Subscript[t, 0]) = 0.

\(v\_y[z_] = z\ t[z] /. sol\_6; 
eq\_0 = v\_0 == v\_y[z\_0]; 

    c\_0 = \(Solve[eq\_0, c\_1\ ]\)[[1
        ]]\)

Then there results a dependence of time t on the parameter z.

τ[z_] = t[z] /. (Subscript[sol, 6]) /. (Subscript[c, 0])

Now we investigate the function τ(z) for different values of the particle's acceleration a. It should be mentioned that if the particle slides (without friction) down an inclined straight line under the influence of gravitational force then the maximum acceleration in horizontal direction is given by (g/2).

But for a curved trajectory we can even obtain values for the acceleration a > (g/2). In this case  τ(z) is defined for any value of z.

In order to determine which value of the parameter z corresponds to the initial instant t = 0 we set e.g. (Subscript[z, 0])=-10 and plot the corresponding graph of τ(z) for different values of (Subscript[v, 0])=-10 resp. -20.

\(\(r\_1 = {g -> 9.8, a -> 10}; \ 
    th\  = \ AbsoluteThickness[2]; 


    Plot[{\(\(τ[z] /. z\_0 -> \(-10\)\) /. \ r\_1\) /. 
          v\_0 -> \(-10\), \ \ \ 
\ \ \ \ \ \ \ \ 
        \(\(τ[z] /. z\_0 -> \(-10\)\) /. r\_1\) /. 
          v\_0 -> \(-20\)}, {z, \(-30\), 30}, 
\ \ \ \ \ \ \ \ \ \ 
      AxesLabel -> {"z", "t"}, 
	\ \ \ \ \ \ \ \ 
      PlotStyle -> {{th, Red}, {th, Magenta}}, 
\ \ \ \ \ \ \ \ \ \ 
      Epilog -> {Text[\*""\!\(v\_0\)=-10"", {35, 0.8}], 
          Text[\*""\!\(v\_0\)=-20"", {35, 1.6}]}\ \ ]; \ \ \)\)

We see that in the case of  a > (g/2) the motion of the particle may be simulated only for a finite time interval depending on the velocity (Subscript[v, 0]) ; the value z = 0 corresponds to the maximum time possible. It means that the final velocity of the particle in y-direction (\(TraditionalForm\`y\& .  = t\ z\)) is zero too and the corresponding tangent to the curve is horizontal.

For the initial time (Subscript[t,  0 ])= 0  the parameter z should tend to ± ∞. The two branches of the function t(z) correspond to two possible signs of \(TraditionalForm\`y\& . \) in eq. (2).

Case (Subscript[z, 0])  -> -∞

If (Subscript[v, 0])<0 then for (Subscript[t,  0]) = 0 we should choose (Subscript[z, 0]) -> -∞ because \(TraditionalForm\`y\& .  = t\ z\) and  t ≥ 0. In this case the velocity (Subscript[v, 0]) will be the initial velocity of the particle in y-direction and we can define the function t(z) as

(Subscript[t, 2])[z_] = 
      Limit[ τ[z], (Subscript[z, 0]) -> -∞] // PowerExpand

Then taking into account the initial condition (Subscript[y, 0]) = 0 and the equation ((d y)/(d t)) = (1/(t ' (z))) ((d y)/(d z)) we define the y-coordinate of the particle as function of the parameter z

\(y\_2[z1_] = 
      ∫\_\(-∞\)\%z1
            \((z\ \((
                    \(t'\)[z] /. 
                      \(Solve[eq\_6, \(t'\)[z]]\)[[1
                        ]])\)\ t[z] /. 
                t[z] -> t\_2[z])\) \[DifferentialD]z\  // 

        		\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
        PowerExpand\)

Now we can visualize the function y(t).

\(\(ParametricPlot[{{\(t\_2[z] /. r\_1\) /. v\_0 -> \(-10\), 
          \(y\_2[z] /. r\_1\) /. v\_0 -> \(-10\)}, 

        \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {
          \(t\_2[z] /. r\_1\) /. v\_0 -> \(-20\), 
          \(y\_2[z] /. r\_1\) /. v\_0 -> \(-20\)}}, 
      
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {z, 
        \(-300\), 0}, 
	\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
      PlotStyle -> {{th, Red}, {th, Blue}}, 

      \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
      AxesLabel -> {"t", "y"}, PlotRange -> All]; \)\)

and finally plot the trajectory y(x) of the particle.

(Subscript[x, 2])[z_] =  (Subscript[x, 1])[(Subscript[t, 2])[z]];

\(\(p\_2 = 
      ParametricPlot[{{\(x\_2[z] /. r\_1\) /. v\_0 -> \(-10\), 
            \(y\_2[z] /. r\_1\) /. v\_0 -> \(-10\)}, \ 

          				\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {
            \(x\_2[z] /. r\_1\) /. v\_\(0\ \) -> \(-20\), 
            \(y\_2[z] /. r\_1\) /. v\_0 -> \(-20\)}}, 

        		\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {z, 
          \(-500\), 0}, AxesLabel -> {"x", "y"}\ , 

        	\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
        PlotStyle -> {{th, Red}, \ {th, Blue}}\ ]; \)\)

Case (Subscript[z, 0])  -> +∞

Choosing (Subscript[z, 0]) -> ∞ for t = 0  we consider the case of (Subscript[v, 0]) > 0 in a similar way.

\(t\_3[z_] = Limit[τ[z], z\_0 -> ∞] // PowerExpand; 

    x\_3[z_] = \ x\_1[t]\  /. t -> t\_3[z]; 

    y\_3[z_] = 
      ∫\_∞\%z
            \((z\ \((
                    \(t'\)[z] /. 
                      \(Solve[eq\_6, \(t'\)[z]]\)[[1
                        ]])\)\ t[z] /. 
                t[z] -> t\_3[z])\) \[DifferentialD]z\  // 

        		\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ PowerExpand\)

\(\(p\_3 = 
      ParametricPlot[{{\(x\_3[z] /. r\_1\) /. v\_0 -> 20, 
            \(y\_3[z] /. r\_1\) /. v\_0 -> 20}, 
          
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {\(x\_3[z] /. r\_1\) /. v\_\(0\ \) -> 30, 
            \(y\_3[z] /. r\_1\) /. v\_0 -> 30}}, 

        		\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {z, 
          0, 500}, \ AxesLabel -> {"x", "y"}, \ \ 

        		\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
        PlotStyle 
          -> {{th, RGBColor[1, 0, 0.7]}, 

            				\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {th, RGBColor[0.7, 0, 1]}}\ ]; 
    \)\)

Thus in the case of  a > (g/2) we have a lot of possible tracks which depend on the value of initial velocity (Subscript[v, 0])≠0 chosen. As \(\(x\_0\)\& . \  = 0\) and  (OverDot[(Subscript[y, 0])])≡(Subscript[v, 0])≠0 the tangent on the curve y =y(x) in the initial point should be a vertical line (y-axis). Several curves are shown below.

Show[(Subscript[p, 2]), (Subscript[p, 3]), PlotRange -> {{0, 8}, {-30, 30}}];

Considering the case of  a<(g/2) we should take into account that t(z) is a real function of z. Hence three possibilities are analyzed when the parameter (Subscript[z, 0]) belongs to one of the following intervals: (i) (Subscript[z, 0]) < (Subscript[z, 1]) , (ii) (Subscript[z, 1]) < (Subscript[z, 0]) < (Subscript[z, 2]) or (iii) (Subscript[z, 2]) < (Subscript[z, 0]) , where (Subscript[z, 1]) and (Subscript[z, 2]) are the two roots of the equation (a^2) + z  (g + z) = 0 . These intervals result from the conditions where the radical under the square root in the denominator of τ(z) vanishes

\(sol\_7 = Solve[z\^2 + g\ z + a\^2 == 0, z]; 
    
{z\_1 = z /. sol\_7[[1]]\ , \ 
      z\_2 = z /. sol\_7[[2]]\ }\)

In order to estimate the behavior of t(z) for (Subscript[z, 0]) < (Subscript[z, 1]) we calculate this function for different z-values.

\(r\_2 = {g -> 9.8, a -> 4}; 

    t\_4[z_] = \(τ[z] /. r\_2\) /. z\_0 -> \((z\_1 /. r\_2)\) - 10; 
    
Table[t\_4[z] /. v\_0 -> \(-10\), {z, \((z\_1 /. r\_2)\) - 5, 
          \((z\_2 /. r\_2)\) + 5, 0.5}] // Chop\)

We see that t(z) is a real function of z only inside the interval z < (Subscript[z, 1]). The corresponding graphs of this function for different values of (Subscript[v, 0]) are shown below.

\(\(Plot[{\ t\_4[z]\  /. v\_0 -> \(-10\), 
		\ \ \ \ \ \ \ \ \ 
        t\_4[z]\  /. v\_0 -> \(-40\)\ \ }, {z, \((z\_1 /. r\_2)\) - 10, 
        \((z\_1 /. r\_2)\)}, 
	\ \ \ \ \ \ \ \ 
      AxesLabel -> {"z", "t"}, \ 
      PlotStyle -> {{th, Blue}, {th, Magenta}}]; \)\)

As in the previous case a > (g/2) the starting point of the particle's motion at (Subscript[t, 0]) = 0 corresponds to (Subscript[z, 0]) -> -∞. But now  t -> ∞ when z -> (Subscript[z, 1]) and hence in contrast to the previous case we can simulate an infinite motion of the particle.

In order to define the function t(z) we again calculate the limit of τ(z) at (Subscript[z, 0]) -> -∞ . But now we cannot apply the function PowerExpand to the result because there the radical 4 (a^2) - (g^2) < 0 under the square root. It should be emphasized that applying  PowerExpand to an expression (Sqrt[(1/Sign)[a]]) (Sqrt[Sign])[a] may turn out to be incorrect.

(Subscript[t, 4])[z_] = Limit[ τ[z], (Subscript[z, 0]) -> -∞]

It is evident that the initial velocity of the particle in y-direction should not be zero in this case and so the tangent on the curve in the initial point is a vertical line. But then the tangent tends to be parallel to the straight line which is determined by (Subscript[sol, 51]) . This follows from the next formula tan α = ((\(y\& . \))/(\(x\& . \))) :

Limit[((z t[z])/(∂_(t[z]))) (Subscript[x, 1])[t[z]], z -> (Subscript[z, 1])]

Now we can plot several curves for different values of the initial velocity (Subscript[v, 0]) of the particle.

\(x\_4[z_] = x\_1[t\_4[z]]; 


    y\_4[z1_] = 
      ∫\_\(-∞\)\%z1
          \((z\ \((
                  \(t'\)[z] /. 
                    \(Solve[eq\_6, \(t'\)[z]]\)[[1
                      ]]\ )\)\ t[z] /. 
              t[z] -> t\_4[z])\) \(\[DifferentialD]z\ \)\)

\(\(p\_4 = 
      ParametricPlot[{{\(x\_4[z] /. r\_2\) /. v\_0 -> \(-10\), 
            \(y\_4[z] /. r\_2\) /. \ v\_0 -> \(-10\)}, \ \ 

          				\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {
            \(x\_4[z] /. r\_2\) /. v\_\(0\ \) -> \(-20\), 
            \(y\_4[z] /. r\_2\) /. v\_0 -> \(-20\)}}, 

        		\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {
          z, \(-500\), \((z\_1 /. r\_2)\)}, \ \ \ \ \ 

        		\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 	\ \ \ \ 
        PlotStyle 
          -> {{th, RGBColor[1, 0, 0.8]}, 
            
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {th, 
              RGBColor[0.8, 0, 1]}}, \ \ \ \ 

        	\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ AxesLabel -> {"x", "y"}]; \)\)

It is clear that if the initial velocity of the particle (Subscript[v, 0])->0 these curves degenerate to the straight line determined by (Subscript[sol, 51]).

Now we consider the case of  (Subscript[z, 2])<(Subscript[z, 0]) in a similar way. First of all we look for the range of the parameter z where the values of the function t(z) are real numbers.

\(t\_5[z_] = \(τ[z] /. r\_2\) /. z\_0 -> \(-1\); 

    Table[t\_5[z] /. v\_0 -> \(-10\), {z, \(-10\), 0, 0.5}] // Chop\)

We see that the parameter z should belong to the interval ((Subscript[z, 2]), ∞). The corresponding graph is shown below.

\(\(Plot[{t\_5[z]\  /. v\_0 -> \(-10\), 
        t\_5[z]\  /. v\_0 -> \(-20\)}, {z, \((z\_2 /. r\_2)\), 20}, 

      	\ \ \ \ \ \ \ \ AxesLabel -> {"z", "t"}, 
      PlotStyle -> {{Blue, th}}]; \)\)

Now t(z) is a restricted function with two branches (as in the case of a > (g/2)). The positive values of z correspond to the case of a positive initial velocity ((Subscript[v, 0]) > 0). And as before, the corresponding curves y = y(x) can be easily found.

\(t\_5[z_] = Limit[τ[z], z\_0 -> ∞]; 

    x\_5[z_] = \ x\_1[t\_5[z]]; 

    y\_5[z_] = 
      ∫\_∞\%z
          \((z\ \((
                  \(t'\)[z] /. 
                    \(Solve[eq\_6, \(t'\)[z]]\)[[1
                      ]])\)\ t[z] /. t[z] -> t\_5[z])
            \) \[DifferentialD]z\ ; 

    p\_5 = ParametricPlot[{{\(x\_5[z] /. r\_2\) /. v\_0 -> \(-20\), 
            \(y\_5[z] /. r\_2\) /. v\_0 -> \(-20\)}, 

          			\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {\(x\_5[z] /. r\_2\) /. v\_\(0\ \) -> \(-40\), 
            \(y\_5[z] /. r\_2\) /. v\_0 -> \(-40\)}}, 

        	\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {
          z, 0, 500}, \ AxesLabel -> {"x", "y"}, \ \ \ \ \ 

        \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 	\ \ \ \ 
        PlotStyle 
          -> {{th, RGBColor[0, 0.8, 1]}, 
            
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {th, 
              RGBColor[0, 0.7, 1]}}\ \ ]; \)

If z is in the interval ((Subscript[z, 2]),  0) then the initial velocity of the particle should be zero (Subscript[(\(y\& . \)), 0])=(Subscript[z, 2]) t((Subscript[z, 2]))  = 0 and the tangent to the curve in this point will coincide with the straight line determined in (Subscript[sol, 52]). Actually, tan α=((\(y\& . \))/(\(x\& . \)))=(z/a) and

(z/a) /. z -> (Subscript[z, 2])

But nevertheless choosing different values of (Subscript[v, 0]) for given (Subscript[z, 0]) from the interval ((Subscript[z, 2]),  0) we get different curves shown below.

\(t\_6[z_] = τ[z]; 

    x\_6[z_] = \ \(x\_1[t\_6[z]] /. z\_0 -> \(-1\)\) /. r\_2; 

    yy = ∫\((
              z\ \((\(t'\)[z] /. 
                      \(Solve[eq\_6, \(t'\)[z]]\)[[1
                        ]])\)\ t[z] /. 
                t[z] -> t\_6[z])\) \[DifferentialD]z\  // \ 

        			\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Simplify; 

    y\_6[z_] = 
      \((\(yy /. z\_0 -> \(-1\)\) /. r\_2)\) - \ 

          			\ \ \ \ \ \ \ \ \ \ 
          Limit[\((\(yy /. z\_0 -> \(-1\)\) /. r\_2)\), 
            z -> \((z\_2 /. r\_2)\), Direction -> \(-1\)] // Chop; 
    \)

\(\(p\_6 = 
      ParametricPlot[{{x\_6[z] /. v\_0 -> \(-10\), \ 
            y\_6[z] /. v\_0 -> \(-10\)}, \ \ 

          				\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {
            x\_6[z] /. v\_\(0\ \) -> \(-12\), \ 
            y\_6[z] /. v\_0 -> \(-12\)}}, 

        		\ \ \ \ 		\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {z, 
          \((z\_2 /. r\_2)\), 0}, \ AxesLabel -> {"x", "y"}, 
        \ \ \ \ 

        	\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
        PlotStyle 
          -> {{th, RGBColor[0.9, 1, 0]}, 
            
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {th, Green}}\ ]; 
    \)\)

Finally, we consider the case where (Subscript[z, 0]) belongs to the interval ((Subscript[z, 1]),  (Subscript[z, 2])). From the table below we see that now the function t(z) is real only inside this interval.

\(t\_7[z_] = \(τ[z] /. r\_2\) /. z\_0 -> \(-4\); 

    Table[t\_7[z] /. v\_0 -> \(-10\), {z, \(-10\), 0, 0.5}] // Chop\)

The corresponding graph is shown below.

\(\(Plot[{t\_7[z]\  /. v\_0 -> \(-10\), 
        t\_7[z]\  /. v\_0 -> \(-20\)}, {z, \((z\_1 /. r\_2)\), 
        \((z\_2 /. r\_2)\)}, 
	\ \ \ \ \ \ \ \ 
      AxesLabel -> {"z", "t"}, PlotStyle -> {{Blue, th}}]; 
    \)\)

As in the previous case the tangent on the curve in the initial point will coincide with the straight line defined by (Subscript[sol, 52]) but now this curve will be infinite and tends to be parallel to another straight line defined by (Subscript[sol, 51]). It should also be mentioned that there are several curves corresponding to different values of (Subscript[v, 0]). Two of them are shown below.

\(x\_7[z_] = \ \(x\_1[τ[z]] /. z\_0 -> \(-4\)\) /. r\_2; 

    y\_7[z_] = 
      \((\(yy /. z\_0 -> \(-4\)\) /. r\_2)\) - \ 

          \ \ \ \ \ \ \ \ \ \ \ \ 		
          Limit[\((\(yy /. z\_0 -> \(-4\)\) /. r\_2)\), \ 
            z -> z\_2 /. r\_2] // Chop\)

\(\(p\_7 = 
      ParametricPlot[{{x\_7[z] /. v\_0 -> \(-10\), 
            y\_7[z] /. v\_0 -> \(-10\)}, \ 

          			\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {x\_7[z] /. v\_\(0\ \) -> \(-30\), y\_7[z] /. v\_0 -> \(-30\)}}, 

        		\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {z, 
          \((z\_1 /. r\_2)\), \((z\_2 /. r\_2)\)}, \ \ \ \ \ 

        	\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
        PlotStyle 
          -> {{th, RGBColor[0.9, 0.6, 0]}, 
            
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {th, 
              RGBColor[0, 0.7, 0]}}, \ \ \ \ 

        	\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
        AxesLabel -> {"x", "y"}]; \)\)

Thus we have found all possible curves which are determined by eq. (2) in the case of an uniformly accelerated motion of the particle. If the initial velocity of the particle in y-direction (Subscript[v, 0])  ≠  0 then eq. (2) has two solutions in accordance with two possible signs of the derivative \(TraditionalForm\`y\& . \) in eq. (2). The following graph shows several curves which were obtained for different values of the initial velocity (Subscript[v, 0]) and acceleration a.

Show[(Subscript[p, 2]), (Subscript[p, 3]), (Subscript[p, 4]), (Subscript[p, 5]), 
      PlotRange -> {{0, 50}, {-100, 75}}];

If initial velocity (Subscript[v, 0])  = 0 then eq. (2) has many solutions. The corresponding curves have a common tangent in the initial point which is a straight line determined by (Subscript[sol, 52]) and which is also a solution of eq. (2). The curves above this straight line are finite and have a horizontal tangent in their final points. The remaining curves are located between two straight lines which are determined by (Subscript[sol, 51]) and (Subscript[sol, 52]) and tend to become parallel to the first one at infinity. Several of these curves are shown in the next graph together with the straight lines which corresponding the solutions (Subscript[sol, 53]) and (Subscript[sol, 54]) (in violet). If the particle slides along any of these curves without friction then its x-coordinate will change according to the  x(t)=((a (t^( 2)))/2).

\(p\_8 = 
      Plot[{\(y /. sol\_51\) /. r\_2, \ \(y /. sol\_52\) /. r\_2}, {x, 0, 
          500}, 
	\ \ \ \ \ \ \ \ \ \ \ \ \ 	
        PlotStyle 
          -> {{th, RGBColor[0.8, 0, 0.8]}, {th, RGBColor[0.8, 0, 1]}}, 
        
\ \ \ \ \ \ \ \ \ \ \ \ \ 
        DisplayFunction -> Identity\ ]; 

    Show[p\_6, p\_7, p\_8, 
	\ \ \ \ \ \ \ \ 
      PlotRange -> {{0, 250}, {\(-300\), 0}}, 
      
\ \ \ \ \ \ \ \ 
      DisplayFunction -> $DisplayFunction]; \)

(2) Parametrization using arc length s(t)

If instead of x(t) we choose the function  s = s(t)  where s is the arc length of the curve then the corresponding functions x(t) and y(t) are determined by the equations (3), (4).

In order to compare the results we solve these equations too in the case of an uniformly accelerated motion. At first we define equations (3) and (4).

eq=x'[t] == s'[t] (Sqrt[1 - (((s''[t]/g))^2)]);
y[t_] := ((s'[0]^2 - s'[t]^2)/(2  g));

If  s=((a (t^2))/2) then equations (3) and (4) become

(Subscript[s, 1])[t_] := ((a (t^2))/2);

(Subscript[eq, 1])=eq/.s->((Subscript[s, 1])[#]&) (Subscript[y, 1])[t_]=y[t]/.s->((Subscript[s, 1])[#]&)

The  solution of the equation (Subscript[eq, 1]) can be easily found with DSolve.

\(sol\_1 = 
      \(DSolve[{eq\_1, x[0] == 0}, x[t], t]\)[[1
        ]]; 
x\_1[t_] = x[t] /. sol\_1\)

Obviously, the corresponding curve is a straight line.

\(\(\(Eliminate[{x\_1[t] == x, \ y\_1[t] == y}, t] // 
          \(Solve[#, y]\)[[1]]&\) // 

        	\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ PowerExpand\) // \ \ Simplify\)

Obviously , the trajectory  y(x)=-((a x)/(Sqrt[(g^2) - (a^2)]))  is  a straight line  because it is linear in x  .

b) Harmonic oscillations

Let us consider the case of  x = a sin(ω t).

(1) Representation in Cartesian coordinates x(t),y(t)

The Cartesian coordinate x of the particle describes a harmonic motion.

Clear["Global`*"];

(Subscript[x, 2])[t_] = a Sin[ω t];

eq=x'[t]^2 + y'[t]^2 + 2 g y[t] == (Subscript[E, 0]);

As before we choose the initial conditions as (Subscript[y, 0]) = 0, y' (0) = 0 so that the equation eq can be rewritten as

initial = {y[0] -> 0, y'[0] -> 0};
(Subscript[eq, 1])=eq/.x->((Subscript[x, 2])[#]&);
(Subscript[eq, 1]) = 
        (Subscript[eq, 1]) /. 
          Solve[(Subscript[eq, 1])  /. t -> 0  /. initial, (Subscript[E, 0])]
              [[1]]

The corresponding equations which determine the function y(t) are given below.

(Subscript[sol, 2]) = Solve[ (Subscript[eq, 1]), y'[t]]  // Simplify

Because no analytial solutions of these equations are found we try numerical solutions. It should be taken into account that the radicant must be positive. Thus at  t = (π/(2  ω))  the derivative y'(t) should be zero. Then it changes its sign and the particle begins to move along the same trajectory. Hence, it is sufficient to calculate the trajectory only for this time intervall.

Case y'(t) >0

At first we consider the case of a positive value of the derivative y'(t) >0.

\(r\_1 = {g -> 9.8, a -> 3, ω -> 1}; 

    eq\_3 = \(y'\)[t] == 
        \((sol\_2[[2, 1, 2]] /. r\_1)
          \)\)

Solving this differential equation the calculations require higher precision in order to avoid that the radicant becomes negative. Thus, the values of PrecisionGoal and WorkingPrecision in DSolve are increased.

\(\(sol\_3 = 
      \(NDSolve[{eq\_3, y[0] == 0}, y, {t, 0. , π/2}, \ 

            \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 		
            PrecisionGoal -> 25, WorkingPrecision -> 20]
            \)[[1]]\  // Chop\ ; \)\)

The trajectory y = y(x) found in this way is shown below. Plotting this curve (which is actually the particle's trajectory in Fig. 1)  we have taken into account that it should be symmetric with respect to the y-axis.

\(th\  = \ AbsoluteThickness[2]; 

    p\_1 = ParametricPlot[{{\ x\_2[t] /. r\_1, \((y[t] /. sol\_3 // Chop)\)}, 
          
		\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {\ \(-x\_2[t]\) /. r\_1, \((y[t] /. sol\_3 // Chop)\)}}, 

        		\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 			\ \ \ \ \ \ \ {t, 0, 
          π/2.01}, 

        \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
        PlotStyle -> {{th, Red}}, AxesLabel -> {"x", "y"}]; 
    \)

Projecting the particle's movement on the horizontal axis there results a harmonic oscillation as wanted.

Case y'(t) < 0

In the case where y'(t)<0 we do not encounter any problems to obtain the solution. The corresponding curve is shown below.

\(r\_1 = {g -> 9.8, a -> 3, ω -> 1}; 

    eq\_4 = \(y'\)[t] == 
        \((sol\_2[[1, 1, 2]]\  /. r\_1)
          \); 
sol\_4 = 
      \(NDSolve[{eq\_4, y[0] == 0}, y, {t, 0. , 50}]\)[[1
        ]]; 

    ParametricPlot[{x\_2[t] /. r\_1, \ 
          y[t] /. sol\_4[[1]]} // 
        Evaluate, {t, 0, 20}, 

      	\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
      PlotStyle -> {{th, Blue}}, AxesLabel -> {"x", "y"}]; 
    \)

In this case the movement of the particle is still oscillatory with respect to the x-axis, but y'(t) is always negative because the particle moves down on the trajectory due to gravitational forces.

(2) Parametrization using arc length s(t)

Now let us consider the case of harmonic oscillations with respect to the arc length; i.e.  s = (Subscript[s, 0]) sin(ω t). At first we define the equations (3) and (4) which determine the trajectory.

eq=x'[t] == s'[t] (Sqrt[1 - (((s''[t]/g))^2)]);
       y[t_] := ((s'[0]^2 - s'[t]^2)/(2 g));

Substituting the corresponding function s(t) into the second equation the y-coordinate of the curve is found as a function of time.

\(s\_2[t_] := \(s\_0\) Sin[ω\ t]; 

    y\_2[t_] = y[t]\  /. s -> \((s\_2[#]&)\) // Simplify\)

In order to find the dependence x(t) the following differential equation must be solved

(Subscript[eq, 1]) = eq  /. s -> ((Subscript[s, 2])[#]&)

Use of  DSolve gives

(Subscript[sol, 1]) = 
      (DSolve[(Subscript[eq, 1]), x[t], t][[1
              ]] //. 
            Cos[2 x_] -> (1 - 2  Sin[x]^2) )  /.  
        C[1] -> (Subscript[c, 1])

The constant (Subscript[c, 1]) is determined from the initial conditions.

(Subscript[c, 11]) =  
      Solve[(
              (Subscript[sol, 1])[[1, 2]] /. 
                t -> 0) == 0, (Subscript[c, 1])][[1
        ]]

Introducing a new parametrization according to the equation  \(TraditionalForm
      \`\(s\_0\ ω\^2\)\/g\ \(Sin(ω\ t)\) = \ 
        Sin(φ)\) we rewrite the solution x(t) in the following form.

\(xφ = 
      \(\(\(\(x[t]\  /. sol\_1\)\  /. 
                Sin[t\ ω] -> 
                  g\ \(Sin[φ]/ω\^2\)/s\_0\  // 

              		\ \ \ \ \ \ \ \ \ \ \ \ Simplify\) // PowerExpand\) // 
          TrigReduce\) // Expand\)

The corresponding function y(φ) is

yφ = 
      (Subscript[y, 2])[t]  /. 
            Sin[t ω] -> 
              g Sin[φ]/(ω^2)/(Subscript[s, 0])  // TrigReduce // 
        Expand

It is easy to test that the parametric curves x(φ) and y(φ)  are oscillating functions and determine a cycloid in parametric form.

\(\(Plot[{
        \(\(x[t] /. sol\_1[[1]]\) /. 
            c\_11\) /. {ω -> 1, s\_0 -> 2, g -> 10}, \ 
        
\ \ \ \ \ \ \ 		
        y\_2[t] /. {ω -> 1, s\_0 -> 2, g -> 10}}, {t, 
        0, 10}, 
	\ \ \ \ \ \ \ \ 
      PlotStyle -> {{th, Magenta}, {th, Blue}}, 
\ \ \ \ \ \ \ \ \ \ 
      Epilog -> {\ Text["x(t)", {9.3, 1.5}], 
          Text[\*""\!\(y\_2\)(t)"", {8. ,  .5}]\ }\ ]; \)\)

The corresponding trajectory is shown below.

\(\(p\_2 = 
      ParametricPlot[{
          \(\(x[t] /. sol\_1[[1]]\) /. 
              c\_11\) /. {ω -> 1, s\_0 -> 2, g -> 10}, 
          
		\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 	
          y\_2[t] /. {ω -> 1, s\_0 -> 2, g -> 10}}, {t, 
          0, 10}, 

        \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 		
        PlotStyle -> {{th, Green}}, 
        AxesLabel -> {"x", "y"}]; \)\)

In order to compare the shape of the two curves x(t) and s(t) which correspond to the oscillaion of the particle in different representations we display them within one graph.

Show[(Subscript[p, 1]), (Subscript[p, 2]),  PlotRange -> {0,  .21}];

The plot shows that the curves are very similar but of course they are (slightly) different curves due to the different approach.

Example 3 : Particle on a Rotating Contour

A wider range of situations can be simulated which is schematically shown in Fig. 2. A body of mass M moves along the horizontal x-axis without friction and at the same time slides on a rotating contour 𝒞. The contour 𝒞  in the xy- plane is shown in Fig.2 (in red) and 𝒞  is given as  r = r(φ).

The contour 𝒞  is rigidly fixed to a cylinder with radius R and a weightless unextensible thread is wrapped around the cylinder with mass m on the other end. Under the influence of gravitational force the mass m  moves down along the vertical y-axis so that the contour 𝒞  rotates. As a result of the rotation of 𝒞  the body M starts moving along the horizontal x-axis according to the equation  x(t) = r(φ(t)). Therefore, if the contour 𝒞  is of an appropriate prescribed shape r = r(φ) then the body M  will move according to the desired law of motion.

It is easy to write the equation of motion for the system under consideration using the fact that the total energy of the system is a constant of motion.

E=(1/2)((Subscript[I, 0])+m (R^2)) ((OverDot[φ])^2)+(1/2)M ((OverDot[x])^2)-mg R φ=const

It is assumed that the contour 𝒞  has no mass. Differentiating (5) we have

((Subscript[I, 0])+m (R^2)) (OverDot[φ]) (Overscript[φ, ¨])+M (OverDot[x]) (Overscript[x, ¨])-mg R (OverDot[φ])=0

If the shape r(φ) of the contour 𝒞  is given then equation (5) determines the motion of the system. But for a given law of motion x=x(t) for the body M the last differential equation offers the possibility to find φ as function of t. So we may find the polar coordinates of the surface r(t) and φ(t) as functions of time t. In comparison to (1) equation (6) has some additional parameters. Thus this equation should have solutions for a wider class of functions x(t).

Solution with Mathematica

a) Uniformly accelerated motion

With this system under consideration we can also simulate an uniformly accelerated motion of the particle. It means that the law of motion should be

Clear["Global`*"];

(Subscript[x, 1])[t_] := (Subscript[x, 0]) + ((a (t^2))/2);

We define the differential equation describing the behaviour of  φ(t) and substitute the function (Subscript[x, 1])(t) into this equation.

(Subscript[eq, 1])=((Subscript[I, 0]) + m (R^2)) φ''[t] == 
            m g R - 
              M ((x'[t] x''[t])/(φ'))[t] /.x->( (Subscript[x, 1])[#]& )

To reduce the order of the differential equation we substitute  φ'-> ω.

(Subscript[eq, 2])=(Subscript[eq, 1])  /.φ -> 
            (∫ω[#] \[DifferentialD]# &)

Obviously, there results a homogeneous differential equation and in order to solve it we use the usual substitution ω-> t·z(t).

(Subscript[eq, 3])=(Subscript[eq, 2])/.ω -> ( # z[#]&);
(Subscript[sol, 3]) = 
        Solve[ (Subscript[eq, 3]), z'[t]][[1
          ]]

As in the previous example there are two solutions corresponding to a constant value of z(t) for which the numerator of the right side of (Subscript[sol, 3]) is zero. In both these cases the contour 𝒞  on the Fig. 2 will rotate with constant angular acceleration (Overscript[φ, ¨]).

\(eq\_4 = sol\_3[[1, 2, 4]] == 0; 

    sol\_4 = Solve[eq\_4, z[t]]\)

The corresponding functions for the rotation angle φ(t) are shown below.

\(φ\_1[t_] = 
      ∫\((
            t\ z[t] /. sol\_4[[1]])\) 
          \[DifferentialD]t; 

    φ\_2[t_] = 
      ∫\((
            t\ z[t] /. sol\_4[[2]])\) 
          \[DifferentialD]t\)

In order to visualize these solutions we introduce some numerical values for the system and define the corresponding functions.

(Subscript[rule, 1]) = {m -> 2,  M -> 0.1,  (Subscript[I, 0]) -> 10, 
        g -> 9.8, R -> 1, a -> 4, (Subscript[x, 0]) -> 2};

\(φ\_11[t_] = φ\_1[t] /. rule\_1; 

    φ\_12[t_] = φ\_2[t] /. rule\_1; 

    x\_11[t_] = x\_1[t] /. rule\_1; \)

\(th\  = \ AbsoluteThickness[2.5]; 

    Plot[{φ\_12[t], x\_1[t] /. rule\_1}, \ {t, 0, 5}, 

      		\ \ \ \ \ \ \ PlotStyle -> {{th, Red}, {th, Blue\ }}, 

      \ \ \ \ \ \ \ \ \ \ AxesLabel\  -> {"t", ""}, 

      \ \ \ \ \ \ \ \ \ \ 
      Epilog -> {Text[\*""\!\(φ\_2\)(t)"", {5.2, 18}], 
          Text[\*""\!\(x\_1\)(t)""\ , {5.2, 50}]}\ ]; \)

The solution (Subscript[φ, 2])(t) (red curve) for given (Subscript[x, 1])(t) = (Subscript[x, 0]) + ((a (t^2))/2) (shown in blue) is the law of motion for the mass M in horizontal direction.

To find another solution of the equation (Subscript[eq, 3]) we rewrite this equation regarding t as a function of z.

(Subscript[eq, 6]) = 
      t'[z] == 
        (1/z'[t] /. (Subscript[sol, 3]) /. z[t] -> z /. t -> t[z])

Then DSolve easily finds the solution which is quite involved.

(Subscript[sol, 6]) = 
      (DSolve[ (Subscript[eq, 6]), t[z], z] // Flatten) /.  C[1] -> (Subscript[c, 1]) // 
        Simplify

The constant (Subscript[c, 1]) is determined by the initial conditions as in the previous example where a similar function t(z) was analyzed in detail. Due to the fact that our system moves under the action of the gravitational force applied to the mass m we only consider the case of positive values of angular velocity  ω= t·z(t). This means that the parameter z should not be negative (z ≥ 0).

On the other side, it is easy to check that in case of m g R  < 2 a (Sqrt[M ((Subscript[I, 0]) + (mR^2))]) the uniformly accelerated motion of the body M can be simulated only for a finite interval of time. So we consider the opposite case where m g R  > 2 a (Sqrt[M ((Subscript[I, 0]) + (mR^2))]).

Then, as in example 2, function t(z) will be real only if the parameter z is inside the interval (Subscript[z, 1])< z < (Subscript[z, 2]) where (Subscript[z, 1]), (Subscript[z, 2]) are the roots of the quadratic equation   \(TraditionalForm
      \`\((\ I\_0 + m\ R\^\(\ 2\))\)\ z\^\(\ 2\) - m\ g\ R\ z + 
          M\ a\^\(\ 2\) = 0\).

\(sol\_7 = 
      Solve[a\^2\ M - g\ m\ R\ z + I\_0\ z\^2 + m\ R\^2\ z\^2 == 0, z]; 
    
{\ z\_1 = 
        \(z /. sol\_7[[1]]\) /. 
          rule\_\(1\ \), 
      z\_2 = \(z /. sol\_7[[2]]\) /. 
          rule\_1\ }\)

Thus, using solution (Subscript[sol, 6]) found then we define the angle φ as function of z.

(Subscript[φ, 3])[z_] = 
      ∫(
            z (t'[z] /. 
                    Solve[(Subscript[eq, 6]), t'[z]][[1
                      ]]) t[z] /. (Subscript[sol, 6])) 
          \[DifferentialD]z

And choosing different values of the constant (Subscript[c, 1]) we shall get a lot of solutions of (Subscript[eq, 6]).
To visualize these solutions we define the corresponding functions.

\(φ\_4[z_]\ \  = 
      \(φ\_3[z]\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \  /. 
            rule\_1\)\  /. c\_1 -> I\  // Chop; 

    x\_12[z_] = 
      \(\(x\_1[t[z]] /. sol\_6\)\  /. rule\_1\) /. 
        c\_1 -> I; \)

From the graphics shown below we can see the shape of the contour 𝒞  for three cases. Of course, the contour 𝒞  may be realized by a wire which is shaped in accordance with the solutions found above.

\(p\_6 = 
      ParametricPlot[{{\ x\_11[t]\ Cos[φ\_11[t]], 
            \(-x\_11[t]\)\ Sin[φ\_11[t]]}, 

          			\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {\ x\_11[t]\ Cos[φ\_12[t]], 
            \(-x\_11[t]\)\ Sin[φ\_12[t]]\ }}, 

        	\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 	{
          t, 0, 5}, 

        		\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
        AspectRatio -> Automatic, 

        	\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ PlotStyle -> {{th, Red}, {th, Blue}}, 

        	\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ DisplayFunction -> Identity]; 

    p\_7 = ParametricPlot[{x\_12[z]\ Cos[φ\_4[z]], 
          \(-x\_12[z]\)\ Sin[φ\_4[z]]}, 

        	\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {
          z, z\_1,  .95\ z\_2}, 

        		\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
        AspectRatio -> Automatic, \ 

        	\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
        PlotStyle -> {{th, Green}}, 

        		\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
        PlotPoints -> 100, DisplayFunction -> Identity\ \ ]; 

    Show[\ p\_6, p\_7, DisplayFunction -> $DisplayFunction]; \)

The red and blue spirals correspond to contours 𝒞  rotating with (Overscript[φ, ¨]) = const. Between these two (limiting) spirals there are additional solutions situated; one example is shown in green.

b) Harmonic oscillations

This system also permits us to simulate harmonic oscillations of the body m. In order to find the shape of the contour 𝒞  we define the corresponding dependence x(t) and solve equation (3) numerically.

Clear["Global`*"];

(Subscript[x, 2])[t_] = 5 + 3 Sin[2 π t];

(Subscript[eq, 8])=((Subscript[I, 0]) + m (R^2)) φ''[t] == 
            m g R - 
              M ((x'[t] x''[t])/(φ'))[t] /.x->((Subscript[x, 2])[#]&)

\(rule\_1 = {m -> 2, \ M -> 0.1, \ I\_0 -> 10, 
        g -> 9.8, R -> 1, a -> 4, x\_0 -> 2}; 

    initial\_1 = {φ[0] == 0, \(φ'\)[0] == 1}; 

    eq\_9 = \ eq\_8 /. rule\_1\ ; 

    sol\_9 = NDSolve[\ Join[{eq\_9}, initial\_1\ ], 
          φ, {t, 0, 5}] // Flatten\)

The solution for the rotation angle φ(t) (in red) found from above in terms of an interpolating function is shown below together with horizontal oscillation x(t) (in blue) of mass M.

\(φ\_M[t_] = φ[t] /. sol\_9; \ \ \ 
    th\  = \ AbsoluteThickness[2.5]; 

    p\_2 = Plot[{φ\_M[t], x\_2[t]}, {t, 0, 3}, 

        	\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
        PlotStyle -> {{th, Red}, {th, Blue}}, 

        	\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ AxesLabel -> {"t", " "}, 
        
		\ \ \ \ \ \ \ \ \ \ \ \ \ 
        Epilog -> {Text["x(t)", {3, 6}], 
            Text[\*""\!\(φ\_M\)(t)"", {3, 10.5}]}\ ]; \)

The shape of the contour 𝒞  which gives rise to oscillations of the mass M is shown below.

\(\(p\_3 = 
      ParametricPlot[{x\_2[t]\ Cos[φ\_M[t]], 
            \(-x\_2[t]\)\ Sin[φ\_M[t]]} // Evaluate, 

        \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 		{t, 
          0, 3.22}, 

        	\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 	
        PlotStyle -> {{th, Red}}, 

        	\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 	
        PlotPoints -> 100, AspectRatio -> Automatic, 

        		\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
        PlotLabel -> "Shape  of Rotating Contour
", 

        		\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
        AxesLabel -> {"x(t)", "y(t)"}]; \)\)

Animation

Finally, we animate the system under investigation using the solution φ(t) obtained above.

At first we should load the package Arrow and make some definitions for graphics.

\(Needs["Graphics`Arrow`"]; 
HB = "Helvetica-Bold"; 
    
$DefaultFont = {HB, 12}; 
th = AbsoluteThickness[2.5]; 

    th1 = AbsoluteThickness[1. ]; \)

Then we prepare the graphics objects to animate our system.

\(x\_M[t1_, t2_] = \ \ 
      x\_2[t1]\ Cos[φ\_M[t1] - φ\_M[t2]]; 

    y\_M[t1_, t2_] = 
      \(-x\_2[t1]\)\ Sin[φ\_M[t1] - φ\_M[t2]]; \ \ \ \ \ 
    
obj[tmax_, t2_, dt_] := 
\ \ 
      Module[\ {y\_m, dn, q1, q2}, 
        
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
        y\_m = \(-2.2\) - φ\_M[t2]; dn = Floor[tmax/dt]; 
        
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
        q2 = ParametricPlot[{x\_M[t1, t2], y\_M[t1, t2]}, {t1, 0, t2 +  .01}, 
            
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ PlotStyle -> {{th, \ Red}}, 
            
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ AspectRatio -> Automatic, \ 
            Ticks -> None, 

            \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
            PlotRange -> {{\(-8\), 9.5}, {\(-15.7\), 8}}, 	

            			\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
            Prolog -> {\ {Green, Disk[{x\_2[t2], 0}, 0.5]}, 
                
\ 		\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {
                  GrayLevel[0. ], th, 
                  
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
                  Line[{{x\_2[t2] - 0.3, 0}, {x\_2[t2] +  .3, 0}}]}, 

                			\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {
                  GrayLevel[0.3], th, 
                  
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
                  Line[{{0, 0}, {x\_M[0, t2], y\_M[0, t2]}}]}, 
                
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
                Table[{GrayLevel[k\ dt/tmax], 
                    
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
                    Line[{{0, 0}, {x\_M[k\ dt, t2], y\_M[k\ dt, t2]}}]}, 
                  
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {k, dn}\ ]

                						\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ }, 

            			\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
            Epilog -> {\ Arrow[{0, 0}, {9.5, 0}], 
                
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {Blue, 
                  Disk[{0, 0}, 1]}, 
                
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {Magenta, 
                  th, Circle[{0, 0}, 1]}, 

                		\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {Magenta, 
                  Rectangle[{\(-1.3\), y\_m - 1}, {\(-0.7\), y\_m}]}, 

                	\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {Magenta, th, 
                  Line[{{\(-1\), y\_m}, {\(-1\), 0}}]}, 
                
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
                Text[FontForm["x", {HB, 12}], {9, 0.7}]

                			\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ }, 
            
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ DisplayFunction -> Identity\ ]; 
        

\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ q1 = ParametricPlot[{x\_M[t1, t2], y\_M[t1, t2]}, {t1, 0, tmax}, 
            
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
            PlotStyle -> {{th1, \ GrayLevel[0.5]}}, 
            
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ AspectRatio -> Automatic, \ 
            Ticks -> None, 

            \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
            PlotRange -> {{\(-8\), 9.5}, {\(-15.7\), 8}}, 
            
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ DisplayFunction -> Identity]; 
		
        
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
        Show[q2, q1, DisplayFunction -> $DisplayFunction]; 
        
\ \ \ \ \ \ \ \ \ \ \ ]; \)

animateSequence = 
      Table[ obj[3.1, (Subscript[t, 2]), 0.05], {(Subscript[t, 2]), 0, 3.1, 0.05}];