CC BY
113
УДК 519.653:521.182
ИССЛЕДОВАНИЕ ЭФФЕКТИВНОСТИ АЛГОРИТМОВ МЕТОДА ЭВЕРХАРТА С ВЫСОКИМ ПОРЯДКОМ АППРОКСИМИРУЮЩИХ ФОРМУЛ
А. А. Заусаев
Самарский государственный технический университет,
443100, Россия, Самара, ул. Молодогвардейская, 244.
E-mail: zausaevaa@mail. ru
Разработан модифицированный алгоритм численного интегрирования уравнений движения небесных тел методом Эверхарта. Проведено исследование эффективности алгоритма для больших порядков аппроксимирующих формул. Показана высокая эффективность метода на примере совместного интегрирования уравнений движения больших планет, Луны, Солнца и малых тел Солнечной системы.
Ключевые слова: метод Эверхарта, численное интегрирование, орбита, малые тела Солнечной системы.
Введение. Метод Эверхарта [1] относится к числу неявных одношаговых методов, что обеспечивает его сходимость и устойчивость. Основным достоинством одношаговых методов является то обстоятельство, что для них разработаны надежные оценки локальной погрешности дискретизации. Алгоритмы численного интегрирования методом Эверхарта разработаны до 27-го порядка, однако использование данных алгоритмов при порядках выше 19-го не приводит к повышению точности вычислений.
В данной работе предлагается модификация алгоритма метода Эверхарта, позволяющая повысить его эффективность при увеличении порядка.
1. Основные уравнения. Рассмотрим основную идею построения метода Эверхарта на примере решения задачи Коши:
х = F(x,t); x(t i) = a?i, x(t\) = x\. (1)
Представим правую часть уравнения (1) в виде временного ряда
х = F(x, t) = F\ + A\t + A2t^ + • • • + Antn, (2)
где F\ = F(x\,ti), Ai — некоторые коэффициенты. Интегрируя (2), получим выражения для определения координат и скоростей:
^2 ^3 -j-n-\-2
x = xl+xlt + Fl- + Al- + --- + An{n + i){n + 2y (3)
t2 t3
x = x\ + F\t + A\— + A2— + • • • + An——. (4)
2 3 n +1
Полиномы (3) и (4) не являются рядами Тейлора, поскольку коэффициенты Ai вычисляются не по известным формулам коэффициентов ряда Тейлора, а
Артём Анатольевич Заусаев (к.ф.-м.н., доц.), доцент, каф. прикладной математики и информатики.
определяются из условия наилучшего приближения X и х с помощью конечных разложений (3) и (4). Данные условия будут введены позднее.
Выразим коэффициенты Ai через разделенные разности. Для этого представим функцию F(x,t) в виде
F(x, t) = Fi + a\t + a2i(i - ¿2) + a3t(t - i2)(i -t3) + ... . (5)
Уравнение (5) усечено по времени tn. В каждый фиксированный момент времени ti имеем
F2 = F\ -\- a\t2,
F3 = Fi + aiU + a2t3(i3 - t2),
Принимая tnj = tn — tj, найдём од через разделённые разности:
«1 = (F2 — Fi)/t2,
«2 = ((F3 — Fi)/ts — ai)/ts2,
аз = ((№ - Fi)/U ~ oli)/U2 ~ a2)/t43,
a4 = ((((-^5 — F\)jtb — CKl)/i52 — «2)7*53 — 013)/t 54,
Приравнивая коэффициенты при одинаковых степенях £ в уравнениях (2) и (5), выразим коэффициенты А* через од:
^4-1 = ОД + (—¿2)02 + (t2tз)o'з + • • • = СпОД + С21ОД + С31ОД + . . . ,
А2 = ОД + (—¿2 — *з)«3 + • • • = С22СК2 + С32ОД + • • • ,
Аз = а,з + ■ ■ ■ = С33ОД + ...,
(6)
Коэффициенты Су определяются из следующих рекуррентных соотношений:
Cij — 1)
Cjl — tiCi—11,
* =
г > 1,
— Cj_ij'_i tiCi— ij, 1 <C j < i.
В частности, для алгоритма интегрирования пятого порядка имеем С41 = — ¿2*3*4, С42 = ¿2*3 + *3*4 + *4*2, С43 = —¿2 — ¿3 “ *4,
где ¿2, *з, *4 ~ корни кубического уравнения
(-¿2*3*4) + (¿2*3 + *3*4 + *4*2)* + (-¿2 - ¿3 - *4)*2 + (I)*3 = 0.
(7)
2. Алгоритм интегрирования. Вопрос нахождения узлов разбиения отрезка [0, Т], равного по длине шагу интегрирования, рассмотрим на примере алгоритма интегрирования пятого порядка.
В начальный момент времени ¿1=0 известны Х\,Х1,Р\. Значения х в моменты времени ¿2) *з,*4 определяются с помощью предсказывающих уравнений:
Х2= Х\+ Х\i2 + Хз=Х\+ Х\¿3 +
Х4 = х\ + xit^ +
m
2
м
2
^1*1
+
A\t^ ¿2*2 ¿3*2
+
+
6
¿1*1 6
¿lit
+
+
12
м
12
+
+
+
¿2*4
12
+
20
м
20
Asi|
и двух исправляющих уравнении для нахождения положения и скорости на конце отрезка [0,Т]:
, ч ^1 т2 АгТ3 А2Та А3Т5
х{Т) -Х1+Х1Т + —+ — + — + —, (9)
^ ^ ^ АгТ2 А2Т3 А3Т4 , Л
х(Т) — Х\ + Р{Т н -----1 ---1 -—. (10)
Эта схема является неявной, так как коэффициенты, стоящие в квадратных скобках (8), неизвестны при первой итерации.
Уравнения (8)—(10) обеспечивают пятый порядок точности относительно Можно увеличить порядок точности в вычислении ж и ж до седьмого порядка путём специального выбора подшагов Ь2, ¿3) ¿4- С этой целью увеличим количество разбиений интервала интегрирования, добавив два дополнительных значения времени ¿5, ¿6- Затем вычислим для ¿5 и ¿6 значения й4 и «5, а также новые значения А!а , А!ъ и А!х, А'2, А!г.
Из уравнения (9) можно найти поправки Аж, улучшающие значения координат:
дж = (¿;1 ~ М)ТЪ (А'2-А2)Т4 (4-А3)Г5 А^ А^_
6 12 20 30 42 ' 1 ]
Выражая в уравнении (7) С51, ..., С54 через С41, С42, С43, а также полагая Л-2 = Ь2/Т, Л,з = £з/Т, ^=¿4/Т,
выражение (11) можно записать в виде
Ах = (0:4 — ¿бС^Т1
С41 , ^42 , ^43 ,
6 12 20 30
+ «5 Т7
С41 , ^42 , ^43 ,
12 20 30 42
Значение Аж в последнем выражении можно обратить в ноль при выполнении следующих условий:
^1^42 ^з1=0 ^41+^42+^43+1=0. (12)
6 12 20 30 ’ 12 20 30 42 у 7
Проводя подобные рассуждения для скорости, приравнивая к нулю Аж, получим третье условие для определения с'А1, С42, С43. Тогда соответствующие данным разбиениям коэффициенты с- будут определяться из системы алгебраических уравнений
^+£?+тЧ=°- 41+х+х+й=°- а«
2345 3456 4567
Первое и второе уравнения (12) являются линейными комбинациями уравнений (13), поэтому ими можно пренебречь.
Из решения системы (13) получим
С41 = -4/35 = —Л.2Л.3Л4, ¿4 2 = 6/7 = Н2Н3 + Л,3Л,4 + Л,2Л.4, (1
С43 = -12/7 = -Н2 - Л-з - Н4. ^ '
Следует отметить, что значения величин Л-2, hs, Л4 являются корнями полинома третьей степени
"3-fft2 + ?'l-é = 0 (15>
и имеют следующие значения:
h2 = t2/T = 0,212340538239 ... , h3 = t3/T = 0,590533135559 ... , h4 = U/T = 0,911141240488 ....
Использование данных узлов позволяет получить решение уравнения (1) с точностью до седьмого порядка для обеих компонент х и х. Полученные по формуле (15) узлы разбиения hi совпадают с узлами квадратурной формулы Гаусса—Радо. Область изменения hi заключена в пределах 0 ^ hi ^ 1.
Порядок метода, определяющий точность интегрирования, зависит от количества разбиений основного шага h на подшаги hi. Порядок данного метода определяется наибольшей степенью временного ряда (3), увеличенной на два. Таким образом, 11-му порядку метода в обозначениях Эверхарта соответствует 9-й порядок, 15-му — 11-й, 19-му — 13-й и 27-му — 17-й порядок точности. В дальнейшем под величиной порядка будем понимать величину, принятую Эверхартом.
В работе Эверхарта [1] приведены узлы разбиения основного шага интегрирования на подшаги, обеспечивающие точность до 15-го порядка включительно. Узлы разбиения отрезка [0,1] на подшаги, вычисленные с помощью алгоритма, основанного на формулах типа (13)—(15), для порядков метода с 9-го по 33-й, приведены в работе [2].
При численном интегрировании уравнений (1) методом Эверхарта используются формулы как открытого, так и замкнутого типов, при этом разбиение шага интегрирования производится по формулам Гаусса—Радо или Гаусса— Лобатто. Алгоритмы численного интегрирования обыкновенных дифференциальных уравнений реализованы Эверхартом для первого случая.
Как известно, при численном решении дифференциальных уравнений повышения точности можно достичь двумя способами: либо путём уменьшения шага интегрирования, либо путём увеличения порядка метода. Недостатком первого способа является увеличение ошибок округления. Второй же способ приводит к усложнению алгоритма и замедлению процесса вычислений. Кроме того, коэффициенты высокого порядка, вычисленные по разностной схеме, часто отягощены большими ошибками, т.е. их значения получаются на уровне шумов. Однако всегда, когда это возможно, следует отдавать предпочтение увеличению порядка.
3. Модификация алгоритма метода Эверхарта. Для пространственного случая коэффициенты Ai вычисляются для каждой из компонент х, у, z. Обозначим соответствующие коэффициенты Aix, Aiy, AiZ.
При порядке метода не выше 15-го, значения величин Aíx, Aiy, Aíz вычисляются с высокой степенью точности, и в модификации алгоритма численного интегрирования нет необходимости.
Неопределённость остается в выборе шага интегрирования для конкретно решаемой задачи. Если Anx,Any,Anz — конечные коэффициенты для нахождения координат х, у, z в соотношении типа (3), то для определения макси-
мальной длины шага интегрирования h может быть использована формула [4] h = {i/R)l^n, R = (\Апх\ + \Апу\ + \Anz\)/(n(n + 1)), (16)
где £ — заданная пользователем абсолютная погрешность.
При использовании методов низкого порядка величина шага интегрирования мала, что не позволяет численно интегрировать, к примеру, уравнения движения небесных тел на интервале времени порядка нескольких столетий. С целью увеличения шага интегрирования применяют методы более высокого порядка.
Однако при использовании метода Эверхарта для решения уравнений движения небесных объектов увеличение порядка метода выше 19-того не приводило к повышению точности и эффективности вычислений. Главная причина заключается в способе нахождения коэффициентов Апх, Апу, Anz.
Вследствие того, что начальные данные задачи Коши имеют ограниченную точность, разделенные разности высоких порядков, вычисленные с помощью этих данных, теряют всякий физический смысл. В силу неопределенности коэффициентов высокого порядка, вычисленных с использованием разностных схем, нахождение решения по формулам (3) и (4) будет приводить к неверным результатам. Следовательно, задача численного интегрирования дифференциальных уравнений методом Эверхарта высокого порядка сводится к оценке коэффициентов Апх, Апу, Anz.
Оценить указанные коэффициенты можно на основании выполнения нижеследующих условий.
1. При возрастании порядка метода величина последнего члена в формуле (3) должна убывать и при п —> оо стремиться к нулю. Выполнение этого условия обеспечивает сходимость временного ряда.
2. Увеличение порядка метода должно приводить к увеличению максимальной длины шага интегрирования. Выполнение этого требования менее очевидно, однако в случае нарушения данного условия увеличение порядка не может привести к более эффективному алгоритму и дальнейшее повышение порядка метода не имеет смысла.
Увеличение эффективности метода при повышении его порядка может быть достигнуто следующим способом. Если Апх, Апу, Anz — конечные коэффициенты для координат х, у, z в формуле (3), то максимальная длина шага определяется из соотношения (16). При этом значение величины R следует задать таким образом, чтобы шаг интегрирования увеличился на величину, заданную пользователем.
Так как система уравнений для нахождения А\, А2, А3, ... записана в неявном виде, с учётом сделанного допущения производится пересчёт коэффициентов Ai при г < п, которые находятся по формулам (6).
Численные расчёты показали, что при высоком порядке метода значение R в формуле (16) следует выбирать достаточно малым независимо от решаемой задачи. При значениях R —> 0 полученное решение приближается к точному решению.
4. Применение модифицированного алгоритма при решении задачи эволюции малых тел Солнечной системы. Рассмотрим изложенный выше алгоритм на примере решения задачи численного интегрирования уравнений движения больших планет, Луны, Солнца и малых тел Солнечной системы.
Дифференциальные уравнения движения небесного тела в гелиоцентрической системе координат с учётом релятивистских эффектов от Солнца имеют следующий вид [5]:
d2X ;2Л1 .X (Xi-X ХЛ
ж—* —тг)+
к2 г к2 X2 (XX)2
+ ^2 [(4 - 2«)^ - (1 + «)^Х + 3«11х + (4
где X — матрица-столбец с элементами ж, у, z\ X¿ — матрица-столбец с элементами Xi, у i, Zi; X — матрица-столбец с элементами ж, у, m, х, у, z — масса и гелиоцентрические координаты возмущаемого тела; m¿, Xi, %i — массы и гелиоцентрические координаты больших планет; г2 = х2 + у2 + z2, А2 = (a;¿ - ж)2 + {yi - у)2 + (z¿ - z)2, г2 = х2 + у2 + z2; к — постоянная Гаусса, с — скорость света, a — параметр, характеризующий выбор системы координат. Случай a = 1 соответствует стандартным координатам, случай a = 0 — гармоническим координатам.
Проведено численное интегрирование уравнений движения (17) как классическим, так и модифицированным методом Эверхарта с целью исследования эволюции орбиты астероида Aten 2003 НТ42.
В табл. 1-7 приведены абсолютные значения разностей коэффициентов \Ai~А1^ для внутренних планет Луны, Солнца и астероида Aten 2003 НТ42 на эпоху 2006 03 06, где Ai — коэффициенты, вычисленные классическим методом, a A!i— коэффициенты, вычисленные модифицированным методом 27-го порядка.
Стоит отметить, что чем выше порядок коэффициента, тем существеннее его изменение в результате модификации алгоритма. При этом каждый из коэффициентов уменьшается по абсолютной величине, что непосредственно отражается на уменьшении абсолютной погрешности на каждом шаге численного интегрирования дифференциальных уравнений движения.
Использование предложенной модификации позволило разработать алгоритм и программный комплекс, позволяющий выполнять численное инте-
2а)Щ±Х
(17)
Таблица 1
Разности коэффициентов | Ai — A!i | для Солнца на эпоху 2006 03 06
i X У Z
1 1,99999589756 ю-24 0 0
2 9,05000036914 Ю-гз 1,32200000513 10 22 2,98800000021 10-23
3 1,70031799998 ю-21 2,48523699999 Ю-21 5,61614999973 10 22
4 1,70905135000 ю-20 2,49800280000 ю-20 5,64499653000 10~21
5 1,04910645049 Ю-19 1,53340673609 Ю-19 3,46519855809 Ю-20
6 4,21800546854 Ю-19 6,16516845862 Ю-19 1,39320718700 кг19
7 1,15414274226 Ю-18 1,68693105886 Ю-18 3,81213342502 кг19
8 2,18758757429 Ю-18 3,19744628453 Ю-18 7,22560165804 кг19
9 2,87350440171 Ю-18 4,20000373052 Ю-18 9,49118490777 кг19
10 2,56738980545 Ю-18 3,75257708120 Ю-18 8,48008841024 кг19
И 1,48882415949 Ю-18 2,17611186542 Ю-18 4,91758613082 кг19
12 5,05313532605 Ю-19 7,38582032710 Ю-19 1,66905057513 кг19
13 7,61883029409 Ю-20 1,11359201811 Ю-19 2,51649961137 Ю-20
Таблица 2
Разности коэффициентов | А* — А!1 | для Меркурия на эпоху 2006 03 06
г X У г
1 0 0 0
2 1,70000041910 ю-17 1,70000174258 ю-18 1,19899999362 КГ17
3 3,24019999941 Ю-16 3,15900002964 КГ17 2,25339999841 Ю-16
4 3,25683610000 ю-15 3,17566899995 Ю-16 2,26499980000 ю-15
5 1,99921891100 ю-14 1,94939459500 ю-15 1,39037716208 ю-14
6 8,03799871400 ю-14 7,83767658130 ю-15 5,59010810600 ю-14
7 2,19938023957 ю-13 2,14456752803 ю-14 1,52958139738 ю-13
8 4,16875374904 ю-13 4,06486052400 ю-14 2,89920227045 ю-13
9 5,47586409274 ю-13 5,33939520565 ю-14 3,80824547714 ю-13
10 4,89251996240 ю-13 4,77058911463 ю-14 3,40255285805 ю-13
И 2,83716243824 ю-13 2,76645498611 ю-14 1,97313352581 10-13
12 9,62945533335 ю-14 9,38947109959 ю-15 6,69690282682 Ю-14
13 1,45187455457 Ю-14 1,41569109558 ю-15 1,00972095227 Ю-14
Таблица 3
Разности коэффициентов |А* — А^\ для Венеры на эпоху 2006 03 06
г X У г
1 0 0 0
2 9,52999994960- ю-18 2,15999998315- ю-18 7,89000023429 • ю-19
3 1,79033330000- Ю-16 4,05678000002 • ю-17 1,48215799996- ю-17
4 1,79952980920- ю-15 4,07762024200 • Ю-16 1,48977196920- Ю-16
5 1,10464693240- ю-14 2,50305978227 • ю-15 9,14501125383- Ю-16
6 4,44130983989- ю-14 1,00637259878- ю-14 3,67681539476- ю-15
7 1,21524392419- ю-13 2,75366554060 • ю-14 1,00606076360- ю-14
8 2,30340010031 • ю-13 5,21935831664- ю-14 1,90690973034- ю-14
9 3,02562987875 • ю-13 6,85588511893- ю-14 2,50482018100- ю-14
10 2,70330934630 • 10-13 6,12552726604- ю-14 2,23798153689- ю-14
И 1,56764362653- 10-13 3,55218088186- ю-14 1,29780097028- ю-14
12 5,32065209809 • Ю-14 1,20562596895- ю-14 4,40479413725 • ю-15
13 8,02217688103- Ю-15 1,81777432482- ю-15 6,64129829239- Ю-16
Табли
Разности коэффициентов | А* — А' | для Земли на эпоху 2006 03 06
г X У г
1 2,00005654670 Ю-19 0 0
2 9,78700000284 Ю-18 1,61299991973 ю-18 1,25000018651 ю-19
3 1,83962900002 Ю-16 3,03079699996 Ю-17 2,34543000053 ю-18
4 1,84907818200 Ю-15 3,04636523000 Ю-16 2,35747724000 Ю-17
5 1,13506235422 Ю-14 1,87002070750 ю-15 1,44714468780 Ю-16
6 4,56359715927 Ю-14 7,51854834800 Ю-15 5,81834589000 Ю-16
7 1,24870453091 10-13 2,05724674213 Ю-14 1,59203247410 Ю-15
8 2,36682207128 10-13 3,89935078694 Ю-14 3,01757341680 Ю-15
9 3,10893777228 10-13 5,12198998658 Ю-14 3,96373182830 Ю-15
10 2,77774244494 10-13 4,57634408603 Ю-14 3,54147523880 Ю-15
И 1,61080722999 10-13 2,65381268670 Ю-14 2,05369433371 Ю-15
12 5,46715128546 Ю-14 9,00715813250 Ю-15 6,97032978707 Ю-16
13 8,24306003073 Ю-15 1,35804811895 Ю-15 1,05094671557 Ю-16
Таблица 5
Разности коэффициентов | Ai — A!i | для Луны на эпоху 2006 03 06
i X У Z
1 0 3,00003188050 ю-18 1,00000180357 ю-18
2 4,28025999988 ю-17 1,12499991144 Ю-16 3,94999984492 кг17
3 8,04510000093 Ю-16 2,11482199995 ю-15 7,42231999943 Ю-16
4 8,08649820000 ю-15 2,12568496000 ю-14 7,46044923000 ю-15
5 4,96392185010 ю-14 1,30485828080 ю-13 4,57961980640 ю-14
6 1,99577931338 ю-13 5,24627349430 ю-13 1,84126800270 ю-13
7 5,46090854263 ю-13 1,43550038584 кг12 5,03813026740 ю-13
8 1,03507263306 ю-12 2,72087904896 кг12 9,54938307640 ю-13
9 1,35961906266 10~12 3,57400911196 10~12 1,25435866548 10~12
10 1,21477876236 Ю-12 3,19326970689 кг12 1,12073176162 кг12
И 7,04447749237 10-13 1,85177064942 кг12 6,49910083556 10-13
12 2,39092695019 10-13 6,28499183382 10-13 2,20582369051 10-13
13 3,60490378814 Ю-14 9,47615353480 Ю-14 3,32581561196 Ю-14
Таблица 6
Разности коэффициентов | Ai — A!i | для Марса на эпоху 2006 03 06
i X У Z
1 0 0 0
2 6,64730000069 ю-18 1,71020000066 ю-18 1,97970000018 ю-18
3 1,24943152074 Ю-16 3,21444537000 ю-17 3,72100184000 ю-17
4 1,25584950947 Ю-15 3,23095709580 Ю-16 3,74011561090 Ю-16
5 7,70907100900 Ю-15 1,98333299420 ю-15 2,29588152160 ю-15
6 3,09948562970 Ю-14 7,97412827970 ю-15 9,23075137750 ю-15
7 8,48090139040 Ю-14 2,18190382838 ю-14 2,52574464107 ю-14
8 1,60748872917 10-13 4,13562857385 ю-14 4,78735202350 ю-14
9 2,11151589690 10-13 5,43235253784 ю-14 6,28842350069 ю-14
10 1,88657598175 10-13 4,85364369615 ю-14 5,61851736803 ю-14
И 1,09402159905 10-13 2,81461817017 ю-14 3,25816686666 ю-14
12 3,71315790008 Ю-14 9,55293908579 ю-15 1,10583630627 ю-14
13 5,59848847704 Ю-15 1,44033786964 ю-15 1,66731714210 ю-15
Таблица 7
Разности коэффициентов |A¿ — А^\ для астероида Aten 2003 НТ42 на эпоху
2006 03 06
i X У Z
1 0 0 0
2 3,30889999678 ю-18 1,29100001396 ю-18 1,43999995430 КГ19
3 6,21942010000 ю-17 2,42642000002 ю-17 2,71540899989 ю-18
4 6,25136755017 Ю-16 2,43888298600 Ю-16 2,72935725000 ю-17
5 3,83742128210 ю-15 1,49711585502 ю-15 1,67542437945 Ю-16
6 1,54286192267 ю-14 6,01925844660 ю-15 6,73616027357 Ю-16
7 4,22162300085 ю-14 1,64700674330 ю-14 1,84316747537 ю-15
8 8,00175721940 ю-14 3,12177285746 ю-14 3,49358022957 ю-15
9 1,05107035994 ю-13 4,10060294385 ю-14 4,58899030392 ю-15
10 9,39099771453 ю-14 3,66376546630 ю-14 4,10012489161 ю-15
И 5,44582059546 ю-14 2,12461018943 ю-14 2,37765413828 ю-15
12 1,84833569868 ю-14 7,21102135125 ю-15 8,06986375306 Ю-16
13 2,78681553255 ю-15 1,08723682184 ю-15 1,21672819871 Ю-16
грирование уравнений движения небесных объектов методом Эверхарта до 27-го порядка включительно. Модифицированный метод Эверхарта позволяет учитывать в аппроксимирующей формуле различное количество членов, в зависимости от требуемой точности.
Применение данного алгоритма к решению систем дифференциальных уравнений движения небесных объектов выявило преимущество увеличения порядка аппроксимирующих формул в алгоритме по сравнению со способами, основанными на процедуре уменьшения шага интегрирования. В частности, уменьшение шага интегрирования вдвое увеличивает время вычислений в 2 раза, а увеличение порядка метода при сохранении величины шага увеличивает время вычислений лишь в 1,2 раза.
Для сопоставления эффективности алгоритмов метода Эверхарта различных порядков проведено совместное численное интегрирование уравнений движения больших планет, Луны, Солнца и астероида Aten 2003 НТ42 (17) на интервале времени с 2006 по 2204 гг. Выбор интервала интегрирования обусловлен необходимостью сопоставления результатов с данными каталога [3]. Выбор самого астероида связан с тем, что на данном интервале времени он не проходит через сферу действия больших планет, что даёт возможность проводить численное интегрирование с постоянным шагом. Для астероидов, проходящих через сферу действия больших планет, следует использовать переменный шаг интегрирования.
В табл. 8 приведены элементы орбиты астероида Aten 2003 НТ42. Здесь Т — момент оскуляции, М — средняя аномалия, a — большая полуось, е — эксцентриситет, со — аргумент перигелия, Q — долгота восходящего узла, г — наклонение, AS — абсолютная величина разности между полученными значениями элементов орбиты и данными каталога [3].
Элементы орбиты в первой строке табл. 8 получены путём совместного интегрирования уравнений движения больших планет, Солнца, Луны и астероида с переменным шагом [3]. В последующих строках табл. 8 приведены элементы орбиты астероида Aten 2003 НТ42, полученные модифицированным методом различных порядков при интегрировании с постоянным шагом.
Видно (см. табл. 8), что для получения верных результатов для 19-го порядка шаг интегрирования не должен превышать 4 дней, для 27-го порядка —
7 дней и для 31-го порядка —6 дней. Отсюда следует, что наиболее эффективно использовать модифицированный алгоритм метода Эверхарта 27-го порядка, поскольку дальнейшее повышение порядка не приводит к увеличению шага интегрирования.
Таблица 8
Элементы орбит астероида Aten 2003 НТ42 на эпоху 2204 12 13
Порядок; шаг М a е ш О i
Каталог [3] 261,7901 0,815063 0,257662 355,3049 33,2163 5,0752
15-ый; 1 день Дй' 261,79 0,0001 0,815063 0 0,257662 0 355,305 0,0001 33,2163 0 5,0752 0
19-ый; 4 дня Д5 261,79 0,0001 0,815063 0 0,257662 0 355,3049 0,0001 33,2163 0 5,0752 0
27-ой; 7 дней Д5 261,79 0,0001 0,815063 0 0,257662 0 355,3049 0,0001 33,2163 0 5,0752 0
31-ый; 6 дней Дй' 261,79 0,0001 0,815063 0 0,257662 0 355,305 0 33,2163 0 5,0752 0
Следует отметить, что верные результаты можно получить, решая данную задачу классическим методом, однако в этом случае для методов 19-го и 27-го порядков шаг интегрирования не должен превышать двух дней.
Таким образом, модификация алгоритма метода Эверхарта позволила более чем в три раза повысить эффективность метода для порядков выше 15-го по сравнению с классическим алгоритмом. Этот факт наиболее важен при исследовании эволюции орбит малых тел Солнечной системы на интервалах времени порядка нескольких сотен и тысяч лет.
Работа выполнена в рамках государственного задания высшим учебным заведениям в части проведения научно-исследовательских работ (проект 2.534.2011).
БИБЛИОГРАФИЧЕСКИЙ СПИСОК
1. Everhart Е. Implicit Single-Sequence Methods for Integrating Orbits// Cel. Mech., 1974. Vol. 10, no. 1. Pp. 35-55.
2. Заусаев А. Ф., Заусаев А. А. Математическое моделирование орбитальной эволюции малых тел Солнечной системы. М.: Машиностроение, 2008. 250 с. [Zausaev A.F., Zausaev A. A. Mathematical modelling of orbital evolution of small bodies of the Solar system. Moscow: Mashinostroenie, 2008. 250 pp.]
3. Заусаев А. Ф., Абрамов В. В., Денисов С. С. Каталог орбитальной эволюции астероидов, сближающихся с Землей с 1800 по 2204 гг.. М.: Машиностроение-1, 2007. 608 с. [Zausaev A. F., Abramov V. V., Denisov S. S. Catalogue of orbital evolution of asteroids approaching to the Earth between 1800 and 2204. Moscow: Mashinostroenie-1, 2007. 608 pp.]
4. Modern numerical methods for ordinary differential equations / eds. G. Hall; J. M. Watt. Oxford: Clarendon Press, 1976. 336 pp.
5. Брумберг В. А. Релятивистская небесная механика. М.: Наука, 1972. 382 с.
[Brumberg V.A. Relativistic celestial mechanics. Moscow: Nauka, 1972. 382 pp.]
Поступила в редакцию 02/V/2012; в окончательном варианте — 21/V/2012.
MSC: 85-08; 70М20, 65L99
RESEARCH OF EFFICIENCY OF ALGORITHMS OF METHOD EVERHART WITH HIGH ORDER OF APPROXIMATING FORMULAS
A. A. Zausaev
Samara State Technical University,
244, Molodogvardeyskaya St., Samara, 443100, Russia.
E-mail: zausaevaa@mail. ru
The modified algorithm for the numerical integration of the equations of celestial motion by the Everhart’s method is developed,. The study of the effectiveness of the algorithm for approximating high-order formulas is carried. High efficiency of the method is shown on the example of joint integration of the equations of motion of major planets, the Moon, the Sun and, the small bodies of Solar system.
Key words: Everhart’s method, numerical integration, orbit, small bodies of Solar system.
Original article submitted 02/V/2012; revision submitted 21/V/2012.
Artem A. Zausaev (Ph.D. (Phys. & Math.)), Associate Professor, Dept, of Applied Mathematics & Computer Science.
