Программа составления уравнений движения

Составим уравнения связей (4) и уравнения Чаплыгина (6) с помощью Mathematica

Текст программы

(* Определение вспомогательных функций *)

DS[T_, q_] := {Table[D[T, q[[i, 1]]], {i, Length[q]}]};

DST[T_, q_] := Transpose[DS[T, q]];

DV[X_, q_] := Table[D[X[[i, 1]], q[[j, 1]]], {i, Length[X]}, {j, Length[q]}];

XX[r_]:={{0,-r[[3,1]],r[[2,1]]},{r[[3,1]],0,-r[[1,1]]},{-r[[2,1]],r[[1,1]],0}};

ASymmQ[A_]:=Simplify[(A+Transpose[A])]==DiagonalMatrix[{0,0,0}];

ReverseXX[A_]:=If[ASymmQ[A],{{A[[3,2]]},{A[[1,3]]},{A[[2,1]]}},Print["Matrix is not ScrewMatrix"]];

XS[r_,p_]:=(Transpose[r].p)[[1,1]];

XV[r_,p_]:=XX[r].p;

OmegaMov[B_] := ReverseXX[Transpose[B].\!\( \*SubscriptBox[\(\[PartialD]\), \(t\)]B\)];

OmegaFix[B_] := ReverseXX[\!\( \*SubscriptBox[\(\[PartialD]\), \(t\)]B\).Transpose[B]];

(* Задание векторов обобщённых координат и обобщённых скоростей q1 , q2 , q1', q2', координаты z центра масс диска *)

q1 = { {ψ [ t ]} , {θ [ t ]} , {φ [ t ]} } ;

q2 = { {x [ t ] },{ y [ t ]} } ;

qv1 = D [ q1 , t ] ;

qv2 = D [ q2 , t ] ;

z = ρ Sin [ θ [ t ] ] ;

(* Введение радиус-вектора и скорости центра масс диска С в неподвижной системе координат OX*Y*Z* *)

rC = {{x[t]}, {y[t]}, {z}};

vC = D[rC, t];

(* Построение матрицы направляющих косинусов осей подвижной системы координат Cxyz относительно осей неподвижной системы координат OX*Y*Z* *)

B = RotationMatrix[ψ[t], {0, 0, 1}].RotationMatrix[θ [ t ], {1, 0, 0}].RotationMatrix[φ[t], {0, 0, 1}];

(*Определение вектора угловой скорости диска Ω в проекциях на оси неподвижной системы координат OX*Y*Z* *)

Ωf = OmegaFix[B] // Simplify;

(* Нахождение вектора угловой скорости диска Ω в проекциях на оси подвижной системы координат Сxyz (формулы Эйлера ) *)

Ωm = OmegaMov[B] // Simplify;

(* Построение вектора CP в проекциях на оси промежуточной системы координат Сx1, Cy1, Cz (Сх1 направлена по горизонтальному диаметру диска, Cz - перпендикулярно плоскости диска) *)

CP = { {0} , - { ρ} ,{0} } ;

(* Определение вектора CP в проекциях на оси неподвижной системы координат OX*Y*Z*)

CPf =RotationMatrix [ ψ [ t ] , { 0 , 0 , 1 } ] . RotationMatrix [ θ [ t ] , { 1 , 0 , 0 } ] . CP ;

(*Нахождение скорости точки P контакта диска с горизонтальной поверхностью (Vp = Vc + [Ω, CP]) с помощью теоремы о распределении скоростей в твердом теле *)

vP = D[rC, t] + XV[Ωf, CPf] // Simplify;

(* Введение уравнений неголономных связей, выражающих условия равенства нулю скорости точки P (отсутствие проскальзывания диска в точке Р ) *)

eq1 = vP[[1, 1]] == 0;

eq2 = vP[[2, 1]] == 0;

(* Разрешение уравнения неголономных связей относительно проекций скоростей центра масс диска *)

sub1 = Solve[{eq1, eq2}, Flatten[qv2]][[1]];

(* Нахождение матрицы A в матричной записи уравнений неголономных связей в виде qv2 = A.qv1*)

A = DV[qv2 /. sub1, qv1] ;

(* Задание тензора инерции диска относительно центра масс С как бесконечно тонкого однородного круглого диска *)

IC = m ρ^ 2 / 4 DiagonalMatrix[{1, 1, 2}];

(* Определение кинетической энергии диска по теореме Кенига как суммы кинетической энергии центра масс диска 1/2 m v²_C и кинетической энергии диска в его движении относительно центра масс С *)

T = 1/2 m XS[vC, vC] + 1/2 XS[ Ωm, IC.Ωm] // Simplify ;

(* Введение силовой функции силы тяжести диска *)

U = - m g z ;

(* Нахождение приведённой кинетической энергии диска *)

Θ = T /. sub1 // Simplify ;

(* Определение вектора обобщённого импульса *)

p2 = DST[T, qv2] /. sub1 // Simplify;

(* Построение левых частей уравнений Чаплыгина d(∂Θ/∂q'₁ )^T/dt - (∂Θ/∂q₁)^T - ( ∂U/∂q₁ )^T - [ d A^T/dt - ( ∂( Aq'₁ )/∂q₁ )^T ] p₂ = 0, *)

EqChaplygin = D[DST[Θ, qv1], t] - DST[Θ, q1] - DST[U, q1] - Transpose[D[A, t] - DV[A.qv1, q1]].p2 // Simplify;

(* Представление полученных уравнений Чаплыгина в традиционной форме записи *)

EqChaplygin == {{0}, {0}, {0}}// TraditionalForm

Результат исполнения программы , уравнения Чаплыгина (6)

В обычной математической нотации полученные уравнения после тождественных преобразований примут вид

( 7 + 5 cos 2θ ) ψ''+ 12 cos θ φ'' - 4 sin θ θ' (φ' + 5 cos θ ψ' ) = 0,

5 θ '' + 5 cos θ sin θ ψ' ² + 6 sin θ φ' ψ' + 4 g/ρ cos θ = 0, (10)

3 (cos θ ψ''+ φ'' ) - 5 sin θ θ' ψ' = 0