CC BY
174
Вычислительные технологии
Том 15, № 4, 2010
Интегратор Гаусса—Эверхарта*
В. А. Авдюшев НИИ прикладной математики и механики Томского государственного университета, Россия e-mail: [email protected]
Представлена новая версия интегратора Эверхарта для решения обыкновенных дифференциальных уравнений первого порядка. Обсуждаются особенности в компьютерной реализации интегратора, исследуются его возможности при решении задач небесной механики.
Ключевые слова: численное интегрирование, обыкновенные дифференциальные уравнения, интегратор Эверхарта, неявные методы Рунге—Кутты, численное моделирование орбит.
Введение
В 1973 г. Э. Эверхарт [1] предложил интегратор, специально разработанный для численного исследования орбит, и продемонстрировал его высокую эффективность в задачах кометной динамики. По-видимому, обнаружив в дальнейшем принадлежность данного интегратора к семейству интеграторов бутчеровекого типа, Эверхарт акцентировал внимание на оригинально реализованный им алгоритм интегрирования и обобщил его для численного решения любых обыкновенных дифференциальных уравнений первого и второго порядка [2, 3], тем самым расширив область применения предложенного интегратора, который, тем не менее, остается одним из самых популярных именно в решении задач небесной механики.
Интегратор Эверхарта (ЕА15) основан на видоизмененных формулах неявных кол-локационных методов Рунге—Кутты бутчеровекого типа, поэтому он наследует все их замечательные свойства [4]. Более того, именно благодаря оригинальному представлению вычислительной схемы интегратор Эверхарта с точки зрения численного интегрирования имеет ряд следующих преимуществ:
1) алгоритм интегрирования универсален для любого порядка;
2) интегратор имеет простой критерий для выбора шага интегрирования;
3) в интеграторе реализован достаточно точный предиктор решения, что позволяет выполнять численное интегрирование всего с двумя итерациями на шаге.
Несмотря на это, программный код Эверхарта RA15 [3] (впрочем, как и любая его модификация типа RAD.M 27). на наш взгляд, довольно существенно ограничивает возможности интегратора и поэтому нуждается в дополнительной редакции.
*Работа выполнена при финансовой поддержке РФФИ (грант № 08-02-00359) и Министерства образования РФ (грант № РНП.2.1.1/2629).
Среди главных недостатков в программной реализации интегратора можно назвать следующие:
1) трудночитаемый и громоздкий код;
2) много констант, связанных с порядком интегратора, что затрудняет обобщение кода на другие порядки;
3) интегратор реализован только для определенных порядков, причем для нечетных с разбиением Гаусса—Радо, хотя известно, что неявные коллокационные методы Рунге—Кутты, построенные на симметричных разбиениях Гаусса—Лобатто и Гаусса Лежандра, обладают геометрическими свойствами [5];
4) алгоритм выбора шага для уравнений первого порядка используется такой же, как и для уравнений второго порядка, поэтому шаг при интегрировании уравнений первого порядка выбирается неверно;
5) стартовый шаг интегрирования в режиме переменного шага выбирается независимо от дифференциальных уравнений, поэтому не всегда оптимально;
6) ограничения на величину выбираемого переменного шага, по-видимому, заданы в интеграторе просто из эмпирических соображений и, кроме того, не зависят от порядка интегратора,
В настоящей работе представлен новый код интегратора Эверхарта СIАI ЯЯ_1">. который позволяет устранить перечисленные выше недостатки:
1) путем использования возможностей Фортран 90 программный код сокращен почти в 2 раза;
2) устранены все константы, связанные с порядком метода (оставлены лишь константы узловых значений на шаге);
3) код позволяет получать решение 2-15 порядка точности (хотя при необходимости код без изменений можно обобщить на любой другой порядок: для этого нужно лишь получить соответствующие узловые значения);
4) исправлен алгоритм выбора переменного шага;
5) стартовый шаг выбирается по оценке интегрирующей схемы второго порядка с учетом поведения правых частей уравнений;
6) накладываются ограничения на выбираемый шаг в соответствии с порядком интегратора.
Кроме того, предложенный интегратор имеет новые возможности:
1) интегрирование на шаге до полной сходимости итерационного процесса;
2) запоминание величины предпоследнего шага после выполнения процедуры интегрирования при многократном использовании программного кода в режиме переменного шага;
3) быстрый выбор стартового шага, требуемый лишь для первого обращения к интегратору (при повторном обращении применяется запоминаемый шаг предыдущего обращения),
В работе представлена общая теория интегратора Эверхарта с внесенными автором коррективами, предлагается новый программный код интегратора, а также показаны результаты тестирования интегратора на примере дифференциальных уравнений задачи двух тел, В дальнейшем ввиду того что интегратор использует гауссовы разбиения, будем называть его интегратором Гаусса—Эверхарта (хотя обычно такие интеграторы называются гауссовыми [5]),
1. Основные формулы
Предположим, что на шаге к мы решаем задачу
х' = ¥(£, х), Хо = х(£о). (1)
Здесь £ — независимая переменная, х — интегрируемые переменные, ¥ — заданная вектор-функция £ и х. Введем перемеиную т = (£ — £0)/к и представим правую часть уравнений (1) в виде полинома степени к:
х' = хТ/к = ¥ = ¥о + А1Т + А2Т2 + Азт3 + ... + Ак тк, (2)
где коэффициенты А^ пока те определены, Интегрируя (2) по т, получаем решение
X = хо + к (ьт + I Агт2 + ^ А2г3 + ^Азг4 + ... + ^ Актк+1^ . (3)
то, т1, . . . , тк
(то = 0):
¥ = + а1т + а2т(т — т1) + азт(т — т1)(т — т2) + ... + акт(т — т1)... (т — тк-1). (4) Из соотношений
¥1 = + а1т1,
¥2 = + а1т2 + а2т2(т2 — т1),
¥з = + а1тз + а2тз(тз — т1) + азтз(тз — п)(тз — т2), ...
а
а1 = (¥1 — ¥о)/т1, а2 = ((¥2 — ¥о)/т2 — а1)/(т2 — т1), аз = (((¥з — ¥о)/тз — а1)/(тз — т^ — а2)/(тз — т2), ...
В свою очередь, сравнивая (2) и (4), будем иметь
А1 = а1 + (—т1 )а2 + (т1т2)аз + ... + (—1)к-1(т1... тк-1)ак, А2 = а2 + (—т1 — т2)аз + ...,
Ак = ак,
или
А1 = сиа1 + С21а2 + ... + Ск1ак, А2 = С22а2 + ... + Ск2ак,
Ак = скк ак.
(7)
Тогда обратный переход от А к а можно представить как
ai = diiAi + d2iA2 + ... + dkiAk, a 2 = d22A2 + ... + dk2Ak,
ak = dkk Ak. (8)
Коэффициенты Cj и dj являются числами Стирлинга и вычисляются по формулам cu = du = 1, cío = dio = 0 (i> 0),
Cij = Cí-ij-1 - Tí—iCí—i,j, dij = dí-ij-i - Tj di-i,j (i>j> 0). (9)
В соответствии с (9) каждый коэффициент cj представляет собой сумму всевозможных произведений i — j величин t1, ..., tí— i со знаком (—1)i—например,
С62 = TiT2T3T4 + TiT2T3T5 + Ti T2 T4 T5 + Ti T3 T4 T5 + T2T3T4T5.
Отсюда согласно теореме Внета значения ti ,..., tí— i будут являться корнями полинома вида
Pi—i = Cii + Cí2T + Cí3T 2 + ... + CííT i—i. (10)
2. Интегрирование на шаге
Величины а определяются по f, которые, в свою очередь, вычисляются по решениям x. Согласно (3) эти решения будем задавать как
XI = х0 + h ^foTi + + i A 2r¡ + ... + ,
xfc = x0 + h (Vfc + 1-KxtI + 1-K2ri + ... + ^AfcT^1^ . (11)
x
ются итерационным способом.
Для получения начального приближения а на следующем шаге h используется информация о коэффициентах A на текущем шare h. Безразмерная независимая переменная следующего шага будет f = (t — th) /h, вде = t0 + h. Отсюда
T = rf+1, (12)
где r = h/h. Согласно (2)
fo + Ait + A2T2 + A3T3 + ... + Ak t k = fo + Aif + A2f2 + A3f3 + ... + Akfk. (13)
Подставляя (12) в (13) и приравнивая коэффициенты при одинаковых степенях т, получаем
fo = eoofo + eioAi + e2oA2 + e3oA3 + ... + ekoAk, Ai = r(eiiAi + e2iA2 + e3iA3 + ... + ekiAk), A2 = r2(e22A2 + e32A3 + ... + ek2Ak),
Ak = rk ekk Ak,
(14)
где eij — числа арифметического треугольника, вычисляемые по рекуррентным формулам
ец = 6io = 1, бij = 6i-i,j-i + 6i-i,j (i > j > 0). (15)
Далее оценка A для h уточняется путем внесения поправки AA, получаемой как разность между значениями A после итераций и оценкой A на текущем шare h. Наконец, пользуясь соотношениями (8), будем иметь начальное приближение а.
Каждая итерация выполняется следующим образом. Сначала определяется решение xi, из которого по первой формуле (5) улучшается значение а1. Далее определяется Х2, по которому улучшаетея а2, и т. д. до xk, Как правило, для получения достаточно хороших а необходимы всего две итерации, очень редко — три.
Как только величины а получены, решение на шаге h (т =1) будет
xfc = х0 + h (f0 + tjAi + ^А2 + ... + • (16)
а
и запускается вышеописанный итерационный процесс. Если начальный шаг достаточно большой, чтобы обеспечить заданную локальную точность, то его следует уменьшить,
а
итерации,
3. Формулы интегратора как одно из представлений неявного метода Рунге—Кутты
Пользуясь (6), (7) и (11), нетрудно показать, что решение (16) предетавимо в виде
k k xh = Хо + h^&iki, где ki = f (to + hTi, xo + h^a^ kj) (i = 0,...,k), (17) i=0 j=0
а коэффициенты ay и bi — постоянные, зависящие только от Ti. Таким образом, интегратор Гаусса—Эверхарта фактически основан на видоизмененных формулах неявного метода Рунге—Кутты, причем эти формулы являются коллокационными [6]. Действительно, в качестве коллокационного полинома в интеграторе выступает приближенное представление решения (3), тогда как условия коллокации формально задаются соотношениями (5),
Произвольный выбор узлов Ti дает метод порядка k +1, Однако, используя свойства неявных методов Рунге—Кутты, Эверхарт выбрал разбиение Гаусса—Радо, что позволило повысить точность метода до порядка 2k + 1,
4. Повышение порядка интегратора
Значения узлов г\,...,гк — свободные параметры и их можно выбирать как угодно, лишь бы они не равнялись между собой, В общем случае интегратор будет иметь порядок к + 1, однако существуют такие узловые значения, которые позволяют значительно повысить точность решения — до порядка 2 к либо 2 к + 1,
Приближенное решение порядка 2к + 1 на сетке Т1,..., т2к будет
X/, = х0 + к ^о + ^ а; + ^А'2 + ... + . (18)
Определим Т1,..., тк так, чтобы решение (16) имел о порядок 2к + 1, Это возможно в случае, если разность решений (16) и (18) будет равна нулю:
, (А1 — А1 А2 — А2 Ак — Ак Ак ,1 А2к \ , ,
Ах/,, = /г —-- + -- + ... + —^-- + —^ + ... + , = 0. 19)
V 2 3 к + 1 к + 2 2к +1 ] К 1
Согласно (6) коэффициенты а1,..., ак в (16) и (18) совпадают. Тогда
А1 - А1 = Ск+1,1ак+1 + ... + С2к,1а2к, А2 - А2 = Ск+1,2ак+1 + ... + С2к,2«2к,
= Ск+1,к ак+1 + ... + с2к,к , Ак+1 = С^+1,к+1ак+1 + ... + С2к,к+1а2к, Ак+2 = Ск+2,к+2ак+2 + ... + с2к,к+2а2к,
А^ = С2к,2к а^. (20)
Используя рекуррентные соотношения (9), выразим в (20) коэффициенты через
Ск+2,1 = —1^+1^+1,ь = — Тк+1Ск+1,^ > 1),
Ск+3,1 = Тк+2Тк+1ск+1,1, = С^+1,^-2 + ( —Тк+1 — Тк+2)Ск+1а-1 + ^+2^+1^+1,^ > 1),
...
Подставляя (21) в (20), получим
А1 - А1 = В1Ск+1,1,
А2 - А2 = В1ск+1,2 + В2ск+1,1,
Ак - Ак = В1ск+1,к + В2ск+1,к-1 + ... + ВкС^+1,1, Ак+1 = В1ск+1,к+1 + В2ск+1,к + ... + Вк ск+1,2, Ак+2 = В2ск+1,к+1 + ... + Вк
(22)
В
В1 = ак+1 + (-Тк+1)ай+2 + (тк+1Тк+2)ак+з + ... + (-1)к-1(тк+1... Т2к-1)а2к,
В2 = ак+2 + (-Тк+1 - Тк+2)ак+3 + . . . ,
Вк = а2к.
Подставляя далее (22) в (19) и уравнивая коэффициенты при одинаковых величинах В, получим следующие уравнения для ск+1а-:
Ск+1,1 Ск+1,2 Ск+1,к 1 _ ф
2 3 "' к + 1 /с + 2 ~ '
ск+1,1 . ск+1,2 . . ск+1,к . 1
к + 1 к + 2 2к 2к + 1
0. (23)
Решения (23) однозначно определяют значения Т1,..., тк, которые вычисляются из уравнения
Ск+1,1 + Ск+1,2Т + Ск+1,эг2 + ... + ск+1уктк-1 + тк = 0. (24)
Узловые значения т», определяемые из (23) и (24), задают (левое) разбиение Гаусса— Радо, Их также можно получить из уравнения
(тк+1(т - 1)к)Тк) = 0. (25)
В (25) левый корень равен нулю (то = 0), Если потребуем, чтобы правый корень был тк = 1
биение Гаусса—Лобатто, узловые значения которого также удовлетворяют уравнению
(тк(т - 1)к)Тк-1) = 0. (26)
2к
5. Выбор шага
В интеграторе Гаусса—Эверхарта контроль шага интегрирования осуществляется по величине последнего члена в (16),
Пусть ||е4о1|| — заданная точность. Потребуем, чтобы на следующем шаге выполнялось равенство
1 ii а ii _ ii ii
к+1
Отсюда, используя последнее соотношение в (14), получим оценку
1
Очевидно, при разбиениях Гаусса—Радо и Гаусса—Лобатто недостаток такой оценки
к
говоря, она не обеспечивает сохранение требуемой локальной точности для решения более высокого порядка.
Во избежание слишком больших (и малых) локальных ошибок на г следует наложить ограничение:
- < гк+1 < а. (28)
а
Чтобы величина последнего члена в (16) была ограничена в пределах одного порядка, значение а должно быть равно л/10. Это следует из /гЦА^Ц ~ гк+1.
Выполнение обоих неравенств проверяется лишь в начале интегрирования при выборе стартового шага: если (28) не выполняется, то интегрирование повторяется с новым шагом к = ^г и т.д., пока не выполнится условие (28), Обычно для получения стартового шага требуется не более четырех итераций, В дальнейшем для ограничения г
г
чение правого предела.
Начальное приближение стартового шага получается из оценки (27) для к =1:
1= у ' ^(¿о + Мо + ^о), (29)
где к — малая величина. Если к настолько мала, что в компьютерной арифметике ^ = £0, то она увеличивается в 10 раз и оценка повторяется снова.
6. Порядок и шаг интегрирования
при компьютерной реализации метода
Теоретически совместное повышение порядка и уменьшение шага метода неограниченно повышает методическую точность численных результатов интегрирования. Однако при компьютерной реализации в арифметике с определенной точностью вследствие ошибок округления существуют такие значения параметров интегрирования, которые дают предельно высокую точность ввиду того, что методические ошибки становятся соизмеримыми с ошибками округления, и в этом случае не имеет смысла предпринимать какие-либо дальнейшие попытки получить более высокую точность численного интегрирования путем повышения порядка и уменьшения шага метода,
В общем случае получить оптимальные параметры интегрирования невозможно, поскольку они непосредственно зависят от специфики решаемой задачи. Впрочем, принимая во внимание, что интегратор Гаусса—Эверхарта используется, как правило, для численного моделирования орбит, найдем оценки оптимальной пары порядок—шаг интегрирования применительно к решению какой-нибудь простой задачи небесной механики, например, круговой двумерной задачи двух тел с уравнениями
г' = V, v' = —j^p-r. (30)
Здесь r = (ri, r2), v = (vi, v2) — векторы положения и скорости cooтвететвенно, ß — гравитационный параметр. Поскольку в круговом случае |r| = а = const, будем полагать, что величина ß/|r|3 = n2 = const, т.е. решение x = (r,v) описывается уравнениями гармонического осциллятора с частотой n и может быть представлено в виде
x1 = r1 = а cos nt, x2 = r2 = а sin nt,
x3 = v1 = —an sin nt, x4 = v2 = an cos nt. (31)
Оценим методическую ошибку |e|M по главному члену погрешности:
_ аУТТ^(пкГ+1 Мм---, (32)
где использована формула = <тор+1л/1 + п2. Здесь р — порядок метода. Ошибку
округления |е|д можно оценить как
|е|д = е|х| = еал/1 + п2, (33)
где е — машинный эпсилон.
Очевидно, что не имеет смысла выбирать такие шаг и порядок интегрирования, при которых методическая ошибка будет меньше ошибки округления. Из условия |е|м = |е|д получим отношение между оптимальными параметрам и интегрирования к и р:
к*"1 = е(р +1)!, (34)
где = пк — шаг то долготе ф = пЪ, соответствующий шагу к. Уравнение (34) дает
к
верхняя граница, задаваемая условием сходимости итерационного процесса для решения нелинейных уравнений в неявных методах [6]:
ь.<
L max¿ ^ | aij |'
где L — постоянная Липшица, которая задает неравенство
|f(t,x) - f(t,y)| < L|x - y| (35)
для любых x и y. Если положить maxi J^j |aij| = 1 (в действительности максимум близок к единице), то получим следующее ограничение на шаг интегрирования: h < 1/L или hv = n/L. Оценим постоянную Липшица L для исследуемой задачи. Рассмотрим отношение
|¿f|2 _ nA\5r\2 + |¿v|2 |¿r|2 + |¿v|2'
где Sx, ir и Sv — всевозможные разности векторов в соответствующих переменных. Принимая | Sr | = £ cos ф, |Sv| = £ sin ф, где £ > 0 и 0 > ф > п/2, будем иметь
= 4-l)cos2V + l.
|Sx|2
Отсюда следует, что все значения отношения лежат между 1 и n4. Следовательно, согласно (35) в качестве постоянной Липшица можно выбрать L = max(1,n2). Тогда получим верхнюю границу шага
n
К <-Тл—2Т - L
max(1, n2)
Таким образом, в лучшем случае, а именно при n =1, когда верхняя граница максимальна, шаг интегрирования hv должен удовлетворять неравенствам
ф +1)! < hp+1 < 1.
Очевидно, условие
ф +1)! > 1 (36)
означает, что порядок метода завышен и использование такого метода при вычислениях в арифметике с точностью е не логично в том смысле, что ту же точность результатов интегрирования можно получить с использованием методов более низких порядков. Оптимальные порядки р неявных методов Рунге—Кутты для различных е, соответствующих одинарной, двойной, расширенной и четверной точности, представлены в табличном виде (следует иметь в виду, что эти порядки получены для задачи с п = 1, в ином случае они могут быть меньше):
е 1.1 • 10"7 2.2 • 10"16 1.1-10"19 1.9 • 10"34
V 9 16 19 29
7. Фортран-код интегратора Гаусса—Эверхарта
Представленный выше фортран-код интегратора Гаусса—Эверхарта С ¡ AI SS_15 (до 15
порядка), основанного на гауссовых разбиениях Лобатто и Радо, а также Лежандра, доступен в интернете по адресу http://astro.tsu.ru/sch/Gauss_15.txt,
GAUSS_15 (X,TS,TF,STEP,ERR,N,NGR,NI,NS,NBS,NF,FUN) — заголовок процедуры интегрирования. Опишем входные и выходные параметры интегратора, помечая их соответственно индексами / и О. Итак, X¿ — интегрируемые переменные (до выполнения процедуры — начальные для значения TS, после — конечные для TF); TS7 и TF7 — начальное и конечное значения независимой переменной соответственно; STEP¿ — величина стартового шага интегрирования (если STEP = 0, то стартовый шаг выбирается по формуле (29)), после выполнения процедуры STEP — величина предпоследнего шага; ERR7 — значен не ||eto11| для выбора переменного шага (27) (условие ERR = 0 задает ре-
=0II интегратора р; NI7 — максимальное число итераций на шare (при N1 < 0 итерационный процесс выполняется до сходимости); NSO — число шагов интегрирования за выполнение процедуры; NBSO — число шагов, на которых итерационный процесс не сходится (несходимость итерационного процесса регистрируется, когда число итераций достигает 100); NFo — число обращений к процедуре правых частей; FUN7 — название процедуры правых частей f.
8. Тестирование интегратора Гаусса—Эверхарта в задаче двух тел
Интегратор тестировался на дифференциальных уравнениях задачи двух тел (30). Начальные условия задачи
ri = 1 - е, г2 = 0, ы = 0, v2 = л/(1 + е)/(1 - е)
соответствуют однопараметричеекому семейству орбит с эксцентриситетом e в качестве параметра и с единичной большой полуосью а. Интегрирование выполнялось на интервале 1000 оборотов для различных эксцентриситетов. При этом исследовались те или иные характеристики интегратора.
8.1. Круговой случай е = 0
Как уже было замечено, шаг в интеграторе выбирается таким образом, чтобы сохранялась величина члена порядка к + 1. Если бы интегратор имел порядок к, то таким способом можно было бы обеспечить сохранение локальной точности на всем интервале
к
На примере круговой задачи (в = 0) исследовалась зависимость реальной локальной точности от задаваемой для выбора переменного шага. В (16) величину к + 1 член а ||ек || можно оценить как
к
ек || =
А
кк+1
х
(к+1)|
к +1"' (к + 1) тогда как ошибка численного решения порядка р будет
,, ,, кр+1
х
(р+1) |
" (Р +1)
х
Цх^Ц = ал/1 + п2п\
Подставляя (39) в (37) и (38), получим
-- -р — к
(алЛ +
11 Р11 (р +1)! В частности, для разбиения Гаусса—Радо (р = 2к + 1)
1
р+1
((А; + 1)!|Ы|) —
lv/ГT
к+1 = п г
п2 к+1+г
г=1
ек |
(37)
(38)
(39)
(40)
<
0.001
1е-005
1е-007
в 1е-009 ^
о
н
| 1е-011 л
(3 1е-013
1е-015 -
1е-017 1е-010
1е-008
1е-006
0.0001
0.01
е^
2
е
Р
Рис. 1. Зависимость локальной точности от параметра ЕП Е
На рис, 1 (см, с, 41) представлены оценки реальной локальной ошибки
||ер|| = || Дх|| = у/Дг2 + Дг2 + At'2 + Дг|
в зависимости от ERR ||ek||) для p = 11. Оценивание выполнялось в режиме STEP = 0 на первом шаге, подбираемом в соответствии с параметром ERR, Видно, что полученные результаты хорошо согласуются с (40). В частности, экспериментальная оценка при ERR > 10-6 (ERR < 10-6 — область доминирования вычислительных ошибок) подтверждает квадратичную зависимость ||ep|| a ||ek||2, отчего главным образом реальная точность численного решения существенно выше задаваемого значения параметра ERR.
8.2. Слабоэксцентричный случай е = 0.1
Дня оценки оптимизации составленного фортран-кода сравнили быстродействия (в затратах процессорного времени) интегратора САи88_15 и широко используемого в небесной механике интегратора ИАБАи_27. При этом в обоих случаях выполнялись все операции, определяемые основными формулами. В то же время расхождение двух численных решений, полученных интеграторами, оказалось меньше методической ошибки па несколько порядков.
На рис. 2 показаны отношения временных затрат, которые потребовались интеграторам дня получения численных решений 7, 11 и 15 порядков. Как видно, интегратор САи88_15 работает быстрее, однако этот выигрыш с увеличением порядка уменьшается. Впрочем, следует заметить, что в задачах со сложной функцией ¥ оперативность интеграторов должна быть одинаковой, поскольку большая часть времени будет затрачиваться па процедуру ее вычисления РОТ и тогда быстродействие главным образом определяется числом обращения к этой процедуре.
1.3
1.2
к н о о
О «
о о с
к
я m
1.1
11
Порядок метода p
15
Рис. 2. Оптимальность интегратора GAUSS_15 относительно R ADAL27
1
7
Далее была оценена ошибка интегрирования |Дг| = л/Аг2 + Аг| в зависимости от величины постоянного шага к для порядков 2-11. Интегрирование выполнялось с шестью итерациями (на шаге) на интервале 1000 оборотов.
В случае к = 1-3 представленные па рис. 3 характеристики показывают, что при достаточно больших к (глобальная) ошибка решения рпорядка, как и ожидалось, с уменьшением шага ведет себя как |Дг| ~ кр, что указывает на соответствие ошибки методу порядка р. Для к = 4 и 5 ситуация несколько иная. Это объясняется тем, что шести итераций недостаточно для сходимости итерационного процесса на шаге и получаемое решение не соответствует порядку интегрирующих формул. Случайное поведение характеристик ниже |Дг| = 10-6 обусловлено влиянием ошибок округления.
Следует заметить, что несмотря на крутой скат характеристики нечетного порядка 2к + 1 в сравнении с характеристикой четного порядка 2к (для к = 1-3) определенная точность для четного порядка достигается при большем шаге, т. е. быстрее. Это связано с тем, что при симметричном разбиении Гаусса—Лобатто интегратор Гаусса—Эверхарта обладает геометрическими свойствами и методическая ошибка вдоль независимой переменной £ эволюционирует медленнее, чем при разбиении Гаусса—Радо.
На рис. 4 показано поведение ошибки с изменением £ для четных и нечетных порядков (к = 3-5). Приведенные результаты были получены в режиме N1 < 0, т.е. до полной сходимости итерационного процесса. При этом на шаге требовалось от 11 до 16 итераций. Как следует из рисунка, при постоянном шаге (к = 2п/16) ошибка интегрирования для четных порядков ведет себя линейно, тогда как для нечетных порядков — квадратичееки. Даже несмотря на то что интегратор порядка 2к + 1 в начале интегриро-
2к
разного роста ошибок первый уже уступает второму в точности почти на два порядка. Очевидно, с увеличением интервала интегрирования это преимущество для каждого интегратора с симметричным разбиением будет возрастать. К сожалению, такие свойства интеграторов четных порядков имеют место только при использовании постоянного
Шаг интегрирования к
Рис. 3. Зависимость точности от величины шага интегрирования
Оборот
Рис. 4. Поведение ошибки для решений четных и нечетных порядков
Быстродействие NF
Рис. 5. Характеристики точность—быстродействие в зависимости от числа итераций на шаге
шага [5]. При автоматическом выборе шага ошибка так же, как и для интеграторов нечетных порядков, ведет себя квадратичным образом.
Далее был исследован вопрос о соотношении точности и быстродействия при увеличении числа итераций на шаге. Путем вариации ERR оценивались точность и быстродействие интегратора 11 порядка с переменным шагом на интервале 1000 оборотов для N1 = 1-6. Результаты приведены на рис. 5. Видно, что уже на второй итерации можно получить довольно хорошее решение, хотя соответствующая характеристика заметно отклоняется от степенной зависимости. Поэтому для уверенного результата необходимо использовать, по крайней мере, три итерации на шаге. Следует также заметить, что увеличение числа итераций не существенно понижает эффективность интегрирования.
8.3. Сильноэксцентричный случай е = 0.9
Очевидно, еилыюэкецентричные орбиты необходимо интегрировать с переменным шагом. Поведение выбираемого шага на одном обороте для е = 0.9 представлено на рис. 6, где также приведены функции, пропорциональные |г|, |г|3/2, |г|2 и |V|—1. Как видно, изменение шага очень близко к поведению | г |3/2. Это означает, что шаг переменной величины по £ соответствует почти постоянному шагу по так называемой эллиптической аномалии.
я
cd я о а
s &
и Ё s
0.1
0.01
0.001
- \-Jri2
1 jr j3'2- Л ~ l r |
_
= / |v j"1
40
80 120 Номер шага
160
200
Рис. 6. Выбор шага интегрирования в случае сильновытянутой орбиты
<
л
н о
1000000
10000000 Быстродействие NF
100000000
Рис. 7. Характеристики точность быстродействие для интеграторов GAUSS_15 и RAT) AU27 в экстраэксцентричном случае
1
8.4. Экстраэксцентричный случай e = 0.999
Наконец, был исследован алгоритм выбора шага при долгосрочном интегрировании очень вытянутой орбиты с эксцентриситетом e = 0.999, Интегрирование выполнялось на интервале 1000 оборотов двумя интеграторами: СIAI SS_1"> и RADAU_27. Результаты показали (рис, 7), что RADAU_27 заметно менее эффективен, нежели СIAI SS_1">.
Так, при одинаковом быстродействии первый интегратор дает решение с точностью ниже на порядок и более. Это связано с тем, что в RADAU_27 заложен неверный алгоритм выбора шага для интегрирования систем с уравнениями первого порядка, а именно, алгоритм не точно соответствует формуле (27): вместо степени 1/(k +1) в RADAU_27 используется 1/(k + 2), Очевидно, эта ошибка легко устраняется.
Таким образом, в работе представлена новая версия интегратора Гаусса—Эверхарта
СIAI SS_1"). На примере плоской задачи двух тел при различных начальных условиях
экспериментально показано, что предложенная программная реализация интегратора значительно эффективнее ранее используемой, а также предоставляет больше возможностей для высокоточного и оперативного численного интегрирования.
Список литературы
[1] Everhart Е. A new method for integrating orbits // Bull. Amer. Astronom. Soc. 1973. Vol. 5. P. 389.
[2] Everhart E. Implicit single sequence methods for integrating orbits // Celest. Mech. 1974. Vol. 10. P. 35-55.
[3] Everhart E. An efficient integrator that uses Gauss^Radau spacings // Dynamics of Comets: Their Origin and Evolution. Proc. of IAU Colloq. 83. Italy, 1984 / Eds. A. Carusi and G.B. Valsecchi. Dordrecht: Reidel, Astrophvs. and Space Science Library Rame, 1985. Vol. 115. P. 185-202.
[4] Butcher J.C. Implicit Runge^Kutta processes // Math. Comput. 1964. Vol. 18. P. 50-64.
[5] Hairer E., Lubich C., Wanner G. Geometric Numerical Integration: Structure-Preserving Algorithms for Ordinary Differential Equations / Springer Series in Comput. Math. Springer, 2002. 536 p.
[61 Hairer E., Norsett S.P., Wanner G. Solving Ordinary Differential Equations. Nonstiff Problems / Springer Series in Comput. Math. Springer, 1993. 544 p.
Поступила в редакцию 27 августа 2009 г., с доработки — 26 сентября 2009 г.
