О новом способе оценки локальной ошибки в методах Гаусса – Эверхарта

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

В небесной динамике наряду с аналитическими [1, 2] и полуаналитическими методами [3,4] широко применяются алгоритмы численного интегрирования. Например, расчет схода с орбиты [5] ввиду сильной нелинейности уравнений невозможен без использования численных методов. Важную роль такие алгоритмы играют при предсказании движения космических аппаратов, участвующих в современных геодезических миссиях. Так, в проектах GRACE и GRACE-FO [6,7] при построении траекторий спутников используется метод интегрирования Коуэлла 7-го порядка с постоянным шагом. Численные методы являются неотъемлемым инструментом для анализа и интерпретации данных лазерной дальномет-рии [8], которые обеспечивают расчет некоторых параметров вращения Земли.

Для численного решения задачи Коши часто используются явные одношаговые методы интегрирования, среди которых можно выделить методы Гунге - Кутты (например, классический метод Гунге - Кутты [9], алгоритмы Дорманда - Принца [10,11] и Гунге -Кутты - Фелберга [12]). Они обладают простой структурой, их легко реализовать, но при этом явные методы имеют ограничения в точности и устойчивости. Это особенно заметно при моделировании орбит с высоким эксцентриситетом или в случае интегрирования траекторий на длительные периоды времени. Также для численного решения задачи Коши применяются многошаговые методы. Эти алгоритмы занимают особое место в задачах интегрирования движения тел солнечной системы [13]. Широко известным представителем многошаговых методов является метод Гаусса - Джексона, который используется центром космических наблюдений США с 60-х годов [14]. Среди его преимуществ можно выделить сдерживание ошибок округления, что подчеркивалось в работах NASA [15]. Кроме того, при интегрировании орбит, близких к круговым, хорошие результаты показывает многошаговый явный симметричный метод интегрирования Куинлэна - Тремейна [16,17]. Тем не менее многошаговые алгоритмы имеют ряд недостатков, среди которых можно выделить, например, ограничение на порядок метода [18]. Отдельного внимания заслуживают неявные одношаговые методы. Их преимущество перед явными одношаговыми методами заключается в том, что для достижения одинакового порядка аппроксимации неявным методом необходимо меньшее число вычислений правой части системы ОДУ. Кроме того, неявные методы, по сравнению с явными, часто обладают лучшей устойчивостью и больше подходят для жестких систем ОДУ. В работе [17] для задач высокоточного моделирования орбитального движения рекомендуются именно неявные методы Гунге - Кутты. Среди них особого внимания заслуживает оригинальный коллокационный интегратор Гаусса -Эверхарта (здесь и далее под ^оригинальными» будем понимать метод интегрирования и способ оценки локальной ошибки, впервые предложенные Эдгаром Эверхартом в [19]). Этот подход хорошо зарекомендовал себя в различных приложениях небесной механики, в особенности в кометной динамике.

Развитие метода Гаусса - Эверхарта не останавливается ввиду его эффективности в задачах механики небесных тел. Преимущества его использования по сравнению с другими методами интегрирования можно найти в [17]. Помимо оригинального метода Гаусса -Эверхарта, существуют его модифицированные версии, которые позволяют повысить точность численного интегрирования. Такие усовершенствования становятся критически важными в задачах, требующих точности. В работе [20] показан отличный от оригинального подход к разложению функций, который позволил авторам увеличить порядок аппроксимации и улучшить результаты моделирования движения астероидов. В [21] также продемонстрирована модификация алгоритма, позволившая повысить быстродействие в несколько раз для порядков аппроксимации выше 15-го. В рамках продолжения развития коллокаци-ohhbix методов был также представлен интегратор Lobbie [22]. Его достоинством является возможноств решения смешанных систем дифференциа.льнв1х уравнений второго и первого порядков, применяющихся в задачах орбитальной динамики. Согласно этой работе, использование стабилизирующих и регуляризирующих преобразований позволило повысить эффективность интегратора Lobbie. Таким образом, неявные методы интегрирования, в особенности метод Гаусса - Эверхарта, продолжают активно развиваться ввиду их значимой роли в задачах небесной механики.

