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

Функции 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. График проекции на вертикаль реакции плоскости в точке контакта диска и плоскости вблизи области нулевых значений проекции реакции