Maple: складання рівнянь Лагранжа 2 роду і метод надлишкових координат

Передмова

За родом професійної та наукової діяльності я механік. Викладаю теоретичну механіку в університеті, пишу докторську дисертацію в галузі динаміки рухомого складу залізниць. Загалом, ця наука поглинає більшу частину мого робочого і навіть вільного часу.

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

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

Зробити все те, що буде запропоновано читачеві під катом, мене завдання принесена ученицею (доводиться займатися ще і репетиторством) зі шкільної олімпіади. Умова завдання таке:
Вантаж, що висить на нитці довжини L = 1,1 м, прив'язаною до цвяха, штовхнули так, що він піднявся, а потім ударився в цвях. Яка його швидкість у момент удару об цвях? Прискорення вільного падіння g = 10 м/с2.
Якщо не чіплятися до некоторонной туманності умови, то задача досить проста, а її рішення, отримане шляхом досить громіздких для школяра викладок, в загальному вигляді дає результат



І ось тут захотілося перевірити рішення, отримане з оглядкою на шкільну програму з фізики незалежним способом, наприклад склавши диференціальні рівняння руху цього маятника, та не просто, а з урахуванням звільнення від зв'язку (у процесі руху нитка, считаемая невагомою, провисає і маятник рухається як вільна точка).

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



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

Що стосується Maple, то його бібліотека для вирішення задач варіаційного числення дає можливість швидко отримати рівняння Ейлера-Лагранжа, рішення яких мінімізує дію з Гамільтону, що застосовується для консервативних систем



де — функція Лагранжа, дорівнює різниці кінетичної і потенційної енергій системи.

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



1. Метод надлишкових координат

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

Облік неудерживающих зв'язків вимагає визначення та аналізу величини їх реакцій, тому необхідно визначити їх величину. Приберемо зазначені зв'язки і введемо додатково r узагальнених координат, висловивши через них кінетичну енергію системи



Складемо s + r рівнянь руху у формі рівнянь Лагранжа 2 роду



містять s+r невідомих координат і r невідомих реакцій зв'язків. Вважаючи зв'язку утримують, доповнюємо дану систему рівнянь зв'язків (для простоти розглядаючи геометричні зв'язку) у вигляді



отримуємо замкнену систему рівнянь, з якої знаходяться значення реакцій



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



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



де 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 повинна містити надлишкові координати в явному вигляді, тобто мати вигляд



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

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 — масив реакцій, щодо яких необхідно розв'язати рівняння руху.

Тепер проілюструємо, як застосовувати описане «господарство» у справі

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]);
# Якщо воно негативно - нитка провисла, реакції немає
as if < 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), v(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 := v(t1);

# Будуємо траєкторію вантажу
plot([x(t), y(t), t=0..t1], view=[-L..L, -L..L]);

В результаті, отримуємо результат, наведений на скріншоті. Швидкість вантажу в момент удару відповідає наведеному в передмові значенням, і видно, що до провисання нитки вантаж рухається по колу, а після провисання нитки рухається як вільна точка під дією сили тяжіння, по параболі.


Зауважу, що похибки попадання у цвях — вимушена міра: у полярних координатах, які були використані, задача має особливість, зрозумілу з рівняння (8). Тому r(t) порівнювалося не з нулем, а з величиною eps достатньо малою, щоб отримати рішення, і досить великий, щоб чисельний вирішувач fsolve() не сходив з розуму. Проте це анітрохи не применшує практичної цінності викладених результатів.

Замість висновку

Можливо, читач дорікне мені, що я стріляю з гармати по горобцях. Однак, хочеться зауважити, що все складне починається з простого, а велика наука — з малих завдань.

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

Дякую за увагу до моєї праці )

Джерело: Хабрахабр

0 коментарів

Тільки зареєстровані та авторизовані користувачі можуть залишати коментарі.