устойчивый предельный цикл в системе О.Д.У.

[LynxGAV]

Для определенных значений параметров $\beta_0, \beta_1, \alpha, \gamma, \mu$ cистема трех О.Д.У. первого порядка

$\left\{ \begin{array}{l} s'(t) = \mu (1-s(t))-\beta_0(1+\beta_1\cos 2\pi t)s(t)i(t),\\ e'(t)=\beta_0(1+\beta_1\cos 2\pi t)s(t)i(t)-(\alpha+\mu)e(t),\\ i'(t)=\alpha e(t)-(\gamma+\mu)i(t). \end{array} \right.$

имеет устойчивый предельный цикл $(\bar s(t), \bar e(t), \bar i(t))$ периода $T=t=1$. Допустим в момент $t=0$ начальные условия $(s(0), e(0), i(0))$ выбраны таким образом, что они принадлежат циклу. Численное интегрирование в течение одного периода (на интервале $t \in [0,1]$) дает замкнутый цикл. Интегрируя в течение одного периода для отрицательных времен (на интервале $t \in [0,-1]$), решение начинает расходится для $t\sim-0.2$. Я предполагала, что проинтегрировав в обратку, тоже получу цикл.. Подскажите, пожалуйста, в чем проблема.

[Gafield]

В обратную сторону цикл может быть неустойчивым, типа траектории от него уже расходятся.

[LynxGAV]

Хмм. то же самое произойдет с уравнением $y'(t)=\cos t$?

[ewert]

Ну тут же цикл не предельный, а просто все решения периодические. Если же некий цикл устойчив на плюс бесконечности, то он, скорее всего, будет неустойчив на минус бесконечности или даже просто в противоположную сторону (как у Вас и наблюдается). Иначе было бы даже как-то странно.

[LynxGAV]

Допустим, чтобы исследовать на устойчивость написанную выше систему нелинейных уравнений, я пользуюсь линейными вариационными уравнениями и смотрю на устойчивость их нулевого решения. В данном случае, их решения в общем виде записываются как $c_1e^{\mu_1 t} \vec p_1(t)+c_2e^{\mu_2 t} \vec p_2(t)+c_3e^{\mu_3 t} \vec p_3(t)$, где $c_1, c_2, c_3$-константы, определяемые из начальных условий, $\vec p_1(t), \vec p_2(t), \vec p_3(t)$ - периодические вектор-функции, $\mu_1, \mu_2, \mu_3$ - экспоненты Флоке. Для интересующих меня параметров, у этих экспонент отрицательная действительная часть, потому решения асимптотически устойчивы. Формально отсюда и видно, что в противоположную сторону решения будут не устойчивы и расти экспоненциально. Типа того?

[Gafield]

Да. Можно еще привести дискретный вариант рассуждений. Если в одну сторону отображение (сдвиг за период) сжимающее в окрестности какой-то точки (лежащей на предельном цикле), то в обратную сторону оно уже будет "расжимающим"

К тому же если показатели отрицательны, то устойчивость (линеаризации, по крайней мере) будет асимптотической. В обратную сторону должна быть неустойчивость. Пример с косинусом сюда не подходит.

[LynxGAV]

Gafield в сообщении #276261 писал(а):

Да. Можно еще привести дискретный вариант рассуждений. Если в одну сторону отображение (сдвиг за период) сжимающее в окрестности какой-то точки (лежащей на предельном цикле), то в обратную сторону оно уже будет "расжимающим"

Это вы наверное имеете ввиду отображение Пуанкаре, да? Ну, если да, то до меня дошло .

[Gafield]

Да.

[LynxGAV]

всем спасибо

[LynxGAV]

У меня возник еще один сопутствующий вопрос по этой теме.

Как я и сказала, чтобы исследовать на устойчивость написанную выше систему нелинейных уравнений, я пользуюсь линейными вариационными уравнениями, которые в общем виде можно записать в следующей форме $x'(t)=A(t) x$ (*), где $x(t)$ - вектор-функция, а $A(t)=A(t+T)$ - периодическая 3х3 матрица, $T$ - период. Для линейных уравнений с периодическими коэффициентами существует теорема Флоке, согласно которой фундаментальная матрица решений данной системы записывается как $X(t)=P(t) e^ {t R}$. Здесь $P(t)=P(t+T)$ - периодическая 3х3 матрица, а $R=\frac{1}{T}\log X(T)$, если $X(0)=E$ единичная матрица. Для данного выбора начальных условий собственные значения матрицы $X(T)$, $\rho_i$, называются характерными мультипликаторами уравнения (*), а характерные показатели $\mu_i$ определяются как $\rho_i=e^{\mu_i T}$.

Моя задача заключается в том, что надо найти матрицу $P(t)=X(t) e^{-t R}$ (и естественно мультипликаторы Флоке), что для меня оказалось непосильным трудом. Два мультипликатора Флоке комплексно-сопряженные, а третий действительный (все по модулю меньше единицы). Третий мультипликатор настолько маленький, что точности double не хватает, чтобы его посчитать. Технически это все выражается в том, что $\det X(T)=\det e^{TR}$ практически ноль и матрица $P(t)$ для значений $0.8<t<T=1$ предствляет шум (только один из ее элементов периодический). Пробовала вычислять $P(t)$ решая дифференциальное уравнение, а также диагонализируя $R$. Короче ничего дельного из этого не вышло по той же самой причине, что теоретически третий мультипликатор Флоке хоть и не ноль, узнать его точное значение не представляется возможным. Идея только одна и абстрактная, cделать замену координат от $x(t)$ к чему-нибудь, чтобы свести трехмерную систему к двухмерной, где третий мультипликатор не будет играть никакой роли. Это можно сделать? Он-то и в трех ее не играет, потому что соответствующий показатель безумно большой отрицательный и система супер сжимается в этом направлении, но вот расчеты сделать невозможно, а без них никак.

[Gafield]

Правая часть $А(t)$ задается аналитически? А то в этом случае математические пакеты могут считать и со ста знаками.

Для сведения к двумерной системе надо знать собственное подпространство для двух других показателей. Типа матрица после линейного преобразования будет $\mathrm{diag}\,(\lambda,A)$, где $A$ - блок два на два. Весь вопрос, с какой точностью. Возможно, в задании подпространства нужна точность порядка обратной к этому большому множителю растяжения.

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

[LynxGAV]

Gafield в сообщении #284736 писал(а):

Правая часть $А(t)$ задается аналитически? А то в этом случае математические пакеты могут считать и со ста знаками.

Исходная система (1):
$\left\{ \begin{array}{l} y_1'(t) = \mu (1-y_1(t))-\beta_0(1+\beta_1\cos 2\pi t)y_1(t)y_3(t),\\ y_2'(t)=\beta_0(1+\beta_1\cos 2\pi t)y_1(t)y_3(t)-(\alpha+\mu)y_2(t),\\ y_3'(t)=\alpha y_2(t)-(\gamma+\mu)y_3(t). \end{array} \right.$

Уравнения линериазации $x'(t)=A(t)x$ (2), где $A=\left ( \begin{array}{ccc} -\mu - \beta_0 (1+\beta_1 \cos 2\pi t) \bar{y_3}(t) & 0 & - \beta_0 (1+\beta_1 \cos 2\pi t) \bar{y_1}(t) \\ \beta_0 (1+\beta_1 \cos 2\pi t) \bar{y_3}(t) & - (\alpha + \mu) & \beta_0 (1+\beta_1 \cos 2\pi t) \bar{y_1}(t) \\ 0 & \alpha & -(\gamma+\mu) \end{array}\right)$. Здесь $\bar{y}(t)=y(t)-x(t)$ - перидиоческая вектор функция, задающая цикл. То есть правая часть аналитически не задается. Более того все начинается с того, что 1. надо знать начальное условие для цикла с точностью х знаков, 2. найти численный метод для решения дифур, который бы выдавал на гора решения хотя бы в течение одного периода с точностью х.

Если положить в системе (1) $\beta_1=0$, то получим автономную нелинейную систему, в которой есть только одна седловая бифуркация. Аналогично из уравнений линеаризации, для тех же значений параметров, я определила, что система асимтотически устойчива. Собственные значения

lambda =
  1.0e+02 *
-0.001287080938828 + 0.027858096249451i
-0.001287080938828 - 0.027858096249451i
-1.359373451767490

То есть первых два и третий отличаются на два порядка и отсуюда растут ноги(!). Насколько я понимаю, система с косинусом будет как бы соответствовать движению этой фиксированной точки со временем. И оценка точности:

exp(lambda(3))
ans =
    9.186733228133013e-60

Matlab сдался без боя, Mathematica просто не интегрирует численно уравнения (Рунге-Кутта 4-5 не хватает и на 35 знаков).. Короче, пришлось мне поставить arbitrary precision libraries, вбить Рунге-Кутта 7-8 и на данный момент я пытаюсь найти более точные начальные условия для цикла, а сходимость надо сказать совсем не быстрая и работать надо уже не с double... Пока что для 35 знаков только один элемент $P(t)=X(t)e^{-tR}$ периодический, в остальных - мусор (в конце периода так точно). Пробовала находить $P(t)$, решая матричное дифференциальное уравнение $P'(t)=A(t)P(t)-P(t)R, P(0)=X(0)=E$ - результат тот же. Пробовала перейти к каноническим переменным, тогда фундаментальная матрица будет $X_0(t)=X(t)S=P_0(t)e^{tR_0}$, где $P_0(t)=P_0(t+T)$, $R_0=\mathrm{diag}(\mu_1,\mu_2,\mu_3)$, $S$ - матрица собственных векторов $X(T)$. В таком представлении элементы в первых двух столбцах $P(t)$ периодические, а третий - мусор.

Gafield в сообщении #284736 писал(а):

Для сведения к двумерной системе надо знать собственное подпространство для двух других показателей. Типа матрица после линейного преобразования будет $\mathrm{diag}\,(\lambda,A)$, где $A$ - блок два на два. Весь вопрос, с какой точностью. Возможно, в задании подпространства нужна точность порядка обратной к этому большому множителю растяжения.

То есть надо перейти от $x$ к новым коодинатам $z$, таким что уравнения станут $z'(t)=Bz(t)$, где $B=\left ( \begin{array}{ccc} b_{11} & b_{12} & 0 \\ b_{21} & b_{22} & 0 \\ 0 & 0 & \mu_{3} \end{array}\right)$, $\mu_3$ - третий показатель? Я не знаю, как это сделать, пока даже "геометрически" представляю приблизительно

Gafield в сообщении #284736 писал(а):

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

А поподробнее? Тоже что-то подобное мне в голову приходило, потому я и начала этот топик, спрашивая по поводу устойчивости цикла в обратную сторону.

[Gafield]

Про обратную сторону имелось в виду следующее. Систему, котоорая имеет переменные коэффициенты, конечно, не преобразуешь. Однако за период все сводится к умножению на некоторую матрицу $A$. Пусть у нее есть один собственный вектор с очень большим собственным значением и собственным двумерным подпространством, в котором вектора не сильно растягиваются при умножении на $A$. Возьмем произвольный вектор $x_0$ единичной длины. Тогда $x_0$ под действием $A$ сильно растянется и вектор $x_1=Ax_0/\|x_0\|$ будет очень похож на собственный (если случайно не окажется из двумерного подпространства, которое почти не растягивается). При желании можно его уточнить, взяв $x_2=Ax_1/\|x_1\|$ и т.д. Отсюда получим с.з. $\lambda$. Дополнив с.в. $f_1$ до базиса $f_1$, $f_2$, $f_3$, можно приближенно найти нерастягивающееся подпространство$\langle e_2,e_3\rangle$: $e_2\approx f_2-Af_2/\lambda$, $e_3\approx f_3-Af_3/\lambda$. Если затем взять в качестве начальных данных ортонормированный базис из $\langle e_2,e_3\rangle$ и посчитать куда он перейдет за период (или его проекцию на $\langle e_2,e_3\rangle$ вдоль с.в., если вычисления приближенные), получится двумерная матрица.

Или, если позволяет точность, можно обойтись и без сведения к двухмерному случаю. Надо просто найти эту матрицу $A$, расчитав траектории с начальными условиями, образующими единичную матрицу, в обратную сторону за период. Тогда требуемая матрица будет $A^{-1}$.

[LynxGAV]

Gafield, я еще осмыслю, что вы написали, но ведь решения расходятся настолько быстро, что я боюсь не проинтегрировать уравнения в течение периода.