Lagrangian Equations

Ales Bardorfer, student of the Faculty of Electrical Engineering, Ljubljana, Slovenia,  May 1998.

Improvement and extension of an item written by Jeff Adams, Wolfram Research Inc.,  October 1992.

Introduction

This notebook illustrates one way of using Mathematica to analyze some advanced theoretical mechanics problems using Lagrangian Mechanics. The example problems are taken from Analytical Mechanics, Grant R. Fowles,  Fourth Edition, Saunders College Publishing, 1986 and Schaum's outline of theory and problems of Lagrangian Dynamics, Dare A. Wells, McGraw-Hill Book Company, 1967.

Functions

The function LagrangianEquations[T,  V,  Q,  Coords] takes four arguments:  the kinetic energy of the system, T, the potential energy of the system, V,  a list of non-conservative forces of the system, Q, and a list of the generalized coordinates of the system, Coords.
The function returns a list of the Lagrangian equations for the particular system.

Note that, as in the development in classical mechanics, it is very important to explicitly declare all dependencies on the generalized coordinates of T, V, and Q, as well as the explicit time dependence of the generalized coordinates themselves. The arguments, Q and Coords, must each be in the form of a list, even when the function is used for only one generalized coordinate. These lists must have the same length so that the non-conservative forces map in a one-to-one fashion to the generalized coordinates.

Recall the definition of MapThread:

The definition of MapThread

Function that generates the Lagrangian Equations of motion:

\(LagrangianEquations::"usage" = 
      "LagrangianEquations[KineticEnergy, PotentialEnergy, {NonConservativeForces}, {GeneralizedCoordinates}] returns a list of the Lagrangian equations for the dynamical system. The alternative form, LagrangianEquations[Lagrangian, {NonConservativeForces}, {GeneralizedCoordinates}], can also be used when the Lagrangian itself is supplied."; 


    LagrangianEquations[T_, V_, Q_List, Coords_List] := 
	
      LagrangianEquations[T - V, Q, Coords]; 


    LagrangianEquations[L_, Q_List, Coords_List] := 
      MapThread[
        ∂\_t\((∂\_\(∂\_t #2\)L)\) - 
              ∂\_#2 L == #1&, {Q, Coords}]; \)

The function KineticEnergy[m,  Coords] returns the kinetic energy of a particle. It takes two arguments: the mass of the particle, m, and a list of coordinates defined by transform equations determining the possible path of a massive particle.

Function that generates the kinetic energy of a particle:

\(KineticEnergy::"usage" = 
      "KineticEnergy[Mass, {Coords}] returns the kinetic energy of a particle of mass "Mass" whose motion depends on the coordinates Coords."; 
    

KineticEnergy[m_, Coords_List] := 
      1\/2\ m\ Plus@@\(\((∂\_t #1&)\)/@Coords\^2\); \)

Example 1 :  Harmonic Oscillator

Coords = {x[t]};
T=(1/2) m ((x'[t])^2);
V = (1/2) k x[t]^2 - m g x[t];
Q = {0};

equations = LagrangianEquations[T,  V,  Q,  Coords]

DSolve[equations,Coords,t]

Example 2 :  Harmonic Oscillator with Velocity Damping

Coords = {x[t]};
T=(1/2) m ((x'[t])^2);
V = (1/2) k x[t]^2;
Q={-f x'[t]};

equations = LagrangianEquations[T,  V,  Q,  Coords]

DSolve[equations, Coords, t]

Example 3 :  Single Particle in a Central Field

Using polar coordinates

Coords = {r[t], θ[t]};
T=(1/2) m (((r'[t])^2)+r[t]^2 ((θ'[t])^2));
V = v[r[t]];
Q = {0, 0};

TableForm[equations = LagrangianEquations[T,  V,  Q,  Coords]]

Example 4 :  Atwood's Machine

A mechanical system known as Atwood's machine consists of two weights of mass m1 and m2 , respectively, connected by a light inextensible cord of length l which passes over a pulley.  The system has one degree of freedom.  Let x represent the configuration of the system, where x is the vertical distance from the pulley to m1.  I is the moment of inertia of the pulley with radius a.

Coords = {x[t]};
T=(1/2) m1 ((x'[t])^2)+(1/2) m2 ((x'[t])^2)+((l ((x'[t])^2))/(2 (a^2)));
V = -m1 q x[t] - m2 q (l - x[t]);
Q = {0};

equations = LagrangianEquations[T,  V,  Q,  Coords]

DSolve[equations, Coords, t]

Example 5 :  Double Atwood  Machine

Replace one of the weights in the Atwood machine by another pulley supporting two weights connected by another cord.  The system now has two degrees of freedom.  We shall specify the coordinates x1 and x2.  For simplicity we shall neglect the masses of the pulleys in this case.

Coords = {x1[t], x2[t]};
T=(1/2) m1 ((x1'[t])^2)+(1/2) m2 (((-x1'[t]+x2'[t]))^2)+(1/2) m3 (((-x1'[t]-x2'[t]))^2);
V = -m1 q x1[t] - m2 q (l1 - x1[t] + x2[t]) - 
          m3 g (l1 - x1[t] + l2 - x2[t]);
Q = {0, 0};

TableForm[equations = LagrangianEquations[T,  V,  Q,  Coords] ]

\(TableForm[
	
      \(Simplify[DSolve[equations, Coords, t]]\)[[1
        ]]
]\)

Example 6 :  Particle Sliding on a Movable Inclined Plane with Velocity Damping

Consider a particle sliding on a smooth inclined plane which itself is free to slide on a smooth horizontal surface. There are two degrees of freedom x1 and x2, the horizontal displacement of the plane from some reference point and the displacement of the particle from some reference point on the plane.  Using the law of cosines for the velocity of the particle, we have:

Coords = {x1[t], x2[t]};
T=(1/2) m (((x1'[t])^2)+((x2'[t])^2)+2 x1'[t] x2'[t] Cos[θ])+(1/2) M ((x1'[t])^2);
V = -m g x2[t] Sin[θ] + κ;
Q={f1 x1'[t],f2 x2'[t]};

\(TableForm[
	equations = LagrangianEquations[T, \ V, \ Q, \ Coords]\ 

      	]\)

TableForm[
	Simplify[Solve[equations,{(x1')'[t],(x2')'[t]}]][[1]]
]

Example 7 :  Spherical Pendulum, or Bar of Soap in a Bowl

Coords = {θ[t], ϕ[t]};
T=(1/2) m (l^2) (((θ'[t])^2)+Sin[θ[t]]^2 ((ϕ'[t])^2));
V = m g l (1 - Cos[θ[t]]);
Q = {0, 0};

\(TableForm[
	equations = LagrangianEquations[T, \ V, \ Q, \ Coords]\ 

      	]\)

Given that ϕ is an ignorable coordinate then m (l^2) Sin[θ[t]]^2 ϕ'[t]=constant.

TableForm[
	equations1=equations/.ϕ'[t]->(h1/(m (l^2) Sin[θ[t]]^2))
]

equations2=Simplify[Solve[equations1,(θ')'[t]]]

θEqn  =   equations2  /.  {h1  ->  h (l^2) m}

In one special case, if ϕ  =  constant then ϕ' is zero and so h = 0 and we recover the simple pendulum case.

Example 8 :  Pendulum Attached to a Cart with Coulomb and Viscous Friction

Consider a pendulum with a mass particle m, attached to a cart with mass M, which can be moved along the X axis on a rough surface (Coulomb friction  -  kf). There are two degrees of freedom:  the cart displacement x1[t], and  the pendulum angle ϕ[t]. Point (x1,  y1) represents the point, at which pendulum is connected to the cart (pendulum rotation center), and point (x2,  y2) represents  the mass particle location.

Clear[x1, y1,  x2,  y2,  kf,  f2,  ϕ,  m,  M,  g,  z];

Transform equations:

y1 = 0;
x2 = x1[t] + l Sin[ϕ[t]];
y2 = y1 + l Cos[ϕ[t]];
Coords = {x1[t], ϕ[t]};
T = Simplify[
          KineticEnergy[M, {x1[t], y1}] + KineticEnergy[m, {x2, y2}]];
V = m g y2;
Q={-kf Sign[x1'[t]],-f2 ϕ'[t]};

\(\(\(TableForm[
	
        equations = Simplify[LagrangianEquations[T, \ V, \ Q, \ Coords]]
	]
        \ \)
\)\)

m = 1;
M = 5;
g = 9.81;
l = 1;
kf = 1;
f2 = 0.3;
tend = 10;

NumSolution=NDSolve[Flatten[{equations,x1[0] == 0,x1'[0]==0,ϕ[0] == 0.5,ϕ'[0]==0}],Coords,{t, 0, tend}];

\(pl1 = 
      Plot[Evaluate[x1[t] /. NumSolution], {t, 0, tend}, 
        FrameLabel -> {t, x1[t]}, DisplayFunction -> Identity, 
        Axes -> False, Frame -> True]; 


    pl2 = Plot[Evaluate[ϕ[t] /. NumSolution], {t, 0, tend}, 
        FrameLabel -> {t, ϕ[t]}, DisplayFunction -> Identity, 
        Axes -> False, Frame -> True]; 


    Show[GraphicsArray[{pl1, pl2}]]; \)

Example 9 :  Two Connected Pendulums with  Velocity Damping

Consider two pendulums attached to each other in series with particles of mass  m1 and m2.  There are two degrees of freedom:  the pendulum angle ϕ[t]  and θ[t]. Point  (x1,  y1) represents the location of the first mass particle and point (x2,  y2) represents the location of  the second massive particle.

Clear[x1, x2, y1, y2, z, r1, r2, ϕ, θ, m1, m2, f1, f2, 
      g];

Transform equations:

x1 = r1 Sin[ϕ[t]];
y1 = r1 Cos[ϕ[t]];
x2 = x1 + r2 Sin[ϕ[t] + θ[t]];
y2 = y1 + r2 Cos[ϕ[t] + θ[t]];
Coords = {ϕ[t], θ[t]};T = Simplify[
          KineticEnergy[m1, {x1, y1}] + KineticEnergy[m2, {x2, y2}]];
V = -m1 g y1 - m2 g y2;
Q={-f1 ϕ'[t],-f2 θ'[t]};

\(\(equations = Simplify[LagrangianEquations[T, \ V, \ Q, \ Coords]]

    \)\)

m1 = 1;
m2 = 2;
r1 = 1;
r2 = 1;
f1 = 1;
f2 = 0.5;
tend = 10;
g = 9.81;

NumSolution=NDSolve[Flatten[{equations,ϕ[0] == 2,ϕ'[0]==0,θ[0] == 1,θ'[0]==0}],Coords,{t, 0, tend}];

\(pl1 = 
      Plot[Evaluate[ϕ[t] /. NumSolution], {t, 0, tend}, 
        FrameLabel -> {t, ϕ[t]}, DisplayFunction -> Identity, 
        Axes -> False, Frame -> True, PlotRange -> All]; 


    pl2 = Plot[Evaluate[θ[t] /. NumSolution], {t, 0, tend}, 
        FrameLabel -> {t, θ[t]}, 
        DisplayFunction -> Identity, Axes -> False, 
        Frame -> True, PlotRange -> All]; 


    Show[GraphicsArray[{pl1, pl2}]]; \)

Example 10 :  Pendulum on a Spring with Velocity Damping

Consider a mass particle on a spring, forming a pendulum with two degrees of freedom:  the pendulum angle, θ[t] and the stretch of a spring, r[t]. Point (x1,  y1) represents the location of the mass particle m.

\(\(
\(Clear[x, y, z, r, θ, m, g, r0, f1, f2, k]; \)\)\)

Transform equations:

x = r[t] Sin[θ[t]];
y = -r[t] Cos[θ[t]];
Coords = {r[t], θ[t]};
T = KineticEnergy[m, {x, y}];V = m g (r0 + y) + (1/2) k (((r[t] - r0))^2);
Q={-f1 r'[t],-f2 θ'[t]};

\(\(\(TableForm[
	
        equations = Simplify[LagrangianEquations[T, \ V, \ Q, \ Coords]]
	]
        
\)\ \)\)

m = 1;
r0 = 1;
g = 9.81;
k = 5;
f1 = 0.1;
f2 = 0.2;
tend = 10;

NumSolution=NDSolve[Flatten[{equations,r[0] == r0,r'[0]==0,θ[0] == 1,θ'[0]==0}],Coords,{t, 0, tend}];

\(pl1 = 
      Plot[Evaluate[r[t] /. NumSolution], {t, 0, tend}, 
        FrameLabel -> {t, r[t]}, DisplayFunction -> Identity, 
        Axes -> False, Frame -> True, PlotRange -> All]; 


    pl2 = Plot[Evaluate[θ[t] /. NumSolution], {t, 0, tend}, 
        FrameLabel -> {t, θ[t]}, 
        DisplayFunction -> Identity, Axes -> False, 
        Frame -> True, PlotRange -> All]; 


    Show[GraphicsArray[{pl1, pl2}]]; \)

Example 11 :  Car Moving on an Inclined Road Through the Air

Consider a car of mass m with an engine, travelling down the slope through the air on an inclined road (α). Assume that the air resistance is proportional to the square of the car's speed. There is only one degree of freedom:  the car displacement, s ,along the road. Point (x,  z) represents the location of the car.

This example shows that not only frictional forces but also external (actuator) forces and torques etc. are non-conservative.

\(\(
\(Clear[x, y, z, α, m, g, f, EngineForce]; \)\)\)

Transform equations:

x = s[t] Cos[α];
z = -s[t] Sin[α];
Coords = {s[t]};
T = KineticEnergy[m, {x, z}];
V = m g z;
Q={EngineForce-f ((s'[t])^2)};

equations =  Simplify[LagrangianEquations[T,  V,  Q,  Coords]]

m = 1000;
g = 9.81;
f = 10;
α = 5 Degree;
EngineForce = 1500;
tend = 25;

NumSolution=NDSolve[Flatten[{equations,s[0] == 0,s'[0]==0}],Coords,{t, 0, tend}];

Plot[Evaluate[D[s[t]  /.  NumSolution, t]], {t, 0, tend},  
      FrameLabel  ->  {t, "s'[t] [m/s]"}, Axes -> False, Frame -> True, 
      PlotRange -> All];

D[s[t]  /.  NumSolution, t]   /.  t -> 25

It is quite obvious that the car practicly stops accelerating after 25 s, travelling at constant speed of approx. 15.3 m/s. This is due to  the air resistance becoming equal to the sum of the gravitational force and engine force.

Example 12 :  Two Link Planar Robot Manipulator in Gravitational Field with a Position Controller

Consider two links of mass m1, m2, length l1 and l2 attached in series by revolute joints, forming a simple 2 DOF planar manipulator. These two degrees of freedom are: joint 1 angle, phi, measured from a direction of the gravitational field and joint 2 angle θ, measured from the first link. Point (x1,y1) represents the end of link1 and (x2, y2) the end of link 2.  In each joint there is an actuator (e.g. a servo motor) providing an external torque. Assume the manipulator is in a gravitational field (along the y axis: (Overscript[g, ⇀]) = (0,  g,  0)). 
Let's make a simple move from the rest position (all joint angles set to zero) to some reference position.
This example shows how controller output signals, which are non-conservative, can be applyed to a list of non-conservative forces, forming a more sophisticated, controller system.

\(\(
\(Clear[x1, x2, y1, y2, xc1, xc2, yc1, yc2, m1, m2, l1, l2, g, 
      ϕ, θ, \ ϕref, \ θref, \ Kp1, \ Kp2, \ Kd1, \ 
      Kd2, NumSolution1, NumSolution2]; \)\)\)

Transform equations:

\(x1 = \(-l1\)\ Sin[ϕ[t]]; 
y1 = \(-l1\)\ Cos[ϕ[t]]; 

    x2 = x1 - l2\ Sin[ϕ[t] + θ[t]]; 

    y2 = y1 - l2\ Cos[ϕ[t] + θ[t]]; \)

Center of mass coordinates for each link:

\(xc1 = x1\/2; 
yc1 = y1\/2; 
xc2 = \(x1 + x2\)\/2; 

    yc2 = \(y1 + y2\)\/2; \)

The kinetic energy for each link consists of a translational movement of the center of mass point and rotating around the same point:

w1=ϕ'[t];
w2=ϕ'[t]+θ'[t];
J1 = ((m1 (l1^2))/12);   
J2 = ((m2 (l2^2))/12);
T = Simplify[
          KineticEnergy[m1, {xc1, yc1}]  +  J1 ((w1^2)/2)  +  
            KineticEnergy[m2, {xc2, yc2}]  +  J2 (( (w2^2))/2)];
V = m1 yc1 g  +  m2 yc2 g;

Coords = {ϕ[t], θ[t]};

A simple PD controller:

Q = {Kp1 (ϕref - ϕ[t])  -  Kd1 ϕ'[t], 
          Kp2 (θref - θ[t])  -  
          Kd2 θ'[t]};

equations = Simplify[LagrangianEquations[T,  V,  Q,  Coords]]

m1 = 1;
m2 = 0.5;
g = 9.81;
l1 = 1;
l2 = 0.5;
tend = 5;
ϕref = 0.5;
θref = -1;
Kp1 = 10;
Kp2 = 10;
Kd1 = 1;
Kd2 = 0.3;

NumSolution1=NDSolve[Flatten[{equations,ϕ[0] == 0,ϕ'[0]==0,θ[0] == 0,θ'[0]==0}],Coords,{t, 0, tend}];

\(pl1 = 
      Plot[Evaluate[ϕ[t] /. NumSolution1], {t, 0, tend}, 
        FrameLabel -> {t, ϕ[t]}, DisplayFunction -> Identity, 
        Axes -> False, Frame -> True, PlotRange -> All]; 


    pl2 = Plot[Evaluate[θ[t] /. NumSolution1], {t, 0, tend}, 
        FrameLabel -> {t, θ[t]}, 
        DisplayFunction -> Identity, Axes -> False, 
        Frame -> True, PlotRange -> All]; 


    Show[GraphicsArray[{pl1, pl2}]]; \)

In order to avoid colisions, of a robot, the damping coefficient for each axis has to be greater than one.  We can see that both links oscilate over the reference angle; therefore both damping coefficients are too small. The static error is quite big too (about -0.2 rad for joint 1 and 0.1 rad for joint 2). By applying some adjustments to a controller parameters, a better response can be obtained:

Kp1 = 1000;
Kp2 = 100;
Kd1 = 55;
Kd2 = 3;
tend = 1;

NumSolution2=NDSolve[Flatten[{equations,ϕ[0] == 0,ϕ'[0]==0,θ[0] == 0,θ'[0]==0}],Coords,{t, 0, tend}];

\(pl1 = 
      Plot[Evaluate[ϕ[t] /. NumSolution2], {t, 0, tend}, 
        FrameLabel -> {t, ϕ[t]}, DisplayFunction -> Identity, 
        Axes -> False, Frame -> True, PlotRange -> All]; 


    pl2 = Plot[Evaluate[θ[t] /. NumSolution2], {t, 0, tend}, 
        FrameLabel -> {t, θ[t]}, 
        DisplayFunction -> Identity, Axes -> False, 
        Frame -> True, PlotRange -> All]; 


    Show[GraphicsArray[{pl1, pl2}]]; \)

Now the robot makes a move near its reference position without oscilating over the reference position in about 0.3s, which is approx. 20 times faster from the example above. The static eroor is reduced to:

\(TableForm[

      	{e1 = \((\(ϕ[t] /. NumSolution2\)\  /. t -> 1)\)\  - ϕref
            \ , e2 = 
          \((\(θ[t] /. \ NumSolution2\) /. t -> 1)\)\ \  - θref}
        
]\)

Theoretically the static error and response time could be minimized by applying bigger proportional and derivative gains, but this is not possible due to the limitations of actuator torques.