Применение численных методов решения задачи Коши нередко предполагает необходимость коррекции шага интегрирования. Эта потребность особенно актуальна в жестких задачах, например при моделировании орбит околоземных тел с высокими эксцентриситетами. Оценка локальной ошибки интегрирования и расчет шага на ее основе играют важную роль в поддержании точности вычислений и предотвращении численной неустойчивости метода. Это, в свою очередь, обеспечивает эффективность алгоритмов. Заметим, что оригинальный метод Гаусса - Эверхарта включает в себя свой собственный механизм оценки локальной ошибки, что является одним из ключевых аспектов работы [19]. Однако такой подход имеет недостаток: в процессе адаптивного выбора шага интегрирования оценка ошибки фактически производится на основе вспомогательного решения существенно меньшего порядка, чем порядок основного метода [17]. В результате полученная оценка может быть значительно завышена, что приводит к излишнему уменьшению шага и, как следствие, к неоправданному увеличению вычислительных затрат. На практике это означает, что даже при высокой точности самого метода интегрирования его эффективность может быть ограничена частыми пересчетами. В задачах, требующих длительного интегрирования траекторий или моделирования большого числа объектов, — например, при исследовании эволюции популяции космического мусора или при прогнозировании движения спутниковых группировок — такие дополнительные вычисления становятся критическим фактором, существенно увеличивающим общее время расчета. Поэтому возникает необходимость в улучшении механизма оценки локальной ошибки, сохраняющего точность исходного алгоритма, но требующего меньших вычислительных ресурсов.

Одним из первых алгоритм оценки предложил Гунге [23]. Его подход заключается в дополнительных вычислениях с уменьшением размера шага, и в сравнении результатов с основным решением. Такой метод требует больших ресурсов, что замедляет процесс решения и снижает эффективность применения алгоритма.

В настоящее время проведены исследования, посвященные разработке и улучшению методов подбора шага интегрирования при решении задач неявными методами. Например, в работе [24] описан подобный алгоритм. Но рассматриваемая в этой статье система ОДУ представлена в упрощенном варианте, что делает бессмысленным перенос идеи [24] на метод Гаусса - Эверхарта, призванного решать задачи высокоточного моделирования. В работе [25] описан подход к оценке ошибки, основанный на таблицах Бутчера. Однако в оригинальном методе Гаусса - Эверхарта не предполагается их составление, и для полноценного применения [25] необходимо предварительное вычисление таблиц Бутчера. Это открывает перспективу для будущих исследований, направленных на разработку алгоритмов, специально адаптированных под метод Гаусса - Эверхарта.

В рамках задачи оценки шага в неявных алгоритмах также можно найти применение автоматического управления и ПИ-регуляторов [26,27]. Например, в работе [28] показано, что использование адаптивного шага по времени с регулятором может снизить вычислительную нагрузку и повысить точность интегрирования. Тем не менее ПИ-регуляторы не лишены ограничений применимости. Они создают препятствия в сценариях расчета методом Гаусса - Эверхарта ввиду своей зависимости от таблицы Бутчера [26,27]. Важно отметить, что теоретически регуляторы могут быть совмещены с алгоритмами вычисления локальной ошибки интегрирования с целью дополнительной оптимизации расчета. Эта комбинация методов может обеспечить более гибкое управление шагом, но этот вопрос тре- бует дополнительных исследований и остается открытым для будущих работ.

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

2.    Кол локационные методы Гаусса — Эверхарта

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

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

У ‘ = f^(t, У), У⁽М = Уо .                                   (!)

Решение у приблизим полиномом g _s степени s на от резке [to,to + h]:

s

У^(t) ~ gsСО = Уо + ^h ^ у ^{т j},                          (2)

1=1 J где т = '—0 — нормализованная переменная времени; коэффициенты а? пока не определены, способ их вычисления будет описан ниже. Представление (2) позволяет описать поведение решения на заданном отрезке [ф,ф + h].

Введем строго возрастающий набор точек коллокации {ti} на от резке [фЦо + h]:

ti = t0 + hci, i = 1,...,s,(3)

где 0 C ci < C2 < ... < cs C 1- Потребуем, чтобы производные полинома gs(t) в точках {ti} удовлетворяли следующему условию:

s k i = g's(ti) = 52 aj Ci-1, i = 1,...,s,(4)

1=1

где {к i } представляют собой значения правой части системы (1) в точках коллокации {ti}:

ki = f(ti, gs(ti)), i = 1,...,s.(5)

