Модифицированная версия интегратора Гаусса-Эверхарта Текст научной статьи по специальности «Математика»

МАТЕМАТИКА

УДК 519.622

МОДИФИЦИРОВАННАЯ ВЕРСИЯ ИНТЕГРАТОРА ГАУССА-ЭВЕРХАРТА

В. Г. Борисов

MODIFIED VERSION OF GAUSS-EVERHART INTEGRATOR

V. G. Borisov

Представлена модифицированная версия интегратора Гаусса-Эверхарта, реализованная в среде Delphi. Для повышения эффективности численного интегрирования были использованы вещественные переменные расширенной арифметики. Добавлены интерфейсные процедуры, позволяющие использовать интерпретатор формул для ввода исходных данных и графического вывода результатов интегрирования. Проведено тестирование модифицированной версии интегратора на примере задачи двух тел.

The modified version of the Gauss-Everchart integrator realized in Delphi is developed. To increase the efficiency of numerical integration, variables of «extended» type are used. The interface procedures allowing to use a parser for input data and graphic output of results are added. Testing of the modified version of integrator on the two-body problem is held.

Ключевые слова: системы ОДУ, численное интегрирование, интегратор Гаусса-Эверхарта, модифицированная версия.

Keywords: ODE systems, numerical integration, Gauss-Everchart integrator, modified version.

В 1973 г. Э. Эверхарт [5] предложил оригинальный алгоритм численного решения задачи Коши для систем обыкновенных дифференциальных уравнений небесной механики. Метод Эверхарта получил широкое распространение при исследовании уравнений небесной механики, что связано с его высокой эффективностью на этом классе задач.

Метод Эверхарта является одношаговым неявным методом высокого порядка типа методов Рунге-Кутты, предназначенным для решения систем уравнений первого и второго порядков, разрешенных относительно старшей производной. Реализация метода Эверхарта для решения задачи Коши вида:

х = / , х), ) = х0 (1)

подробно описана в [2]. Кратко идея метода состоит в следующем. На каждом шаге интегрирования, длину которого обозначим через И, дифференциальное уравнение (1) в локальной безразмерной независимой переменной 0 < т < 1 аппроксимируется формулой:

^ = Ь ■ (0 + ±Л,т'),

ат г=1 (2)

где к - некоторое натуральное число, 1° - известное значение правой части уравнения (1) при т = 0, Л1 - коэффициенты, подлежащие определению. Интегрирование равенства (2) по т, дает следующее выражение для искомой функции на рассматриваемом интервале:

х = х0 + Ь ■ (т / + 1-Л-тм),

/=1 I + 1 (3)

где х - известное значение искомой функции при т = 0. Далее на безразмерном интервале 0 < т < 1 задается разбиение Гаусса-Радо с узловыми точками т0 = 0, ть..., тк, и второе слагаемое в скобках правой части (2) представляется в виде интерполяционного многочлена Ньютона на данном разбиении с неизвестными коэффициентами а1, 1 = 1,..,к: к

^ ЛТ = а1т+а2т(т-т1) +

1=1

+а3т(т-т1)(т-т2) +... +акт(т-т)...(т-Т-1). (4)

Значения искомой функции х-1, ) = 1,..,к в узловых точках разбиения, согласно (3) будут выражаться формулами

к Л

xJ

x0 + h ■ (т

■ f0 +1 t! i +1

т+1).

(5)

Формулы (4), (5), а также само уравнение (1) дают неявные уравнения, связывающие величины xJ, а, и a,, i, j = 0,..,k. Эти уравнения на каждом шаге интегрирования решаются итерационным методом, при этом исходными данными для итераций являются значения величин xj и а,, вычисленные на предыдущем шаге. Алгоритм метода и теоретическое обоснование эффективного алгоритма выбора переменного шага для данной реализации метода детально изложены в [2]. Там же приведены результаты тестирования на примере задачи двух тел программы Gauss_15 на языке Fortran (названной автором интегратором Гаусса-Эверхарта), реализующей предложенный алгоритм.

В настоящей работе описана модифицированная версия интегратора Гаусса-Эверхарта. В качестве исходного кода для модификации был использован код бета-версии интегратора Gauss_32_mod (аналогичного Gauss_15) на языке Fortran [1]. Основными входными параметрами как исходной, так и модифицированной версии интегратора, являются размерность системы N, начальная TS, конечная TF, точки интервала интегрирования, начальное значение искомой функции и следующие параметры алгоритма:

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

- STEP - стартовый шаг интегрирования, при этом, если все элементы массива ERR равны нулю, то интегрирование на всем интервале проводится с постоянным шагом STEP. Если существуют ненулевые (положительные) элементы массива ERR и STEP=0, то интегрирование проводится с переменным шагом и величина стартового шага определяется автоматически с помощью итерационной процедуры, предваряющей процесс интегрирования.

