Variational Calculus in Classical Mechanics and the Euler-Lagrange Equation

Abstract:
In this class, we will review the derivation of the Euler-Lagrange equation of Analytical Mechanics through the use of variational calculus techniques, and from this, we will show in detail its application in solving the Brachistochrone problem.

Learning Objectives:
By the end of this class, the student will be able to:

1. Understand Hamilton’s principle of least action
2. Demonstrate the Euler-Lagrange equation
3. Solve the Brachistochrone problem using the Euler-Lagrange equation.

TABLE OF CONTENTS:
WHY VARIATIONAL CALCULUS IN CLASSICAL MECHANICS
FORMULATION OF THE VARIATIONAL PROBLEM
THE EULER-LAGRANGE EQUATION
THE BRACHISTOCHRONE PROBLEM
GITHUB REPOSITORY WITH WOLFRAM ALGORITHM

Why Variational Calculus in Classical Mechanics

Newtonian physics presents numerous problems that can be more effectively addressed using variational calculus. This approach is fundamental in the Lagrange equations and Hamilton’s principle of least action. Essentially, this method involves finding trajectories that maximize or minimize a certain quantity. For example, one can seek the trajectory between two points that minimizes the distance traveled or the travel time. An example of this approach is Fermat’s principle, which states that light always follows the path that minimizes the travel time, which in turn leads to Snell’s law of light refraction.

Variational calculus has multiple advantages in classical mechanics. For instance, it allows obtaining exact analytical solutions for systems with symmetry, and approximate solutions through variational perturbation theory for more complex systems. Moreover, in situations where it is difficult to express the forces in terms of differential equations, the principle of least action provides a more efficient method for solving problems in classical mechanics. In summary, variational calculus is a fundamental tool that offers an alternative formulation of Newton’s laws, a unification of the laws of physics, greater efficiency in problem-solving, and higher precision in predicting experimental results.

Formulation of the Variational Problem

Variational calculus focuses on finding the function y(x) that extremizes the value of the functional:

J(x,y(x))=\displaystyle \int_{x_1}^{x_2} f\left(x,y(x),\frac{dy(x)}{dx}\right)dx,

with the aim of finding its maximum or minimum value. In this equation, the functional J depends on the function y(x) and its derivative dy(x)/dx, while the limits of integration remain fixed. To extremize the integral, variations are applied to the function y(x), seeking to obtain the function that makes the value of the functional an extremum. For example, if the integral reaches its minimum value, any function within its neighborhood, no matter how close it is to y(x), will increase the value of the functional.

To establish the concept of a neighboring function, we can assign a parametric representation y(\alpha,x) to all possible functions y, so that if \alpha=0, then y(0,x)=y(x) is the function that extremizes J. This can be expressed as follows:

y(\alpha, x) = y(x) + \alpha \eta(x),

where \eta(x) is some function of class \mathcal{C}^1 that vanishes at x_1 and x_2, so that the function y(\alpha,x) that includes this variation is identical to y(x) at the initial and final points of the integration path.

By substituting the function y(\alpha,x) that includes the variation \eta(x) instead of y(x) in the integral that defines the functional J, we obtain a new functional that depends on the parameter \alpha:

J(x,y(\alpha, x)) = \displaystyle \int_{x_1}^{x_2} f\left(x,y(\alpha,x), \dfrac{d}{dx}y(\alpha,x)\right)dx

For local extrema to exist, it is necessary that the condition holds:

\left.\dfrac{\partial J(x,y(\alpha,x))}{\partial \alpha}\right|_{\alpha=0} = 0

for any function \eta(x).

The Euler-Lagrange Equation

By analyzing the derivative \partial J(x,y(\alpha,x))/\partial \alpha, we obtain:

\begin{array}{rll} {}\dfrac{\partial J(x,y(\alpha,x))}{\partial \alpha} &=&\dfrac{\partial}{\partial \alpha} \displaystyle \int_{x_1}^{x_2} f\left(x,y(\alpha,x),\dfrac{dy(\alpha, x)}{dx}\right)dx \\ \\ &=&\displaystyle \int_{x_1}^{x_2} \left(\dfrac{\partial f}{\partial x}\dfrac{\partial x}{\partial \alpha} + \dfrac{\partial f}{\partial y(\alpha, x)}\dfrac{\partial y(\alpha, x)}{\partial \alpha} + \dfrac{\partial f }{ \partial \frac{dy(\alpha,x)}{dx}} \dfrac{\partial \frac{dy(\alpha,x)}{dx}}{\partial \alpha}\right)dx \\ \end{array}

From this point, it’s important to note that:

\begin{array}{rll} \dfrac{\partial x}{\partial \alpha} &=& 0 \\ \\ \dfrac{\partial y(\alpha,x)}{\partial \alpha} &=& \dfrac{\partial}{\partial \alpha} \left(y(x) + \alpha \eta(x) \right) = \eta(x) \\ \\ \dfrac{\partial}{\partial \alpha}\left( \dfrac{dy(\alpha, x)}{dx} \right)&=& \dfrac{\partial}{\partial \alpha} \left(\dfrac{dy(x)}{dx} + \alpha\dfrac{d\eta(x)}{dx} \right) = \dfrac{d\eta}{dx} \end{array}

Therefore, the expression reduces as shown below:

\begin{array} {} \dfrac{\partial J(x,y(\alpha,x))}{\partial \alpha} &=& \displaystyle \int_{x_1}^{x_2} \left(\dfrac{\partial f}{\partial y(\alpha,x)}\eta(x) + \dfrac{\partial f}{\partial \frac{dy(\alpha,x)}{dx}} \dfrac{d\eta(x)}{dx} \right)dx \\ \\ &=&\displaystyle \int_{x_1}^{x_2} \dfrac{\partial f}{\partial y(\alpha,x)}\eta(x) dx + \int_{x_1}^{x_2} \dfrac{\partial f}{\partial \frac{dy(\alpha,x)}{dx}} \dfrac{d\eta(x)}{dx} dx \end{array}

Then, if we observe the second integral, we will see that it can be simplified using integration by parts:

\begin{array}{rll} \displaystyle \int_{x_1}^{x_2} \dfrac{\partial f}{\partial \frac{dy(\alpha,x)}{dx}} \dfrac{d\eta}{dx} dx &=& \left. \dfrac{\partial f}{\partial \frac{dy(\alpha,x)}{dx}} \eta(x)\right|_{x_1}^{x_2} - \displaystyle \int_{x_1}^{x_2}\eta(x) \dfrac{d}{dx}\left( \dfrac{\partial f}{\partial \frac{dy(\alpha, x)}{dx}} \right) dx\\ \\ &=& - \displaystyle \int_{x_1}^{x_2}\eta(x) \dfrac{d}{dx}\left( \dfrac{\partial f}{\partial \frac{dy(\alpha, x)}{dx}} \right)dx \end{array}

And therefore

\begin{array}{rll} {} \dfrac{\partial J(x,y(\alpha,x))}{\partial \alpha} &=& \displaystyle \int_{x_1}^{x_2} \left[ \eta(x) \dfrac{\partial f}{\partial y(\alpha, x)} - \eta(x) \dfrac{d}{dx}\left( \dfrac{\partial f}{\partial \frac{dy(\alpha,x)}{dx}} \right) \right]dx \\ \\ &=& \displaystyle \int_{x_1}^{x_2} \left[ \dfrac{\partial f}{\partial y(\alpha, x)} - \dfrac{d}{dx}\left( \dfrac{\partial f}{\partial \frac{dy(\alpha,x)}{dx}} \right) \right] \eta(x) dx \end{array}

So, by the condition that \left.\dfrac{\partial J (x,y(\alpha, x))}{\partial \alpha}\right|_{\alpha=0} = 0, and since \eta(x) is any function subject only to the condition of vanishing at x_1 and x_2, we have:

\dfrac{\partial f}{\partial y(0, x)} - \dfrac{d}{dx}\left( \dfrac{\partial f}{\partial \frac{dy(0,x)}{dx}}\right) = \dfrac{\partial f}{\partial y(x)} - \dfrac{d}{dx}\left( \dfrac{\partial f}{\partial \frac{dy(x)}{dx}}\right) = 0.

Finally, “unloading the notation” in this last expression, we arrive at what is known as the Euler-Lagrange Equation:

\boxed{\dfrac{\partial f}{\partial y}= \dfrac{d}{dx}\left( \dfrac{\partial f}{\partial y^\prime} \right)},

and this represents in a much simpler way the necessary condition for the functional J to reach an extreme value.

The Brachistochrone Problem

Problem Formulation

The brachistochrone problem is a classic in mechanical physics that is solved using variational calculus. The situation presented is as follows: Suppose we have a material object moving under the effect of a constant force field, traveling from an initial point (x_1,y_1) to a final point (x_2,y_2), where the initial point is at a higher elevation than the final point. The question posed is: What is the trajectory the particle should follow to reach the final point in the shortest possible time?

Solution Formulation

To solve the brachistochrone problem, it’s useful to consider the situation in a simple form. Therefore, we can set the starting point (x_1, y_1) at the origin of coordinates, while the endpoint (x_2,y_2) is to the right of the origin and below the \hat{x} axis.

In this situation, we can consider a force field acting downward (in the -\hat{y} direction) generated by gravity, and assume the motion occurs without friction. In this context, the particle is constrained to follow different trajectories connecting the departure and arrival points to find which minimizes the travel time.

Examining the Energy

To solve this problem, we can utilize the conservation of energy in the gravitational system. The total energy of the system will remain constant, considering both the kinetic energy E_{kin}=\frac{1}{2}mv^2 and the gravitational potential energy E_{pot,g}, where m is the particle’s mass and v is its velocity. For the potential energy, the origin is taken as the reference, so E_{pot,g}(y=0)=0, while at any other height y, E_{pot,g}(y)=mgy.

Since the particle starts from rest at the origin, its total energy equals zero. Therefore, we have:

E_{kin} + E_{pot,g}=0

As the particle falls below the reference point, its potential energy will be negative, and its kinetic energy will be positive. Thus, we can solve for the velocity v using the energy conservation equation to obtain:

\begin{array}{rl} {} &\dfrac{1}{2}mv^2 + (-mgy) = 0 \\ \\ \vdash &\dfrac{1}{2}mv^2 = mgy \\ \\ \vdash &v^2 = 2gy \\ \\ \vdash &v = \sqrt{2gy} \end{array}

In this way, we can calculate the particle’s velocity at any point along its trajectory as a function of the height y at that point.

Examining the Travel Time

Once we have obtained the speed of movement, we can construct the travel time element using the displacement element ds=\sqrt{dx^2 + dy^2} as follows:

\begin{array}{rl} {} dt &= \dfrac{ds}{v} = \dfrac{\sqrt{dx^2 + dy^2}}{\sqrt{2gy}}\\ \\ &= \sqrt{\dfrac{dx^2 + dy^2}{2gy} } \end{array}

Therefore, the travel time between points (x_1,y_1) and (x_2,y_2) can be obtained by integrating

\begin{array}{rl} {} t &= \displaystyle \int_{(x_1,y_1)}^{(x_2,y_2)} dt \\ \\ &= \displaystyle \int_{(x_1,y_1)}^{(x_2,y_2)} \sqrt{\dfrac{dx^2 + dy^2}{2gy}} \\ \\ &= \displaystyle \dfrac{1}{\sqrt{2g}}\int_{y_1}^{y_2} \sqrt{\dfrac{1+ \left(\dfrac{dx}{dy}\right)^2 }{y}}dy \\ \\ \end{array}

Variational Problem Formulation

With this last expression, we have managed to express the time as a functional of the form

{}t = J(y,x(y)) = \displaystyle \int_{y_1}^{y_2} f\left(y,x(y),\dfrac{dx(y)}{dy} \right) dy

where

f\left(y,x(y), \dfrac{dx(y)}{dx}\right) = \sqrt{\dfrac{1+ \left(\dfrac{dx(y)}{dy} \right)^2}{y}}

At this point, we can overlook the factor \sqrt{2g}, because optimizing J is the same as optimizing \sqrt{2g}J.

With the above, we can now construct the Euler-Lagrange equation following the same procedure used earlier, finally arriving at:

