Предисловие
По роду профессиональной и научной деятельности я механик. Преподаю теоретическую механику в университете, пишу докторскую диссертацию в области динамики подвижного состава железных дорог. В общем, эта наука поглощает большую часть моего рабочего и даже свободного времени.
С Maple (на кафедре была 6-я версия, а у лоточников домой была куплена 8-я) познакомился ещё студентом, когда начинал работать над будущей кандидатской под крылом моего первого (ныше покойного) научного руководителя. Были и добрые люди, что помогли на самом первом этапе разобраться с пакетом и начать работать.
И вот так постепенно на его плечи была переложена большая часть вычислительной работы по подготовке диссертации. Диссертация была защищена, а Maple навсегда остался надёжным помошником в научном труде. Часто бывает необходимо быстро оценить какую-нибудь задачу, составить уравнения, исследовать их аналитически, быстро получить численное решение, построить графики. В этом отношении Maple просто незаменим для меня (ни в коем разе не хочу обидеть приверженцев других пакетов).
Сделать всё то, что будет предложено читателю под катом, меня подвигла задача принесенная ученицей (приходится ещё заниматься и репетиторством) со школьной олимпиады. Условие задачи таково:
Груз, висящий на нити длины L = 1,1 м, привязанной к гвоздю, толкнули так, что он поднялся, а затем ударился в гвоздь. Какова его скорость в момент удара о гвоздь? Ускорение свободного падения g = 10 м/с2.
Если не придираться к некоторонной туманности условия, то задача достаточно проста, а её решение, полученное путем довольно громоздких для школьника выкладок, в общем виде дает результат
![v = \sqrt{gL\sqrt{3}}](https://habrastorage.org/getpro/habr/post_images/528/490/18f/52849018f722c050f2d3458a9820c088.gif)
И вот тут захотелось проверить решение, полученное с оглядкой на школьную программу по физике независимым способом, например составив дифференциальные уравнения движения этого маятника, да не просто, а с учетом освобождения от связи (в процессе движения нить, считаемая невесомой, провисает и маятник движется как свободная точка).
Это послужило катализатором для того, чтобы взять да и откопать свои старые задумки, накопленные ещё со времен работы в оргкомитете Всероссийской Олимпиады студентов по теоретической механике — три года подряд занимался там подготовкой задач компьютерного конкурса. Задумки касались автоматизации построения уравнений движений для механических систем с неудерживающими связями и трением, используя известные всем уравнения Лагранжа 2 рода
![\frac{d}{dt}\left(\frac{\partial T}{\partial \dot{q}_i} \right ) - \frac{\partial T}{\partial q_i} = Q_i, \quad i=\overline{1,s}](https://habrastorage.org/getpro/habr/post_images/49e/930/2f7/49e9302f7d3c41f9a99d7a4a5d2689b8.gif)
поборов стереотип многих преподавателей о том, что уравнения эти неприменимы к системам с неудерживающими связями и трением.
Что касается Maple, то его библиотека для решения задач вариационного исчисления дает возможность быстро получить уравнения Эйлера-Лагранжа, решение которых минимизирует действие по Гамильтону, что применимо для консервативных систем
![\frac{d}{dt}\left(\frac{\partial L}{\partial \dot{q}_i} \right ) - \frac{\partial L}{\partial q_i} = 0, \quad i=\overline{1,s}](https://habrastorage.org/getpro/habr/post_images/fb5/1c3/696/fb51c369678e831f7d8e55b742c941c8.gif)
где
![L = T - \Pi](https://habrastorage.org/getpro/habr/post_images/fd7/ae7/5d0/fd7ae75d01f259d1dc8a9f025df7d048.gif)
Так как расматриваемые задачи не относятся к классу консервативных, то автором была предпринята попытка самостоятельно реализовать автоматизацию построения и анализа уравнений движений. Что из этого вышло, изложено под катом
1. Метод избыточных координат
Рассматриваем механическую систему, имеющую s степеней свободы, положение которой описывается вектором обобщенных координат
![\vec{q} = \left[q_1\left(t \right ), q_2\left(t \right ),...,q_s\left(t \right ) \right ]^T](https://habrastorage.org/getpro/habr/post_images/65e/5a8/ca3/65e5a8ca3b73816cd2186a5dc4ffd618.gif)
Учет неудерживающих связей требует от нас определения и анализа величины их реакций, поэтому необходимо так же определить их величину. Уберем указанные связи и введем дополнительно r обобщенных координат, выразив через них кинетическую энергию системы
![T = T\left(q_1,...,q_s, q_{s+1},...,q_{s+r}, \dot{q}_1,...,\dot{q}_s, \dot{q}_{s+1},..., \dot{q}_{s+r} \right )](https://habrastorage.org/getpro/habr/post_images/032/d72/e09/032d72e09daf20e235544d2e8455c1ef.gif)
Составим s + r уравнений движения в форме уравнений Лагранжа 2 рода
![\frac{d}{dt}\left(\frac{\partial T}{\partial \dot{q}_i} \right ) - \frac{\partial T}{\partial q_i} = Q_i, \quad i=\overline{1,s+r}\quad\quad(1)](https://habrastorage.org/getpro/habr/post_images/3e0/c56/ab4/3e0c56ab4447cf3cd64f71600ab5be1c.gif)
содержащие s+r неизвестных координат и r неизвестных реакций связей. Считая связи удерживающими, дополняем данную систему уравнениями связей (для простоты рассматривая геометрические связи) в виде
![f_j\left(q_1,...,q_{s+r} \right ) = 0, \quad j=\overline{1,r}](https://habrastorage.org/getpro/habr/post_images/5cf/7d6/7ba/5cf7d67bae6f3809b6880fe75d9aa3a0.gif)
получаем замкнутую систему уравнений, из которой находятся значения реакций
![R_j = R_j\left(q_1,...,q_{s},\dot{q}_1,...\dot{q}_s \right ),\quad j=\overline{1,r}\quad\quad(2)](https://habrastorage.org/getpro/habr/post_images/9c4/c22/2ff/9c4c222ff256f764166975ab40d8e878.gif)
являющиеся функциями первых s (независимых) обобщенных координат и скоростей и они могут быть расчитаны на любом шаге интегрирования уравнений движения (1). Для удерживающих связей типа «нить/поверхность» уравнения (1) и (2) надо дополнить условием освобождения от связи
![R_j = \left\{\begin{matrix} 0,\quad R_j \le 0 \\ R_j,\quad R_j > 0 \end{matrix}\right.\quad \quad (3)](https://habrastorage.org/getpro/habr/post_images/a92/995/239/a9299523939bfe79bf5bf9a45f31d50c.gif)
а для связей с сухим трением вида
![F_j\vec{\tau}+N_j\vec{n}](https://habrastorage.org/getpro/habr/post_images/4d9/4f5/071/4d94f50710ea5e35283329dab1914952.gif)
![F_j = \left\{\begin{matrix} &F_j, \quad \left|F_j\right| \le fN_j\\ &-fN_j \cfrac{v_j}{\left|v_j\right|}, \quad \left|F_j\right| > fN_j \end{matrix}\right. \quad (4)](https://habrastorage.org/getpro/habr/post_images/f85/007/107/f85007107b61d6623e0aa09454f051c4.gif)
где Fj и Nj соответственно касательная и нормальная составляющая реакции; vj — проекция скорости относительного проскальзывани точки приложения реакции.
Таким образом, уравнения (1) — (4) представляют собой полную математическую модель движения рассматриваемой механической системы.
Засим с теорией можно покончить и перейти к практике
2. Maple-функции построения и анализа уравнений Лагранжа
Для решения этой задачи была написана 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 должна содержать избыточные координаты в явном виде, то есть иметь вид
![q_{s+j} = g\left(q_1,..q_s \right ) \quad j=\overline{1,r}](https://habrastorage.org/getpro/habr/post_images/574/c98/8c6/574c988c6e7746faecab55aaf7f967ae.gif)
в противном случае функции библиотеки не смогут обработать уравнения правильно. Для тестирования возможностей библиотеки сойдет и так, но в дальнейшем этот момент будет переработан: просто пока неясно, будет ли гарантированно разрешена система уравнений связи относительно угловых избыточных координат.
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(R);
# Формируем вектор переменных, относительно которых будем решать уравнения движения
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 — массив реакций, относительно которых необходимо разрешить уравнения движения.
Теперь проиллюстрируем, как применять описанное «хозяйство» в деле
3. Задача о маятнике на тонкой нерастяжимой нити
Расчетная схема будет такой. В качестве обобщенной координаты выбираем угол
![\varphi](https://habrastorage.org/getpro/habr/post_images/171/69a/674/17169a6741a1bb27b86a1280ebc0586a.gif)
![](https://habrastorage.org/files/116/312/00c/11631200c37c4b47823cc2b27f975f4d.png)
Поскольку нить — неудерживающая связь, нас будет интересовать её реакция, а значит введем дополнительную, избыточную координату 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-нотации)
![T = 1/2\,m \left( \left( r \left( t \right) \right) ^{2} \left( {\frac { d}{dt}}\phi \left( t \right) \right) ^{2}+ \left( {\frac {d}{dt}}r \left( t \right) \right) ^{2} \right)](https://habrastorage.org/getpro/habr/post_images/606/7f1/e23/6067f1e23de923b9c8e49f86491867ec.gif)
Формируем массив сил и массив координат точек их приложения
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):
получая на выходе уравнения движения
![m{\frac {d^{2}}{d{t}^{2}}}r \left( t \right) -mr \left( t \right) \left( {\frac {d}{dt}}\phi \left( t \right) \right) ^{2}=mg\cos \left( \phi \left( t \right) \right) -S \quad \quad (5)](https://habrastorage.org/getpro/habr/post_images/7b7/6ad/58a/7b76ad58a593f7252c64786f74ba2a08.gif)
![2\,mr \left( t \right) \left( {\frac {d}{dt}}\phi \left( t \right) \right) {\frac {d}{dt}}r \left( t \right) +m \left( r \left( t \right) \right) ^{2}{\frac {d^{2}}{d{t}^{2}}}\phi \left( t \right) = -mgr \left( t \right) \sin \left( \phi \left( t \right) \right) \quad \quad (6)](https://habrastorage.org/getpro/habr/post_images/ad4/809/3f7/ad48093f7530af6c07986b2f7d155cf9.gif)
Нетрудно убедится, что функция отработала нормально — для иллюстрации специально выбрана не слишком громоздкая задача.
Далее задаем уравнение связи — пока нить натянута, справедливо условие
![r(t) = L](https://habrastorage.org/getpro/habr/post_images/bb5/40b/0d1/bb540b0d1c271dd2594db8b4301b6f40.gif)
преобразуем систему с учетом этого условия и находим реакцию связи
link_eqs := {r(t) = L};
simple_eqs := ReduceSystem(EQs, link_eqs, q);
solv1 := SolveAccelsReacts(simple_eqs, [phi(t)], [S]);
Сила натяжения нити равна
![S=m \left( {\frac {d}{dt}}\phi \left( t \right) \right) ^{2}L+mg\cos \left( \phi \left( t \right) \right) \quad \quad (7)](https://habrastorage.org/getpro/habr/post_images/4a6/111/430/4a61114307b36e9b026543fe14b55a12.gif)
Система (5) — (7) является полной системой уравнений движения груза, с учетом возможности провисания нити. Теперь подготовим её к численному интегрированию. Для начала разрешим её относительно ускорений, передав в SolveAccelsReacts() уравнения (5) и (6), вектор обобщенных координат и пустой массив реакций
EQs2 := SolveAccelsReacts(EQs, q,[]);
получая на выходе
![{\frac {d^{2}}{d{t}^{2}}}\phi \left( t \right) =-{\frac {2\, \left( { \frac {d}{dt}}r \left( t \right) \right) {\frac {d}{dt}}\phi \left( t \right) +\sin \left( \phi \left( t \right) \right) g}{r \left( t \right) }} \quad \quad (8)](https://habrastorage.org/getpro/habr/post_images/897/e54/52e/897e5452eb0f71d1e662ad71f257d5db.gif)
![{\frac {d^{2}}{d{t}^{2}}}r \left( t \right) ={\frac {mr \left( t \right) \left( {\frac {d}{dt}}\phi \left( t \right) \right) ^{2}+mg \cos \left( \phi \left( t \right) \right) -S}{m}} \quad \quad (9)](https://habrastorage.org/getpro/habr/post_images/d89/bca/5ca/d89bca5ca6191cc1a80118083edc08a5.gif)
Для численного моделирования, хоть это и не спортивно, напишем отдельный код, дабы не забивать голову читателя длительной обработкой полученной системы напильником. Тем более что моделирование будет иметь свои особенности.
Готовим исходные данные и систему уравнений движения
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]);
В итоге, получаем результат, приведенный на скриншоте. Скорость груза в момент удара соответствует приведенному в предисловии значению, и видно, что до провисания нити груз движется по окружности, а после провисания нити движется как свободная точка под действием силы тяжести, по параболе.
![](https://habrastorage.org/files/026/5d8/e6c/0265d8e6c32740a691339366ccf95eab.png)
Замечу, что погрешности попадания в гвоздь — вынужденная мера: в полярных координатах, которые были использованы, задача имеет особенность, понятную из уравнения (8). Поэтому r(t) сравнивалось не с нулем, а с величиной eps достаточно малой, чтобы получить решение, и достаточно большой, чтобы численный решатель fsolve() не сходил с ума. Однако это нисколько не умаляет практической ценности изложенных результатов.
Вместо заключения
Возможно, читатель упрекнет меня, что я стреляю из пушки по воробьям. Однако, хочется заметить, что всё сложное начинается с простого, а большая наука — с малых задач.
Тестовую версию библиотеки можно качнуть тут
Благодарю за внимание к моему труду )