Методы анализа нелинейных динамических моделей. - Холодниок М.
Скачать (прямая ссылка):
мультипликаторами) соответствуют собственным числам линейной части
отображения Пуанкаре и определяют устойчивость периодического решения.
При этом периодическое решение является устойчивым (точнее, орбитально
асимптотически устойчивым), если выполняются неравенства
\Xt\<l, i = 2, ...,". (5.8.18)
Периодическое решение является неустойчивым, если хотя бы для одного
мультипликатора имеет место условие |А,г|> 1. Если
5.8. Вычисление периодических решений в автономном случае
217
|Я.|| ^ Ю-4-10-5, можно говорить о сильно неустойчивом периодически
решении. В п. 5.8.5 мы рассмотрим различные случаи изменения устойчивости
периодических решений при изменении параметров задачи.
5.8.4. Продолжение периодических решений по параметру
Займемся теперь исследованием решения задачи (5.8.1),
(5.8.2) (или, что то же самое, (5.8.4), (5.8.5)) в зависимости от
параметра а. Опишем алгоритм для нахождения указанной зависимости
периодических решений, который основывается на алгоритме DERPAR из § 5.2,
использовавшемся для продолжения решений нелинейных (алгебраических)
уравнений.
При фиксированном r\k (значение индекса k также не изменяется)
соотношения (5.8.10) принимают вид
Л(Л> Ч*-ь Чн1 Чп, Т, а) = 0, г = 1,2,..., п. (5.8.19)
Решение этой системы п уравнений относительно п + 1 неизвестных t]i, ...,
r)fe_i, Цй+ь ..., цп, Т, а можно исследовать с помощью алгоритма DERPAR.
Для процесса продолжения нам потребуется вычислять частные производные
функций Fi по tj, Т и а. Вычисление dFi/dr\j, dFi/dt описано в п. 5.8.2.
(Непосредственно для самого процесса продолжения производные dFi/дцк не
нужны; для определения устойчивости по матрице В (5.8.17) они
необходимы.) Покажем, как вычислить производную dF/да. Обозначим r?- =
dxi/dа; дифференцируя уравнение (5.8.4) по а, получим уравнения в
вариациях
r'>=Ttwr.+T% '='-г л <5-8-20>
s= 1
с начальными условиями
г,(0) = 0. (5.8.21)
К формулам (5.8.15) и (5.8.16) следует теперь добавить соотношения
-Ц- = г?(1), г = 1, 2, ...,". (5.8.22)
Успешно продолжать решение системы (5.8.19) (т. е. двигаться вдоль
(связной компоненты) кривой (5.8.19).- Ред.) можно до тех пор, пока r\k
будет лежать на профиле xk(z) (см. рис. 5.24). Если же fjfe "покинет"
профиль xk(z), т. е. при данном выборе fjfe периодическое решение не
будет существовать,
218
Глава 5
то данный алгоритм не сработает. Поэтому в процессе продолжения мы будем
адаптивно "приспосабливать" значение r\k к функции Xk \z).
Рассмотрим этот вопрос более подробно. Предположим, что нам известно
некоторое решение системы (5.8.19) и найдено соответствующее
периодическое решение системы (5.8.4). Опишем, как выбирается значение
fjft при следующем значении а (см. рис. 5.24).
1. Возьмем монотонный участок периодического профиля на промежутке [0,
Zi](J [z2, 1] (где лежит значение rjft). Через х^ях (или соответственно
х(tm)1п) обозначим максимум (или соответственно минимум) на этом участке.
Введем следующие обозначения (0 < coi < 1):
xt =Т К1 ~ т0 х?п + 0 + (5.8.23а)
** =Т[0 +(r),)*Гп+(1 -(r)!)*ГХ]. (5.8.23Ь)
При этом для rjft потребуем выполнения неравенства
хк 1* < xt- (5.8.24)
Это условие означает, что fjft располагается достаточно далеко от
локального максимума и локального минимума. Поэтому можно надеяться, что
для следующего значения параметра а значение % "не пропадает" с профиля
Хк.
2. Обозначим через Ятах максимальную разность между
максимумом и минимумом на монотонном участке профиля
xk(z), если этот монотонный участок выбирать из всего проме-
жутка ге[ 0, 1]:
Ящах == (^"*") %k (2 )¦
Потребуем теперь, чтобы
хГ-хГп>(1~<о2)Нтах. (5.8.25)
Это; второе условие означает, что v\k лежит на достаточно широком
монотонном участке.
Опыт показывает, что обычно разумно выбирать значения сох и (о2 в
диапазоне [0.5, 0.8].
Если хотя бы одно из условий (5.8.24) и (5.8.25) не выполняется, то
значение г\к выбирается заново по формуле
jjHOBoe = [Хк (2+) + (г-)]. (5.8.26)
1) Речь идет о монотонности функции хк (г).-Прим. ред.
5.8. Вычисление периодических решений в автономном случае
219
Остальные составляющие Tjft-i, т)*+1, ть должны
пересчитываться с использованием профилей х\ (z), х"(г),
которые были найдены для последнего значения параметра а. Для этого
пересчета выберем координату г* между zr и z+ таким образом, чтобы
выполнялось условие
** (2*) ~ ЛГвое- (5.8.27)
Тогда значения т^овое = хг (z*), ?=1, ..., п, будут представлять собой
стартовые значения для дальнейшего процесса продолжения. При этом
многошаговая схема интегрирования (например, схема Адамса-Бэшфорта),
используемая в качестве предиктора в алгоритме DERPAR, вновь начинается
со схемы первого порядка (т. е. схемы Эйлера). При такой организации
вычислений должны также пересчитываться и некоторые управляющие параметры
в алгоритме DERPAR. Прежде всего это относится к параметрам направления
Nh поскольку нам необходимо сохранять направление движения вдоль ветви
периодических решений. Все детали указанного подхода читатель может найти
