Вариационное исчисление в классической механике и уравнение Эйлера-Лагранжа

Резюме:
В этом уроке мы рассмотрим вывод уравнения Эйлера-Лагранжа из аналитической механики с использованием методов вариационного исчисления, и на основе этого подробно покажем его применение для решения проблемы брахистохроны.

Цели обучения:
По завершении этого урока студент сможет:

1. Понять принцип наименьшего действия Гамильтона
2. Продемонстрировать уравнение Эйлера-Лагранжа
3. Решить проблему брахистохроны, используя уравнение Эйлера-Лагранжа.

ОГЛАВЛЕНИЕ:
ЗАЧЕМ НУЖНО ВАРИАЦИОННОЕ ИСЧИСЛЕНИЕ В КЛАССИЧЕСКОЙ МЕХАНИКЕ
ФОРМУЛИРОВКА ВАРИАЦИОННОЙ ЗАДАЧИ
УРАНЕНИЕ ЭЙЛЕРА-ЛАГРАНЖА
ПРОБЛЕМА БРАХИСТОХРОНЫ
РЕПОЗИТОРИЙ GITHUB С АЛГОРИТМОМ WOLFRAM

Зачем нужно вариационное исчисление в классической механике

Ньютоновская физика представляет множество задач, которые могут быть более эффективно решены с использованием вариационного исчисления. Этот подход является фундаментальным для уравнений Лагранжа и принципа наименьшего действия Гамильтона. По сути, этот метод заключается в нахождении траекторий, которые максимизируют или минимизируют определенную величину. Например, можно найти траекторию между двумя точками, которая минимизирует пройденное расстояние или время путешествия. Примером такого подхода является принцип Ферма, который гласит, что свет всегда следует траектории, минимизирующей время прохождения, что, в свою очередь, приводит к закону преломления света Снеллиуса.

Вариационное исчисление имеет множество преимуществ в классической механике. Например, оно позволяет получить точные аналитические решения для систем с симметрией и приближенные решения с помощью теории вариационных возмущений для более сложных систем. Кроме того, в ситуациях, когда сложно выразить силы в терминах дифференциальных уравнений, принцип наименьшего действия предоставляет более эффективный метод для решения задач классической механики. В общем, вариационное исчисление — это фундаментальный инструмент, который предлагает альтернативную формулировку законов Ньютона, унификацию законов физики, повышенную эффективность в решении задач и большую точность в предсказании экспериментальных результатов.

Формулировка вариационной задачи

Вариационное исчисление направлено на поиск функции y(x), которая экстремизирует значение функционала:

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

с целью нахождения его максимума или минимума. В этом уравнении функционал J зависит от функции y(x) и ее производной dy(x)/dx, в то время как пределы интегрирования остаются фиксированными. Чтобы экстремизировать интеграл, на функцию y(x) накладываются вариации, с целью найти функцию, которая делает значение функционала экстремальным. Например, если удается добиться минимального значения интеграла, любая функция в его окрестности, независимо от того, насколько она близка к y(x), увеличит значение функционала.

Для введения понятия соседней функции мы можем присвоить параметрическое представление y=(\alpha,x) всем возможным функциям y, так что если \alpha=0, то y(0,x)=y(x) является функцией, которая экстремизирует J. Это можно выразить следующим образом:

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

где \eta(x) — некоторая функция класса \mathcal{C}^1, которая равна нулю в x_1 и x_2, так что функция y(\alpha,x), включающая это вариацию, идентична y(x) в начальных и конечных точках траектории интегрирования.

Подставив функцию y(\alpha,x), включающую вариацию \eta(x) вместо y(x) в интеграл, определяющий функционал J, мы получаем новый функционал, зависящий от параметра \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

Для существования локальных экстремумов необходимо выполнение условия:

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

для любой функции \eta(x).

Уравнение Эйлера-Лагранжа

Анализируя производную \partial J(x,y(\alpha,x))/\partial \alpha, получаем:

\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}

С этого момента важно заметить, что:

\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}

Таким образом, выражение сокращается, как показано ниже:

\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}

