Оценка точности метода Эверхарта при решении уравнений движения больших планет на интервале времени 10 000 лет Текст научной статьи по специальности «Математика»

Численные методы

УДК 521.1, 521.4

А. Ф. Заусаев, А. А. Заусаев, А. Г. Ольхин

ОЦЕНКА ТОЧНОСТИ МЕТОДА ЭВЕРХАРТА ПРИ РЕШЕНИИ УРАВНЕНИЙ ДВИЖЕНИЯ

БОЛЬШИХ ПЛАНЕТ НА ИНТЕРВАЛЕ ВРЕМЕНИ 10 000 ЛЕТ

Исследована возможность применения метода Эверхарта для решения уравнений движения планет на интервале времени 10 000 лет. Оценка точности метода была проведена на основе алгоритма численного интегрирования 23 порядка с шагом 3 дня.

Целью данной работы является решение двух задач: а) исследование точности метода Эверхарта при интегрировании уравнений движения больших планет (Меркурий-Плутон) на интервале времени порядка 10 000 лет, путем решения задачи п тел; б) разработка аппроксимирующих формул высокого порядка.

Создаваемые ранее базы данных координат больших планет охватывали период телескопических наблюдений, начиная с 1660 г. [1]. Для исследования эволюции орбит малых тел Солнечной системы этот интервал времени мал, так как период изменения угловых элементов орбит составляет десятки и сотни тысяч лет.

Применяемые явные многошаговые методы Коуэлла не позволяли увеличить интервал интегрирования более чем на 500 лет, из-за быстрых накоплений ошибок округления [2]. Разработанный Эверхартом метод численного интегрирования обыкновенных дифференциальных уравнений относится к числу неявных одношаговых методов [3].

Основным достоинством одношаговых методов является то обстоятельство, что для них разработаны надежные оценки локальной погрешности дискретизации. Существуют три основных метода для оценки полной погрешности: экстраполяция, вложенные методы, многошаговые оценки погрешности [4]. Экстраполяция и вложенные методы являются наиболее распространенными для оценок полной погрешности метода интегрирования. Экстраполяция основана на следующей идее: интегрирование производится дважды, с шагом И и И/2 . Тогда погрешность при шаге И / 2 приблизительно равна разности решений с шагом И и И / 2 , умноженной на 2 р [4].

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

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

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

1. Метод Эверхарта. Метод Эверхарта является одной из разновидностей методов Рунге-

Кутты. Он относится к неявным одношаговым методам, что обеспечивает его сходимость и ус-

тойчивость [3, 6]. Рассмотрим основную идею построения метода Эверхарта на примере решения уравнения вида

x = F(x,t). (1)

Представим правую часть (1) в виде временного ряда

x = F(x,t) = F1 + A1t + A2t2 +... + Antn. (2)

Интегрируя (2), получим выражения для определения координат и скоростей:

t2 t3 tn+2

x = Xj + X-, t + F--+ A-, —+... + An----------------, (3)

1 1 1 2 1 6 n (n + 1)(n + 2)

12 t3 tn+1

x = x1 + F1t + A1 — + A2 — + ... + An - . (4)

2 6 n +1

Полиномы (3) и (4) не являются рядами Тейлора, а коэффициенты Ai вычисляются из условия наилучшего приближения x и x с помощью конечных разложений (3) и (4). Для связи A-значений с F-значениями воспользуемся вспомогательным выражением

F = F1 + a1t + a2t(t -12) + a3t(t -12)(t -13) +... (5)

Уравнение (5) усечено по времени tn . В каждый фиксированный момент времени tt имеем

F = F + « t 1 2 J 1 ~ кл\12'>

F3 = F1 + a1t3 + a2t3(t3 -12), (6)

Принимая tnj = tn - tj, найдем ai через разделенные разности:

«1 = (F2 - F1)/t2, a 2 = ((F3 - F1)/t3 - «1)/132 ,

«3 = (((F4 - F1)/14 - «1)/142 - «2)/143 , (7)

«4 = ((((F5 - F1)/ 15 - «1)/ 152 - «2)/ 153 - «3)/ t54 ,

Приравнивая коэффициенты при одинаковых степенях t в уравнениях (2) и (5), выразим коэффициенты Д через а.:

Д = а1 + (—^)а2 + (t2 ^)а3 +... = с11а1 + с 21а2 + с31а3 +

Д2 = а 2 + (—12 — tз )«з +... = с22а 2 + съ2аъ +...,

A3 = «3 +... = С33«3 +...,

Коэффициенты с. определяются из следующих рекуррентных соотношений:

