Кинематическая интерпретация системы ДУ

Рассмотрим ДС, моделируемые конечным числом обыкновенных дифференциальных уравнений (ОДУ). Применительно к таким системам сохранились представления и терминология, первоначально возникшие в механике. В рассматриваемом случае для определения ДС необходимо указать объект, допускающий описание состояния заданием величин x₁, x₂, …, x_n в некоторый момент времени t = t₀. Величины x_i могут принимать произвольные значения, причем двум различным наборам величин x_i и x_i’ отвечают два разных состояния. Закон эволюции динамической системы во времени описывается системой ОДУ

Если рассматривать величины х₁, х₂,..., x_N как координаты точки x в N-мерном пространстве, то получается наглядное геометрическое представление состояния динамической системы в виде этой точки. Последнюю называют изображающей, а чаще фазовой точкой, а пространство состояний фазовым пространством динамической системы. Изменению состояния системы во времени отвечает движение фазовой точки вдоль некоторой линии, называемой фазовой траекторией. В фазовом пространстве системы уравнениями (1.5) определяется векторное поле скоростей, сопоставляющее каждой точке x выходящий из нее вектор скорости F(x), компоненты которого даются правыми частями уравнений (1.5):

Динамическая система (1.5) может быть записана в векторной форме:

где F(x) вектор-функция размерности N. Необходимо уточнить взаимосвязь понятий числа степеней свободы и размерности фазового пространства динамической системы. Под числом степеней свободы понимается наименьшее число независимых координат, необходимых для однозначного определения состояния системы. Под координатами первоначально понимались именно пространственные переменные, характеризующие взаиморасположение тел и объектов. В то же время для однозначного решения соответствующих уравнений движения необходимо, помимо координат, задать соответствующие начальные значения импульсов или скоростей. В связи с этим система с n степенями свободы характеризуется фазовым пространством размерности в два раза большей (Ν = 2n).

2.4. Эволюция ДС. Качественные особенности эволюции ДС проявляются в характере фазовых траекторий. Например, состоянию равновесия отвечает вырожденная траектория точка в ФП, периодическому движению замкнутая траектория. Траектория квазипериодического движения с т несоизмеримыми частотами w_i (т.е. такими, что не существует отличных от нуля целых чисел k_i, удовлетворяющих равенству ) сколь угодно близко проходит около любой точки m-мерного тора (всюду плотна на нём). Вообще, для стационарного режима (установившегося движения системы) характерны траектории, плотные в некотором подмножестве фазового пространства, а для переходного процесса траектории, не возвращающиеся в окрестность своих начальных точек.

2.5. Уравнения движения маятника

Рассмотрим движение маятника без трения (Higham, Higham 2005 (177-179)). Уравнение движения в этом случае имеет вид:

здесь θ – угол отклонения маятника от вертикали.

Определив y₁(t) = θ(t) и y₂(t) = dθ(t)/dt, можно ОДУ 2-го порядка записать как систему из двух ОДУ первого порядка:

Эта система уравнений может быть запрограммирована следующей m–функцией:

function yprime=pendd(t,y)

yprime=[y(2);-sin(y(1))];

Следующий m-файл дает возможность вычислить движение маятника на временном интервале 0 < t < 8 для трех различных начальных условий yazero=[1;1], ybzero=[-5;2], yczero=[5;-2].

Первая строка это название m-файла. Третий оператор снизу дает возможность поместить на график с помощью мыши название соответствующим графикам. Оператор xlim([-5 5]) ограничивает отображение интервала колебаний маятника.

%ppend

tspan=[0 8];

yazero=[1;1];

ybzero=[-5;2];

yczero=[5;-2];

[ta,ya]=ode45(@pendd,tspan,yazero);

[tb,yb]=ode45(@pendd,tspan,ybzero);

[tc,yc]=ode45(@pendd,tspan,yczero);

hold on

plot(ya(:,1),ya(:,2))

plot(yb(:,1),yb(:,2))

plot(yc(:,1),yc(:,2))

xlim([-5 5])

xlabel('угол'); ylabel('скорость');

gtext('колебание');gtext('вращение');

hold off

grid on

Результат вычислений можно видеть на рисунке 1.

Рис. 1.

Предположим, что величина сил трения, которые, в конечном счете, приводят к остановке маятника, пропорциональна скорости маятника (Hunt et al. 2008 (сс. 216-220)). Предположим также, что длина маятника равна 1 м, масса на конце маятника составляет 1 кг, а коэффициент трения примем равным 0.5. В таком случае, уравнения движения для маятника будут следующими.

(1)

где t представляет время в секундах, х обозначает угол отклонения маятника от вертикальной линии в радианах (то есть х = 0 – это исходное положение), y выражает угловую скорость маятника в радианах в секунду, а постоянная 9.81 – это приблизительное ускорение, вызванное действием силы тяжести, в метрах в секунду в квадрате.

На переменную g в Command Window (CW) запрограммирована система уравнений (1):

g = @(t,x)[x(2);-0.5*x(2)- 9.81*sin(x(1))];

Возьмем в качестве начального положения х(0) = 0 и начальной скоростью у(0) = 5 и введем в CW следующую последовательность операторов:

[t, xa] = ode45(g, [0:0.01:20], [0 5]);

plot(xa(:, 1), xa(:, 2))

grid on

На рис. 2 показан график х и у как функции от t на промежутке времени 0 < t <20 для этих начальных данных. (Чтобы использовать инструмент среды MATLAB для численного решения дифференциальных уравнений ode45, мы объединим х и у в один вектор х – смотрите оперативную справку для команды ode45.)

Начав в точке (0; 5), по мере увеличения t, мы следуем за кривой, в то время как она закручивается по часовой стрелке к точке (0; 0). Маятник колеблется назад и вперед, но с каждым колебанием угол отклонения становится меньше, пока маятник совсем не вернется в состояние покоя ко времени t = 20. Одновременно скорость также периодически изменяется, достигая своих наибольших значений в течение каждого колебания, когда маятник находится в середине своего колебания (угол близок к нулю), и достигает нуля, когда маятник находится в конце своего колебания.

Рис. 2.

Увеличим начальную скорость до 10.

[t, xa] = ode45(g, [0:0.01:20], [0 10]);

Рис. 3.

На этот раз (Рис. 3.) угол увеличивается до значения, превышающего 14 радиан, перед тем, как кривая подходит к точке примерно (12.5; 0). Если быть более точными, кривая закручивается к точке (4π; 0), потому что 4π радиан представляет то же самое положение маятника, что и 0 радиан. Маятник выполнил два полных оборота перед началом своих затухающих колебаний. Скорость сначала уменьшается, но затем повышается после того, как угол проходит через π, поскольку маятник проходит вертикальное положение и получает импульс. Импульса маятника будет достаточно только чтобы еще раз пройти через вертикальное положение в угле 3π.