Если затем рассмотреть второй интеграл, мы увидим, что его можно упростить, используя интегрирование по частям:

\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}

Таким образом

\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}

Таким образом, по условию, что \left.\dfrac{\partial J (x,y(\alpha, x))}{\partial \alpha}\right|_{\alpha=0} = 0, и так как \eta(x) — это произвольная функция, подчиненная единственному условию быть равной нулю в x_1 и x_2,, имеем:

\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.

Наконец, «сбросив обозначения» в этом последнем выражении, мы приходим к тому, что известно как уравнение Эйлера-Лагранжа:

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

и это представляет гораздо более простым способом необходимое условие для того, чтобы функционал J достигал экстремального значения.

Задача о брахистохроне

Постановка задачи

Задача о брахистохроне — это классическая задача механической физики, решаемая с помощью вариационного исчисления. Рассмотрим следующую ситуацию: предположим, что у нас есть материальный объект, который движется под действием постоянного силового поля и перемещается из начальной точки (x_1,y_1) в конечную точку (x_2,y_2), где начальная точка находится выше конечной. Вопрос заключается в следующем: какой траекторией должна следовать частица, чтобы достичь конечной точки за минимальное время?

Постановка решения

Для решения задачи о брахистохроне полезно рассмотреть ситуацию в упрощенной форме. Таким образом, можно установить начальную точку (x_1, y_1) в начале координат, а конечная точка (x_2,y_2) находится справа от начала и ниже оси \hat{x}.

В этой ситуации можно рассмотреть силовое поле, действующее вниз (в направлении -\hat{y}), создаваемое гравитацией, и предположить, что движение происходит без трения. В этом контексте частица ограничена различными траекториями, соединяющими точки выхода и прибытия, с целью найти ту, которая минимизирует время путешествия.

Исследование энергии

Для решения этой задачи мы можем воспользоваться законом сохранения энергии гравитационной системы. Полная энергия системы будет оставаться постоянной, учитывая как кинетическую энергию E_{cin}=\frac{1}{2}mv^2, так и потенциальную гравитационную энергию E_{pot,g}, где m — масса частицы, а v — ее скорость. Для потенциальной энергии был взят за основу источник, таким образом, что E_{pot,g}(y=0)=0, тогда как на любой другой высоте y имеем E_{pot,g}(y)=mgy.

Поскольку частица начинает движение из покоя, ее полная энергия равна нулю. Таким образом, имеем:

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

Так как частица падает ниже точки отсчета, ее потенциальная энергия будет отрицательной, а кинетическая — положительной. Таким образом, мы можем выразить скорость v из уравнения сохранения энергии и получить:

\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}

Таким образом, мы можем вычислить скорость частицы в любой точке ее траектории в зависимости от высоты y, на которой она находится.

Исследование времени прохождения

Как только мы получили скорость движения, мы можем построить элемент времени прохождения, используя элемент смещения ds=\sqrt{dx^2 + dy^2} следующим образом:

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

Таким образом, время прохождения между точками (x_1,y_1) и (x_2,y_2) можно получить, интегрируя

\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}

Постановка вариационной задачи

С помощью этого последнего выражения мы смогли выразить время как функционал в следующей форме

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

где

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

На этом этапе мы можем опустить множитель \sqrt{2g}, поскольку оптимизация J идентична оптимизации \sqrt{2g}J.

Таким образом, мы можем теперь построить уравнение Эйлера-Лагранжа, следуя той же процедуре, которую мы использовали ранее, и в итоге получить:

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

Однако здесь мы видим, что \dfrac{\partial f}{\partial x} = 0,, следовательно, имеем

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

или, другими словами

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

где a — произвольная постоянная, записанная таким образом, потому что это «удобно» для последующих разработок.

Решение вариационной задачи

Подставив функцию f в это последнее выражение, получаем:

\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}

Чтобы решить этот интеграл, можно рассмотреть следующее подстановку

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

С этим имеем:

\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}

Можно заметить, что кривая брахистохрона может быть выражена как параметрическая кривая в полярных координатах, которая совпадает с циклойдой, начинающейся в начале координат.

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

