Построение бифуркационной картины аттрактора Рёсслера

[optimden]

Здравствуйте. Для изучения влияния параметров уравнений аттрактора хочу построить бифуркационные картины, чтобы посмотреть, при каких значениях параметров какие состояния системы возникают. То есть по сути фиксирую все параметры, и изменяю только один. При этом уравнения системы зависят от времени, что я учитываю во внутреннем цикле за счёт переменной $dt$. Для каждого значения изменяемой переменной я заполняю массив значениями $zd$ для каждой итерации:

Используется синтаксис Matlab M

N = 1000;
xd=zeros(N,1);
yd=zeros(N,1);
zd=zeros(N,1);
xd(1) = 5;
yd(1) = 10;
zd(1) = 20;
a = 0.2;
c = 5.7;
dt = 0.01;
alpha = 6;
hold on;
for p=0:0.01:10
    for i = 2:N
        xd(i,1) = xd(i-1)-dt*(alpha*(yd(i-1)+zd(i-1)));
        yd(i,1) = yd(i-1)+dt*(alpha*(xd(i-1)+a*yd(i-1)));
        zd(i,1) = zd(i-1)+dt*(alpha*(p+zd(i-1)*(xd(i-1)-c)));
    end
zd(1:N) = [];
t(:,1)=p;
plot(t,zd,'k.','MarkerSize',3);
end

hold off;

Однако в итоге матлаб выдаёт пустой график. Подскажите, пожалуйста, в чём может быть ошибка?

[Pphantom]

Вы бы определились, являются ли у вас массивы xd, yd и zd одномерными или двумерными. А то вы меняете это как попало, в итоге скрипт просто падает с ошибкой и до рисования графиков дело не доходит.

-- 02.12.2019, 21:55 --

А, да, удалять данные перед построением по ним графика - это тоже очень, очень хорошая идея.

-- 02.12.2019, 21:57 --

Собрался писать про третью ошибку, после чего подумал, что надо остановиться. Удалите этот ужас и сделайте все заново. Настоятельно рекомендую писать комментарии к коду с пояснениями, что именно вы надеялись сделать в каждой строке.

[optimden]

Pphantom в сообщении #1428638 писал(а):

Вы бы определились, являются ли у вас массивы xd, yd и zd одномерными или двумерными. А то вы меняете это как попало, в итоге скрипт просто падает с ошибкой и до рисования графиков дело не доходит.

-- 02.12.2019, 21:55 --

А, да, удалять данные перед построением по ним графика - это тоже очень, очень хорошая идея.

-- 02.12.2019, 21:57 --

Собрался писать про третью ошибку, после чего подумал, что надо остановиться. Удалите этот ужас и сделайте все заново. Настоятельно рекомендую писать комментарии к коду с пояснениями, что именно вы надеялись сделать в каждой строке.

Да, вы правы. Я этот код писал на основании кода, найденного в интернете. Он был для более простой функции, но график строился нормально. И перед построением данных стирались первые $N-120$ значений переменной. Так и не смог понять, почему именно так.
Написал свой код:

Используется синтаксис Matlab M

%Бифуркационная диаграмма для аттрактора Ресслера
%для p
%Задаем параметры
N = 100; %количество итераций по переменным x, y, z
a = 0.2;
c = 5.7;
dt = 0.01;
alpha = 6;
pstep = 0.01; %шаг по переменной p
pmax = 10; %максимальное значение переменной p
Np = 1000; %количество итераций по переменной p
% задание пустых массивов
xd = [];
yd = [];
zd = [];
p = [];
% задание начальных условий для переменных
p(1:N,1)=0;
xd(1,1:Np) = 5;
yd(1,1:Np) = 10;
zd(1,1:Np) = 20;

for j = 2:1000
    for i = 2:N
        xd(i,j-1) = xd(i-1,j-1)-dt*(alpha*(yd(i-1,j-1)+zd(i-1,j-1)));
        yd(i,j-1) = yd(i-1,j-1)+dt*(alpha*(xd(i-1,j-1)+a*yd(i-1,j-1)));
        zd(i,j-1) = zd(i-1,j-1)+dt*(alpha*(p(i-1,j-1)+zd(i-1,j-1)*(xd(i-1,j-1)-c)));
    end
p(:,j) = p(:,j-1)+pstep; % изменение переменной p на величину шага по этой переменной
end
plot(p,xd,'k.','MarkerSize',3); %построение графика
xlabel('p'); ylabel('x'); %подпись осей

График строится, но он сильно отличается от той бифуркационной диаграммы, которая приведена для параметра $p$ аттрактора Рёсслера в Википедии.
То, что получилось у меня:

То, что должно получиться в теории:

[Pphantom]

А почему, собственно, вы решили, что строите бифуркационную диаграмму? По-видимому, вам надо сначала посмотреть, что это такое.

[optimden]

Pphantom в сообщении #1428666 писал(а):

А почему, собственно, вы решили, что строите бифуркационную диаграмму? По-видимому, вам надо сначала посмотреть, что это такое.

Насколько я понимаю, бифуркационная диаграмма отображает возможные состояния динамической системы при разных значениях бифуркационного параметра. Рассматриваемая мной система меняется со временем. Соответственно для каждого значения бифуркационного параметра я рассматриваю все возможные состояния системы, которые она может приобретать в разные моменты времени. При критических значениях бифуркационного параметра будут точки притяжения и система будет устойчива, не будет изменяться со временем.
Я правильно понимаю?
Также прочитал, что некоторое количество первых итераций по переменной пропускается, а бифуркационная диаграмма строится по результатам последующих итераций. Не очень понятно, для чего необходим этот пропуск.

