Численное интегрирование уравнений движения
Функции x, y, ψ, θ, φ, ψ', θ', φ' могут быть найдены путём численного интегрирования полученных уравнений связей (4) и уравнений динамики (10) с помощью представленной далее программы. При использовании NDSolve, стандартной функции Mathematica, интегрируемую систему уравнений можно не приводить к нормальной форме Коши. Предполагается, что предварительно выполненапрограмма составления уравнений. Возможен также независимый ввод уравнений (4),(10) в явном виде внутри программы.
Числовые значения параметров задачи заданы в единицах СИ.
Программа численного интегрирования уравнений движения диска
(*Преобразование уравнения Чаплыгина к одномерному списку (уравнения получены ранее)*)
EqChaplygin1 = (# == 0) & /@ (EqChaplygin[[All, 1]]/m ) // Simplify;
(*Построение системы уравнений кинематических связей (4) (уравнения получены ранее) и уравнений Чаплыгина (10) *)
Eq = Join[{eq1}, {eq2}, EqChaplygin1];
(*Задание начальных условий движения, интервала интегрирования tk , числовых значений массы m и радиуса ρ диска *)
InitialConditions = { x[0] == 0,y[0] == 0,ψ[0] == 0,θ[0] == 1.918,φ[0] ==0,ψ'[0] == 2.569,θ'[0] == 3.0743, φ'[0] == 0.1229};
tk = 25;
Num = {m -> 0.8, ρ -> 0.3, g -> 9.8};
(*Нахождение численного решения уравнений движения диска *)
sol = NDSolve[Join[Eq, InitialConditions] /. Num, {x,y,ψ,θ,φ,ψ',θ',φ'}, {t, 0, tk}, MaxSteps -> 9000000];
(*Построение графиков обобщённых координат и скоростей диска*)
VarsPlots = Table[Plot[Evaluate[f[t]] /. sol, {t, 0, tk}, LabelStyle -> FontSize -> 14, PlotStyle -> Thick,
PlotRange -> All, PlotLabel -> Style[f[t], 16]], {f, {x,y,ψ,θ,φ,ψ',θ',φ'}}]
(*Построение траектории центра масс диска *)
CenterMassTrajectory = ParametricPlot[{x[t], y[t]} /. sol, {t, 0, tk}, PlotRange -> All, PlotStyle -> Thick, AspectRatio -> 1, LabelStyle -> FontSize -> 14, AxesLabel -> {"X", "Y"}]
(*Построение графика проекции на ось OZ* нормальной реакции плоскости в точке контакта P диска и плоскости (7) *)
NormalReaction = Plot[m (ρ Cos[θ[t]] θ''[t] - ρ Sin[θ[t]] θ'[t]^2 + g) /. sol /. Num, {t, 0., tk}, PlotRange -> All,
AxesLabel -> {" t (s) ", "Rn (N)"}, LabelStyle -> FontSize -> 14]
(*Построение графика проекции на ось OZ* нормальной реакции плоскости в точке контакта P диска и плоскости (7) вблизи области нулевых значений реакции*)
NormalReaction0 = Plot[m (ρ Cos[θ[t]] θ''[t] - ρ Sin[θ[t]] θ'[t]^2 + g) /. sol /. Num, {t, 0., tk},
PlotRange -> {{0, tk/5}, {-5, 20}}, AxesLabel -> {" t (s) ", "Rn (N)"}, LabelStyle -> FontSize -> 14]
Рис. 2. Графики обобщённых координат и скоростей диска
Рис. 3. Траектория центра масс диска
Рис. 4. График проекции на вертикаль реакции плоскости в точке контакта диска и плоскости
Рис. 5. График проекции на вертикаль реакции плоскости в точке контакта диска и плоскости вблизи области нулевых значений проекции реакции