- NOR - порядок интегратора. Этот параметр определяет величину k в (2) по формуле k = [NOR/2] ([] -целая часть числа). Здесь мы используем только нечетные значения параметра NOR.

- NI - число итераций на шаге для вычисления xJ, Ai и Oj. Практика показывает, что итерации сходятся достаточно быстро, и значения параметра NI большее 3 обычно не приводит к увеличению эффективности вычислений.

Выходными параметрами интегратора являются значение искомой функции в конечной точке интервала интегрирования и NF - суммарное число обращений к процедуре вычисления правой части системы уравнений (1) в течение всего процесса интегрирования. Показатель NF вместе с точностью полученного решения позволяет судить об эффективности процедуры интегрирования. Более эффективная процедура позволяет получить решение с заданной точностью при меньшем значении показателя NF.

Модификация интегратора, предпринятая автором при активном участии И.Е. Палехова [3], которому автор выражает благодарность, заключалась в нескольких моментах, речь о которых пойдет ниже.

Исходный код интегратора Gauss_32_mod.for был переписан на язык Pascal. Выбор языка был обоснован несколькими причинами, главной из которых является возможность использования типа данных «extended», что по сравнению с использованием переменных типа «Real*8», в перваначальной версии программы, дает повышение точности вычислений. Кроме этого, в модифицированной версии интегратора предполагалось использование графического вывода результатов вычислений и интерпретатора формул, функционирующих в программе Diff5 [3], разработанной в среде Delphi.

Для подтверждения корректности переноса кода программы на другой язык, были проведены тесты, повторяющие тесты, приведенные в [2]. Тестирование

проводилось для исходного кода на языке Fortran и двух версий модифицированного кода, скомпилированных в Lazarus и Delphi. Результаты интегрирования тестовых задач разными версиями программы сравнивались между собой и с результатами [2]. Интегрировалась задача двух тел

r=v, v = -// • r • r| r 1(0) = 1 - e, r2(0) = 0,

v1 (0) = 0, v2 (0) = V(1 + e)/(1 -e), (6) где r = (r1, r2) - координаты тела, v = (v1, v2) - вектор его скорости, 0 < e <1 - эксцентриситет орбиты, д > 0 -гравитационная постоянная,

И = у1 (г 1)2 + (r2)2 .

Известно, что решение задачи Коши (6) периодично с периодом 2п/д. В связи с этим, погрешность численного решения в конечной точке интервала интегрирования, длина которого кратна периоду, можно определить, сравнивая значения решения в начальной и конечной точках интервала. Ниже для оценки погрешности решения используется величина

Ar = V (r1(0) - r\TF ))2 + (r 2(0) - r \TF ))2,

где через TF обозначена конечная точка интервала интегрирования. Некоторые результаты тестирования модифицированной версии опубликованы в [4]. В качестве примера, на рис. 1 представлены диаграммы точность - быстродействие для Fortran (треугольные маркеры) и Delphi (квадратные маркеры) версий интегратора. Интегрировалась задача (6) c д = 1 и e = 0,999. Интегрирование проводилось на интервале [0, 2000п] (1000 оборотов) с параметрами NOR = 13, NI = 3, STEP = 0. Все элементы массива ERR полагались равными между собой и величина ERR варьировалась в интервале от 1e-5 до 1e-10 c мультипликативным декрементом 3,16. Каждому значению ERR на рис. 1 соответствует точка графика с ординатой, равной точности полученного решения Ar и абсциссой, равной NF.

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

Ax = [(и (0) - И (Tf ))2 + (V (0) - V (Tf ))2 ]

дают кривые, подобные приведенным на рис. 1. При использовании компилятора Lazarus существенных изменений в точности получаемых решений и величины NF обнаружено не было, лишь увеличивалось в 3 -4 раза физическое время интегрирования задачи. Физическое время счета для исходной версии приблизительно на 15 % лучше, чем время счета в тех же условиях той же задачи с помощью модифицированной версии, скомпилированной в Delphi.

1.E-Ü1 1.E-QZ 1.Е-03 1.E-Ö4 1.E-Q5 1.Е-05 1.Е-07 1.E-DB 1.Е-0Э

\

\

Ч л— к

К X/

—а— NOR=13 Fortran —■— NOR=13 Delphi —^NOFt=13 Delphi ERR \

f

1.00Е+06

NF

1.00Е+07

Рис. 1. Диаграммы точность - быстродействие для Fortran и Delphi версий интегратора

Результаты тестирования соответствуют результатам [2] с той лишь разницей, что достижимая точность вычислений модифицированной версии оказалась на 2 - 3 порядка выше, чем исходной, что связано с большей точностью представления чисел с плавающей точкой в расширенной арифметике.

