Maple: составление уравнений Лагранжа 2 рода и метод избыточных координат
По роду профессиональной и научной деятельности я механик. Преподаю теоретическую механику в университете, пишу докторскую диссертацию в области динамики подвижного состава железных дорог. В общем, эта наука поглощает большую часть моего рабочего и даже свободного времени.
С Maple (на кафедре была 6-я версия, а у лоточников домой была куплена 8-я) познакомился ещё студентом, когда начинал работать над будущей кандидатской под крылом моего первого (ныше покойного) научного руководителя. Были и добрые люди, что помогли на самом первом этапе разобраться с пакетом и начать работать.
И вот так постепенно на его плечи была переложена большая часть вычислительной работы по подготовке диссертации. Диссертация была защищена, а Maple навсегда остался надёжным помошником в научном труде. Часто бывает необходимо быстро оценить какую-нибудь задачу, составить уравнения, исследовать их аналитически, быстро получить численное решение, построить графики. В этом отношении Maple просто незаменим для меня (ни в коем разе не хочу обидеть приверженцев других пакетов).
Сделать всё то, что будет предложено читателю под катом, меня задача принесенная ученицей (приходится ещё заниматься и репетиторством) со школьной олимпиады. Условие задачи таково:
Груз, висящий на нити длины L = 1,1 м, привязанной к гвоздю, толкнули так, что он поднялся, а затем ударился в гвоздь. Какова его скорость в момент удара о гвоздь? Ускорение свободного падения g = 10 м/с2.
Если не придираться к некоторонной туманности условия, то задача достаточно проста, а её решение, полученное путем довольно громоздких для школьника выкладок, в общем виде дает результат
И вот тут захотелось проверить решение, полученное с оглядкой на школьную программу по физике независимым способом, например составив дифференциальные уравнения движения этого маятника, да не просто, а с учетом освобождения от связи (в процессе движения нить, считаемая невесомой, провисает и маятник движется как свободная точка).
Это послужило катализатором для того, чтобы взять да и откопать свои старые задумки, накопленные ещё со времен работы в оргкомитете Всероссийской Олимпиады студентов по теоретической механике — три года подряд занимался там подготовкой задач компьютерного конкурса. Задумки касались автоматизации построения уравнений движений для механических систем с неудерживающими связями и трением, используя известные всем уравнения Лагранжа 2 рода
поборов стереотип многих преподавателей о том, что уравнения эти неприменимы к системам с неудерживающими связями и трением.
Что касается Maple, то его библиотека для решения задач вариационного исчисления дает возможность быстро получить уравнения Эйлера-Лагранжа, решение которых минимизирует действие по Гамильтону, что применимо для консервативных систем
где — функция Лагранжа, равная разности кинетической и потенциальной энергий системы.
Так как расматриваемые задачи не относятся к классу консервативных, то автором была предпринята попытка самостоятельно реализовать автоматизацию построения и анализа уравнений движений. Что из этого вышло, изложено под катом
Рассматриваем механическую систему, имеющую s степеней свободы, положение которой описывается вектором обобщенных координат . Пусть также имеется r неудерживающих связей, к числу реакций которых можно причислить и трение покоя, при превышении предельного значения переходящее в активную силу трения скольжения, направление которой противоположно направлению относительной скорости скольжения.
Учет неудерживающих связей требует от нас определения и анализа величины их реакций, поэтому необходимо так же определить их величину. Уберем указанные связи и введем дополнительно r обобщенных координат, выразив через них кинетическую энергию системы
Составим s + r уравнений движения в форме уравнений Лагранжа 2 рода
содержащие s+r неизвестных координат и r неизвестных реакций связей. Считая связи удерживающими, дополняем данную систему уравнениями связей (для простоты рассматривая геометрические связи) в виде
получаем замкнутую систему уравнений, из которой находятся значения реакций
являющиеся функциями первых s (независимых) обобщенных координат и скоростей и они могут быть расчитаны на любом шаге интегрирования уравнений движения (1). Для удерживающих связей типа «нить/поверхность» уравнения (1) и (2) надо дополнить условием освобождения от связи
а для свзяей с сухим трением вида
где Fj и Nj соответственно касательная и нормальная составляющая реакции; vj — проекция скорости относительного проскальзывани точки приложения реакции.
Таким образом, уравнения (1) — (4) представляют собой полную математическую модель движения рассматриваемой механической системы.
Засим с теорией можно покончить и перейти к практике
Для решения этой задачи была написана Maple-библиотека lagrange, содержащая четыре функции
LagrangeEQs — построение уравнений движения в форме Лагранжа 2 рода
LagrangeEQs:= proc (T, q, r, F)
local s:= numelems (q); local n:= numelems (rk); local i, k; local T1, dT1dv; local dTdv, dTdvdt; local T2, dT2dq; local dTdq; local left_part; local Q; local summa; local r1, dr1dq, drdq;
# Получение левой части уравнений движения for i from 1 to s do # Дифференцируем кинетическую энергию по обобщенным скоростям и времени T1[i] := subs (diff (q[i], t) = v[i], T); dT1dv[i] := diff (T1[i], v[i]); dTdv[i] := subs (v[i] = diff (q[i], t), dT1dv[i]); dTdvdt[i] := diff (dTdv[i], t);
# Дифференцируем кинетическую энергию по обобщенным координатам T2[i] := subs (q[i] = q1[i], T); dT2dq[i] := diff (T2[i], q1[i]); dTdq[i] := subs (q1[i] = q[i], dT2dq[i]); # Формируем левую часть уравнения движения left_part[i] := expand (simplify (dTdvdt[i] — dTdq[i])); end do;
VectorCalculus[BasisFormat](false);
# Вычисляем обобщенные силы (правая часть уравнений движения) for i from 1 to s do summa:= 0; for k from 1 to n do
# Дифференцируем радиус-ректор точки приложения k-й силы по i-й обобщенной координате r1[k] := subs (q[i] = q1[i], r[k]); dr1dq[k] := VectorCalculus[diff](r1[k], q1[i]); drdq[k] := subs (q1[i] = q[i], dr1dq[k]); # Скалярно перемножаем вектор силы на производную от радиус-вектора по обобщенной координате # и накапливаем результат summa:= summa + LinearAlgebra:-DotProduct (F[k], drdq[k], conjugate = false); end do;
Q[i] := expand (simplify (summa)); end do;
# Окончательно формируем уравнения и возвращаем результатq return {seq (left_part[i] = Q[i], i=1…s)};
end proc: В качестве входных параметров функция принимает выражение кинетической энергии T как функцию обобщенных координат и обобщенных скоростей; массив обобщенных координат q; массив радиус-векторов точек приложения сил r и массив векторов сил F.
LinksEQs — получение уравнений дифференциальных связей из уравнений геометрических связей
LinksEQs:= proc (eqs)
local Eq1, Eq2; local i; local r:= numelems (eqs); # Дважды дифференцируем уравнения связей по времени for i from 1 to r do Eq1[i] := diff (lhs (eqs[i]), t) = diff (rhs (eqs[i]), t); Eq2[i] := diff (diff (lhs (eqs[i]), t), t) = diff (diff (rhs (eqs[i]), t), t); end do;
# и возвращаем результат return {seq (eqs[i], i=1…r), seq (Eq1[i], i=1…r), seq (Eq2[i], i=1…r)};
end proc: Здесь надо отметить, что система уравнений геометрических связей eqs должна содержать избыточные координаты в явном виде, то есть иметь вид
в противном случае функции библиотеки не смогут обработать уравнения правильно. Для тестирования возможностей библиотеки сойдет и так, но в дальнейшем этот момент будет переработан: просто пока неясно, будет ли гарантированно разрешена система уравнений связи относительно угловых избыточных координат.
ReduceSystem — преобразование уравнений движения с учетом уравнений связей
ReduceSystem:= proc (eqs, links, q)
local i, j, k; local links_eqs:= LinksEQs (links);
local r:= numelems (links_eqs); local s:= numelems (q);
local eq:= [seq (eqs[i], i=1…s)];
for i from 1 to s do for j from 1 to r do eq[i] := simplify (algsubs (links_eqs[j], eq[i])); end do: end do:
return {seq (eq[i], i=1…s)};
end proc: Данный код в подробных пояснениях не нуждается — тут выполняется подстановка избыточных обобщенных координат, скоростей и ускорений, выражаемых уравнениями геометрических и дифференциальных связей в уравнения движения, с целью приведения их к виду, пригодному для вычисления реакций неудерживающих связей
SolveAccelsReacts — решение уравнений движения относительно реакций и обобщенных ускорений
SolveAccelsReacts:= proc (eqs, q, R)
local s:= numelems (q); local r:= numelems®;
# Формируем вектор переменных, относительно которых будем решать уравнения движения local vars:= [seq (diff (diff (q[i], t), t), i=1…s), seq (R[i], i=1…r)]; local eq:= [seq (eqs[i], i=1…numelems (eqs))]; local i, j; local x; local solv; # Вводим подстановку — заменяем «иксами» все искомые переменные for i from 1 to numelems (eqs) do for j from 1 to s + r do eq[i] := subs (vars[j] = x[j], eq[i]); end do: end do;
# ищем «иксы» (система всегда линейна относительно них) solv:= solve ({seq (eq[i], i=1…numelems (eq))}, {seq (x[i], i=1…s+r)});
# Связываем иксы с найденными значениями assign (solv);
# Возвращаем уравнения, решенные относительно обобщенных ускорений и реакций return {seq (vars[i] = x[i], i=1…s+r)};
end proc: Данная функция принимает на вход систему уравнений движения eqs, преобразованную с учетом уравнений связей. Она линейна относительно вторых производных независимых координат и реакций связей. Другие входные параметры: q — вектор независимых координат; R — массив реакций, относительно которых необходимо разрешить уравнения движения.
Теперь проиллюстрируем, как применять описанное «хозяйство» в деле
Расчетная схема будет такой. В качетве обобщенной координаты выбираем угол наклона нити к вертикали.
Поскольку нить — неудерживающая связь, нас будет интересовать её реакция, а значит введем дополнительную, избыточную координату r (t).
Приступаем. Чистим память и подключаем библиотеку линейной алгебры
restart; with (LinearAlgebra): Подключаем библиотеку lagrange read `/home/maisvendoo/work/maplelibs/mechanics/lagrange.m`; Определяем вектор обобщенных координат, вычисляем координаты и скорость груза, а так же кинетическую энергию системы q:= [r (t), phi (t)];
xM:= q[1]*sin (q[2]); yM:= -q[1]*cos (q[2]);
vMx:= diff (xM, t); vMy:= diff (yM, t);
T:= simplify (m*(vMx^2 + vMy^2)/2);
На выходе получаем выражение для кинетической энергии (для вставки сюда использована функция latex (), генерирующая результат в LaTeX-нотации)Формируем массив сил и массив координат точек их приложения
Mg:= Vector ([0, -m*g]);
React:= Vector ([-S*sin (q[2]), S*cos (q[2])]);
rM:= Vector ([xM, yM]);
Fk:= [Mg, React];
rk:= [rM, rM];
Скармливаем всё функции LagrangeEQs ()
EQs:= LagrangeEQs (T, q, rk, Fk):
получая на выходе уравнения движения
Нетрудно убедится, что функция отработала нормально — для иллюстрации специально выбрана не слишком громоздкая задача.
Далее задаем уравнение связи — пока нить натянута, справедливо условие
преобразуем систему с учетом этого условия и находим реакцию связи
link_eqs:= {r (t) = L};
simple_eqs:= ReduceSystem (EQs, link_eqs, q);
solv1:= SolveAccelsReacts (simple_eqs, [phi (t)], [S]);
Сила натяжения нити равнаСистема (5) — (7) является полной системой уравнений движения груза, с учетом возможности провисания нити. Теперь подготовим её к численному интегрированию. Для начала разрешим её относительно ускорений, передав в SolveAccelsReacts () уравнения (5) и (6), вектор обобщенных координат и пустой массив реакций
EQs2:= SolveAccelsReacts (EQs, q,[]);
получая на выходе
Для численного моделирования, хоть это и не спортивно, напишем отдельный код, дабы не забивать голову читателя длительной обработкой полученной системы напильником. Тем более что моделирование будет иметь свои особенности.
Готовим исходные данные и систему уравнений движения
L:= 1.1: g:= 10.0:
# Функция вычисляет производные фазовых координат EQs_func:= proc (N, t, Y, dYdt) # Ускорение силы натяжения нити (as = S/m) local as:= 0; # Если нить уже провисла, то реакции нет if Y[1] < L then as := 0; else # Если нить натянута, вычисляем ускорение её реакции as := L*Y[4]^2 + g*cos(Y[2]); # Если оно отрицательно - нить провисла, реакции нет if as < 0 then as := 0; end if; end if;
# Собственно система уравнений в форме Коши # Y[1] → r (t) — расстояние от груза до гвоздя # Y[2] → phi (t) — угол радиус-вектора груза к вертикали # Y[3] → vr (t) — радиальная скорость груза # Y[4] → omega (t) — угловая скорость поворота радиус-вектора dYdt[1] := Y[3]; dYdt[2] := Y[4]; dYdt[3] := Y[1]*Y[4]^2 + g*cos (Y[2]) — as; dYdt[4] := -(2*Y[3]*Y[4] + g*sin (Y[2]))/Y[1]; end proc: Строим функцию вычисления состояния системы, при заданной горизонтальной начальной скорости груза
sys_pos:= proc (v0)
# Формируем начальные условия local initc:= Array ([L, 0, 0, v0/L]); # Задаем функции, которые ищем local q:= [r (t), phi (t), vr (t), omega (t)];
# Численно решаем систему ОДУ движения local dsolv:= dsolve (numeric, number = 4, procedure = EQs_func, start = 0, initial = initc, procvars = q, output=listprocedure);
# Выделяем из решения полученные функции local R:= eval (r (t), dsolv); local Phi:= eval (phi (t), dsolv); local Vr:= eval (vr (t), dsolv); local Omega:= eval (omega (t), dsolv);
return [R, Phi, Vr, Omega];
end proc: Теперь проверяем «школьное» решение задачи
# Такая начальная скорость должна быть, согласно школьному решению задачи v0:= evalf (sqrt (g*L*(2 + sqrt (3)))):
# Погрешность попадания груза в гвоздь eps:= 1e-5:
# Интегрируем уравнения и получаем решение r:= sys_pos (v0)[1]: phi:= sys_pos (v0)[2]: vr:= sys_pos (v0)[3]:
# Строим декартовы координаты груза x:= t→r (t)*sin (phi (t)): y:= t→-r (t)*cos (phi (t)):
# Определяем момент удара о гвоздь t1:= fsolve (r (t) = eps, t=0…10.0):
# Вычисляем скорость в момент удара v:= vr (t1);
# Строим траекторию груза
plot ([x (t), y (t), t=0…t1], view=[-L…L, -L…L]);
В итоге, получаем результат, приведенный на скриншоте. Скорость груза в момент удара соответствует приведенному в предисловии значению, и видно, что до провисания нити груз движется по окружности, а после провисания нити движется как свободная точка под действием силы тяжести, по параболе.
Замечу, что погрешности попадания в гвоздь — вынужденная мера: в полярных координатах, которые были использованы, задача имеет особенность, понятную из уравнения (8). Поэтому r (t) сравнивалось не с нулем, а с величиной eps достаточно малой, чтобы получить решение, и достаточно большой, чтобы численный решатель fsolve () не сходил с ума. Однако это нисколько не умаляет практической ценности изложенных результатов.
Возможно, читатель упрекнет меня, что я стреляю из пушки по воробьям. Однако, хочется заметить, что всё сложное начинается с простого, а большая наука — с малых задач.
Тестовую версию библиотеки можно качнуть тут
Благодарю за внимание к моему труду)