Константа интегрирования C была равна нулю, чтобы удовлетворить начальное условие, что траектория начинается в начале координат. Кроме того, можно заметить, что существует пара уравнений, которые предоставляют возможные решения задачи, где константу a можно настроить так, чтобы кривая проходила через точку (x_2,y_2) в конце траектории. Эти уравнения:

Вариант 1:\boxed{\begin{array} {} & x(\theta) &=& a(\theta - \sin(\theta)) \\ & y(\theta) &=& a(1-\cos(\theta)) \end{array}}

Вариант 2:\boxed{\begin{array} {} & x(\theta) &=& - a(\theta - \sin(\theta)) \\ & y(\theta) &=& a(1-\cos(\theta)) \end{array}}

Жизнеспособное решение для этой задачи дается вторым вариантом, и, настроив константу a как отрицательное значение, мы получаем кривую, которая удовлетворяет необходимым условиям для решения.

Окончательная настройка решения

После последних настроек кривая брахистохрона имеет параметрическую форму:

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

Здесь было подставлено a=-b, где 0\lt b. Кривая имеет период 2b\pi и должна удовлетворять условию x_2 \in ]0,2b\pi[ и y_2 \in ]-2b,0[. Последнее особенно важно, так как требует, чтобы кривая брахистохрона представлялась как одна дуга циклоида, поскольку решение перестанет быть действительным, если частица снова вернется в состояние покоя при возвращении к точке нулевой высоты.

Чтобы настроить эти уравнения для задачи, нам необходимо найти значения \theta и b, которые удовлетворяют системе:

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

Эта нелинейная система, похоже, не имеет аналитических решений, поэтому мы используем численные методы в Wolfram Mathematica. Далее приведены шаги для решения задачи:

Шаг 1: Установить систему

Установить уравнения, составляющие систему для решения

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

Шаг 2: Определить конечную точку

Определить точку, в которую частица прибудет в конце своего пути. В данном случае установим ее в (x_2,y_2)=(1,-2). Эти значения можно изменить, чтобы попробовать другие аналогичные конфигурации.

x2val = 1; y2val = -2;

Шаг 3: Численно вычислить искомые значения

Использовать функцию «FindRoot» для численного решения задачи

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

Здесь были использованы значения b=1 и \theta=1 в качестве начальной точки для численного приближения решения. Таким образом, получаем решение b\approx 2.4056 и \theta \approx 1.40138

Шаг 4: Проверка результатов

Напомним, что для того чтобы эти ответы имели физический смысл, необходимо, чтобы x_2 \in ]0,2b\pi[ и y_2 \in ]-2b, 0[. Мы можем быстро проверить, что это так, с помощью следующего процесса

Сначала извлекаем значения b и \theta, полученные в качестве решения

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

Затем даем команду на выполнение подтверждения

If[0 < x2val < 2*Pi*bval && -2*bval < y2val < 0 "Допустимые значения", "Недопустимые значения"]

Если все прошло хорошо, мы должны получить "Допустимые значения" в выводе. Этот фрагмент кода поможет вам проверить, правильно ли смоделирована физическая ситуация.

С помощью этих процедур у нас наконец-то полностью настроена наша кривая решения, которая соединяет точки (x_1,y_1)=(0,0) и (x_2,y_2)=(1,-2). Получившаяся кривая:

\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]

Графически это выглядит так:

Репозиторий GitHub с алгоритмом Wolfram

Полный код решения задачи о брахистохроне, включая алгоритм, разработанный в Wolfram Mathematica, доступен для скачивания и ознакомления в моем репозитории GitHub. Этот репозиторий включает в себя файл `.nb` с кодом в формате интерактивного блокнота, а также версию в виде простого текста `.m` для тех, кто предпочитает просматривать код напрямую.

Вы можете скачать репозиторий с GitHub здесь.

Помимо кода, репозиторий содержит файл "README" с подробными инструкциями по использованию и пониманию алгоритма, а также пошаговое объяснение решения задачи о брахистохроне. Надеюсь, это будет полезно!