maisvendoo December 4, 2014 at 23:29

Maple: compiling Lagrange equations of the 2nd kind and the method of excess coordinates

Foreword

By the nature of my professional and scientific activity, I am a mechanic. I teach theoretical mechanics at the university, I am writing a doctoral dissertation in the field of dynamics of rolling stock of railways. In general, this science absorbs most of my work and even free time.

I became a student when I started working on a future candidate under the wing of my first (now deceased) scientific adviser with Maple (the department had the 6th version, and the 8th was bought from the hawkers home). There were also kind people who helped at the very first stage to deal with the package and start working.

And so, gradually, most of the computational work on the dissertation was shifted onto his shoulders. The dissertation was defended, and Maple forever remained a reliable assistant in scientific work. It is often necessary to quickly evaluate a problem, draw up equations, analyze them analytically, quickly get a numerical solution, and build graphs. In this regard, Maple is simply irreplaceable for me (I never want to offend the supporters of other packages).

To do everything that will be offered to the reader under the cut, I was inspired by the task brought by the student (I also have to deal with tutoring) from the school Olympiad. The condition of the problem is as follows:

A load hanging on a thread of length L = 1.1 m tied to a nail was pushed so that it lifted, and then hit the nail. What is its speed when it hits a nail? Acceleration of gravity g = 10 m / s ² .

If you don’t find fault with a certain nebula of conditions, then the problem is quite simple, and its solution, obtained by calculations that are rather cumbersome for the student, in general gives the result. And here I wanted to check the solution obtained with an eye on the school curriculum in physics in an independent way, for example by writing the differential equations of motion of this pendulum, and not just, but taking into account the release from communication (in the process of movement the thread, considered weightless, sags and the pendulum moves as a free point).

$v = \ sqrt {gL \ sqrt {3}}$

This served as a catalyst in order to dig up and dig out your old ideas, accumulated from the time of work in the organizing committee of the All-Russian Olympiad of students in theoretical mechanics - for three years in a row he was engaged in the preparation of the tasks of a computer competition. The ideas concerned automation of the construction of equations of motion for mechanical systems with non-holding bonds and friction, using the well-known Lagrange equations of the second kind of requisitions, the stereotype of many teachers that these equations are not applicable to systems with non-holding bonds and friction. As for Maple, its library for solving variational calculus problems makes it possible to quickly obtain Euler-Lagrange equations, the solution of which minimizes the Hamilton action, which is applicable for conservative systems where

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

$\ frac {d} {dt} \ left (\ frac {\ partial L} {\ partial \ dot {q} _i} \ right) - \ frac {\ partial L} {\ partial q_i} = 0, \ quad i = \ overline {1, s}$

$L = T - \ Pi$ Is the Lagrange function equal to the difference of the kinetic and potential energies of the system.

Since the problems under consideration are not conservative, the author made an attempt to independently automate the construction and analysis of equations of motion. What came of it is stated under the cut

1. The method of excess coordinates

We consider a mechanical system having s degrees of freedom, the position of which is described by a vector of generalized coordinates $\ vec {q} = \ left [q_1 \ left (t \ right), q_2 \ left (t \ right), ..., q_s \ left (t \ right) \ right] ^ T$ . Let there also be r non-retaining bonds, the number of reactions of which can also include rest friction, which, when the limit value is exceeded, passes into the active force of sliding friction, the direction of which is opposite to the direction of the relative sliding speed.

Taking into account non-retaining bonds requires us to determine and analyze the magnitude of their reactions, therefore it is also necessary to determine their magnitude. We remove the indicated connections and introduce an additional r generalized coordinates, expressing through them the kinetic energy of the system. Compose s + 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)$

equations of motion in the form of Lagrange equations of the 2nd kind containing s + r unknown coordinates and r unknown coupling reactions. Considering the constraints as constraints, we supplement the given system with constraint equations (for simplicity, considering geometric constraints) in the form we obtain a closed system of equations from which the reaction values are functions of the first s (independent) generalized coordinates and velocities and they can be calculated at any step of integration of the equations of motion (1). For the “thread / surface” type of holding bonds, equations (1) and (2) must be supplemented with the condition of release from the bond, and for bonds with dry friction of the form where F _j and N

$\ 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)$

$f_j \ left (q_1, ..., q_ {s + r} \ right) = 0, \ quad j = \ overline {1, r}$

$R_j = R_j \ left (q_1, ..., q_ {s}, \ dot {q} _1, ... \ dot {q} _s \ right), \ quad j = \ overline {1, r} \ quad \ quad (2)$

$R_j = \ left \ {\ begin {matrix} 0, \ quad R_j \ le 0 \\ R_j, \ quad R_j> 0 \ end {matrix} \ right. \ Quad \ quad (3)$

$F_j \ vec {\ tau} + N_j \ vec {n}$

$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)$

_j respectively the tangent and normal component of the reaction; v _j is the projection of the relative slippage rate of the reaction application point.

Thus, equations (1) - (4) are a complete mathematical model of the motion of the mechanical system under consideration.

