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

 

Рассмотрим ДС, моделируемые конечным числом обыкновенных дифференциальных уравнений (ОДУ). Применительно к таким системам сохранились представления и терминология, первоначально возникшие в механике. В рассматриваемом случае для определения ДС необходимо указать объект, допускающий описание состояния заданием величин 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 степенями свободы характеризуется фазовым пространством размерности в два раза большей (Ν = 2 n).

 

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π.