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

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



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

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

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

    Similar posts

    AdBlock has stolen the banner, but banners are not teeth — they will be back

    More
    Ads

    Comments 13

      +3
      Прошу прощения, но условие действительно туманно (хотя, возможно, только мне это сложно представить):
      Груз, висящий на нити длины L = 1,1 м, привязанной к гвоздю, толкнули так, что он поднялся, а затем ударился в гвоздь. Какова его скорость в момент удара о гвоздь? Ускорение свободного падения g = 10 м/с2.

      Могу себе только представить, что груз толкнули вверх (строго по направлению нити), чтобы тот ударился о гвоздь. Если его толкнуть в бок, то в идеальном случае он опишет окружность вокруг гвоздя. Не в идеальном — нить намотается на гвоздь.
      Как все-таки толкнули груз, что он ударился в гвоздь? По какой траектории?
        +2
        Дополню: Если все-таки толкнули вверх, то могли толкнуть так, что в точке соприкосновения с гвоздем скорость может иметь значения от нуля (при минимальной силе, достаточной для того, чтобы груз достиг гвоздя) и выше. В этом случае, условие неполно.
          +1
          нуля (при минимальной силе, достаточной для того, чтобы груз достиг гвоздя)

          Не хочу показаться занудой, но я бы поправил Вас, сказав о минимальной начальной скорости направленной вертикально ().

          Сила которая сообщила грузу такую скорость перестала дейсвовать в момент начала полета и в рассмотрении не участвует. Хотя, опосредованно, о минимальной силе, разумеется можно говорить :)

          P.S.: Почувствовал себя Споком и Шелдоном одновременно…
            0
            Да, в данном случае о скорости говорить разумнее, т.к. от масса, которая в задаче не сообщается, нас интересовать не будет.
          +3
          Последний рисунок в статье- траектория решения. Гвоздь- в центре координат.
            +3
            Как говорил Док Браун: С пространственным воображением у тебя проблемы, Марти!
            Спасибо, все встало на места
            0
            Могу себе только представить, что груз толкнули вверх (строго по направлению нити), чтобы тот ударился о гвоздь.


            Согласен с Вами — условие не просто туманно, оно сформулировано так, что допускает разное толкование процесса. Я бы сформулировал задачу по другому.

            Мы с ученицей посчитали предложенный Вами вариант тривиальным и сразу его отмели :)

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

            P.S.: В статью надо вставить свой вариант условия.
              +1
              Да не то что хотят запутать — это прямой признак неграмотности составителей задачи.
              С озвученными условиями задача имеет бесконечное множество решений — ведь неизвестно, с какой силой изначально толкнули груз. Вдруг он развил первую космическую и подбил американский спутник-шпион??!
            +1
            ведь неизвестно, с какой силой изначально толкнули груз


            Дело в том, что попадание в гвозь возможно только при вполне определенной горизонтальной начальной скорости , а положение нити при котором она провисает вообще от её длины не зависит

            Единственное что авторам надо было указать — что груз толкают горизонтально
              0
              При том что задача решается арифметически в несколько строк (и в этом ее соль) может не стоило городить такое громоздкое численное решение?
                0
                Суть статьи не в том, что решается школьная задача. Её решение лишь илюстрация предлагаемой методики — метод избыточных координат для систем с освобождающими связями и сухим трением (статья о сухом трении скоро будет). Речь идет о применимости формализма Лагранжа к задачам подобного рода и более сложным. Пример выбран так, чтобы любой читатель мог при желании быстро проверить его не вдаваясь в громоздкие выкладки.

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

                Методика изложена в литературе (А. И. Огурцов. Аналитическая механика), но внимания к ней довольно мало
                0
                Что интересно, в системе Wolfram Mathematica 10 получение уравнений Лагранжа выглядит гораздо компактнее
                LagrangeEQs[T_, q_, r_, F_]:= Module[{ret=Range[Length[q]]},
                	For[i = 1, i <= Length[q], i++,
                		ret[[i]] = D[D[T, D[q[[i]], t]], t] - D[T, q[[i]]] == Sum[F[[k]].D[r[[k]], q[[i]]], {k, 1, Length[F]}]
                	];
                	ret
                ]
                
                  0
                  А всё потому что вольфрам умеет брать производную по функции и по производной от функции, а Maple до сих пор не понимает конструкций вида
                  diff(T, diff(q(t), t))
                  

                  и проиходится вводить подстановки

                Only users with full accounts can post comments. Log in, please.