\dfrac{\partial f}{\partial x} = \dfrac{d}{dy} \dfrac{\partial f}{\partial x^\prime}

However, here we can see that \dfrac{\partial f}{\partial x} = 0, so it will be

\dfrac{d}{dy}\dfrac{\partial f}{\partial x^\prime} = 0,

or in other words

\dfrac{\partial f}{\partial x^\prime} = \dfrac{1}{\sqrt{2a}},

where a is an arbitrary constant written that way because it is “convenient” for later developments.

Resolution of the Variational Problem

Substituting the function f into this last expression, we have:

\begin{array}{rl} {} &\dfrac{\partial }{\partial x^\prime} \sqrt{\dfrac{1+ x^{\prime 2}}{y}} = \dfrac{1}{\sqrt{2a}} \\ \\ \vdash & \dfrac{1}{2}\left( \dfrac{1 + x^{\prime 2} }{y} \right)^{-1/2} \left(\dfrac{2x^\prime}{y} \right) = \dfrac{1}{\sqrt{2a}} \\ \\ \vdash & \dfrac{1}{2}\sqrt{\dfrac{y}{1 + x^{\prime 2}}} \left(\dfrac{2x^\prime}{y} \right) = \dfrac{1}{\sqrt{2a}} \\ \\ \vdash & \sqrt{\dfrac{4x^{\prime 2} y}{4y^2 (1 + x^{\prime 2})} } = \sqrt{\dfrac{1}{2a}} \\ \\ \vdash & \dfrac{y x^{\prime 2} }{y^2 (1 + x^{\prime 2})} = \dfrac{1}{ 2a} \\ \\ \vdash & 2ayx^{\prime 2} = y^2 + y^2 x^{\prime 2} \\ \\ \vdash & x^{\prime 2} (2ay - y^2) = y^2 \\ \\ \vdash & \left(\dfrac{dx}{dy}\right)^2 = \dfrac{y^2}{2ay - y^2} \\ \\ \vdash & \dfrac{dx}{dy} = \pm \sqrt{\dfrac{y^2}{2ay - y^2}} \\ \\ \vdash & dx = \pm \dfrac{ydy}{\sqrt{2ay - y^2}} \\ \\ \vdash & x = \displaystyle \pm \int \dfrac{y}{\sqrt{2ay - y^2}}dy \end{array}

To solve this integral, one option to consider is performing the following substitution

\begin{array} {} y &=& a[1-\cos(\theta)] \\ dy &=& a\sin(\theta) d\theta \end{array}

With this we have:

\begin{array}{rl} {} x= & \pm \displaystyle \int \dfrac{y}{\sqrt{2ay - y^2}}dy = \displaystyle \int \dfrac{a[1-\cos(\theta)]a\sin(\theta)}{\sqrt{2a^2[1-\cos(\theta)] - a^2[1-\cos(\theta)]^2 }}d\theta \\ \\ & {} = \pm \displaystyle \int \dfrac{a^2[1-\cos(\theta)]\sin(\theta)}{\sqrt{a^2[1-\cos(\theta)]\left\{ 2 - [1-\cos(\theta)] \right\} }}d\theta \\ \\ & {} = \pm \displaystyle \int \dfrac{a[1-\cos(\theta)]\sin(\theta)}{\sqrt{[1-\cos(\theta)] [1 + \cos(\theta)] }}d\theta \\ \\ & {} = \pm \displaystyle \int \dfrac{a[1-\cos(\theta)]\sin(\theta)}{\sqrt{ 1-\cos^2(\theta)}}d\theta \\ \\ & {} = \pm \displaystyle \int \dfrac{a[1-\cos(\theta)]\sin(\theta)}{\sin(\theta)}d\theta \\ \\ & {} = \pm \displaystyle \int a[1-\cos(\theta)] d\theta \\ \\ & {} = \pm a(\theta - \sin(\theta)) + C \end{array}

We can see that the brachistochrone curve can be expressed as a parametric curve in polar coordinates, which coincides with a cycloid that starts at the origin.

\begin{array} {} & x(\theta) &=& \pm a(\theta - \sin(\theta)) \\ & y(\theta) &=& a(1-\cos(\theta)) \end{array}