Для модифицированной версии исследовалась зависимость от величины ERR суммарного количества несошедшихся итераций на всем интервале интегрирования, обозначенного ниже через I*. В таблице 1 приведена зависимость величин Ar, NF и I* от ERR для задачи (6) с д = 1, e = 0,99 и входными параметрами NOR = 15, NI = 3, STEP = 0, TF = 2000п. Все элементы массива ERR полагались равными между собой.

Таблица1

Зависимость Ar, NF и I* от ERR для задачи (6)

ERR Ar NF I*

1,00e-09 4,09e-09 4,78e+06 1,05e+05

3,16e-10 4,51e-10 5,52e+06 2,31e+04

1,00e-10 1,91e-10 6,38e+06 5,66e+03

3,16e-11 4,51e-10 7,38e+06 7,18e+02

1,00e-11 1,28e-11 8,52e+06 1,73e+02

3,16e-12 2,50e-10 9,84e+06 4,00e+01

1,00e-12 3,08e-10 1,14e+07 1,00e+01

Данные таблицы 1 показывают монотонное возрастание NF и убывание I* при уменьшении величины задаваемой точности ERR от 1e-9 до 1e-12. При этом, величина Ar, убывая в начале интервала вариации ERR, затем переходит к упомянутым выше случайным колебаниям, вызванным влиянием ошибок округления. Таким образом, большое значение величины I* соответствует интегрированию с точностью, не достигающей предельно возможного значения и с относительно малым значением NF, то есть с меньшим временем счета. Слишком малое значение величины I* соответствует зоне случайных колебаний Ar и неоправданно большому времени счета.

В связи с этим было принято решение ввести в модифицированную версию интегратора механизм, авто-

матически изменяющии величины элементов массива ERR в соответствии с характером сходимости итераций для соответствующей компоненты решения. Действие механизма состоит в том, что в ходе последней итерации на каждом шаге формируются множители, изменяющие величины элементов массива ERR, используемых на следующем шаге. Проводились эксперименты с различными вариантами механизма. Расчеты, проведенные на задаче (6) и других задачах, показывают, что введение такого рода отрицательной обратной связи приводит к тому, что величина NF практически перестает зависеть от входных значений параметра ERR, а величина достигнутой точности колеблется в некотором диапазоне, в связи с влиянием ошибок округления. На рис. 1 положениями круглых маркеров изображены соотношения точность -быстродействие при решении задачи (6) c д = 1 и e = 0,999 на интервале [0,2000п] с параметрами интегрирования NOR = 13, NI = 3, STEP = 0, при использовании одного из вариантов механизма автоматического изменения элементов массива ERR.

Механизм автоматического изменения ERR, очевидно, имеет смысл применять только при интегрировании с переменным шагом, при этом его действие накладывается на действие механизма изменения величины шага, исходно существующего в программе. Расчеты, проведенные на задаче (6) с разными значениями эксцентриситета показали, что изменение величины ERR не препятствует выбору оптимальной величины переменного шага в ходе интегрирования. На рис. 2 приведен график зависимости величины шага H от номера шага на одном обороте для задачи (6) c д = 1, e = 0,99 и входными параметрами NOR = 15, NI = 3, STEP = 0. На этом же рисунке приведен график зависимости от номера шага, величины, пропорциональной |r|15. Тот факт, что эти два графика близки, означает, что при интегрировании данной задачи используются величины шагов, близкие к теоретически обоснованным в [2] оптимальным величинам. На этом же рисунке пунктирной линией приведен график зависимости величины ERR[2]*1e10 от номера шага для одного из вариантов механизма автоматического изменения ERR.

Рис. 2. Зависимость величины шага интегрирования Н от номера шага на одном обороте

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

ri vi =-Mi ■ ri ■ |rj| , t = 1,..,4, (7)

i о ^ i о ^

где ri=(ri, ri, ri), vi=(vi, vi, vi) - координаты и векторы скорости тел, |r| = ^(r1 )2 + (r2 )2 + (r3 )2 . Начальные условия для системы (7) и параметры ^ выбирались так, чтобы движения происходили в ортогональных координатных плоскостях с различными периодами Ti и эксцентриситетами ei. Выбранные значения эксцентриситетов ei, параметров ^ и соответствующих им периодов Ti приведены в таблице 2.

Таблица 2

Значения ei9 T

В силу выбора в качестве ^ чисел, обратных к простым, минимальный период решения системы (7) с указанными параметрами и начальными условиями равняется 546п.

На рис. 3 изображены диаграммы точность - быстродействие для задачи (7). Линии на графиках соответствуют значениям погрешностей

^ = (0) - (ГР ))2 ,

