Лаб6 / 7408_MMOiSU_LR_6_Brig_3
.docxМИНОБРНАУКИ РОССИИ
Санкт-Петербургский государственный
электротехнический университет
«ЛЭТИ» им. В.И. Ульянова (Ленина)
Кафедра КСУ
отчет
по лабораторной работе №6
по дисциплине «Математическое моделирование объектов и
систем управления»
Тема: Исследование хаотических систем
Вариант 3
Студент гр. 7408 |
|
Лебедев Р. А. Трусова Е. С. |
Преподаватель |
|
Шпекторов А. Г. |
Санкт-Петербург
2021
Цель работы:
Ознакомиться с классом хаотических систем, освоить средства MATLAB для моделирования «странных аттракторов», изучить практические способы определения параметров хаотических систем.
Задание:
1. Написать программу решения системы уравнений в соответствии с вариантом;
2. Построить график фазовой траектории странного аттрактора;
3. Определить характеристики хаотической системы (определение старшего показателя Ляпунова).
Исходные данные:
Вариант №3. Аттрактор Лотки-Вольтерра
Определение старшего показателя Ляпунова хаотической системы
1 способ. Анализ уравнений хаотической системы.
1) Построим график фазовой
траектории странного аттрактора, на
котором выберем начальную точку
clc; clear;
a = 2.9851; b = 3; c = 2;
%Старший показатель Ляпунова
%Первый способ определения СПЛ
tspan = [0 1000]; % диапазон интегрирования
x0 = [1 1 1]; % начальная точка
[t, x] = ode45(@(t,x) odefcn(t,x,a,b,c), tspan, x0); % решение СНДУ
figure (1);
plot3 (x(:,1),x(:,2),x(:,3));
grid on; xlabel('x1'); ylabel('x2'); zlabel('x3');
Рис. 1 – Фазовая траектория странного аттрактора
2) Определим опорную траекторию
как решение системы
при выбранном начальном условии
x0 = x(end,:); % начальная точка опорной траектории на аттракторе
[t,x1] = ode45(@(t,x) odefcn(t,x,a,b,c),tspan,x0); % решение СНДУ
figure (2);
plot3 (x1(:,1),x1(:,2),x1(:,3)); grid on;
grid on; xlabel('x1'); ylabel('x2'); zlabel('x3');
Рис. 2 – Опорная траектория
3) Определим линеаризованную матрицу системы.
4) Исследуем собственные
числа матрицы А на интервале интегрирования
и построим график зависимости СПЛ от
времени. СПЛ определяется как максимальное
значение собственного числа матрицы
.
SPL1 = [];
for i = 1:1:length(t)
A = [1-x1(i,2)+2*c*x1(i,1)-2*a*x1(i,1)*x1(i,3) -x1(i,1) -a*x1(i,1)^2;
x1(i,2) -1+x1(i,1) 0;
2*a*x1(i,1)*x1(i,3) 0 -b+a*x1(i,1)^2]; % матрица линеаризации A
SPL1(i) = max (real(eig(A))); % вычисление СПЛ по 1-ому методу
end
figure (3);
plot(t, SPL1); xlim ([0 50]);grid on;xlabel 't'; ylabel 'СПЛ';
Рис. 3 – Зависимость СПЛ от времени (1 способ)
2 способ. Анализ решений системы уравнений.
1) Зададим два начальных условия на поверхности странного аттрактора, разлиающеся на малую величину ε.
%Второй способ определения СПЛ
e1 = [0.0001; 0; 0]; % отклонения по координате
e2 = [0; 0.0001; 0]; % - // -
e3 = [0; 0; 0.0001]; % - // -
% x0op = x1(end,:); % первая начальная точка
spl2 = []; % для определения среднего значения по выборке
SPL2 = []; % СПЛ для различных моментов времени
T = [];
2) Найдем решения для выбранных
условий и определим процесс рассогласования
траекторий
и
рассчитаем сигнал вида
.
Определим СПЛ как максимальный коэффициент
наклона сигнала
m=1;
for i = 1:10:length(t)/100
x0op = x1(i,:);
for j = 1:1:3 % вычислние 3-х СПЛ для вычисления среднего значения
x0 = [x0op(1)+ e1(j) x0op(2)+ e2(j) x0op(3)+ e3(j)]; % вторая начальная точка
e = x0 - x0op; % вектор разности координат начальных точек
[t1,xop] = ode45(@(t1,x) odefcn(t1,x,a,b,c), tspan, x0op);
[t1,x] = ode45(@(t1,x) odefcn(t1,x,a,b,c), t1, x0);
deltax = x - xop; % вектор рассогласования траекторий
z = [];
for k = 2:1:length(t1) % вычисление Функции z
z(k) = log(norm(deltax(k,:))/norm(e));
end
p = polyfit(t1, z', 10); % аппроксимация функции z полиномом p
f = polyval(p, t1); % вычисление значений полинома
dp = polyder(p); % вычисление производной полинома p
spl2(j) = max(dp); % вычисление СПЛ
end
i
SPL2(m) = mean(spl2);
T(m) = t(i);
m = m+1;
end
figure (4)
plot (T,SPL2); xlim([0 10]); grid on; xlabel 't'; ylabel 'СПЛ';
Рис. 4 – Зависимость СПЛ от времени (2 способ)
Вывод.
В ходе работы были изучены два способа нахождения старшего показателя Ляпунова хаотической системы.
Для обоих способов были получены графики зависимости СПЛ от времени для разных начальных условий. На графике зависимости СПЛ от времени, полученном в результате вычисления первым способом, имеются отрицательные интервалы. Для второго способа на всей траектории наблюдаются положительные значения СПЛ, что может свидетельствовать о хаотичности системы.
Стоит отметить, что вычисление СПЛ вторым способом в пакете MatLab заняло на порядок больше времени.
ПРИЛОЖЕНИЕ
Функция, вычисляющая правую часть системы.
%функция вычисляющая правую часть
function dxdt = odefcn(t,x,a,b,c)
dxdt = zeros (3,1);
dxdt(1) = x(1)-x(1)*x(2)+c*x(1)^2-a*(x(1)^2)*x(3);
dxdt(2) = -x(2)+x(1)*x(2);
dxdt(3) = -b*x(3)+a*(x(1)^2)*x(3);
end
Код программы.
clc; clear;
a = 2.9851; b = 3; c = 2;
%Старший показатель Ляпунова
%Первый способ определения СПЛ
tspan = [0 300]; % диапазон интегрирования
x0 = [1 1 1]; % начальная точка
[t, x] = ode45(@(t,x) odefcn(t,x,a,b,c), tspan, x0); % решение СНДУ
figure (1);
plot3 (x(:,1),x(:,2),x(:,3));
grid on; xlabel('x1'); ylabel('x2'); zlabel('x3');
x0 = x(end,:); % начальная точка опорной траектории на аттракторе
[t,x1] = ode45(@(t,x) odefcn(t,x,a,b,c),tspan,x0); % решение СНДУ
figure (2);
plot3 (x1(:,1),x1(:,2),x1(:,3)); grid on;
grid on; xlabel('x1'); ylabel('x2'); zlabel('x3');
SPL1 = [];
for i = 1:1:length(t)
A = [1-x1(i,2)+2*c*x1(i,1)-2*a*x1(i,1)*x1(i,3) -x1(i,1) -a*x1(i,1)^2;
x1(i,2) -1+x1(i,1) 0;
2*a*x1(i,1)*x1(i,3) 0 -b+a*x1(i,1)^2]; % матрица линеаризации A
SPL1(i) = max (real(eig(A))); % вычисление СПЛ по 1-ому методу
end
figure (3);
plot(t, SPL1);grid on;xlabel 't'; ylabel 'СПЛ';
%Второй способ определения СПЛ
e1 = [0.00001; 0; 0]; % отклонения по координате
e2 = [0; 0.00001; 0]; % - // -
e3 = [0; 0; 0.00001]; % - // -
% x0op = x1(end,:); % первая начальная точка
spl2 = []; % для определения среднего значения по выборке
SPL2 = []; % СПЛ для различных моментов времени
T = [];
m=1;
for i = 1:10:length(t)
x0op = x1(i,:);
for j = 1:1:3 % вычислние 3-х СПЛ для вычисления среднего значения
x0 = [x0op(1)+ e1(j) x0op(2)+ e2(j) x0op(3)+ e3(j)]; % вторая начальная точка
e = x0 - x0op; % вектор разности координат начальных точек
[t1,xop] = ode45(@(t1,x) odefcn(t1,x,a,b,c), tspan, x0op);
[t1,x] = ode45(@(t1,x) odefcn(t1,x,a,b,c), t1, x0);
deltax = x - xop; % вектор рассогласования траекторий
z = [];
for k = 1:1:length(t1) % вычисление Функции z
z(k) = log(norm(deltax(k,:))/norm(e));
end
p = polyfit(t1, z', 23); % аппроксимация функции z полиномом p
f = polyval(p, t1); % вычисление значений полинома
dp = polyder(p); % вычисление производной полинома p
spl2(j) = max(dp); % вычисление СПЛ
end
i
SPL2(m) = mean(spl2);
T(m) = t(i);
m = m+1;
end
figure (4)
plot (T,SPL2); grid on; xlabel 't'; ylabel 'СПЛ';