Запишем полином g{ в виде интерполянта Ньютона:

sj g‘s(to + hT) = «1 + 52“? П(т - Cm).(6)

j=2

Равенства (4) позволяют найти набор из s разделенных разностей {a_i}:

« 1 = ki,

« 2 = (к2 - « 1)/(C2 - C1),

« 3 = ((к3 - « 1)/(сз - C1) - « 2)/(сз - C2),

При этом коэффициенты {«i} и {ai} связаны соотношением:

s аj = 52 ^j«i, J = 1,..., s, i=j

где X i,j — числа Стирлинга, значения которых определяются выражениями:

^ i,j — ^ i - 1 ,j - 1     ^c i - 1 ^ i - 1 ,j ,

X i,i — 1 ,   Ai, 0 — 0 ( i> 0) .

Далее, для нахождения приближенного решения z i = g_s(ti), i — 1,...,s, в точках коллокации {t i }, нужно решитв нелинейную систему уравнений:

z — z (a (a (k (z)))) ,

(Ю)

где отсутствие индексов обозначает набор векторов, соответствующих точкам коллокации. Решение системы (10) осуществляется аналогом метода Гаусса - Зейделя для нелинейных систем.

После решения системы приближенное решение в точке to + h вычисляется с использованием формулы (2) при т — 1:

S

У1 ^— Уо + ^h ^ —.                            (11)

3 =1 ²

Более подробное описание коллокационного метода можно найти в [17]. В этой монографии показано, что коллокационный метод обеспечивает порядок аппроксимации не менее s. Также заметим, что можно повысить порядок аппроксимации при BBi6ope специфического разбиения {ci}. Например, при применении разбиения Гаусса - Лежандра, порядок составляет 2s, а при исполвзовании разбиения Гаусса - Радо (с ci — 0, что удобно для интегрирования задачи Коши) порядок равен 2s — 1.

3.    Оценка локальной ошибки и выбор шага интегрирования

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

3.1.    Оригинальный способ

В оригинальной работе [19] оценка локальной ошибки и выбор шага основаны на предположении, что абсолютная ошибка интегрирования совпадает с последним слагаемым полинома в выражении (2) при т — 1

Z abs — ^h

Ы* s

где II • II2 - евклидова норма. Соответствующая относительная ошибка определяется выражением:

^ abs

e —      ’ где У1 ^ приближенное решение (11).

Для того чтобы новый шаг интегрирования h_new соответствовал допустимой относительной локальной точности tol, предлагается выбирать шаг по следующему правилу:

^h new

— h () ^1/s.

Проведем аналогию оригинального подхода с оценкой ошибки во вложенных явных методах. Пусть у 1 — решение с порядком аппроксимации l, полученное с использованием формулы (11). При выборе в качестве коллокационных точек разбиения Гаусса - Лежандра имеем I = 2s, а для левого разбиения Гаусса - Радо — I = 2s — 1.3а у^ обозначим вспомо-гателвное решение, определяемое формулой:

S —1

У? = Уо + h £  .

з =1 ^J

Здесь индекс р обозначает порядок аппроксимации вспомогательного решения. В таком случае относительная ошибка £rei может быть представлена как отношение разности двух решений к одному из них:

_   _ Ну 1 ^— У1 Н 2

^£rel =     Ну ? Н 2

Такой метод оценки локальной ошибки дает более точные результаты, когда порядки аппроксимации I и р близки друг к другу.

Решение у 1 можно формально рассматривать как полученное после интегрирования полинома в виде Ньютона gД1, который совпадает с полиномом (6) за исключением последнего слагаемого:

g 3 -1

з       3 - 1

«1 +    a j Ц (т — С т ) —a_sr ^3-1

3 =2    т =1

Здесь учтено, что, согласно (8) и (9),

O' S — ^s,s ^a s — ^a s .

Полином (15) не удовлетворяет условиям (4) и в общем случае не является интерполяционным. Значит, нельзя гарантировать, что порядок аппроксимации р будет удовлетворять неравенству р ^ s — 1. Таким образ ом, порядки I и р могут существенно различаться, что, в свою очередь, способно привести к завышенной оценке локальной ошибки (12).

3.2.    Авторский способ

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

Как и в оригинальной работе Эверхарта [19], в настоящем исследовании используются узлы левого разбиения Гаусса-Гадо. При этом, после вычисления решения у[, мы предлагаем исключить из множества (3) некоторую точку, кроме первой. Следует подчеркнуть, что удаление точки с г — 1 (т.е. С1 — 0) недопустимо по соображениям вычислительной эффективности. Дело в том, что С1 соответствует g'_S(to) для любого порядка схемы, а, значит, решение, полученное на предыдущем шаге, может быть использовано напрямую без пересчета в итерациях метода Гаусса - Зейделя при решении (10).

После исключения одной из точек коллокации из исходного набора (3) сформируем новый набор точек ДД на от резке [to, to + h]:

ti — to + hci, i — 1,... ,s — 1, где {ci} образовано множеством {ci} (3) исключением одного элемента: {ci} — {ci}\ Cj для некоторого j G {2,..., s}.

Соответствующие значения правой части системы (1) обозначим через {k i}. Они получаются из исходного набора {к i } (5) путем удаления элемента, связанного с исключенной точкой Cj. Построим набор из s — 1 разделенных разностей {«i} аналогично (7):

«1 = ki,

«2 = (к2 - «1)/(с2 - С1),

«3 = ((к3 - «1)/(сз - С1) - «2)/(сз - С2),

Теперь, используя числа Стирлинга X i,j для нового коллокационного набора точек {с;}, вычислим s - 1 значение коэффициентов {а^}:

в—1

a j = ^A_i,j «■id = 1,...,s - 1.                            (17)

^; = ^j

Искомое вспомогательное решение определяется формулой:

в—1 _

/1 = л + Е • j=1 J

Далее, оценку локальной ошибки и выбор шага интегрирования можно проводить при помощи формул (13) ,(14).

Авторский способ оценки согласуется с принципами коллокационных методов. Он гарантирует, что у 1 имеет порядок аппроксимации не менее s - 1. Заметим, что при исключении последней точки из набора (3) выполнять вычисления (16) не нужно. В таком случае набор коэффициентов {«;} может быть получен путем отбрасывания последнего элемента в наборе {а;}. Обратим внимание на то, что авторский способ не включает в себя дополнительные вычисления правой части системы дифференциальных уравнений {кД, а следовательно существенно не увеличивает ресурсоемкость и время расчетов.

4.    Анализ эффективности авторского способа

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

4.1.    Орбита Аренсторфа

Начнем тестирование авторского способа с построения орбиты Аренсторфа [29]. Орбита Аренсторфа — плоская периодическая траектория объекта в поле тяжести двух тел, находящихся в точках (0, 0) и (1, 0). Отношение их масс у и у‘ соответствует отношению масс Луны и Земли, причем у = 0.012277471, а сумма у + у‘ равна 1.

Координаты тела г(t) = (x(t'),y(t')')^T на орбите Аренсторфа являются решением задачи Коши системы обыкновенных дифференциальных уравнений второго порядка:

Г .        /х + у x -у‘

Х = х + 2у -у—у——

^1

О.       / УУ у = У - 2х - у— -у—,

D1

где

D 1 = ((х + у)² + у²)³/², D2 = ((х -у ‘)² + у²)^3/2.

Начальные условия для системы (18) имеют вид

0.994

Г (^о) =10 Г

V(*0) = (

-2.00158510637908252240537862224

где v(1) = (x(t),y(t))^T.

В рамках работы было проведено тестирование методик оценки локальной ошибки на примере построения орбиты Аренсторфа. Результаты расчетов представлены на рис. 1. Задача Коши (18) - (19) решалась методом Гаусса - Эверхарта 15-го порядка с левым разбиением Гаусса - Радо и с допустимой локальной точностью tol = 10^-9. Локальная ошибка оценивалась как с помощью оригинального алгоритма Эверхарта, так и с применением авторского способа с исключением второй, пятой и восьмой точек из (3). Красными точками обозначены притягивающие центры: Земля и Луна.

Рис. 1. Построение орбиты Аренсторфа с использованием авторского и оригинального способов оценки ошибки. Цветовая шкала иллюстрирует значение шага интегрирования

Заметим, что графики на рис. 1 обладают схожим распределением величины шага. Вблизи Луны отношение шага интегрирования к периоду уменьшается до значений порядка 10-4, что свидетельствует о высокой чувствительности алгоритмов к резким изменениям траектории в этой области. При отдалении от Луны и Земли шаг увеличивается до 10-2-10-3 периода, обеспечивая эффективное интегрирование на больших расстояниях. Авторский способ оценки реагирует на особенности динамики задачи подобно оригинальному подходу, что подтверждает работоспособность авторского способа.

В авторском подходе максимальный шаг интегрирования оказывается в 2-2.5 раза больше в зависимости от номера исключаемой точки. Это наблюдение объясняется тем, что авторский способ допускает более крупные шаги интегрирования при сохранении заданной допустимой локальной точности. Увеличение шага приводит к уменьшению числа обращений к правой части системы ОДУ (1), что на практике сокращает общее время расчёта.

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

4.2.    Движение в потенциальном поле точечной массы

Для сравнения описанных способов оценки локальной ошибки также рассматривалась задача ограниченного движения материальной точки в центральном гравитационном поле:

Г = V, V =

г

-ц й,

Г (to) = Го, v(to) = V0, где г и v — векторы положения и скорости тела соответственно, а ц — гравитационный параметр. Такой выбор задачи обуславливался наличием аналитического решения, а также возможностью управлять кривизной траектории за счет варьирования эксцентриситета орбиты при помощи начальных условий.

Для исследования поведения алгоритмов оценки ошибки при построении орбит решалась задача (20) для различных значений эксцентриситета (е = 0, 0.1, 0.4 и 0.8) и при единичных значениях большой полуоси и гравитационного параметра ц. Интегрирование осуществлялось методом Гаусса - Эверхарта 15-го порядка с левым разбиением Гаусса -Радо.

На рисунках 2-5 представлены результаты в виде точек, полученных при численном интегрировании. Цвет каждой точки соответствует отношению длины шага интегрирования к периоду движения; положение притягивающего центра отмечено красной точкой. Такая интерпретация позволяет визуально отследить реакцию метода на локальные изменения траектории.

Координата х

ф

з

ф

ф

Рис. 2. Решение задачи (20) для различных способов оценки локальной ошибки при е = 0

Координата х

Цветовая шкала вдоль траектории на рис. 2 однородна, что соответствует равномерному движению по окружности. Это подтверждает корректность авторского способа в рассматриваемом случае. Вместе с тем, при сравнении значений шага можно сделать ряд важных замечаний:

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

• •    Исключение точек из разбиения приводит к увеличению шага по сравнению с оригинальным методом и, следовательно, к уменьшению общего числа вычислений. Особенно заметен рост шага при исключении точек с меньшими номерами.

  

  Рис. 3. Решение задачи (20) для различных способов оценки локальной ошибки при е = 0.1

  
  
  
  
  

  Координата х

  
  
  

На рисунках 3-5 продемонстрировано, как изменяется поведение метода при увеличении эксцентриситета орбиты. В отличие от случая с круговой орбитой, траектория характеризуется неравномерностью движения. Это согласуется с аналитическим решением и необходимостью более частых пересчётов в области максимального возмущения. Такое поведение чётко прослеживается на всех графиках — цветовая шкала демонстрирует падение шага при приближении к перицентру и увеличение при удалении от него. Помимо всего прочего, можно провести параллель с результатами на круговой орбите:

• •    Оригинальный подход Эверхарта снова демонстрирует наименьший шаг интегрирования вдоль всей траектории. При этом в некоторых случаях максимальная величина шага оказывается до 3.5 раз меньше, чем у авторского способа.

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

  

  Координата х

  
  
  

  Координата х

  
  
  
  
  
  

  Координата х

  
  
  

  Рис. 4. Решение задачи (20) для различных способов оценки локальной ошибки при е = 0.4

  
  
  

  Рис. 5. Решение задачи (20) для различных способов оценки локальной ошибки при е = 0.8

  
  

Описанные выше эксперименты позволили рассмотреть влияние различных способов оценки локальной ошибки на поведение метода Гаусса - Эверхарта как в случае круговой орбиты, так и при ненулевом эксцентриситете. Было показано, что выбор подхода к оценке локальной ошибки влияет на эффективность интегрирования, в частности — на средний шаг при фиксированной допустимой погрешности. Далее, проводится более строгое количественное сравнение этих схем, позволяющее оценить различия между ними при решении задачи (20) методом Гаусса - Эверхарта 11-го порядка.

В качестве входных параметров расчета задавалась величина допустимой локальной погрешности tol и эксцентриситет орбиты. При этом полагалось равенство единице большой полуоси и гравитационного параметра ц. В результате численных экспериментов были получены следующие характеристики:

• 1)    Средний шаг интегрирования h_or i _g, соответствующий оригинальному способу оценки локальной ошибки.

• 2)    Средние шаги интегрирования h j , j = 2,..., 6, соответствующие авторскому способу оценки локальной ошибки за исключением одной из пяти коллокационных точек (индекс j указывает на номер исключаемой точки C j из миожества {c_i }).

• 3)    Шаг Н — гипотетически идеальный шаг интегрирования, при котором фактическая локальная ошибка равна заданной допустимой точности tol. Эта зависимость была получена путем сравнения численного решения с аналитическим решением для движения по круговой и эллиптической орбитам.

Для каждой кривой усреднение шагов проводилось по 100 периодам обращения.

Рис. 6. Шаги интегрирования для различных способов оценки локальной ошибки. Нулевой эксцентриситет орбиты

На рисунках 6-8 показаны зависимости шага от заданного допустимого значения ошибки tol для различных эксцентриситетов. По горизонтальной оси в логарифмическом масштабе отложена величина tol, по вертикальной — средний шаг, нормированный на период обращения. Штриховой линией обозначен шаг h_orig’, штрихпунктирными линиями — шаги h j , j = 2,..., 6, соответствующие схемам с исключением отдельной точки; сплошной линией — гипотетически оптимальный шаг Н(tol) В левом верхнем углу рис. 6-8 приведен увеличенный фрагмент, иллюстрирующий различия между авторскими схемами с исключенными коллокационными точками C j , j = 2,..., 6; номера рядом с кривыми соответствуют значениям индекса j.

Анализ графиков показывает, что авторский способ с исключением одной из точек позволяет более точно и реалистично оценивать локальную ошибку интегрирования. Это, в свою очередь, обеспечивает заметное увеличение среднего шага по сравнению с шагом, определяемым оригинальным подходом Эверхарта. При е = 0 и tol = 10-⁸ оригинальный способ выбирает шаг около 0.021 периода, тогда как авторские модифицированные схемы позволяют использовать шаги от 0.04 до 0.05 периода. Это соответствует увеличению шага на 90-140% в зависимости от исключаемой точки, при этом шаг остается под кривой Н(tol)

Аналогичный результат в 90-140% наблюдается и при усреднении величины роста шага по всем значениям tol. Для орбит се = 0.1 рост шага в среднем по tol составляет 120-170%, а для е = 0.3 достигает 130-180%. Это означает, что ошибка контролируется не менее надёжно и оценивается более адекватно, что позволяет избежать избыточно малого шага интегрирования. Именно это — более точная и реалистичная оценка ошибки — и является главным преимуществом авторского подхода. В результате при его использовании снижается число шагов и, как следствие, сокращается общее время расчёта. Заметим, что среди рассмотренных примеров наибольшую эффективность демонстрирует исключение крайней левой (второй) точки.

Рис. 7. Шаги интегрирования для различных способов оценки локальной ошибки. Эксцентриситет 0.1

Рис. 8. Шаги интегрирования для различных способов оценки локальной ошибки. Эксцентриситет 0.3

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

5.    Заключение

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

В ходе тестирования авторский способ оценки ошибки продемонстрировал свою работоспособность. Преимущества авторского способа подтверждаются сравнительным анализом на модельных задачах движения материальной точки. В задаче кругового движения авторский способ позволяет выбирать шаг, превышающий шаг оригинального подхода на 90-140%, а при увеличении эксцентриситета орбиты до е = 0.3 прирост достигает 180%. В то же время важным достижением является сохранение его вычислительной эффективности ввиду отсутствия дополнительного расчета правой части дифференциального уравнения.

В дальнейших работах предлагается обратить внимание на перспективы использования авторского способа оценки ошибки в комбинации с ПИ-регуляторами. Возможно, такое сочетание сможет обеспечить более качественное управление шагом. Кроме того, целесообразно вычислить таблицу Бутчера метода Гаусса - Эверхарта и проанализировать ее применение в рамках других способов оценки локальной ошибки.

6.    Вклад авторов

• •    Фукин И.И. Реализация численной схемы метода Гаусса-Эверхарта 15-го порядка с адаптивным шагом интегрирования. Внедрение авторского способа оценки локальной ошибки и обеспечение корректности его вычислительной реализации.

• •    Кузнецов А.А. Предложение концепции авторского способа оценки локальной ошибки в коллокационном методе Гаусса - Эверхарта. Теоретическое обоснование исключения отдельных коллокационных точек и направления численного эксперимента.

• •    Носырев А.Н. Формулировка и подготовка модельных задач для тестирования авторского способа, включая построение орбиты Аренсторфа и движение тела в потенциальном поле точечной массы. Проведение расчётов для различных значений эксцентриситета и параметров допустимой погрешности.