где TF = 546п. Интегрирование проводилось на интервале [0,546п] с входными параметрами NOR = 15, NI = 3, STEP = 0. Величина ERR варьировалась от 1e-6 до 1e-12 с мультипликативным декрементом 10.

Линии на рис. 3 соответствуют погрешностям Ari, i=1,..,4 при интегрировании без вариации ERR, значки большого размера, серого цвета - величинам Ari, полученным при интегрировании с одним из вариантов механизма автоматического изменения ERR.

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

Была добавлена опциональная возможность задания системы уравнений, параметров и начальных данных во внешнем текстовом файле определенного формата. При этом оставлена возможность задания системы дифференциальных уравнений в коде программы. Выбор этой опции управляется дополнительным входным параметром. Модифицированная версия интегратора представляет собой исполняемый файл, который с помощью встроенного в него интерпретатора формул обращается к внешнему текстовому файлу, содержащему систему уравнений, начальные условия, параметры интегрирования и другие параметры. Работа с использованием интерпретатора формул имеет существенное преимущество, заключающееся в том, что для изменения системы нет необходимости перекомпиляции кода. Расчеты, проведенные с системой, заданной в коде программы и заданной во внешнем файле, дают идентичные результаты по точности решения и величине NF, при этом время счета с интерпретатором формул оказывается значительно большим.

i Hi ei Ti

1 1 0,9 2п

2 1/3 0,1 6п

3 1/7 0,99 14п

4 1/13 0,999 26п

Рис. 3. Диаграммы точность - быстродействие для компонент решения задачи (7)

Была добавлена возможность выдачи результатов интегрирования в табличном и графическом видах. Для этого был использован интерфейс приложения Diff5. Для реализации возможности выдачи результатов интегрирования не только в конечной точке ТБ, как это сделано в исходной версии, но и в промежуточных точках, в интегратор был добавлен входной параметр 51ерои1 Выдача промежуточных результатов с использованием параметра Б1ерои1 не отражается на точности и быстродействии вычислений, поскольку дополнительных вычислений при выводе промежуточных результатов не производится. Процедура вывода промежуточных результатов состоит в том, что в ходе вычислений во вспомогательный массив выборочно записываются вычисленные значения решения на шаге, а именно, на каждом первом шаге, превышающем Т8 + п*Б1ерои1 для всех 0 < п < (ТБ -ТЗЭМерои! Массив становится доступным для обработки с целью выдачи результата в табличном виде

либо в виде 2D- и 3D-графиков выбранных компонент решения по завершению процесса интегрирования.

Таким образом, тестирование модифицированной версии подтвердило адекватность переноса кода на другой язык и увеличение эффективности вычислений по сравнению с исходной версией, что связано с использованием расширенной арифметики. Тестирование различных вариантов механизма автоматического изменения величин ERR.li], введенного в код программы, показало несущественную зависимость величины КБ и точности вычислений от входных значений этих величин для исследуемого класса задач. В то же время универсального механизма, дающего увеличение эффективности вычислений на всем классе задач, обнаружено не было. Подтверждена идентичность результатов интегрирования и количества обращений к процедуре вычисления правых частей КБ с использованием интерпретатора формул и без него.

Литература

1. Авдюшев В. А. Интегратор Гаусса-Эверхарта. Новый фортран-код // Фундаментальные и прикладные проблемы современной механики: материалы Всероссийской конференции. Томск, 3 - 5 октября 2006 г. Томск: Изд-во ТГУ, 2006 с. С. 411 - 412. Режим доступа: http://www.scharmn.narod.ru/AVD/Software.htm

2. Авдюшев В. А. Интегратор Гаусса-Эверхарта // Вычислительные технологии. 2010. Т. 15. № 4. С. 31 - 47.

3. Каков Р. Н., Ганеев Д. Р. Программа Diff4 для численного и качественного анализа решений обыкновенных дифференциальных уравнений // Материалы V Международной научно-практической конференции студентов, аспирантов и молодых ученых. Кемерово, 2010. T. 2. C. 93 - 95.

4. Палехов И. Е., Борисов В. Г. Тестирование модифицированной версии интегратора Гаусса-Эверхарта // Информации в технологиях и образовании: материалы VII Международной научной конференции. Белово, 2014. T. 2. С. 78 - 81.

5. Everhart E. A New Method for Integrating Orbits // Bulletin of the American Astronomical Society. 1973. V. 5. P. 389.

Информация об авторе:

Борисов Владимир Геральдович - кандидат физико-математических наук, доцент кафедры фундаментальной математики КемГУ, [email protected].

Vladimir G. Borisov - Candidate of Physics and Mathematics, Assistant Professor at the Department of Fundamental Mathematics, Kemerovo State University.

Статья поступила в редколлегию 23.04.2015 г.