(8)

С,. = I i =j,

ij

сп =—tгCг-^^, 1 >1> (9)

с.. = с , . , — ^с. , . 1<7 <1.

V 1 —1,} —1 г г—1,7 ^

Для алгоритма интегрирования пятого порядка имеем:

с41 12t зt 4, с42 1213 + 13t 4’ с43 12 13 14’ (10)

где t2 ,t3, t4 являются корнями кубического уравнения

(—t2^зt 4 ) + (t2tз + tзt4 + t4t2 ) + (—t2 — tз — t4 )t + (1) = 0. (11)

Таким образом, нахождение решения уравнения (1) сводится к нахождению узлов разбие-

ния ti шага И .

2. Нахождение узлов разбиения шага интегрирования. Нахождение узлов разбиения шага И=[0,Т] рассмотрим на примере алгоритма интегрирования пятого порядка.

В начальный момент времени tl = 0 известны х; ,Х; . Значения х в моменты времени

t2,t3,t4 определяются с помощью трех предсказывающих уравнений:

Х2 = х, + Х, t2 + 12 /2 + [А,12 /6 + А2t2 /12 + ^ / 20]; (12)

Х3 = х, + Х^3 + 3 /2 + А]tз /6 + [Д2tз /12 + Aзtз / 20]; (13)

Х4 = х, + Х, t4 + 14/2 + .А-! 14/6 + Д214 /12 + [А314 / 20]; (14)

и двух исправляющих уравнений для нахождения положения и скорости на конце шага И:

х(Т) = Х1 + ХТ + ^Т2/2 + А1Т 3/6 + А2Т4/12 + А3Т5/20; (15)

Х(Т) = х, + ^Т + А1Т2 / 2 + А2Т3 / 3 + А3Т4 / 4. (16)

Эта схема является неявной, так как коэффициенты, стоящие в квадратных скобках (12-14), неизвестны при первой итерации.

Уравнения (12)-(16) обеспечивают пятый порядок точности относительно t. Можно увеличить порядок точности в вычислении Х и Х до седьмого порядка путем специального выбора подшагов 12, 13, t4. С этой целью увеличим количество разбиений интервала интегрирования, добавив два дополнительных времени t5,16. Затем вычислим для t5 и 16 значения а4 и а5, а также новые значения А4, А5 и А-, А2, А^.

Из уравнения (15) можно найти поправки Дх , улучшающие значения координат:

Ах = (а; — А1)Т3 /6 + (а; — А2)Т4 /12 + (А3' — А3)Т5 /20 + А'аТ6 /30 + А5Т7 /42. (17)

Выражая в уравнении (9) с51,...,с54 через с41, с42, с43, а также полагая

И2 = 12/ Т, И3 = t3 / Т и И4 = t4 / Т, (18)

формула (17) может быть записана в виде

Дх = (а4 — ^а5)Т6[с4, /6 + с42/12 + с43 /20 +1/30] + а5Т7[с4,/12 + с42 /30 +1/42].

Значение Дх в последнем выражении можно обратить в ноль при выполнении следующих условий:

с4, /6 + с42 /12 + с43 /20 +1/30 = 0 ,

41 42 43 (19)

с4, /12 + с42 /20 + с43 /30 +1/42 = 0 .

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

с4;/2 + с 42/3 + с43/4 +1/5 = 0,

< с4, /3 + с42 /4 + с43 / 5 +1/6 = 0, (20)

с4,/4 + с 42/5 + с43/6 +1/7 = 0.

Из решения этой системы

с4, = —4/35 = —И2И3И4; с42 = 6/7 = И2И3 + И3И4 + И2И4; с43 = —12/7 = —И2 — И3 — И4 (21)

следует, что значения величин И2 ,И3, И4 являются корнями следующего полинома третьей сте-

пени:

И3 + (—12 / 7 )И2 + (6 / 7 )И — 4 / 35 = 0. (22)

которые имеют следующие значения:

И2 = t2/ Т = 0.212340538239...; И3 = ^ / Т = 0.590533135559...;

И4 = ^/ Т = 0.91141240488... (23)

Использование этих узлов позволяет получить решение уравнения (1) с точностью до седьмого порядка для обеих компонент х и Х . Полученные по формуле (22) узлы разбиения И совпадают с узлами квадратурной формулы Гаусса-Радо. Область изменения И заключена в пределах 0 < И < 1.

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

На данный метод Э. Эверхартом разработаны алгоритмы и программа на языке Фортран. Программа Эверхарта адаптирована на ЭВМ БЭСМ-6 С.В. Тарасевием в Институте теоретической астрономии АН СССР.

В работе Э. Эверхарта [3] приведены узлы разбиения основного шага интегрирования на подшаги, обеспечивающие точность до 15 порядка включительно.

t

В табл. 1 приведены узлы разбиения отрезка [0, 7] на подшаги И = ^ для 19 и 23 поряд-

Известно, что положение небесного объекта в заданной системе координат определяется шестью элементами орбит: М, а,е, ш, 0,1, где М - средняя аномалия, а -большая полуось, е -

эксцентриситет, ш -аргумент перигелия, О -долгота восходящего узла, 1 -наклонение.

Положение планеты на орбите определяется двумя величинами: радиус-вектором г и истинной аномалией V . Средняя аномалия представляет собой дугу круга, которую описывала бы планета за время (V — 10), если бы она двигалась равномерно по окружности радиуса а со

средней угловой скоростью п . Связь средней аномалии с эксцентрической определяется уравнением Кеплера

Таким образом, зная M - среднюю аномалию, можно легко определить V - истинную аномалию.

Для оценки интервала интегрирования, на котором можно эффективно использовать метод Эверхарта для получения координат больших планет, нами проведено численное интегрирование уравнений движения больших планет [2] с шагом 0.5 дня и 1, 2, 3 дня с различным порядком точности, 19 и 23 соответственно.

Начальные данные масс, координат и скоростей планет (Меркурий-Плутон) на момент времени 2433280.5 JED взяты из Справочного руководства по небесной механике и астродинамике [7].

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

Численные результаты элементов орбит планет, полученные методом Эверхарта с различным шагом интегрирования и различным порядком, приведены в табл. 2, 3.

Данные эксперимента показывают, что на момент 2798530.5 JED результаты вычислений, полученные методом Эверхарта 19 и 23 порядками точности и шагом интегрирования 1 и 3 дня (соответственно) совпадают с точностью до шестого знака после запятой для всех элементов орбит, включая среднюю аномалию.

В табл. 3 приведены значения максимальной разности погрешности средней аномалии Меркурия на конце интервала интегрирования для различных вариантов метода Эверхарта.

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

Оценка полной погрешности метода Эверхарта экстраполяционным и вложенным методами позволяет сделать следующий вывод: максимальная погрешность в средних аномалиях планет не превосходит 0.1 угловых секунд дуги.

Полученный результат позволяет применять метод Эверхарта с порядком не менее 19 и с шагом интегрирования не более 1 дня для создания банка данных координат больших планет (Меркурий-Плутон) на интервале времени порядка 10 000 лет.

Основными вопросами при создании банка данных являются вопросы эффективности, точности и устойчивости используемых методов.

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

ков метода.

E - e sin E = n(t - T) = M .

Истинная аномалия выражается через эксцентрическую с помощью формулы

кого порядка, вычисленные с помощью разностной схемы, теряют всякий физический смысл. С другой стороны, уменьшение шага вдвое, как правило, удваивает время вычислений, в то время как повышение порядка метода увеличивает работу в 1.3 раза. Сопоставляя эффективность различных модификаций метода Эверхарта, следует отметить, что метод Эверхарта с использованием 23 порядка и шагом интегрирования 3 дня является более эффективным, чем метод Эверхарта 19 порядка с шагом 1 день.

Как показывают расчеты, наибольшая погрешность наблюдается в средней аномалии Меркурия. Это обстоятельство накладывает ограничение на интервал времени, на котором координаты Меркурия могут быть вычислены точно. Для остальных планет (Венера-Плутон) точные координаты с ошибкой в средней аномалии менее 1 секунды дуги могут быть вычислены методом Эверхарта на значительно большем интервале времени, чем 10 000 лет.

В связи с тем, что основной вклад в эволюционные процессы малых тел Солнечной системы вносят внешние планеты: Юпитер, Сатурн, Уран, Нептун, их координаты с помощью рассмотренных алгоритмов и программ можно вычислять с высокой степенью точности (ДМ < 1"’) на значительно больших интервалах времени.

Т а б л и ц а 1

Разбиение интервала (0,1) на подшаги Г аусса-Радо

Порядок К

19 0.036257S12SS320946094

0.11S07S97S7S999S70019

0.2371769S4S149603S531

0.3S1SS276530470597536

0.53S02959S91S9S906511

0.6903324200723621S294

0.S23SS3343S3700471S14

0.925612610290S0395536

0.9S55S759035112345137

23 0.025273б20397520349419925

0.0S304161344740514574191S

0.1691751003771S1424343219

0.27779б715109032072344951

0.401502720232S60S14519170

0.531S623S6910415955S04065

0.659991S420S5334S10022770

0.77715939295б1б2143241701

0.S753S0774S55556925520646

0.94796454SS72S19447093136

0.9S99S171953S319594093396

Т а б л и ц а 2

Элементы орбит планет (Меркурий-Плутон) 2950 г. T=2798530.5 JED (Метод Эверхарта 23-го порядка, шаг 1 день)

Планеты M a e W Ш i

Меркурий Венера Земля + Луна Марс Юпитер Сатурн Уран Нептун Плутон 315.бб0бб3 12б.530241 34б.539534 4б.953004 4S.576655 40.S3S023 25б.544033 16S.3769S0 314.320135 0.3S709S 0.72333б 1.000002 1.523734 5.2074S4 9.56S937 19.2991S2 29.990290 З9.б14951 0.205S33 0.00б357 0.01б27б 0.0941SS 0.050057 0.054025 0.04703б 0.0126S1 0.24S6S0 31.66S529 57.3б4557 293.225930 293.40б555 274.755700 347.71S653 93.079095 279.34S63S 113.123774 4б.47б440 73.434SS2 172.025772 46.1667S2 101.794149 110.б3073б 74.50бб45 131.213339 109.455772 6.94434s 3.3S23S6 0.129775 1.7бб314 1.290321 2.50S759 0.75SS59 1.776997 17.119466

Элементы орбит планет (Меркурий-Плутон) 11950 г. T=6085780.5 JED ____________(Метод Эверхарта 23-го порядка, шаг 1 день)________

Планеты M a e W Ш i

Меркурий Венера Земля + Луна Марс Юпитер Сатурн Уран Нептун Плутон 3.5Q3Q77 3Q8.947234 258.9бЗ782 35.348Q6Q 288.251259 153.31б517 28Q.651233 44.428б58 бЗ^З9б7 Q.387Q99 Q.72333Q 1.QQQQQ1 1.523бб8 5.2Q8511 9.58933Q 19.141б97 3Q.Q3Q115 З9.б2З8б7 Q.2Q7237 Q.QQ35Q4 Q.Q1147Q Q.1Q127Q Q.Q59721 Q.Q15922 Q.Q43851 Q.Q1Q463 Q.246783 57.Q51115 65.142485 345.277957 12.288292 28Q.864Q8Q 9Q.326583 1QQ.926398 267.792744 115.648776 34.19Q516 46.7Q9965 15Q.643196 5.864932 11б.9б2бЗЗ 87.26Q732 84.783361 13Q.5Q7739 1Q8.7Q5552 6.451584 3.QQ1259 1.18Q695 Q.8Q445Q 1.4448Q7 2.297616 Q.648967 1.799445 17.128889

Т а б л и ц а 3

Разности в средних аномалиях для Меркурия в различных вариантах метода Эверхарта.

T=6085780.5 JED

DM (23п2d и23п1d ) DM (19пQ.5d и 23п1d ) DM (19 п1d и 23 п1d ) DM (23п3d и23п1d )

Q."Q612 Q."Q36 Q."Q18 Q."Q648

БИБЛИОГРАФИЧЕСКИЙ СПИСОК.

1. Беляев Н.А. Эволюция орбиты кометы Даниэля 19Q9IY за 4QQ лет (166Q-2Q6Q гг.). Предварительное исследование // Л.: Бюлл.ИТА. 19бб. Т.Ю. № 1Q(157). С.б9б-7Ю.

2. Чеботарев Г.А. Аналитические и численные методы небесной механики. М.-Л.: Наука, 19б5. Зб8 с.

3. EverhartE. Implicit single methods for integrating orbits. // Celestial mechanics. 1974. №.1Q. Р.З5-55.

4. Холл Д., Уатт Д.. Современные численные методы решения обыкновенных дифференциальных уравнений. М.: Мир, 1979. 312 с.

5. Бордовицина Т.В. Современные численные методы в задачах небесной механики. М.: Наука, 1984. 13б с.

6. Борунов В.П., Иванов В.А., Миронов С.В. Сравнение эффективности методов Фелберга, Дормана-Принса и Эверхарта численного интегрирования систем обыкновенных дифференциальных уравнений. Применение систем Mathematica и Maple в научных исследованиях. М.: ВЦ РАН, 2QQ1, С.17-б2.

7. Абалакин В.К., Аксенов Е.П., Гребенников Е.А., Рябов Ю.А. Справочное руководство по небесной механике. М.: Наука, 197б. 8б2 с.

Поступила 29.08.2004 г.