 |
Нахождение начальной скорости, которая заставляет маятник
|
|
|
|
Совершать полные вращения
Предположим, что мы хотим с точностью до 0.001 найти наименьшее значение начальной скорости, которая требуется, чтобы заставить маятник, начинающий движение из своего исходного положения, выполнить полное вращение один раз. Будет полезно показать решения, которые соответствуют нескольким различным начальным скоростям на одном графике.
Сначала мы рассмотрим целые значения скорости в промежутке от 5 до 10.
>> hold on
>> for а = 5:10
[t, xa] = ode45(g, [0:0.01:20], [0 а]);
plot(xa(:, 1), xa(:, 2))
end
>> hold off
Рис. 4.
Начальные скорости 5, 6 и 7 не являются достаточно большими (рис. 4), чтобы увеличить угол более π, но начальные скорости 8, 9 и 10 достаточны, чтобы заставить маятник совершать полный оборот. Посмотрим, что происходит на промежутке между 7 и 8.
Второй шаг
>> hold on
>> for а = 7.0:0.2:8.0
[t, ха] = ode45(g, [0:0.01:20], [0 а]);
plot(xa(:, 1), xa(:, 2))
end
>> hold off
Рис. 5.
Можно заметить (рис. 6), что ответ находится где-то между 7.2 и 7.4. Давайте выполним еще одно уточнение.
Третий шаг
>> hold on
>> for а = 7.2:0.05:7.4
[t, ха] = ode45(g, [0:0.01:20], [0 а]);
plot(xa(:, 1), ха(:, 2))
end
» hold off
Рис. 6.
Следует сделать вывод, что наименьшая необходимая скорость с точностью 0,01 находится где-то между 7.25 и 7.3 (рис. 7 и 8).
Четвертый шаг
for a = 7.25:0.01:7.3
[t, xa] = ode45(g, [0:0.01:20], [0 a]);
plot(xa(:, 1), xa(:, 2))
end
Рис. 7.
Для более точного анализа увеличим область графика, где происходит смена режима колебания.
Рис. 8.
Видно, что наименьшая необходимая скорость находится где-то между 7.29 и 7.3.
Следует продолжить нахождение более точного значения скорости смены режима колебания.
Динамика осциллятора Ван дер Поля при w2 = 2 и c = 1
|
|
|
Предельный цикл – устойчивый режим периодических колебаний в нелинейных системах после завершения переходных процессов
Файл-функция, для этих параметров, имеет следующий вид (см. Андреев 2013, с. 117):
function dydt = vdp1(t,y)
dydt = zeros(2,1); dydt(1) = y(2);
dydt(2) = 1*(1-y(1).^2).*y(2)-2*y(1);
При начальном положении на предельном цикле (х0 = 2, v0= 0) вызов файл-функции имеет следующий вид:
[t,y] = ode23(@vdp1,[0 25],[2;0]);
Следующий оператор дает возможность получить зависимость х и v от времени
plot(t,y(:,1),t,y(:,2)), grid on
Результат показан на рис. 1.
Рис. 1.
Следующий оператор дает возможность получить фазовый портрет системы (рис. 2):
plot(y(:,1),y(:,2)), grid on
Рис. 2.
Если взять начальные данные вне цикла (х0 = -0.5, v0= 5) вызов файл-функции имеет следующий вид:
[t,y] = ode23(@vdp1,[0 25],[-0.5;5]);
Результат показан на рис. 3.
Рис. 3.
Фазовый портрет показан на рис. 4.
Рис. 4.
Если взять начальные данные внутри цикла (х0 = -0.05, v0= -0.05) вызов файл-функции имеет следующий вид:
[t,y] = ode23(@vdp1,[0 50],[-0.05; -0.05]);
Результат показан на рис. 5.
Рис. 5.
Фазовый портрет показан на рис. 6.
Рис. 6.
Система Ван дер Поля с периодическим возмущением
Файл-функция для этой системы при значениях A 0 = 2 и W = 2 имеет следующий вид:
function dydt = vdp1s(t,y)
dydt = zeros(2,1); dydt(1) = y(2);
dydt(2) = 1*(1-y(1).^2).*y(2)-2*y(1)+2*sin(2*t);
Фазовый портрет показан на рис. 7.
Рис. 7.
Лабораторная работа № 6 – 7
Качественный анализ линейных ДС
ЦЕЛЬ РАБОТЫ
Ознакомление с методами качественного анализа систем ОДУ
Рассмотрим линейную однородную систему с постоянными коэффициентами, которая представляет собой линейную динамическую систему (ЛДС):
| (1) |
Координатную плоскость xOy называют ее фазовой плоскостью (ФП). Через любую точку плоскости проходит одна и только одна фазовая кривая (траектория) (ФТ). В системе (1) возможны три типа ФТ:
|
|
|
· точка,
· замкнутая кривая,
· незамкнутая кривая.
Точка на ФП соответствует стационарному решению (положению равновесия, точке покоя) системы (1), замкнутая кривая периодическому решению, а незамкнутая непериодическому.
Положения равновесия ДС
Положения равновесия системы (1) найдем, решая систему:
| (2) |
Система (1) имеет единственное нулевое положение равновесия, если определитель матрицы системы:
|
Если же det A = 0, то, кроме нулевого положения равновесия, есть и другие, так как в этом случае система (2) имеет бесконечное множество решений.
Качественное поведение ФТ (тип положения равновесия) определяется собственными числами матрицы системы.
Классификация точек покоя
Собственные числа матрицы системы найдем, решая уравнение:
| (3) |
Заметим, что a + d = tr A (след матрицы) и ad – bc = det A.
В СКМ MATLAB определять собственные числа и собственные вектора системы ОДУ удобно с помощью функции eig(A).
Собственные числа l i и собственные векторы ui ≠ 0 квадратной матрицы А удовлетворяют равенствам А ui = l i ui. Функция eig, вызванная с входным аргументом матрицей, находит все собственные числа матрицы и записывает их в выходной аргумент – вектор:
>> A = [2. 3; 3 5];
>> lam = eig(A)
lam =
0.1459
6.8541
>>
Для одновременного вычисления всех собственных векторов и чисел следует вызвать eig с двумя выходными аргументами.
>> [U, Lam] = eig(A);
Первый выходной аргумент и представляет собой матрицу, столбцы которой являются собственными векторами. Для доступа, например, к первому собственному вектору следует использовать индексацию при помощи двоеточия
>> u1=U(:,1)
u1 =
-0.8507
0.5257
Вторым выходным аргументом Lam возвращается диагональная матрица, содержащая собственные числа исходной матрицы.
>> Lam
Lam =
0.1459 0
0 6.8541
Проверьте, правильно ли найдены, например, второе собственное число и соответствующий ему собственный вектор. Воспользуйтесь определением:
>> A*U(:, 2) - Lam(2, 2)*U(:, 2)
ans =
1.0e-015 *
0.4441
-0.8882
Классификация точек покоя в случае, когда det A ≠ 0, приведена в таблице:
|
|
|
Устойчивость точек покоя
Собственные значения матрицы системы (1) однозначно определяют характер устойчивости положений равновесия:
Фазовые портреты
Бесконечное множество точек покоя
Если det A = 0, то система (1) имеет бесконечное множество положений равновесия. При этом возможны три случая:
Во втором случае любая точка покоя устойчива по Ляпунову. В первом же случае только, если l2 < 0.
Направление на фазовой кривой указывает направление движения фазовой точки по кривой при возрастании t.
|
|
|