[Pphantom]

optimden в сообщении #1428669 писал(а):

Я правильно понимаю?

Нет, неправильно. Посмотрите в википедии, там в статье, из которой вы брали картинку, есть ссылка.

[optimden]

Pphantom в сообщении #1428672 писал(а):

optimden в сообщении #1428669 писал(а):

Я правильно понимаю?

Нет, неправильно. Посмотрите в википедии, там в статье, из которой вы брали картинку, есть ссылка.

Вроде более менее разобрался. Получается, бифуркационная диаграмма состоит из точек, для которых выполняется условие:
$x_n=f\left(x_{n-1},c\right)$
и
$\left\lvert\dot{f}\left(x\right)\right\rvert<1$
где $c$ - бифуркационный параметр.
Тот код, на который я ориентировался изначально, предполагал разбиение диапазона изменения переменной, допустим $x$ на некоторое число значений. Затем для каждого из значений в цикле вычисляется $f\left(x_{n-1},c\right)$. Полученные значения образуют орбиту точки, т.е. $x, f\left(x\right), f\left(f\left(x\right)\right)$. Предполагается, что если взять большое число итераций, и откинуть большую часть значений, начиная с первой итерации, то последние значения будут соответствовать стабильной периодической орбите.
Иначе надо вводить условие для модуля производной.
Теперь я в верном направлении двигаюсь?

[Pphantom]

optimden в сообщении #1428675 писал(а):

Теперь я в верном направлении двигаюсь?

Примерно да, с поправкой на то, что бифуркационные диаграммы бывают разными. Только подобными методами строить бифуркационную диаграмму вы будете до второго пришествия, поэтому, если она все же нужна, стоит подумать об использовании более подходящих инструментов (хотя бы даже в рамках MatLAB'а).

[optimden]

Pphantom в сообщении #1428676 писал(а):

optimden в сообщении #1428675 писал(а):

Теперь я в верном направлении двигаюсь?

Примерно да, с поправкой на то, что бифуркационные диаграммы бывают разными. Только подобными методами строить бифуркационную диаграмму вы будете до второго пришествия, поэтому, если она все же нужна, стоит подумать об использовании более подходящих инструментов (хотя бы даже в рамках MatLAB'а).

Например? Я просто никогда раньше не работал в матлабе с динамическими системами.

[Pphantom]

Там есть встроенные интеграторы. Но если вы с ними не работали и MatLAB'ом толком пользоваться не умеете, может быть, проще взять инструмент, с которым опыт работы есть?

[optimden]

Pphantom в сообщении #1428702 писал(а):

Там есть встроенные интеграторы. Но если вы с ними не работали и MatLAB'ом толком пользоваться не умеете, может быть, проще взять инструмент, с которым опыт работы есть?

Mathcad, LabView, Wolfram Mathematica. Но ни в одном из данных программных пакетов я не работал с динамическими системами. Какой в целом, более быстрый и простой подход вы предлагаете?

[Pphantom]

Кхм... задался вопросом, что же вы все же знаете, полез смотреть ваши старые темы. Нашел как минимум две, где обсуждался ode45 из MatLAB'а, окончательно перестал понимать, зачем надо было организовывать закат солнца вручную.

В общем, изложите-ка для начала мотивировку собственной деятельности. В то, что решение этой задачи вам необходимо для какой-то конкретной работы, верится очень слабо, а во всех прочих случаях надо бы сначала освоить основы программирования и вычислительной математики.

[optimden]

Pphantom в сообщении #1428704 писал(а):

Кхм... задался вопросом, что же вы все же знаете, полез смотреть ваши старые темы. Нашел как минимум две, где обсуждался ode45 из MatLAB'а, окончательно перестал понимать, зачем надо было организовывать закат солнца вручную.

В общем, изложите-ка для начала мотивировку собственной деятельности. В то, что решение этой задачи вам необходимо для какой-то конкретной работы, верится очень слабо, а во всех прочих случаях надо бы сначала освоить основы программирования и вычислительной математики.

Надо исследовать влияние параметров уравнений и коэффициентов связи на фазовые портреты связанной системы аттракторов Рёсслера и Лоренца.
Я решил сделать это, построив бифуркационные диаграммы, чтобы посмотреть, при каких значениях параметров происходит изменение состояния системы (становится хаотичной или периодической).

Ode45 не стал использовать просто потому, что наткнулся на способ, описанный выше, где откидывается большое количество начальных значений орбиты. Для простой функции он работал.

Хотя наверное проверить каждую точку в цикле на выполнение условия $\left\lvert{\dot{f}\left(x\right)\right\rvert}<1$ не сильно сложнее.

[Pphantom]

optimden в сообщении #1428705 писал(а):

Надо исследовать

Что значит "надо"? Исследовать что-то при отсутствии знакомства с предметом исследования не надо.

[optimden]

Pphantom в сообщении #1428706 писал(а):

optimden в сообщении #1428705 писал(а):

Надо исследовать

Что значит "надо"? Исследовать что-то при отсутствии знакомства с предметом исследования не надо.

Надо, потому что это задание на курсовую. А исследование чего-то, видимо, предполагает и знакомство с предметом исследования.