Вычисление периодических решений для нелинейных автономных систем обыкновенных дифференциальных уравнений

Человеческое время не обращается, но кругу, а бежит по прямой. И в этом причина, по которой человек не может быть счастлив, ибо счастье есть жажда повторения.

Милан Кундера.

Невыносимая легкость бытия Как правило, при исследовании периодического решения автономной системы ОДУ сам период неизвестен и подлежит определению одновременно с решением системы ОДУ. Рассмотрим вкратце основные идеи построения численных методов, позволяющие найти периодическое решение нелинейной задачи. Как и ранее, рассмотрим автономную систему вида.

с условием u (0) = u(Т), причем период Т подлежит определению. Сделаем в системе (8.14) замену независимой переменной: t = 7 т. Тогда система (8.14) превращается в автономную систему вида.

с условием периодичности, принимающим вид u (0) = и (1). Теперь значение периода, вошедшее в правую часть системы ОДУ, можно объявить одной из переменных задачи. Получить разностные формулы можно следующим образом. Заменим выражение (8.15) на разностное соотношение, например, по формуле трапеций:

Из условий периодичности, кроме того, имеем u₀ = u_v. Уравнения (8.16) вместе с граничным условием (периодичности) образуют систему: имеется К х (N + 1) нелинейных уравнений вида (8.16) относительно К х (JV + 1) + 1 неизвестных. Неизвестными системы являются Г, u₀, Uj, …, Uv_iЗдесь К — число уравнений в дифференциальной системе (8.14). Возможные подходы к решению такой системы нелинейных уравнений таковы. Представим (8.16) с учетом граничных условий в виде.

Зафиксируем значение одной из переменных задачи. Такой подход связан с некоторым риском, так как не всегда априори бывает ясно, соответствует ли такое значение какому-либо периодическому решению. Иногда существование такого фиксированного значения на границе периода бывает очевидным из физических соображений. Так, производная периодической функции на периоде обязана обращаться в нуль и т. п. Иногда выбор такого фиксированного значения сопряжен с известным риском — выбранное значение не будет соответствовать какому-либо периодическому решению задачи. Как правило, поиску периодического решения соответствует рассмотрение свойств задачи Коши, и некоторая априорная информация у исследователя уже имеется.

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

Без ограничения общности можно положить =0. Тогда число неизвестных, подлежащих определению, равно числу уравнений. Дальнейшее решение системы связано с методом квазилинеаризации (по Ньютону)¹ и решением линейной системы. Отметим, что вид матрицы Якоби похож на тот, что получается при применении описанного в предыдущем параграфе метода прогонки. Единственное отличие — матрица линеаризованной системы получается почти блочно-диагональной (от нулевых отличаются лишь блоки на главной диагонали и над ней, блок в левом нижнем углу матрицы и заполненный последний столбец). Применению алгоритмов работы с блочными матрицами препятствует структура блоков в первом столбце: их размер получается (К — 1) х К. Тем нс менее для матриц подобной структуры можно построить эффективно работающие алгоритмы, основанные на методе Гаусса^[1][2][3][4].

В связи с быстрым развитием вычислительной техники многие ограничения на размеры памяти и быстродействие компьютеров становятся несущественными и можно применить более «честный» вариант решения задачи. Задаваясь каким-либо начальным приближением (обычно его получают из решения задачи Коши), запишем формулы итерационного процесса:

Вектор х составлен из величин Т, u₀, и_и…, u_v__t. Ввиду того что размеры векторов х и F не совпадают, матрица Якоби системы является прямоугольной, и обратная матрица к ней, вообще говоря, не существует. Последняя формула тем не менее имеет право на существование и является некоторым обобщением лметода Пыотона, если понимать матрицу в приведенной выше формуле как псевдообратную^{3 4}. Для переопределенных линейных систем пссвдообратная матрица дает наилучшее решение в смысле наименьших квадратов.

Приведем алгоритм вычисления псевдообратной матрицы — метод Гревилля*. Пусть А — (т х /г)-магрица, а_А, — ее k-i столбец. Введем обозначение A_fi = (a_jt …, а/_г) — матрица, образованная первыми k столбцами матрицы А. Пусть пссвдообратная к, А матрица есть А⁺, а Ь_к, — последняя строка в матрице А?. Индекс k в приведенных формулах пробегает все значения от 1 до п. В соответствии с введенными обозначениями Aj = а_1? А_п = А. Тогда если А, = a_t = 0, то А[ = 0, иначе.

Для определения строки используется следующий алгоритм:

• 1) если c_k = a_k — А_к_^А_кФ 0, то b_fi = с? =(а^ - A^_₁d_/(,)⁺, знак «+» (верхний индекс в последней формуле) означает псевдообращение (см. формулу для aj^-);
• 2) если с_к = а_к- A*_, d* = 0, то b* =(l + d[d_A)d[A^.

Метод Гревилля не требует вычисления обратных матриц, поэтому может быть применен и для обращения квадратных матриц.

Рассмотрим альтернативные подходы к решению периодической краевой задачи (8.15). Простейшим примером такой альтернативы может служить метод пристрелки (стрельбы) для решения краевых задач. К сожалению, в большинстве учебников по численным методам содержится лишь краткое описание основной идеи метода. Исключение здесь составляет книга Р. П. Федоренко¹, с которой необходимо ознакомиться перед реализацией метода пристрелки для крупной вычислительной задачи. Рассмотрим дальнейшее развитие идеи метода пристрелки^[5][6].

Метод параллельной стрельбы. Как и ранее, нас интересуют периодические решения автономной нелинейной системы (8.14). Вообще говоря, рассмотрим чуть более сложный случай системы, зависящей от параметров:

После замены переменных, как и ранее, имеем.

Делим отрезок [ 0; 1 ] на несколько частей точками т_к 0 = т₀ < т_{ < … < х_м = = 1. В каждой точке задаем вектор начальных условий и (Т/) = V/ = (v_ib а_/2, …, v_jN) и решаем численно задачу Коши с помощью какого-либо метода, описанного в гл. 5, на отрезке |т^; т*₊₁]. Потребуем, чтобы решение было непрерывным. Тогда должно выполняться равенство.

В равенстве (8.17) через и (Т/ | V/ _v Г, а) обозначено значение численного решения задачи Коши на правом конце отрезка [Т/ _t; т/]. Это равенство выполняется при 1 < /< М — 1. Кроме того, в силу периодичности имеем.

Равенства (8.17) и (8.18) образуют нелинейную систему, включающую N х М скалярных уравнений для определения М х N + 1 неизвестного.

Напомним, что период нам пока неизвестен и подлежит определению одновременно с решением.

Дальнейшая процедура поиска решения похожа на рассмотренную выше. Либо при т = 0 фиксируется одна из координат вектора v₀, например полагается нулем из физических соображений, либо применяется обобщение метода Ньютона с алгоритмом Гревилля.

Преимущество описанного вкратце выше метода перед методами линеаризации — меньшая размерность получающейся нелинейной системы уравнений. В качества платы за упрощение несколько возрастает объем вычислительной работы. Действительно, каким бы способом ни находилось решение, нам надо знать значения всех частных производных, входящих в матрицу Якоби,.

Эти матрицы могут быть вычислены путем интегрирования задачи Коши для уравнения в вариациях.

где Y — квадратная матрица размером Nx N. Выбирая в качестве начальных данных для решения задачи Коши для уравнения в вариациях Y (i/) = Е, на другом конце отрезка в результате численного решения задачи получаем нужное значение матрицы Якоби. Па практике система в вариациях обычно решается (численно) одновременно с рассматриваемой нелинейной системой ОДУ.

Трудности при реализации методов могут возникнуть, когда при данном периоде или данном значении параметра решение задачи не единственно.

• [1] Ко. иаггщ Л. Указ, соч.; Беллман Р., Калаба Р. Квазилинеаризация и нелинейные уравнения. М.: Мир, 1968.
• [2] Такие методы рассматриваются в работе: Холодлиок М. и др. Указ. соч.
• [3] См.: Гонтмахер Ф. Р. Указ. соч.
• [4] См.: Там же.
• [5] Федоренко Р. П.

  Введение

  в вычислительную физику. Долгопрудный: Интеллект, 2009.

• [6] Подробно оно описано в работе: Холодниок М. и др. Указ. соч. Также подробное описание включено в книгу: Карпов В. E. t Лобанов А. И. Численные методы, алгоритмы и программы.

  Введение

  в распараллеливание. М.: Физматкнига, 2014.