Предисловие

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

С Maple (на кафедре была 6-я версия, а у лоточников домой была куплена 8-я) познакомился ещё студентом, когда начинал работать над будущей кандидатской под крылом моего первого (ныше покойного) научного руководителя. Были и добрые люди, что помогли на самом первом этапе разобраться с пакетом и начать работать.

И вот так постепенно на его плечи была переложена большая часть вычислительной работы по подготовке диссертации. Диссертация была защищена, а Maple навсегда остался надёжным помошником в научном труде. Часто бывает необходимо быстро оценить какую-нибудь задачу, составить уравнения, исследовать их аналитически, быстро получить численное решение, построить графики. В этом отношении Maple просто незаменим для меня (ни в коем разе не хочу обидеть приверженцев других пакетов).

Сделать всё то, что будет предложено читателю под катом, меня подвигла задача принесенная ученицей (приходится ещё заниматься и репетиторством) со школьной олимпиады. Условие задачи таково:

Груз, висящий на нити длины L = 1,1 м, привязанной к гвоздю, толкнули так, что он поднялся, а затем ударился в гвоздь. Какова его скорость в момент удара о гвоздь? Ускорение свободного падения g = 10 м/с².

Если не придираться к некоторонной туманности условия, то задача достаточно проста, а её решение, полученное путем довольно громоздких для школьника выкладок, в общем виде дает результат

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

Это послужило катализатором для того, чтобы взять да и откопать свои старые задумки, накопленные ещё со времен работы в оргкомитете Всероссийской Олимпиады студентов по теоретической механике — три года подряд занимался там подготовкой задач компьютерного конкурса. Задумки касались автоматизации построения уравнений движений для механических систем с неудерживающими связями и трением, используя известные всем уравнения Лагранжа 2 рода

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

Что касается Maple, то его библиотека для решения задач вариационного исчисления дает возможность быстро получить уравнения Эйлера-Лагранжа, решение которых минимизирует действие по Гамильтону, что применимо для консервативных систем

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

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

1. Метод избыточных координат

Рассматриваем механическую систему, имеющую s степеней свободы, положение которой описывается вектором обобщенных координат . Пусть также имеется r неудерживающих связей, к числу реакций которых можно причислить и трение покоя, при превышении предельного значения переходящее в активную силу трения скольжения, направление которой противоположно направлению относительной скорости скольжения.

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

Составим s + r уравнений движения в форме уравнений Лагранжа 2 рода

содержащие s+r неизвестных координат и r неизвестных реакций связей. Считая связи удерживающими, дополняем данную систему уравнениями связей (для простоты рассматривая геометрические связи) в виде

получаем замкнутую систему уравнений, из которой находятся значения реакций

являющиеся функциями первых s (независимых) обобщенных координат и скоростей и они могут быть расчитаны на любом шаге интегрирования уравнений движения (1). Для удерживающих связей типа «нить/поверхность» уравнения (1) и (2) надо дополнить условием освобождения от связи

а для свзяей с сухим трением вида

где F_j и N_j соответственно касательная и нормальная составляющая реакции; v_j — проекция скорости относительного проскальзывани точки приложения реакции.

Таким образом, уравнения (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 должна содержать избыточные координаты в явном виде, то есть иметь вид

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

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. Задача о маятнике на тонкой нерастяжимой нити

Расчетная схема будет такой. В качетве обобщенной координаты выбираем угол наклона нити к вертикали.

Поскольку нить — неудерживающая связь, нас будет интересовать её реакция, а значит введем дополнительную, избыточную координату 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() не сходил с ума. Однако это нисколько не умаляет практической ценности изложенных результатов.

Вместо заключения

Возможно, читатель упрекнет меня, что я стреляю из пушки по воробьям. Однако, хочется заметить, что всё сложное начинается с простого, а большая наука — с малых задач.

Тестовую версию библиотеки можно качнуть тут

Благодарю за внимание к моему труду )