Then you can end theory and move on to practice.

2. Maple-functions for constructing and analyzing Lagrange equations

To solve this problem was written Maple-library lagrange , containing four functions

LagrangeEQs - construction of the Lagrange equations of motion in the form of 2 genus

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:

As input parameters, the function takes the expression of the kinetic energy T as a function of generalized coordinates and generalized velocities; array of generalized coordinates q ; an array of radius vectors of forces application points r and the array of force vectors F .

LinksEQs - Obtaining differential link equations from geometric link equations

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:

It should be noted here that the system of equations of geometric relationships eqs must contain redundant coordinates in an explicit form, that is, have the form otherwise the library functions will not be able to process the equations correctly. To test the capabilities of the library, it will do so, but in the future this moment will be revised: it is simply not yet clear whether the system of equations of coupling with respect to angular redundant coordinates is guaranteed to be resolved. ReduceSystem - transformation of equations of motion taking into account the equations of relations

$q_ {s + j} = g \ left (q_1, .. q_s \ right) \ quad j = \ overline {1, r}$

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:

This code does not need detailed explanations - here the substitution of redundant generalized coordinates, velocities and accelerations, expressed by equations of geometric and differential constraints, into the equations of motion is carried out in order to bring them to a form suitable for calculating reactions of unstable bonds

SolveAccelsReacts - solving equations of motion with respect to reactions and generalized accelerations

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:

This function accepts the input system of equations of motion eqs, transformed taking into account the equations of constraints. It is linear with respect to the second derivatives of independent coordinates and bond reactions. Other input parameters: q - vector of independent coordinates; R is an array of reactions for which it is necessary to solve the equations of motion.

Now we illustrate how to apply the described "economy" in business

3. The problem of a pendulum on a thin inextensible thread

The design scheme will be like this. As a generalized coordinate, select the angle of $\ varphi$ inclination of the thread to the vertical.

Since the thread is an unrestraining connection, we will be interested in its reaction, and therefore we introduce an additional, excess coordinate r ( t ).

Getting down. We clean the memory and connect the linear algebra library

restart;
with(LinearAlgebra):

Connect the lagrange library

read `/home/maisvendoo/work/maplelibs/mechanics/lagrange.m`;

We determine the vector of generalized coordinates, calculate the coordinates and speed of the load, as well as the kinetic energy of the system

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);

At the output, we get the expression for kinetic energy (the latex () function was used to insert here, generating the result in LaTeX notation). We form an array of forces and an array of coordinates of the points of their application
$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)$

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

We feed all functions of LagrangeEQs ()

EQs := LagrangeEQs(T, q, rk, Fk):

receiving the equations of motion at the output It is easy to verify that the function worked fine - a not too cumbersome task was specially chosen for illustration. Next, we set the coupling equation - while the thread is stretched, the condition is true, we transform the system taking this condition into account and find the coupling reaction
$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)$

$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)$

link_eqs := {r(t) = L};
simple_eqs := ReduceSystem(EQs, link_eqs, q); 
solv1 := SolveAccelsReacts(simple_eqs, [phi(t)], [S]);

The force of tension of the thread is equal. System (5) - (7) is a complete system of equations of motion of the load, taking into account the possibility of sagging of the thread. Now prepare it for numerical integration. To begin with, let us solve it with respect to accelerations, passing equations (5) and (6), a vector of generalized coordinates, and an empty reaction array to SolveAccelsReacts ()
$S = m \ left ({\ frac {d} {dt}} \ phi \ left (t \ right) \ right) ^ {2} L + mg \ cos \ left (\ phi \ left (t \ right) \ right) \ quad \ quad (7)$

EQs2 := SolveAccelsReacts(EQs, q,[]);

getting the output For numerical modeling, although it’s not sports, we will write a separate code so as not to clog the reader’s head with a lengthy file processing of the resulting system. Moreover, modeling will have its own characteristics. We prepare the initial data and the system of equations of motion

${\ 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)$

${\ 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)$

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:

We build the function of calculating the state of the system, for a given horizontal initial load speed

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:

Now we check the "school" solution to the problem

# Такая начальная скорость должна быть, согласно школьному решению задачи
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]);

As a result, we get the result shown in the screenshot. The speed of the load at the moment of impact corresponds to the value given in the preface, and it can be seen that until the sag of the thread, the load moves in a circle, and after the sag of the thread moves as a free point under the action of gravity, along a parabola.

I note that the error in hitting the nail is a necessary measure: in the polar coordinates that were used, the problem has a feature that is clear from equation (8). Therefore, r ( t ) was compared not with zero, but with eps small enough to get a solution, and big enough so that the numerical solver fsolve () would not go crazy. However, this does not detract from the practical value of the results presented.

Instead of a conclusion

Perhaps the reader will reproach me for shooting sparrows from a cannon. However, I would like to note that everything complicated begins with the simple, and big science - with small tasks.

The test version of the library can be downloaded here

Thank you for your attention to my work)

Tags: