Новый коллокационный интегратор для решения задач динамики. I. Теоретические основы
Предлагается новый коллокационный интегратор на разбиении Лобатто для численного решения смешанных систем дифференциальных уравнений динамики первого и второго порядков. Излагается общая теория коллокационных интеграторов, из которой выводятся основные формулы нового интегратора.
New collocation integrator for solving dynamic problems.pdf Введение Коллокационные интеграторы [1-4] представляют собой простой, изящный, удобный и вместе с тем мощный инструментарий для решения дифференциальных уравнений динамики. Идея коллокационных интеграторов весьма прозрачна и для ее понимания совершенно не требуется знаний других интеграторов [3], например, Рунге - Кутты, Грэгга - Булирша - Штера или Адамса. Достаточно лишь иметь представление, что такое дифференциальные уравнения, интегрирование и интерполирование. Поэтому вполне обоснованно рассматривать коллокационные интеграторы в целом как самостоятельный класс, из которого наиболее известными представителями являются неявные коллокационные интеграторы Рунге - Кутты. Примечательной особенностью коллокационных интеграторов является то, что их теоретическая основа, как и программная реализация, универсальна для любого порядка [5]. Практически порядок определяется разбиением на шаге, а именно количеством и спецификой распределения узловых значений, через которые выражаются все остальные константы интегратора. Кроме того, на разбиениях гауссовых квадратур Лежандра и Лобатто коллокационные интеграторы становятся геометрическими [4]: симметричными и орбитально устойчивыми**, а на разбиении Лежандра еще и симплектическими. Следует также отметить, что, в отличие от других интеграторов, коллокационные позволяют на каждом шаге легко конструировать приближенное аналитическое решение, чем удобно пользоваться для частого вывода результатов на плотной временной сетке. В настоящей работе предлагается новый коллокационный интегратор Lobbie на разбиении Лобатто. Фактически его прототипом стал широко используемый в динамической астрономии интегратор Эверхарта [6, 7]. Точнее говоря, он является результатом кардинальной редакции своего предшественника, хотя в итоге теория, алгоритмизация и программный код нового интегратора могут лишь только концептуально напоминать об авторской версии прославленного интегратора Эверхарта. Кратко излагается общая теория коллокационных интеграторов применительно к решению дифференциальных уравнений первого и второго порядков. Приводятся частные примеры, в том числе интегратор Эверхарта. Далее выводятся основные формулы интегратора Lobbie, а также отмечаются особенности в их программной реализации. Описывается программная процедура Lobbie. 1. Коллокационные методы 1.1. Дифференциальные уравнения первого порядка Пусть динамическое состояние описывается во времени векторным дифференциальным уравнением первого порядка (1) при известном начальном динамическом состоянии на момент : (2) Здесь штрих обозначает полную производную по времени; - известная вектор-функция времени и динамического состояния. Необходимо определить динамическое состояние системы на момент : (3) где - малый параметр (длина шага интегрирования). Представим приближенно решение уравнения (1) в виде полинома (4) который должен точно удовлетворять уравнению на некоторые промежуточные моменты времени (точки коллокаций) (5) а также начальному условию (2) (6) Тогда решение (3) в соответствии с (4) приближенно определяется как (7) Геометрический смысл соотношений (5) состоит в том, что касательные к полиному в точках коллокаций должны совпадать (коллоцировать) с направлениями векторного поля, создаваемого функцией дифференциального уравнения (рис. 1). При этом значения самого полинома могут заметно отличаться от точного решения. Рис. 1. Условия коллокаций для дифференциального уравнения при начальном условии . Стрелками указаны направления векторного поля в точках , где . Черная кривая - коллокационный полином (приближенное решение) на разбиении Лобатто (черные точки) при ; пунктирная кривая - точное решение ; точки коллокаций на плоскости выделены белым цветом Условия коллокаций (5) можно рассматривать как условия Лагранжа, налагаемые на производную от полинома (4): (8) которая в таком случае выступает в роли полиномиального интерполянта функции по безразмерной переменной . Согласно (5), условия Лагранжа на узловых значениях безразмерной переменной можно представить в виде (9) Формируя из условий (9) полиномиальный интерполянт и затем интегрируя соотношение (8) по на отрезке , учитывая (6) и (7): получаем приближенное решение Интерполянт (8) конструируется из промежуточных приближенных решений (9), т.е. . Каждое -е из этих решений также определяется путем интегрирования соотношения (8) по , но на отрезке . Таким образом, коллокационный метод для решения дифференциального уравнения (1) можно представить как сводку формул (10) Промежуточные решения в (10) задаются неявным образом через нелинейные уравнения. По этой причине все коллокационные методы - неявные . Уравнения решаются, как правило, методом простых итераций в модификации Зейделя, т.е. поочередно уточняя промежуточные решения на каждой итерации. Фактически итерационное решение нелинейных уравнений является ядром любого коллокационного интегратора, и его эффективность во многом зависит от того, насколько удачно организован итерационный процесс. Порядок метода обусловлен количеством точек коллокаций , а также спецификой их распределения. Принцип коллокаций позволяет получить практически любой порядок. Так, для произвольного распределения (например, равномерного), по меньшей мере, [3]. Однако при использовании узловых значений гауссовых квадратур Лежандра, Радау и Лобатто порядок можно повысить до , и соответственно [3, 4, 8, 9]. Эти узловые значения являются решениями алгебраических уравнений (11) Коллокационные методы примечательны тем, что фактически они дают аналитическое решение на шаге в форме коллокационного полинома (а не только решение на конце шага и несколько промежуточных ), что очень удобно для вывода приближенных решений на плотной временной сетке. Впрочем, следует помнить, что порядок точности внутри шага понижается до [3]. Кроме того, наличие интерполянта позволяет путем экстраполирования получать достаточно хорошие начальные приближения функции на следующем шаге: для дальнейшего итерационного определения промежуточных решений (10). Следует особо отметить, что на разбиениях Лежандра и Лобатто (11) коллокационные методы становятся геометрическими [4]: симметричными и орбитально устойчивыми, а на разбиении Лежандра еще и симплектическими. Казалось бы, разбиение Лежандра предпочтительней, поскольку при одном и том же количестве точек коллокаций его порядок на два выше. Однако, несмотря на то, что порядок метода на разбиении Лобатто ниже, он работает немного быстрее. Так, на каждом шаге при итерациях для определения промежуточных решений количество вычислений функции правых частей составляет , тогда как на разбиении Лежандра - . Кроме того, интерполянт правой части уравнения строится на всем отрезке интегрирования , поскольку и , в отличие от разбиения Лежандра, где все узловые значения находятся внутри отрезка. Следовательно, предиктор метода на разбиении Лобатто будет лучше, что очень важно для его программной реализации. В качестве простых примеров коллокационных методов можно привести хорошо и давно известные методы Рунге - Кутты: явный и неявный Эйлера (Radau I & II: , ); средней точки (Legendre: , ); трапеции (Lobatto: , ); Симпсона (Lobatto: , ). Если в качестве интерполянта функции дифференциального уравнения принять полином Лагранжа то коллокационный метод (10) принимает классическую форму метода Рунге - Кутты [1-4]: (12) где постоянные выражаются через интегралы от базисных функций Лагранжа: Вскоре после первых работ по коллокационным методам Рунге - Кутты [1, 2] Э. Эверхарт [6, 7] предложил в качестве интерполянта использовать полином в каноническом виде Здесь . Простая форма интерполяции позволяет получить достаточно простую форму приближенного решения (13) Однако, чтобы выразить коэффициенты полинома через коллокационные значения функции , автор прибегает к разделенным разностям интерполяционного полинома Ньютона , которые непосредственно определяются из коллокационных значений. Между тем коэффициенты канонического полинома выражаются через разделенные разности посредством линейных соотношений где постоянные вычисляются из узловых значений разбиения по рекуррентным формулам 1.2. Дифференциальные уравнения второго порядка Предположим теперь, что динамическая система описывается векторным дифференциальным уравнением второго порядка (14) при известном начальном динамическом состоянии на момент : (15) Необходимо определить динамическое состояние системы на момент : (16) Здесь и будем рассматривать как векторы координат и скоростей соответственно. Представим приближенно решения для координат и скоростей в виде полиномов (17) которые должны удовлетворять уравнению (14) на моменты коллокаций : (18) а также начальному условию (15): (19) Тогда решения (16) в соответствии с (17) приближенно определяются как (20) Согласно (18), условия Лагранжа для интерполянта функции (21) можно представить в виде (22) Формируя из условий (22) полиномиальный интерполянт и затем интегрируя соотношение (21) по на отрезках и , получаем сводку формул для коллокационного метода применительно к дифференциальному уравнению (14): (23) где подынтегральный полином зависит также от промежуточных решений и . 1.3. Смешанные системы дифференциальных уравнений первого и второго порядков Задачу (14) и (15) можно представить в виде (1) и (2): (24) Несмотря на то, что обе задачи описывают одну и ту же динамическую систему, их численные решения в соответствии с принципом коллокаций будут принципиально разными. Так, согласно (10), для альтернативной задачи будем иметь приближенные решения (25) где полином интерполирует правую часть уравнения для вектора координат. Точность решений (25) (как для координат, так и для скоростей) - одного порядка . Однако в (23) решение для координат вследствие двойного интегрирования имеет порядок , иначе говоря, координаты в (23) определяются точнее, нежели в (25). Тем не менее пользователь порой вынужден обращаться к представлению динамической системы в виде (24) и тем самым прибегать к интегратору для уравнений первого порядка (25), вполне осознавая, что это неизбежно влечет значительное понижение точности приближенного решения. Такая необходимость возникает, когда уравнение динамической системы (14) дополняется уравнениями первого порядка для каких-либо вспомогательных динамических величин. Например, смешанные системы применяются для исследования динамического хаоса [10], а также для линеаризации, регуляризации и стабилизации уравнений динамики [11-14]. Так, чтобы привести все уравнения смешанной системы к единому порядку, вместо (14) используются уравнения (24). Альтернативный вариант - это дифференцирование дополнительных уравнений, с тем чтобы вся система в целом имела вид (14). Хотя для сложных динамических систем взятие производных от функций дополнительных уравнений часто практически невозможно. Естественный и эффективный подход к решению смешанной системы уравнений - это применение гибридного интегратора. Если для моделирования динамической системы необходимо совместно с (14) решать векторное уравнение первого порядка для вспомогательных динамических величин : (26) при начальном условии , сводку (23) следует дополнить формулами . (27) Здесь интерполянт функции конструируется по аналогии с интерполянтом из промежуточных решений , и . Таким образом, гибридный коллокационный метод для уравнений (14) и (26) будет иметь вид (23) и (27). 2. Интегратор Lobbie Применительно к смешанной системе дифференциальных уравнений (14) и (26): (28) примем в качестве интерполянтов полиномы Ньютона (29) на разбиении Лобатто . Здесь . Разделенные разности в (29) определяются из узловых значений функций и по рекуррентным формулам (30) Подставляя интерполянты (29) в (23) и (27), а также учитывая, что и , получим коллокационный метод для системы (28) в виде (31) где (32) Обозначим -кратный интеграл от -й базисной функции Ньютона как . (33) Если величины рассматривать как элементы некой матрицы размера , то константы коллокационного метода - интегралы от базисных функций Ньютона (32) - будут составлять первые ее два столбца: Между тем элементы матрицы для произвольного значения вычисляются построчно через рекуррентные соотношения: (34) Нелинейные уравнения (31) относительно , и решаются на каждом шаге методом простых итераций в модификации Зейделя. До итераций известны решения , и , а также разделенные разности на первой точке коллокаций . В начале итерационного процесса на второй точке коллокаций определяется группа решений , и , а из них - и , по которым уточняются разделенные разности и в соответствии с рекуррентными формулами (30). Затем таким же образом уточняются разделенные разности и на третьей точке коллокаций и т.д. После поочередного уточнения всех разделенных разностей на шаге итерация повторяется. Итерационный процесс продолжается до тех пор, пока не выполнится неравенство (35) Здесь - решение на предыдущей итерации; - малая величина, которая задает точность сходимости. Впрочем, количество итераций на шаге можно ограничить заданным числом, не учитывая выполнение условия (35). В таком случае нужно иметь в виду, что получаемое решение уже не будет соответствовать задаваемому порядку, а метод лишится геометрических свойств. В качестве начальных приближений разделенных разностей и принимаются их оценки, получаемые из рекуррентных формул (30) по экстраполированным значениям функций и : (36) На первом шаге при отсутствии таких оценок итерационный процесс начинается буквально с нуля, т.е. при нулевых значениях разделенных разностей, и продолжается до обязательного выполнения условия (35). При многократном выводе приближенных решений на плотной временной сетке удобно и целесообразно пользоваться подручным коллокационным полиномом, накрывающим на определенном шаге временные моменты сетки, вместо того чтобы выполнять пошаговое интегрирование на каждый из этих моментов с очень мелким шагом. Коллокационные полиномы для системы (28) на произвольный момент времени можно представить в виде (37) где величины и вычисляются по формулам (34). Альтернативный способ получения временных рядов приближенных решений - это полиномиальное интерполирование на промежуточных решениях, например, типа (29): (38) Однако точность интерполяционных полиномов (38) между точками коллокаций оказывается значительно ниже, чем точность коллокационных полиномов (37), получаемых прямым интегрированием полиномов (29). Длина шага как параметр интегратора задается пользователем. Однако возможен режим автоматического выбора длины шага в процессе пошагового интегрирования. Величина выбирается таким образом, чтобы сохранялась приближенная оценка члена ряда Тейлора -го порядка для вектора скорости [5]: (39) т.е. если рассматривать член ряда как главный член погрешности, длина шага будет определяться как для метода порядка , хотя порядок на разбиениях гауссовых квадратур значительно больше (не менее чем в 2 раза). К сожалению, получить на шаге оценку более высокого порядка практически не удается. Предположим, оценка (39) должна быть равна постоянной величине , задаваемой пользователем. Поскольку то требуемую для обеспечения величины длину шага можно оценить как (40) После получения на текущем шаге оценки (40) интегрирование не повторяется с новой длиной шага , но она используется для следующего шага по формуле (41) где и - номера текущего и следующего шагов соответственно. Заметим, что автоматический выбор длины шага предполагает модификацию предиктора (36), а именно (42) Алгоритм (41) должен быть эффективным, если и для любого различаются несущественно. Заметные изменения в последовательности оценок (41) происходят, когда функция ведет себя нерегулярно, например, при интегрировании вблизи ее сингулярностей. Чтобы избежать большого различия между и , к отношению следует применять демпфирование, т.е. наложить на него ограничения: (43) Величина - это диапазон допустимого изменения значения . Для изменения в пределах одного порядка . При невыполнении левого неравенства (43) шаг повторяется. Начальный размер шага определяется из оценки (41) для [5]: (44) Здесь - некоторая малая величина. Если мала настолько, что в компьютерной арифметике , то она увеличивается на порядок и оценка (44) выполняется снова. Оценка стартового размера шага (44) в действительности адекватна только для метода Эйлера первого порядка. Поэтому для метода любого другого порядка стартовый шаг будет значительно меньше ожидаемого. Тем не менее пошаговый процесс интегрирования начинается с оценки (44), но в дальнейшем с использованием алгоритма (41) и с учетом ограничений (43) размер шага постепенно будет выходить на должную величину штатного режима работы интегратора с автоматическим выбором длины шага. Последний шаг выявляется из условия , (45) где - длина всего интервала интегрирования; - момент времени на -м шаге. Тогда, чтобы выйти на конечный момент времени , размер последнего шага задается принудительно как и переопределяется его отношение к длине текущего шага : . 3. Процедура Lobbie Интегратор реализован на процедурном языке Фортран до 32-го порядка для компьютерной арифметики с двойной и четверной точностью (double & quadruple precision). Вызов программной процедуры интегратора Lobbie выполняется командой call lobbie(x, y, z, ts, tf, step, etol, nxy, nz, ns, ni, nst, ncf, fun). Здесь x, y, z - массивы интегрируемых переменных соответственно: на входе значения на начальный момент времени (ts), на выходе значения на конечный момент времени (tf); step - начальный размер шага интегрирования : при автоматическом выборе (41) на выходе значение (размер предпоследнего шага), при нулевом значении step величина задается интегратором, согласно оценке (44); etol - : при нулевом значении - режим постоянного шага интегрирования; nxy и nz - размерности массивов x, y ( ) и z ( ) соответственно; ns - количество узловых значений ; ni - максимальное число итераций на шаге для решения нелинейных уравнений (31) относительно , и ; nst и ncf - количество выполненных шагов и обращений к процедуре fun для вычисления функций и на всем интервале интегрирования. Процедура fun задается как subroutine fun(t, x, y, z, f), где t - текущий момент времени ; x, y, z - массивы интегрируемых переменных со значениями на момент ; f - выходной массив значений функций и размерности . Вычислительный процесс внутри процедуры Lobbie выполняется в соответствии с блок-схемой, представленной на рис. 2. Опишем коротко основные этапы пошагового интегрирования в случае переменного шага с автоматическим выбором его стартовой величины, а также при возрастающем изменении времени . I. Из прилагаемого к процедуре блока данных считывается массив узловых значений и рекуррентно вычисляются константы интегратора (32) с использованием (34). Оценивается стартовая величина шага (44). Из выполнения условия устанавливается, что стартовый шаг является последним, и тогда задается величина . Рис. 2. Блок-схема программной процедуры Lobbie II. Итерационно определяются промежуточные решения , и (31), а также функции и (28) вместе с разделенными разностями и . Итерации завершаются при выполнении условия (35) или по достижении максимально допустимого количества итераций . III. После итерационного процесса формируются решения на шаге (31) с учетом уточненных разделенных разностей и на последней итерации. Оценивается масштабирующий множитель (41) с учетом демпфирующих условий (43). IV. Если шаг последний, процедура завершается. Иначе определяется величина следующего шага (41). Выполнение условия (45) устанавливает последний шаг интегрирования и тогда его величина переопределяется так, чтобы обеспечить выход на конечный момент времени . V. Экстраполируются значения разделенных разностей (30) с использованием значений функций дифференциальных уравнений (42) и вычислительный процесс повторяется от этапа II, где полученные решения принимаются за начальные . Заключение Таким образом, в статье изложены теоретические основы предлагаемого нового коллокационного интегратора для решения смешанных систем дифференциальных уравнений динамики первого и второго порядков. Акцентированы особенности в практической реализации интегратора, а также описана его программная процедура Lobbie. В дальнейших работах будут представлены результаты тестирования нового интегратора на примере простых динамических систем, а также показана его эффективность в сравнении с другими широко используемыми интеграторами.
Скачать электронную версию публикации
Загружен, раз: 19
Ключевые слова
динамические системы, обыкновенные дифференциальные уравнения, коллокационные методы, численные интеграторыАвторы
ФИО | Организация | Дополнительно | |
Авдюшев Виктор Анатольевич | Национальный исследовательский Томский государственный университет | д.ф.-м.н., профессор НИ ТГУ | scharmn@mail.ru |
Ссылки
Шефер В.А. // Астрон. журн. - 1991. - Т. 68. - С. 197-205.
Baumgarte J. // Comp. Math. Appl. Mech. Eng. - 1972. - V. 1. - P. 1-16.
Kustaanheimo P. and Stiefel E. // J. Reine Angew. Math. - 1965. - V. 218. - P. 204-219.
Cincotta P.M., Giordano C.M., and Simó C. // Physica D. - 2003. - V. 182. - P. 151-178.
Burdet C.A. // Z. Angew. Math. Phys. - 1968. - V. 19. - P. 345-368.
Kuntzmann J. // Z. Angew. Math. Mech. - 1961. - V. 41. - P. 28-31.
Butcher J.C. // Math. Comput. - 1964. - V. 18. - P. 50-64.
Everhart E. // Celest. Mech. - 1974. - V. 10. - P. 35-55.
Everhart E. // Dynamics of Comets: Their Origin and Evolution. Proceedings of IAU Colloq. 83, held in Rome, Italy, June 11-15, 1984 / ed. by A. Carusi and G. B. Valsecchi. - Dordrecht: Reidel, Astrophysics and Space Science Library, 1985. - V. 115. - P. 185-202.
Авдюшев В.А. Численное моделирование орбит небесных тел. - Томск: Издат. Дом Томского государственного университета, 2015. - 336 c.
Hairer E., Lubich C., Wanner G. Geometric Numerical Integration. Structure-Preserving Algorithms for Ordinary Differential Equations. - Springer, 2006. - 644 p.
Guillou A. and Soule J.L. // Rev. Francaise Informat. Recherche Oprationnelle 3. - 1969. - Ser. R-3. - P. 17-44.
Wright K. // BIT. - 1970. - V. 10. - P. 217-227.
Hairer E., Norsett S.P., Wanner G. Solving Ordinary Differential Equations. Nonstiff Problems. - Springer, 2008. - 528 p.
