CC BY
90
ВЫЧИСЛЕНИЕ ПОКАЗАТЕЛЕЙ ЛЯПУНОВА ДЛЯ РАСПРЕДЕЛЁННЫХ СИСТЕМ: ПРЕИМУЩЕСТВА И НЕДОСТАТКИ РАЗЛИЧНЫХ ЧИСЛЕННЫХ МЕТОДОВ
П.В. Купцов
При вычислении показателей Ляпунова для распределённых систем возникают специфические сложности, обусловленные природой этих систем. В этой статье обсуждается точность разных алгоритмов ортогонализации применительно к возникающим в ходе расчётов плохо обусловленным матрицам большого размера. Также исследуется паразитное возбуждение коротковолновых пространственных гармоник, которое, как было обнаружено, может происходить при использовании для решения уравнений метода конечных разностей и приводит к грубым ошибкам вычисления показателей. На основе результатов выполненного анализа формулируются практические рекомендации по выбору оптимальных численных методов.
Ключевые слова: Показатели Ляпунова, пространственно-временной хаос, метод конечных разностей, комплексное уравнение Гинзбурга-Ландау, рЯ-разложение.
Введение
Показатели Ляпунова характеризуют устойчивость движения динамических систем в фазовом пространстве. Математическое обоснование существования этих показателей для сосредоточенных систем даёт мультипликативная эргодическая теорема Оселедца [1,2]. Распределённые системы имеют бесконечное число степеней свободы, и поэтому теорема Оселедца напрямую к ним не применима. В работе [3] обсуждаются обобщения теоремы для случая бесконечномерного фазового пространства.
Рассуждая нестрого, можно сказать, что если протяжённость распределённой системы конечна, то для неё можно построить спектральное представление, получив при этом вместо уравнения в частных производных бесконечную систему обыкновенных дифференциальных уравнений. Далее можно перейти к приближенному описанию с сохранением только конечного числа уравнений. Для такой приближённой модели показатели Ляпунова уже будут хорошо определены.
Если рассматривается диссипативная распределённая система, которая в ходе своей эволюции осуществляет сжатие фазового объёма, то её аттрактор может
быть вложенным в так называемое инерциальное многообразие, имеющее конечную размерность [4]. Это означает, что динамику такой системы можно, по крайней мере, в принципе описать при помощи обыкновенных дифференциальных уравнений. И снова это означает возможность вычисления показателей Ляпунова для распределённой системы.
Пока не известна практически реализуемая процедура нахождения инерци-альных многообразий. Однако в недавней работе [5] было показано, что в касательном пространстве аттракторов диссипативных систем можно при помощи численной процедуры выделить конечное число активных направлений, взаимное расположение которых, по всей видимости, в значительно степени определяет динамику системы на аттракторе. Все прочие направления (их имеется бесконечно количество) участвуют только в переходных процессах и представляют собой собственные моды однородного состояния равновесия [6]. В работе [5] высказано предположение, что обнаруженное конечное множество активных направлений есть не что иное как локальный базис инерциального многообразия.
Вычисление показателей Ляпунова для распределённых систем кажется на первый взгляд вполне «понятной» процедурой. Действительно, ставший фактически стандартным для сосредоточенных систем алгоритм Бенетти [7] естественным образом переносится на случай систем, заданных уравнениями в частных производных. Тем не менее, аккуратный анализ показывает, что при этом требуется соблюдать определённую осторожность. Цель настоящей работы состоит в том, чтобы выявить, какие проблемы могут возникнуть при воспроизведении алгоритма Бенетти применительно к распределённым системам.
В разделе 1 приводится определение показателей Ляпунова и описывается стандартный алгоритм Бенетти, используемый для их вычисления. Раздел 2 посвя-щён рассмотрению преимуществ и недостатков различных процедур ортогонализа-ции, применяемых в процессе вычисления показателей. В разделе 3 анализируется паразитное возбуждение коротковолновых пространственных гармоник, возникновение которого вызывает грубые ошибки при вычислении показателей Ляпунова. Наконец, в разделе 3 приводятся практические рекомендации, выработанные по результатам выполненного исследования.
1. Показатели Ляпунова
Пусть имеется система, заданная обыкновенным дифференциальным уравнением
и = f (u,t), (1)
где u = u(t) Е Rm - m-мерный вектор, описывающий состояние системы в момент времени t, а f (u, t) Е Km - нелинейная векторная функция. Рассмотрим возмущения, которые можно придать траектории этой системы. Если их амплитуда всё время остаётся бесконечно малой по отношению к исходному фазовому пространству, то такие возмущения называются инфинитезимальными и описываются линейным уравнением вида
v = J(u, t)v, (2)
где v = v(t) Е Km - m-мерный вектор возмущения, а J(u, t) - якобиан, то есть матрица, построенная из производных компонент векторной функции f (u, t) по компонентам вектора u. Такие возмущения называют также вариациями траектории, а
само уравнение (2) - вариационным уравнением. Векторы V представляют собой касательные к траекториям исходной системы. Всевозможные касательные векторы образуют пространство, которое называется тангенциальным или касательным и имеет размерность, равную размерности фазового пространства т. Говорят, что уравнение (2) описывает динамику в касательном пространстве.
В зависимости от направления приложения возмущения и свойств системы, вектор V нарастает или затухает, причём в силу линейности уравнения (2), это происходит в среднем по экспоненциальному закону. Как гласит мультипликативная эр-годическая теорема Оселедца [1], существует набор чисел А^ ^ Х2 ^ ... ^ Хт, количество которых равно размерности касательного пространства, таких что для любого начального возмущений у(0) существует показатель
где Х(г/(0)) принимает значения из набора Х в зависимости от выбора г7(0). Числа Х называют показателями Ляпунова. Также из мультипликативной эргодической теоремы следует, что сумма первых к показателей Ляпунова представляет собой средний показатель экспоненциального сжатия или растяжения к-мерного фазового объёма.
Старший показатель Ляпунова А.1 вычисляется следующим образом. Мы должны в течение достаточно большого промежутка времени решать совместно уравнения (1) и (2), периодически выполняя перенормировки V и накапливая логарифмы норм, а затем усреднить накопленные значения за время счёта. Однако эта процедура не годится, когда нужно вычислить два или более показателей. Если попробовать решать сразу несколько уравнений вида (2), то независимо от выбора начальных условий все решения после некоторого времени счёта будут вести себя одинаковым образом и выстраиваться вдоль одного и того же направления, соответствующего Хь Это связано с тем, что старший показатель доминирует над всеми остальными. Чтобы этого избежать, в работе [7] было предложено периодически выполнять ортого-нализацию системы решений таким способом, чтобы минимизировать это доминирование.
Алгоритм вычисления показателей Ляпунова выглядит следующим образом [7,8].
1. В течение интервала времени Т решаем уравнение (1) совместно с п уравнениями для инфинитезимальных возмущений (2), где п - количество показателей Ляпунова, которые мы хотим найти, п ^ т. Из полученных решений уравнений (2) строим матрицу V.
2. Выполняем рЯ-разложение [9] матрицы V, представляя её в виде произведения V = QR. Матрица Q является ортогональной, причём к-й столбец V является линейной комбинацией первых к столбцов Q. Матрица R имеет верхнюю треугольную структуру, её элементы суть проекции столбцов V на столбцы Q.
3. Произведение первых к диагональных элементов матрицы R равно объёму к-мерного параллелепипеда в касательном пространстве, в который трансформируется к-мерный единичный куб за время Т. Поэтому для вычисления показателей Ляпунова мы накапливаем логарифмы диагональных элементов Я.
4. Используем матрицу Q в качестве начального состояния уравнений для ин-финитезимальных возмущений, решаем их в течение времени Т и повторяем процедуру ортогонализации.
(3)
5. После достаточно большого числа повторов усредняем накопленные логарифмы за время счёта и получаем показатели Ляпунова.
Для того чтобы инициализировать эту процедуру, можно сгенерировать случайную матрицу V и выполнить её QR-разложение. Полученную матрицу R нужно отбросить, а Q использовать как стартовое значение для итераций.
Описанная процедура как таковая не зависит ни от природы решаемых уравнений, ни от метода поиска решений. Единственное требование состоит в том, что мы должны использовать столбцы матрицы Q как начальные состояния уравнений для инфинитезимальных возмущений, а найдя решения этих уравнений за время T, представить их в виде матрицы V. Поэтому, если система задана уравнением в частных производных и её решение ищется численно, то показатели Ляпунова для такой системы могут быть найдены при помощи данного алгоритма.
2. Выбор метода вычисления QR-разложения
Основная проблема при вычислении показателей Ляпунова состоит в том, что из-за доминирования старших показателей над младшими столбцы матрицы V стремятся выстроиться вдоль одного и того же направления. Для предотвращения этого должна быть использована процедура ортогонализации, которая хорошо справляется с матрицами, имеющими столбцы почти параллельные друг другу. Широко применяются следующие методы [9]:
• преобразование Хаусхолдера (HQR);
• модифицированный метод Грама-Шмидта (MGS);
• классический метод Грама-Шмидта (CGS).
Известно, что наиболее высокую точность имеет метод на основе преобразования Хаусхолдера [9]. Именно его рекомендуется использовать при вычислении показателей Ляпунова [10]. Однако исторически сложилось так, что чаще всего применяются методы Грама-Шмидта, которые, как известно из учебников по численным методам, могут иметь большую погрешность [9]. Поэтому далее мы протестируем три упомянутых метода применительно к системе с высокой размерностью фазового пространства, чтобы выяснить, в каких случаях методы Грама-Шмидта работают хорошо, а когда их использование нежелательно.
Идея классического метода Грама-Шмидта такова. Первый вектор-столбец щ матрицы V мы только нормируем на единицу, получая ql. Второй вектор V2 мы поворачиваем в плоскости векторов ql и V2 таким образом, чтобы он стал перпендикулярным ql, а затем нормируем его на единицу и получаем q2. Тем самым мы сохраняем подпространство, которому принадлежали векторы vi и v2 до ортогонализации, и поэтому V2 может быть выражен как линейная комбинация ql и q2. Далее, третий вектор точно также поворачиваем в пространстве векторов ql, q2 и V3 так, чтобы он стал перпендикулярным qi и q2, нормируем его и получаем q3. С остальными поступаем аналогичным образом. После того как процедура выполнена, мы получаем матрицу Q, столбцами которой являются qi. Одновременно мы находим проекции щ на qi, которые являются диагональными элементами матрицы R. Все прочие элементы этой матрицы не используются при вычислении показателей Ляпунова.
Модифицированный метод Грама-Шмидта работает точно так же, за тем ис-
ключением, что на k-м шаге процедуры мы одновременно поворачиваем сразу все пока не обработанные векторы Vj, j = k,k + 1,... так, чтобы Vk занял правильное положение.
Существенной частью процедур Грама-Шмидта является вычисление скалярных произведений, что фактически означает нахождение косинусов углов между векторами. Так как производная косинуса в нуле обращается в нуль, то в случае малых углов их косинусы достаточно быстро становятся неразличимыми. Поэтому, если углы между столбцами исходной матрицы малы, то процедура Грама-Шмидта даёт неверный результат. Частично эту проблему решает модифицированный метод Грама-Шмидта, в котором удаётся получить более высокую точность за счёт изменения порядка вычисления скалярных произведений.
Вращение Хаусхолдера состоит в том, что для вектора v ищется такая ортогональная матрица H (она называется матрицей Хаусхолдера), что HV = ||v||ei, где ei = (1, 0,...). Иными словами, умножение на соответствующую матрицу Хаусхолдера приводит к тому, что вектор v ориентируется вдоль первой координатной оси. Существенно, что процедура нахождения H очень проста и выполняется достаточно быстро и с высокой точностью [9]. В общем случае матрица вращения кроме косинусов угла поворота содержит также и синусы, которые хорошо определены для малых углов. По этой причине преобразование Хаусхолдера хорошо работает даже когда поворот осуществляется на малый угол.
При вычислении QR-разложения сначала ищется матрица Хаусхолдера для первого столбца исходной матрицы V. После умножения её на V получается матрица, у которой в первом столбце все элементы, кроме первого, равны нулю. Далее строится новая матрица вращения, которая оставляет без изменения первый столбец, а остальные преобразует таким образом, что во втором столбце ниже главной диагонали появляются нули. Затем строится матрица, поворачивающая все столбцы, кроме первых двух, и в результате в третьем столбце ниже главной диагонали также появляются нули. В конечном итоге вместо исходной матрицы получается верхняя треугольная матрица R, а перемноженные матрицы Хаусхолдера образуют матрицу QT [9].
Для проверки качества работы методов QR-разложения вычисляем погрешность ортогонализации по следующей формуле:
Sort = ||QTQ - 11|max, (4)
где ||A||max = maxjj laj | - норма матрицы. Вычислим также число обусловленности ортогонализируемых матриц по формуле
C = ||A||max||A~ || max. (5)
Большое значение C в нашем случае означает малость углов между столбцами матрицы. При вычислении C следует принимать во внимание, что инверсия плохо обусловленной матрицы может быть вычислена с большой погрешностью. Мы будем искать обратные матрицы, используя LU-разложение, как описано в книге [11], и проверять точность выполненного обращения, вычисляя величину Sinv = || A_1A—I ||max. Максимальной приемлемой погрешностью будем считать Sinv = 10"3.
Рассмотрим комплексное уравнение Гинзбурга-Ландау
dtu = u — (1 + ic)|u|2u + (1 + ib)dX u, (6)
где u = u(x, t) - комплексная динамическая переменная, а c и b - управляющие параметры. В зависимости от выбора значений c и b это уравнение может демонстрировать большое количество различных видов поведения. Мы возьмём c = 3, b = -2, что соответствует одному из хаотических режимов, а именно режиму амплитудного хаоса [12]. Протяжённость системы в пространстве возьмём равной X = 50, на границах положим
9xu\x=0,x = °. (7)
Будем решать это уравнение численно с шагом дискретизации по времени At = 0.01 и шагом по координате Ax = X/(M — 1), где M - число точек пространственной сетки. Так как переменная u комплексная, то размерность фазового пространства численной модели равна m = 2M.
Для решения уравнения (6) будем применять метод конечных разностей и использовать нелинейную смешанную схему. Идея состоит в том, что для перехода от временного слоя t к слою t + At строится система из m алгебраических уравнений выражающих решение на новом слое через m уже известных значений переменной u на предыдущем слое. При этом, в силу нелинейности уравнения (6), система алгебраических уравнений также получается нелинейной, и для её решения применяется метод Ньютона.
Совместно с основным уравнением будем решать n уравнений для инфините-зимальных возмущений v = v(x,t):
dtv = v — 2(1 + ic)\u\2v — (1 + ic)uV + (1 + \b)ôXv- (8)
Для вычисления решения мы будем применять неявную схему, описание которой приведено ниже, в разделе 3.4. Решения образуют столбцы матрицы V, которая, для предотвращения расхождения и выстраивания, после каждых T единиц времени подвергается QR-разложению одним из описанных выше методов. Чем больше
интервал T, тем сильнее проявляется выстраивание и тем меньше могут быть углы между столбцами V. Иными словами, число обусловленности этой матрицы растёт с ростом T.
На рис. 1 показана зависимость погрешности ортогонализации, выполненной различными методами, от размерности фазового пространства m и количества уравнений n. Интервал между QR-процедурами T равен 0.1. При таком T максимальное число обусловленности матрицы V при n = m = 250 не превышает нескольких единиц. Из рисунка видно, что в данном случае все три QR-метода имеют очень низкую погрешность.
Если увеличить T до 0.6, то число обусловленности при n = m = 50
Рис. 1. Погрешность ортогонализации матрицы решений n уравнений для инфинитезимальных возмущений (8); m - размерность фазового пространства численной модели; HQR - преобразование Хаусхолдера, MGS и CGS - модифицированный и классический методы Грама-Шмидта, соответственно. Интервал между QR-процедурами T = 0.1
составляет несколько единиц, при п = е, = т = 150 оно достигает нескольких десятков, а при п = т = 250 его значение имеет порядок 106. Как видно из рис. 2, в последнем случае погрешность методов Грама-Шмидта быстро нарастает с ростом п. При этом погрешность классического метода достигает неприемлемо высоких значений. Иными словами, методы Грама-Шмидта работают в данном случае хорошо, только если число искомых показателей п много меньше их полного количества т.
Метод на основе преобразования Хаусхолдера хорошо работает даже когда число обусловленности ортогонали-зируемой матрицы велико. Его погрешность находится на уровне 10_15. Это значение близко к машинному эпсилон*, которое для чисел двойной точности приблизительно равно 10_16. Такая точность, по всей вероятности, является максимально возможной.
Для выяснения границ применимости метода Хаусхолдера на рис. 3 построены значения диагональных элементов матриц И,, полученных в результате выполнения нескольких РЯ-шагов. Вычисления производились при
орт Ю-14
ю-15
Ю-16
Ю-13 Ю-14
Ю-15
10"3
ю-9
Ю-15
.да=50; — HQR : ----- MGS : CGS 1 1 1
^ да 150 A.../V'""""' ' ■ i i i i
:m=25Ö ..............
i—1 * i i i i
50
100
150 200
Рис. 2. То же, что на рис. 1, при T = 0.6
Рис. 3. Диагональные элементы матриц R, полученных после нескольких QR-шагов, выполненных с использованием вращения Хаусхолдера при различных T: 1 - 1.0; 2 - 2.0; 3 - 5.0. По вертикальной оси использован логарифмический масштаб. Горизонтальная пунктирная линия проведена на уровне машинного эпсилон 10~16
Ах = 0.5. Зависимости 1 и 2 представляют по 10 наборов значений т, вычисленных при Т = 1.0 и Т = 2.0, соответственно. Так как зависимости 1 и 2 получены при разных Т, то они имеют разные масштабы. Однако их структура одинакова и, если выполнить усреднение их логарифмов за время счёта, получим одни и те же значения показателей Ляпунова. Диагональные элементы тц, рассчитанные при Т = 5.0, выглядят иначе: примерно при г = 100 зависимость претерпевает излом (см. кривую 3), которого нет на двух других кривых. Это происходит, когда тц становится меньше машинного эпсилон, отмеченного на рисунке горизонтальной пунктирной линией. Так как диагональные элементы И суть составляющие векторов щ вдоль (¡г, то можно сделать вывод, что ортогонализация на основе преобразования Хаус-холдера перестаёт корректно работать, когда угол между векторами настолько мал, что ортогональные их направлениям проекции имеют величину меньше машинного эпсилон.
Таким образом, можно заключить, что наилучшим является метод ортогона-лизации на основе преобразования Хаусхолдера. Однако, чтобы избежать ошибки
* Напомним, что машинное эпсилон - это величина шага дискретизации машинного представления вещественных чисел в окрестности единицы.
при слишком большом Т, нужно в процессе накопления диагональных элементов И, проверять их значения - они должны быть хотя бы на несколько порядков больше машинного эпсилон.
Основным недостатком методов Грама-Шмидта является рост погрешности с увеличением интервала Т. Особенно существенно это для классического метода. Применение методов Грама-Шмидта для вычисления показателей Ляпунова допустимо, если количество искомых показателей много меньше размерности фазового пространства дискретного представления системы. В случае, когда требуется найти все или почти все показатели, нужно выполнять ортогонализации достаточно часто, так чтобы число обусловленности матрицы решений уравнений для инфинитези-мальных возмущений оставалось низким.
3. Паразитное возбуждение коротковолновых пространственных гармоник
Вычисляя показатели Ляпунова, мы должны решать систему уравнений для инфинитезимальных возмущений. Если для этого используется метод конечных разностей, то при неправильном выборе шага по времени по сравнению с шагом по координате, в решениях могут быть возбуждены коротковолновые гармоники, что приводит к вычислению некорректных значений показателей. В этом разделе мы найдём условия, которые должны быть наложены на величины шагов, чтобы избежать паразитного возбуждения.
3.1. Множитель роста пространственных гармоник. Рассмотрим распределённую систему, заданную уравнением в частных производных
дги = /(и, х, Ь) + Бд2й, (9)
где и = и(х, Ь) - вектор динамических переменных, /(и,х,Ь) - нелинейная векторная функция, и Б - матрица с неотрицательными вещественными частями собственных чисел. Если матрица Б диагональная, то уравнение (9) описывает класс реакционно-диффузионных систем. Эти системы широко используются в качестве моделей активных сред и могут демонстрировать как различные виды структурообразования, так и хаотическую динамику. Комплексное уравнение Гинзбурга-Ландау (6), представленное в виде пары уравнений для вещественной и мнимой частей, также имеет структуру, задаваемую уравнением (9).
Уравнение для инфинитезимальных возмущений V = ь(х, Ь) имеет вид:
д^ = 3(и,х,Ь)у + (10)
где 3(и, Ь) - матрица Якоби, построенная из производных компонент векторной функции /(и,х,Ь) по компонентам вектора и. Это уравнение можно интерпретировать как линейную систему под внешним параметрическим воздействием, за включение которого отвечает матрица Якоби. Выбирая разные начальные распределения, мы можем наблюдать экспоненциальный рост или затухание решения, причём усреднённый показатель экспоненты равен одному из показателей Ляпунова системы (9).
Рассмотрим диффузионную подсистему системы (10), отбросив слагаемое, пропорциональное матрице Якоби
дгУ = Бд^у. (11)
Разложим решение этого уравнения по пространственным гармоникам вида V = = V. Здесь V - некоторый постоянный вектор, в общем случае имеющий комплексные компоненты, ю и к - вещественные частота и волновое число, соответственно, а а > 0 - вещественный множитель роста гармоники. Подставив выражение для гармоники в уравнение (11), получим а = е-к Ке где й - одно из собственных чисел матрицы Б. Если Ие й < 0, то все пространственные гармоники неустойчивы, что не представляет для нас интереса, и далее этот случай рассматриваться не будет. Когда Ие й = 0, все гармоники находятся в состоянии безразличного равновесия. При Ие й > 0 имеет место диффузионный процесс, когда исходное распределение затухает и стремится к однородному состоянию. При этом коротковолновые гармоники затухают быстрее, так что длинноволновые гармоники преобладают.
Вычисляя показатели Ляпунова, мы строим конечно-разностное уравнение, приближённо соответствующее уравнению (10). Для того чтобы такое приближение было адекватным, необходимо, в частности, добиться, чтобы диффузионная подсистема численной модели вела себя качественно подобно диффузионной подсистеме (11), полученной для непрерывной системы. А именно, если решение соответствует собственному числу с Ие й > 0, то в его структуре должны преобладать длинноволновые гармоники, а если Ие й = 0, то все пространственные гармоники должны быть безразлично устойчивыми.
Пусть протяжённость системы равна X. Введём пространственную сетку из М точек. Шаг сетки равен Ах = Х/(М — 1). Шаг дискретизации по времени обозначим как АЬ. Пусть в есть номер точки в пространстве и ] - номер шага по времени. Запишем конечно-разностное представление уравнения (11) в следующем виде:
— $8,] = АЬ [(1 — а)БЛу8 ] + аБЛй^+1], (12)
где у8] - значение вектора V в точке сетки в на шаге времени а Лу8] = = ($8+1— 2у8] + у8-1 ])/Ах2 есть кончено-разностная вторая производная по координате. Параметр а позволяет выбирать вариант численной схемы: а = 0 - явная схема Эйлера; а = 1 - неявная схема; а = 1/2 - смешанная схема Кранка-Николсона [11,13].
По аналогии с рассмотренным выше непрерывным случаем, запишем пространственную гармонику для уравнения (12) как у8] = Уа]е1(ю]+^8), где волновое число обозначено буквой к. Подставляя выражение для гармоники в уравнение, получим
2 = 1 — у(1 — а)(8 + 8*) + у2(1 — а)288*
а 1+ уа(8 + 8*) + у2а288* ' (3)
где «*» обозначает комплексное сопряжение,
V = 8ш2(к/2), 8 = 4йАЬ/Ах2. (14)
Комплексный, в общем случае, параметр 8 характеризует пространственно-временную
решётку, на которой ищется решение. Его значение определяется значениями шагов
по времени и по координате, а также одним из собственных чисел d матрицы Б. Разным d соответствуют разные ветви решения и разные значения множителя роста.
Величина V зависит от волнового числа к и всегда принимает значения в диапазоне от 0 до 1. Её удобно использовать в качестве номера гармоники. Обсудим, как соотносятся значения V с длинами волн гармоник. Если на границах либо само решение, либо его первые производные по координате обращаются в нуль, то волновое число непрерывной системы будет принимать значения к = 1п/Х, где I = 0,1, 2,... - номер гармоники. Переходя к дискретной системе, получим: кх = (1п/Х)(зХ/(М — 1)) = кз, где к = 1п/(М — 1). Следовательно, волновое число к пробегает дискретный ряд значений от 0 до п. При этом V меняется монотонно от 0 до 1. Это значит, что близкие к нулю V соответствуют длинноволновым гармоникам с малым к, а V вблизи единицы отвечают коротковолновым гармоникам. Если на систему наложены периодические граничные условия, то, рассуждая аналогично, получим, что волновое число к пробегает значения от 0 до 2п. При этом V сначала растёт, достигает максимума V = 1 при к = п, а потом убывает и обращается снова в нуль при к = 2п. Заметим, однако, что гармоники, расположенные симметрично относительно точки к = п, с одной стороны, имеют одинаковые длины волн, так как е1(2п-&)« = а с другой - им соответствуют одинаковые значения V. Таким об-
разом, интерпретация значений V сохраняется - малые V отвечают длинноволновым гармоникам, а V вблизи единицы соответствуют коротковолновой части спектра.
Рассмотрим поведение множителя роста для различных численных схем, соответствующих разным значениям о.
3.2. Смешанная схема Кранка-Николсона. Пусть о = 1/2, что даёт смешанную схему Кранка-Николсона. Квадрат множителя роста (13) при таком о преобразуется к виду
2 4 — 2v(6 + 6*)+ v266*
а2 =------(15)
4 + 2v(6 + 6*)+ v266*' ^ '
Смешанная схема имеет второй порядок точности локальной аппроксимации производных по времени и по координате и, кроме того, она абсолютно устойчива, так как независимо от 6 множитель роста меньше или равен единице для всех пространственных гармоник. В силу этого именно смешанная схема традиционно считается предпочтительной для решения уравнений вида (9).
Предположим сначала, что Ие I = 0. Тогда а = 1 независимо от V. Это в точности соответствует тому, что мы получили для множителя роста гармоник непрерывной системы (11).
Если Ие I > 0, то, в отличие от непрерывного случая, функция а2^) (15) немонотонно зависит от V. Она равна единице при V = 0, убывает и достигает минимума в точке
V0 = 2/|6|. (16)
Отметим, что функция а2^) имеет также особенности в точках, где её знаменатель обращается в нуль. Однако в случае комплексного I в этих точках V также
комплексное, а если Щ вещественное, то эта особенность находится на отрицательной полуоси из-за того, что мы рассматриваем только d > 0. Поэтому имеет значение только положение минимума. Если |8| < 2, то минимум находится правее единицы, и, следовательно, в интересующей нас области изменения V £ [0,1] множитель роста монотонно убывает, что качественно соответствует поведению множителя роста гармоник а для аналитического решения уравнения (11). Принимая во внимание выражение для 8 (14), можно записать условие монотонности vo > 1 в следующей форме:
Ах2
V! У0
Рис. 4. Множитель роста пространственных гармоник для различных 8: 1 -1.0, 2 - 6.0; смешанная схема при а = 1/2
At <
2\d\
(17)
Зтг/4
Till
к/4
0
■■c.irt <r
у*
fx
у У
у
к-
1
50 100 150
Рис. 5. Фурье-спектры столбцов матрицы Q для уравнения Гинзбурга-Ландау. Смешанная схема. Длина системы X = 25, параметры с = 3, Ь = -2, граничные условия (7). Ах = 0.25, АЬ = 0.013. Условие монотонности (17) выполняется
Это условие накладывает достаточно сильное ограничение на допустимые значения АЬ и Ах. Если Ах = 0.1 и Щ принять равным единице, то Дt < 0.005.
На рис. 4 построены графики зависимости а2^) при различных значениях 8, которые для простоты считаются вещественными. Кривая 1 иллюстрирует случай, когда |8| мало и условие монотонности (17) выполняется.
Для иллюстрации случая монотонного поведения множителя роста обратимся к уравнению Гинзбурга-Ландау (6) с граничными условиями (7). Как и раньше, будем решать это уравнение при помощи смешанного нелинейного метода, идея которого описана в разделе 2. Значения Ах и Дt подберём так, чтобы выполнялось условие монотонности (17). Чтобы найти Щ, разложим уравнение на вещественную и мнимую части и, построив матрицу Б в явном виде, получим: Щ = 1 ± 16. Вместе с основным уравнением будем решать систему уравнений вида (8), используя смешанную схему Кранка-Николсона и выполняя рЯ-разложения матрицы решений по методу Хаусхолдера после каждых Т единиц времени. На рис. 5 изображены фурье-спектры столбцов матрицы Q, полученной после очередной рЯ-процедуры. По горизонтальной оси отложены номера столбцов, которые соответствуют номерам показателей Ляпунова, а вдоль вертикальной оси отложены значения волнового числа. Оттенки серого показывают значения фурье-амплитуд. Для вычисления спектров из каждого столбца Q отбирались элементы, соответствующие только вещественным частям переменной и.
Из рисунка видно, что спектры имеют чётко выраженное основное волновое число и его значение пропорционально номеру столбца, то есть номеру показателя Ляпунова. Устройство спектров обусловлено разбиением касательного пространства уравнения Гинзбурга-Ландау на активные, или физические, и изолированные моды. Сгущение в левом нижнем углу плоскости образовано активными модам. Есть
основания полагать, что наблюдаемая динамика системы в значительной степени обусловлена взаимным расположением этих мод. Диагональная структура в центре представляет изолированные моды, которые являются касательными к траекториям, отвечающим за переходные процессы [5,6]. Сгущение в правом верхнем углу, по всей вероятности, также соответствует активным модам, которые задействованы, когда динамика рассматривается в обратном времени [6].
Предположим теперь, что условие монотонности (17) нарушено. В этом случае существует такая гармоника VI < vo, множитель роста которой совпадает с множителем роста коротковолновой гармоники V = 1. Это иллюстрирует кривая 2 на рис. 4. Из уравнения а2^) = а2(1) можно найти значение VI:
VI = 4/|6|2. (18)
Все гармоники с номерами в диапазоне от VI до vo, имеют коротковолновых «двойников» с одинаковыми множителями роста. Когда решается одиночное уравнение и шаги пространственно-временной сетки выбраны «разумно», так что точка VI расположена достаточно далеко от нуля, то наличие пар пространственных гармоник с одинаковыми а не играет существенной роли. Они все сравнительно быстро затухают, и в решении преобладают длинноволновые гармоники с более высокими значениями множителя роста. Именно это делает возможным использование смешанной схемы несмотря на аномальное поведение а.
Существенно иная ситуация возникает, когда для нахождения показателей Ляпунова численно решается система уравнений вида (10) и периодически применяется ортогонализация с целью исключить доминирование старших показателей над младшими. Представляется естественным предположение о том, что значение множителя роста связано с величиной показателя Ляпунова. Это значит, что в процессе ортогонализации должно происходить исключение гармоник с высокими множителями роста. Для нахождения первого показателя мы решаем первое уравнение без какого-либо вмешательства в структуру решения; в решении второго уравнения, связанного со вторым показателем, несколько первых гармоник должны быть удалены; в решении третьего уравнения извлекаются ещё несколько длинноволновых гармоник и так далее.
Последовательное удаление из структуры решений длинноволновых гармоник ведёт к тому, что решение с номером щ, связанное с вычислением показателя щ, не будет иметь быстрорастущих гармоник с волновыми числами V < VI. Следовательно, это и все последующие решения с номерами п ^ п1 будут содержать пары гармоник - с каждой гармоникой из нижней части спектра будет сосуществовать её коротковолновый «двойник» с одинаковым множителем роста. Такое паразитное возбуждение коротковолновых гармоник искажает значения соответствующих показателей.
На рис. 6 построена плоскость фурье-спектров столбцов матрицы Q уравнения Гинзбурга-Ландау для случая, когда условие монотонности (17) нарушено. Параметры и граничные условия те же, что и на рис. 5. Видно, что столбцы из левой половины рисунка имеют «нормальную» фурье-структуру, а именно: доминирующее волновое число пропорционально номеру столбца. В центре плоскости мы видим
Зтг/4
я/2
я/4
*>
к0
к1
/ У О. 4 к
1
50 100 150
очень четкую границу, правее которой возникает паразитное возбуждение гармоник из коротковолновой части спектра. Волновое число к\, при котором это происходит, можно вычислить, используя формулу (18) и выражение для V (14). Правее границы спектры содержат по два набора гармоник, которые приближаются друг к другу с ростом номера г и сливаются при ко в точке минимума множителя роста, который можно найти из (16). Найденные таким образом к0 и к1 на рисунке обозначены горизонтальными пунктирными линиями. Видно, что они очень хорошо соответствуют моменту начала паразитного возбуждения и точке слияния двух ветвей спектра.
Из рис. 5 видно, что волновое число к гармоники, преобладающей в решении с номером г, пропорционально этому номеру. Такую закономерность можно записать как
г = сМк/п, (19)
где с - количество компонент системы уравнений, в нашем случае равное двум, а М - число точек пространственной сетки. Иными словами, произведение сМ есть максимальное значение, которое может принимать номер г, а п в знаменателе появляется из-за того, что мы рассматриваем граничные условия (7), при которых максимальное значение к есть п. Этим соотношением можно воспользоваться, чтобы получить оценку количества достоверных показателей Ляпунова, вычисляемых по смешанной схеме при заданных значениях Ах и АЬ. Сначала из (14) и (18) выразим к1, а затем воспользуемся нашим феноменологическим соотношением (19). В итоге получим
2сМ Ах2
Рис. 6. То же, что на рис. 5 при АЬ = 0.02. Условие монотонности (17) нарушено. Горизонтальные пунктирные линии отмечают возникновение паразитного возбуждения при к1 ~ 1.55 и точку минимума множителя роста при к0 ~ 1.98. Вертикальная пунктирная линия отмечает количество (¿1 = 99) достоверных показателей Ляпунова, вычисленное по формуле (20)
г1
аГ081П -
п
(20)
2ЩАЬ_
где |_ _| обозначает операцию нахождения целой части числа. Если при постоянном Ах менять значение АЬ, то с ростом АЬ ¿1 стремится к нулю. Если АЬ уменьшается, то г1 растет и достигает максимального значения, когда выражение под знаком агташ становится равным единице. Это соответствует попаданию минимума множителя роста в точку V = 1. При еще меньших значениях АЬ выполняется условие монотонности (17) и, следовательно, все вычисляемые показатели Ляпунова являются достоверными.
На рис. 6 значение г1 отмечено вертикальной пунктирной линией. Видно, что, хотя паразитное возбуждение возникает немного раньше, чем предсказывает формула (20), тем не менее оценка работает достаточно хорошо. Отметим, что формула (20) справедлива только для граничных условий вида (7). Однако аналогичные оценки можно получить и для других граничных условий.
Формула (20) получена в предположении, что в спектре решения уравнения возмущения (10), связанного с вычислением г-го показателя, есть чётко выраженное ведущее волновое число и оно пропорционально г, формула (19). Основанием для такого предположения служит уже упоминавшееся разбиение касательного пространства на активные и изолированные моды [5]. Действительно, изолированные моды в центральной части спектра образуют структуру, которая описывается формулой (19) (см. рис. 5 и 6). Однако из тех же рисунков видно, что если попадает в область активных мод в левом нижнем углу, то соотношение (19) и, следовательно, оценка (20) работают уже не так хорошо. Кроме того в работе [6] показано, что если, решая уравнение Гинзбурга-Ландау, взять достаточно большой шаг в пространстве, то изолированные моды вообще исчезают. Это значит, что оценка (20) вообще перестаёт быть адекватной. Отсюда следует, что этой оценкой следует пользоваться с осторожностью, так как существуют случаи, когда она не работает.
На рис. 7 показаны спектры показателей Ляпунова, вычисленные для уравнения Гинзбурга-Ландау, когда выполняется (кривая 1) и не выполняется (кривая 2) условие монотонности (17). Левые части обеих кривых совпадают, а в точке возникновения паразитного возбуждения, отмеченной вертикальной линией (20), на кривой 2 происходит излом. Правее точки излома она идёт практически горизонтально над кривой 1. Это значит, что соответствующие показатели Ляпунова оказываются существенно завышенными. Такой излом - типичное проявление паразитного возбуждения. Появление этой особенности на кривой спектра показателей Ляпунова может служить качественным критерием недостоверности найденных значений правее точки излома.
Таким образом, используя смешанную схему для вычисления полного спектра показателей Ляпунова, необходимо подобрать значения Ах и Дt так, чтобы выполнялось условие монотонности (17). Однако это условие требует, чтобы значение шага по времени было достаточно малым, что может существенно замедлить вычисления. В случае, когда полный спектр не нужен, можно допустить нарушение условия монотонности, выбрав большее значение АЬ.
Если среди собственных чисел матрицы Б есть такие, что И,е Щ = 0, то соответствующие ветви решения, полученные при помощи смешанной схемы, ведут себя правильно: все пространственные гармоники таких решений имеют множитель роста а = 1.
3.3. Явная схема Эйлера. Рассмотрим явную схему, которая возникает при а = 0. Эта схема имеет первый порядок точности локальной аппроксимации производной по времени и второй порядок - по координате. Квадрат множителя роста для явной схемы имеет вид
а2 = 1 - v(8 + 8*)+ V288*. (21)
Рис. 7. Показатели Ляпунова уравнения Гинзбурга-Ландау, вычисленные при помощи смешанной схемы с различными значениями АЬ: 1 - 0.013, 2 - 0.020. Значения параметров соответствуют рис. 5 и 6.
Если И,е d = 0, то все пространственные гармоники, кроме нулевой, неустойчивы. При этом а растёт с ростом V, так что максимальный множитель роста имеет коротковолновая гармоника с V = 1.
Пусть И,е d > 0. Тогда при V = 0 функция а2^) = 1, убывает до точки минимума vo
8 + 8*
vo =
288*
затем растёт и в точке V2
V2
8 + 8* 88*
(22)
(23)
снова проходит через единицу.
Когда V2 < 1, коротковолновые гармоники с номерами V > V2, неустойчивы и решение расходится. Из требования V2 > 1 можно найти стандартное условие устойчивости явной схемы [11,13]
Аг < Ах2
И,е d 2Й2
(24)
1.6 1.2 0.8 0.4 0
х........... 3 2 "
0
VI
0.2
0.4
V Уо
0.6
0.8
На рис. 8 кривая 1 иллюстрирует поведение а2^), когда условие устойчивости (24) нарушено. Вертикальная пунктирная линия отмечает точку V = V2, в которой множитель роста проходит через единицу.
Когда ищется решение одного уравнения, неустойчивость явной схемы проявляется очевидным образом: решение расходится, и программа перестаёт нормально работать. В случае же поиска показателей Ляпунова расходимость будет завуалирована вследствие периодических перенормировок решений. Поэтому вычисления будут происходить в нормальном режиме, даже когда условие устойчивости не выполняется или когда в спектре собственных чисел Б есть нули или чисто мнимые значение. При этом вместо ожидаемого для диффузионных задач доминирования длинноволновых структур, значения старших показателей Ляпунова будут определяться скоростью роста гармоник с наименьшими длинами волн.
На рис. 9 показано, как выглядят фурье-спектры столбцов матрицы Q, вычисленной для уравнения Гинзбурга-Ландау на одном из шагов процедуры поиска показателей Ляпунова при нарушенном условии устойчивости (24). Видно, что паразитное возбуждение коротковолновых гармоник возникает уже
Рис. 8. Множитель роста пространственных гармоник для различных 8: 1 - 2.3, 2 - 1.8, 3 - 0.8; явная схема при а = 0
Рис. 9. Фурье-спектры столбцов матрицы р. Явная схема. Параметры как на рис. 5, за исключением АЬ = 0.007. Условие устойчивости (24) нарушено
при г = 1, то есть ни один из вычисленных таким образом показателей Ляпунова не является достоверным.
Как обсуждалось в разделе 3.2, для поиска показателей Ляпунова недостаточно одной только устойчивости численной схемы. Необходимо, чтобы множитель роста убывал монотонно на интервале от нуля до единицы. Это имеет место, когда точка минимума функции а2^) находится справа от единицы, то есть vo > 1. Отсюда следует условие монотонности
АЬ< Ах2( (25)
На рис. 8 кривая 2 демонстрирует поведение функции а2^) когда числовая схема устойчива, а условие монотонности нарушено. Видно, что, хотя множитель роста всюду меньше единицы, начиная с точки VI для каждой гармоники существует пара с одинаковым а. При вычислении показателей Ляпунова это приводит к паразитному возбуждению коротковолновых гармоник, что, как обсуждалось в разделе 3.2, делает недостоверным значения младших показателей. Кривая 3 на этом рисунке представляет случай, когда условие монотонности (25) выполняется, так что все показатели Ляпунова являются достоверными.
Найдём количество вычисляемых достоверно показателей Ляпунова в случае, когда условие устойчивости выполнено, а условие монотонности нарушено. Воспроизводя рассуждения раздела 3.2, найдём номер гармоники VI, имеющей такой же показатель роста, что и коротковолновая гармоника при V = 1
V = 8+8* - 1 (26)
Далее, используя формулы (14) и (19), получим, что количество достоверных показателей Ляпунова можно оценить как
г\
2сМ Ах2Яе Щ -аггаш д —, , 1|0— 1
п V 2АЬЩ2
(27)
Зафиксируем Ах и будем увеличивать АЬ. При этом выражение под знаком корня может обратиться в нуль и, следовательно, тоже станет равным нулю. Как видно из неравенства (24), при дальнейшем увеличении шага численная схема становится неустойчивой и вычисление показателей Ляпунова невозможно. Если уменьшать АЬ, то подкоренное выражение достигнет единицы. При этом значение будет максимальным. Неравенство (25) показывает, что при дальнейшем уменьшении шага выполняется условие монотонности и все вычисляемые показатели являются достоверными. Отметим, что замечания, сделанные в разделе 3.2 относительно границ применимости оценки (20), справедливы и для формулы (27).
На рис. 10 построены спектры показателей Ляпунова, вычисленные для уравнения Гинзбурга-Ландау с использованием явной схемы. Кривая 1 изображает «правильный» ляпуновский спектр, вычисленный при АЬ, удовлетворяющем условию монотонности (25). Если немного увеличить АЬ (кривая 2), то условие монотон-
ности перестаёт выполняться, а условие устойчивости (24) выполняется. В точке %\ кривая спектра претерпевает излом по причине паразитного возбуждения коротковолновых гармоник. Видно, что паразитное возбуждение проявляется точно также, как и в смешанной схеме (ср. с кривой 2 на рис. 7). В этом случае достаточно хорошей оценкой количества достоверных показателей Ляпунова является величина %\, вычисляемая по формуле (27). При дальнейшем увеличении ДЬ нарушается условие устойчивости (24). Как уже обсуждалось выше, в этом случае все найденные показатели Ляпунова оказываются недостоверными. Иллюстрирующая этот случая кривая 3 проходит значительно выше двух других кривых.
Таким образом, если среди собственных чисел матрицы О есть такие, что Ие d = 0, то явную схему использовать нельзя, так как показатели Ляпунова, соответствующие этой ветке решения, будут вычислены некорректно - их значения будут завышены. Если же Ие d > 0, то, применяя явную схему, необходимо обязательно добиться выполнения условия устойчивости (24). Для вычисления полного спектра требуется также выполнение условия монотонности.
3.4. Неявная схема. Пусть теперь о =1. Получающаяся в этом случае схема называется неявной. Она, так же как и явная схема, имеет первый порядок точности локальной аппроксимации производных по времени и второй порядок - по координате. Множитель роста для этой схемы имеет вид
а2 = (1 + у(6 + б*) + у2бб* )-1. (28)
При Ие d ^ 0 множитель роста в точке V = 0 равен 1 и убывает при увеличении V. Все особые точки этой функции находятся либо на комплексной плоскости, либо на отрицательной полуоси. Таким образом, при любом б множитель роста убывает монотонно на отрезке от нуля до единицы, что иллюстрирует рис. 11.
Следовательно, неявная схема позволяет находить показатели Ляпунова при любой комбинации Дх и ДЬ. Тем не менее не следует забывать, что точность вычисления существенным образом зависит от значения шагов численной сетки. Так как схема имеет первый порядок по Ь, высокую точность можно достичь, если взять ДЬ < Дх2, что по сути дела накладывает на ДЬ такое же сильное ограничение, что и условия монотонности (17) и (25). Отчасти это снижает преимущества неявной схемы, так как мы снова вынуждены брать достаточно малый шаг по времени. Ещё один недостаток
-5.0 -15.0 -25.0 -35.0
-
-......... 3
1
50
100
<1
150
Рис. 10. Показатели Ляпунова уравнения Гинзбурга-Ландау, вычисленные при помощи явной схемы с различными значениями Д!: 1 -0.003, 2 - 0.004, 3 - 0.007. Значения параметров соответствуют рис. 9
а2 0.8 0.6 0.4 0.2
2
\ —
------------------
0
0.2
0.4
0.6
0.8
Рис. 11. Множитель роста пространственных гармоник, неявная схема при о =1 и различных значениях б: 1 - 4.0, 2 - 1.0
неявной схемы - «неправильное» поведение множителя роста при И,е I = 0. Он убывает, и, следовательно, только несколько старших показателей могут будут вычислены в этом случае достаточно хорошо, тогда как все другие будут иметь заниженные значения.
Заключение
Проанализированы численные методы, используемые для нахождения показателей Ляпунова для распределённых систем. По результатам выполненного анализа можно сформулировать следующие практические рекомендации.
Для вычисления показателей Ляпунова требуется периодически выполнять ор-тогонализацию матрицы решений уравнений для инфинитезимальных возмущений. При этом предпочтительнее всего использовать метод на основе преобразования Ха-усхолдера, который имеет очень высокую точность даже когда интервал между орто-гонализациями выбран достаточно большим и матрицы плохо обусловлены. Однако при этом нужно обязательно проверять значения накапливаемых для вычисления показателей величин. Они должны быть хотя бы на несколько порядков больше значения машинного эпсилон.
Методы Грама-Шмидта, классический и модифицированный, могут работать с высокой точностью при условии, что матрицы хорошо обусловлены. Этого можно добиться, выполняя процедуры ортогонализации достаточно часто. Кроме того, эти методы хорошо работают, когда количество вычисляемых показателей много меньше размерности фазового пространства. Тем не менее, чтобы получить гарантированно хороший результат, применяя эти методы, весьма желательно проверить погрешность ортогонализации, по крайней мере, на нескольких первых шага счёта.
Если показатели Ляпунова ищутся для уравнения реакционно-диффузионного типа, и уравнения для возмущений решаются методом конечных разностей, то выбор слишком большого шага по времени по сравнению с шагом по координате может вызвать паразитное возбуждение коротковолновых пространственных гармоник, что, в свою очередь, приводит к грубым ошибкам при вычислении показателей. Надёжный способ избежать этого - использовать чистую неявную схему, которая не подвержена этому эффекту. Чтобы воспользоваться смешанной или явной схемами для вычисления полного спектра показателей, требуется подобрать шаги по времени и по координате таким образом, чтобы выполнялись условия монотонности, разные для разных схем. Если же полный спектр не требуется, то можно допустить нарушение условия монотонности, взяв шаг по времени немного большим.
Однако, если в спектре собственных чисел матрицы диффузии системы имеются нули или чисто мнимые значения, то наилучшим методом является смешанная схема, так как только она обеспечивает безразличное равновесие для соответствующих пространственных гармоник. Сформулированные выше требования отсутствия влияния паразитного возбуждения должны также выполняться, чтобы обеспечить правильное поведение гармоник, соответствующих другим веткам решения.
Анализ паразитного возбуждения был выполнен для реакционно-диффузионной системы общего вида, решаемой при помощи метода конечных разностей. При этом
паразитное возбуждение возникало при использовании явной или частично явной (то есть смешанной) схемы. По всей вероятности, подобный эффект может иметь место и в других случаях, когда решение ищется с применением явной схемы. В частности, его следует ожидать при использовании метода Рунге-Кутты для решения системы большого числа связанных обыкновенных дифференциальных уравнений.
Автор выражает признательность С.П. Кузнецову за полезное обсуждение этой работы.
Работа выполнена при поддержке гранта РФФИ-ННИО № 08-02-91963. Библиографический список
1. Оселедец В.И. Мультипликативная эргодическая теорема. характеристические показатели Ляпунова динамических систем // Труды Моск. матем. об-ва. 1968. T. 19. С. 197.
2. Eckmann J.P., Ruelle D. Ergodic theory of chaos and strange attractors// Rev. Mod. Phys. 1985. Jul. Vol. 57, № 3. P 617.
3. Schaumlöffel K.-U. Multiplicative ergodic theorems in infinite dimensions // Lyapu-nov Exponents. Springer Berlin / Heidelberg, 1991. T. 1486/1991. Lecture Notes in Mathematics. C. 187.
4. Robinson J.C. Finite dimensional behavior in dissipative partial differential equations // Chaos. 1995. Vol. 5. P. 330.
5. Yang H.-L., Takeuchi K.A., Ginelli F., Chate H., Radons G. Hyperbolicity and the effective dimension of spatially-extended dissipative systems // Phys. Rev. Lett. 2009. Vol. 102. P. 074102.
6. Kuptsov P. V., Parlitz U. Strict and fussy modes splitting in the tangent space of the Ginzburg-Landau equation // Phys. Rev. E. 2010. Vol. 81. P. 036214.
7. Benettin G., Galgani L., Giorgilli A., Strelcyn J.M.Lyapunov characteristic exponents for smooth dynamical systems and for hamiltonian systems: A method for computing all of them. Part I: Theory. Part II: Numerical application // Meccanica. 1980. Vol. 15. P. 9.
8. Parker T.S., Chua L.O. Practical numerical algorithms for chaotic systems. SpringerVerlag, 1989. P. 348.
9. Golub G.H., van Loan C.F. Matrix computations. Third Edition. The Johns Hopkins University Press, Baltimore, MD. 1996. P. 694.
10. Geist K., Parlitz U., Lauterborn W. Comparision of different methods for computing Lyapunov exponents // Prog. Theor. Phys. 1990. Vol. 83, № 5. P. 875.
11. Numerical recipes in C / W.H. Press, S.A. Teukolsky, W.T. Vettering, B.P. Flannery. Cambridge University Press, 1992. P. 994.
12. Aranson I.S., Kramer L. The world of the complex Ginzburg-Landau equation // Rev. Mod. Phys. 2002. Vol. 74. P. 99.
13. Калиткин Н.Н. Численные методы: Учебное пособие. М.: Наука, 1978. C. 512.
Саратовский государственный Поступила в редакцию 26.02.2010
технический университет После доработки 14.05.2010
COMPUTATION OF LYAPUNOV EXPONENTS FOR SPATIALLY EXTENDED SYSTEMS: ADVANTAGES AND LIMITATIONS OF VARIOUS NUMERICAL METHODS
P. V. Kuptsov
Problems emerging in computations of Lyapunov exponents for spatially extended systems are considered. We concentrate on the incorrect orthogonalization of high sized ill conditioned matrices appearing in course of the computation, and on large errors emerging under certain conditions if the finite difference numerical method is applied to solve equations. The practical guidelines helping to avoid the mentioned problems are represented.
Keywords: Lyapunov exponents, spatio-temporal chaos, finite difference method, complex Ginzburg-Landau equation, QR-decomposition.
Павел Владимирович Купцов - родился в Саратове (1972), окончил физический факультет Саратовского государственного университета (1994), защитил кандидатскую диссертацию (1998). В настоящее время - доцент кафедры технической кибернетики и информатики Саратовского государственного технического университета, а также сотрудник научной группы теоретической нелинейной динамики в Саратовском филиале Института радиотехники и электроники РАН. Область научных интересов - сложные колебания в распределённых системах, хаотическая динамика высокой размерности, численные методы.
410054 Саратов, ул. Политехническая, 77 Саратовский государственный технический университет E-mail: p.kuptsov@rambler.ru