The integration constant C has been nullified to satisfy the initial condition that the trajectory starts at the origin. Moreover, we can observe that there are a couple of equations that provide possible solutions to the problem, where the constant a can be adjusted so that the curve passes through the point (x_2,y_2) at the end of the path. These equations are:

Option 1:\boxed{\begin{array} {} & x(\theta) &=& a(\theta - \sin(\theta)) \\ & y(\theta) &=& a(1-\cos(\theta)) \end{array}}

Option 2:\boxed{\begin{array} {} & x(\theta) &=& - a(\theta - \sin(\theta)) \\ & y(\theta) &=& a(1-\cos(\theta)) \end{array}}

The viable solution for this problem is given by the second option, and by adjusting the constant a as a negative value, we obtain a curve that meets the necessary conditions to be the solution.

Final Adjustment of the Solution

After the latest adjustments, the brachistochrone curve has the parametric form:

\begin{array} {} x(\theta) &= b(\theta - \sin(\theta)) \\ y(\theta) &= -b(1-\cos(\theta)) \end{array}

It was substituted a=-b, where 0\lt b. The curve has a period 2b\pi and must satisfy the condition x_2 \in ]0,2b\pi[ and y_2 \in ]-2b,0[. This last part is crucial, as it requires that the brachistochrone curve be represented as a single arc of a cycloid, since the solution will cease to be valid if the particle returns to rest when returning to a point of zero height.

To adjust these equations to the problem, we need to find the values of \theta and b that satisfy the system:

\begin{array} {} x_2 &= b(\theta - \sin(\theta))\\ y_2 &= - b(1-\cos(\theta)) \end{array}

This nonlinear system does not seem to have analytical solutions, so we will use numerical methods in Wolfram Mathematica. Below is a series of steps to solve the problem:

Step 1: Set up the system

Set the equations that form the system to be solved

eq1 = x2 == b*(theta - Sin[theta])
eq2 = y2 == -b*(1 - Cos[theta])

Step 2: Define the arrival point

Set the point the particle will reach at the end of its path. In this case, we will set it at (x_2,y_2)=(1,-2). You can modify these values to test other similar configurations.

x2val = 1; y2val = -2;

Step 3: Calculate the sought values numerically

Use the “FindRoot” function to numerically compute the solution to the problem

sol = FindRoot[{eq1, eq2} /. {x2 -> x2val, y2 -> y2val}, {{b,1}, {theta, 1}}]

Here, the values b=1 and \theta=1 have been used as a starting point for the numerical approximation of the solution. This gives us the solution b\approx 2.4056 and \theta \approx 1.40138

Step 4: Verification of Results

Remember that for these answers to make physical sense, it is necessary that x_2 \in ]0,2b\pi[ and y_2 \in ]-2b, 0[. We can quickly check if this is the case using the following procedure

First, extract the values of b and \theta obtained as the solution

bval = sol[[1, 2]]; thetaval = sol[[2, 2]];

Then, command the confirmation

If[0 < x2val < 2*Pi*bval && -2*bval < y2val < 0 "Valid Values", "Invalid Values"]

If everything went well, we should get "Valid Values" in the output. This piece of code will help you check if the physical situation is correctly modeled.

With these procedures, we finally have our completely adjusted solution curve that connects the points (x_1,y_1)=(0,0) and (x_2,y_2)=(1,-2). The resulting curve is:

\begin{array} {} x(\theta) &\approx 2.4056(\theta - \sin(\theta)) \\ y(\theta) &\approx -2.4056(1-\cos(\theta)) \end{array}\;\;;\theta\in [0, 1.40138]

Which graphically looks like this:

GitHub Repository with Wolfram Algorithm

The complete code for the solution to the brachistochrone problem, including the algorithm developed in Wolfram Mathematica, is available for download and reference in my GitHub repository. This repository includes a `.nb` file with the code in an interactive notebook format, as well as a plain text version `.m` for those who prefer to view the code directly.

You can download the repository from GitHub here.

In addition to the code, the repository contains a "README" file with detailed instructions on how to use and understand the algorithm, as well as a step-by-step explanation of the solution to the brachistochrone problem. I hope you find it useful!