Оценка локальных ошибок дискретизации в численном интегрировании движения больших планет (Меркурий-Плутон) методом Эверхарта

УДК 521.1

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

ОЦЕНКА ЛОКАЛЬНЫХ ОШИБОК ДИСКРЕТИЗАЦИИ В ЧИСЛЕННОМ ИНТЕГРИРОВАНИИ УРАВНЕНИЙ ДВИЖЕНИЯ БОЛЬШИХ ПЛАНЕТ (МЕРКУРИЙ-ПЛУТОН) МЕТОДОМ ЭВЕРХАРТА

Проведено исследование локальной погрешности в численном интегрировании уравнений больших планет (Меркурий-Плутон) методом Эверхарта. Показано, что для внешних планет (Юпитер, Сатурн, Уран, Нептун) координаты и скорость планет могут быть вычислены численным методом Эверхарта с высокой степенью точности на интервале времени порядка 10-20 тысяч лет.

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

Для решения этой задачи на интервале времени с 1950 по 965 г.г. проведено численное интегрирование уравнений движения больших планет (Меркурий-Плутон) методом Эверхарта [1] с переменным шагом и порядком. На основании анализа полученных результатов оценены локальные ошибки дискретизации на исследуемом интервале времени для каждой планеты в отдельности.

Метод Эверхарта является одним из самых эффективных по точности и быстродействию численным методом для решения задач небесной механики [2].

Рассмотрим основную идею построения алгоритмов Эверхарта на примере решения дифференциальных уравнений вида

^2 х ч /1Ч

—т = Р

ёг2

Представим правую часть (1) в виде ряда по степеням X в окрестности 11=0.

х = Р

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

х = Х1 + Р^ + А^2/2 + А2 г 3/з +... + Апгп+1 /п +1,

х = х1 + х^ + 2/2 + А^ 3/б +... + Апгп+2 /

Важно отметить, что полином, стоящий в правой части (2) не является отрезком ряда Тейлора. Коэффициенты А! вычисляются из условий наилучшего приближения х и х в момент времени

Т, соответствующий значению решения на конце шага Ь по формулам (3) и (4). Для связи А-

значений с Б-значениями используется выражение

Р = Р1 + а^ + а 2 г

В каждый момент времени X! имеем Р2 = Р1 + а1г 2 ,

Р3 = Р1 + а1^3 + а 2 г3 (г3 — ^2 ) ,

Принимая ХщИп-Х , получим

«1 = (Р2 — Р1)/г 2, а2 = ((Р3 — Р1) / г3 — «1) / t32,

а 3 — (((Р4 Р1) / г 4 «1) / г 42 « 2) / ^32 ,

«4 =

Связь А-значений с а-значениями находится из следующих соотношений:

А1 = а1 + (—г 2 ^)а 2 + (г 2 ^)а 3 = Сца1 + с 21а1а 2 + с31а 3 +... ,

Aз —а з + ••• — Cззa з + •••

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

Сц — 1 , г — ],

Сц —-^С,-1,1 г > 1>

Сц — С,-1,Ц-1 - t1C1 -1,; , 1 < ц < /•

В дальнейшем нахождение решения уравнения (1) сводится к нахождению узлов разбиения 1! шага Ь Покажем нахождение узлов разбиения шага Ь=[0,Т] на примере алгоритма интегрирования пятого порядка (п=5) В начальный момент времени 1^=0 нам известны х1 , Х1 , Р1 • Значения X в момент времени 12, 13 и 14 вычисляются по формулам

Х2 — Х1 + JC1t2 + 2 /2 + [2 /2• з + A2t2 /з• 4 + Aзt2 /4*5],

Хз — Х1 + Х1 tз + Fjtз /2 + Л11з /2 * з + [Л2tз /з * 4 + Лзtз /4 * 5],

Х4 — Х1 + Х^4 + 4 / 2 + Л114 /2 * з + Л214 / з * 4 + [Лзt4 /4 * 5] •

Выражения (9) являются предсказывающими уравнениями для определения коэффициентов ЛР С помощью следующих двух исправляющих уравнений находятся значения решения на конце шага Ь:

х(Т) — Х1 + Х1Т + Р1Т2 /2 + Л1Т3 /2 * з + Л2Т4/з * 4 + Л3Т5 /4 * 5,

Х

Значения 12 , 1з , 14 определяются из условий, которые позволили бы с помощью формул (10) получить решения с точностью до седьмого порядка, тх^ разности решений пятого и седьмого порядков должны обратиться в нуль:

Ах — 0 , Ах — 0

Вводя безразмерные величины к, — ^, можно получить систему алгебраических уравнений

для соответствующих им коэффициентов вида

с41/2 + с42/з + с 4з/4 +1/5 — 0, с41 /з + с42 /4 + с4з /5 + 1/6 — 0,

С41 /4 + с^ /5 + с4з / 6 + 1/7 — 0^

Решение этой системы

с441 — -4/з5 — -к2Нзк4 ; с^2 — 6/7 — Н2Нз + Нзк4 + Н2к4 ; с^3 — -12/7 — -Н2 -Нз - к4

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

кз +

Из (14) следует, что полученные узлы разбиения шага Ь совпадают с квадратурной формулой

Гаусса-Радо, а величины Ь2, Ьз, Ь4 равны корням полинома Лежандра Р3(2Ь-1)=0^

Метод Эверхарта является неявным одношаговым методом^ Для интегрирования обыкновенных дифференциальных уравнений самим Эверхартом разработана универсальная программа КЛБЛ19[1], написанная на Фортране и адаптированная нами для С++^

Для оценки локальной погрешности численного метода Эверхарта нами решались дифференциальные уравнения движения, которые в гелиоцентрической системе координат имеют следующий вид:

—-к «о+„,) 4+^ к 2 щ (^ _ х.),

dt2 гз ; Аз г3

£ — -к’-(1 + т) У + > к’- щ (^ - 4),

dt2 г3 т а:; г,1

d 2 ,2/л \ 2 V" ; 2 ' 2 - 2 2

' + т) — + > к т

— -к 2

dt2 Г ; а: “

где т, х, у, ъ - масса и гелиоцентрические координаты возмущаемой планеты; ш1, х1, у1, ъ - массы и гелиоцентрические координаты возмущающих планет; г, А1, г1 -расстояния, вычисляемые по формулам

2 2 I 2,2*2/ \ 2  / \ 2  / \2 2 2  2.2

Г — Х + у + 2 , А, — (Х - Х,) + (у - уг) + (г - ) , г — Х, + уг + ;

к- постоянная Гаусса; масса солнца ш0 принята за единицу •

Начальные значения координат и компонент скорости в эпоху 1949дек 30.0ET=JED 2433280•5 (ЕТ-эфемеридное время, JED-юлианская дата в эфемеридном времени) приведены в таблице 1 [3], где х, у, ъ выражены в астрономических единицах, Х, у, г - в астрономических единицах в эфемеридные сутки •

Т а б л и ц а 1

Экваториальные прямоугольные координаты на эпоху 1949дек 30.0 ЕТ^ЕБ2433280.5

Меркурий Сатурн

х= 0.3439262826 х'= -0.0084663252 х= -8.9724944563 х'= -0.0018582755

У= 0.0456120867 У'= 0.0256147777 У= 2.2797122950 У'= -0.0049838514

7= -0.0109252378 ^= 0.0145867628 7= 1.3303683212 7'= -0.0019802491

Венера Уран

х= 0.1429705295 х'= -0.0198937931 х= -1.0029132696 х'= -0.0039552572

У= 0.6470040699 У'= 0.0031132345 У= 17.3234919405 У'= -0.0003759007

7= 0.2824815251 7= 0.0026594965 7= 7.6048360565 7'= -0.0001088437

Земля+Луна Нептун

х= -0.1363613301 х'= -0.0173200155 х= -29.1942677427 х'= 0.0008207898

У= 0.8933977569 У'= -0.0022442677 У= -7.7192336013 У'= -0.0027720922

7= 0.3874584316 7= -0.0009733419 7= -2.4272492809 7'= -0.0011561171

Марс Плутон

х= -1.3698312567 х'= -0.0073845893 х= -26.2317561603 х'= -0.0013157218

У= 0.8431371618 У'= -0.0094773446 У= 20.5615297781 У'= -0.0026198366

7= 0.4238360949 7= -0.0041516538 7= 14.4436902948 7'= -0.0004270513

Юпитер

х= 3.3493566446 х'= 0.0055856550

У= -3.4737659939 У'= 0.0049622487

7= -1.5721612031 7= 0.0019922661

Данным значениям координат и скоростей планет соответствуют значения их обратных масс, представленные в табл^ 2^

Оценка локальной погрешности дискретизации численного интегрирования больших планет (Меркурий -Плутон) методом Эверхарта проводилась нами двумя методами: экстраполяцией и вложенным методом [4]

Экстраполяция является наиболее распространенным методом оценки погрешности Интегрирование на исследуемом интервале времени (1-10) выполнялось дважды: один раз с шагом Ь, другой раз с шагом Ы2. Тогда

dn —-----------------------, где ап - локальная ошибка

п к

дискретизации; р-порядок метода; 1 -начальный момент интегрирования 1949^дек 30^0 ED=2443280•5JED; 10- конечный момент интегрирования - 2073280.5JED.

При использовании вложенного метода решение уравнений (16) на интервале времени (1-10) также проводилось дважды: методом порядка р и методом порядка р+1- Разность полученных значений дает оценку величины

В табл^ 3 и 4 ниже приведены координаты и скорости планет на конечный момент интегрирования, а также локальные погрешности в их значениях, вычисленные методом экстраполяции и вложенным методом •

Т а б л и ц а 2 Обратные массы планет

Планеты шо/ш

Меркурий 6000000

Венера 408000

Земля+Луна 329390

Марс 3093500

Юпитер 1047.355

Сатурн 3501.6

Уран 22869

Нептун 19314

Плутон 360000

Т а б л и ц а 3

Экваториальные прямоугольные координаты на момент 2073280.5№Б

Меркурий Сатурн

х= -0.067450226 х'= 0.022247753 х= 9.46735285 х'= 0.000367267

У= -0.40951509 У'= -0.001326453 У= -0.93364407 у'= 0.005148022

7= -0.211509717 7'= -0.003056474 7= -0.780356583 7'= 0.002104081

Венера Уран

х= 0.618793222 х'= -0.010552375 х= -16.82062061 х'= 0.001604489

У= 0.356552709 У'= 0.015464962 У= -7.108241065 У'= -0.00342724

7= 0.119980039 7'= 0.007591038 7= -2.875325429 7'= -0.00152506

Земля+Луна Нептун

х= -0.509051005 х'= 0.014583928 х= -27.29215056 х'= 0.001348825

У= -0.802897292 У'= -0.007965873 У= -12.42017229 У'= -0.00258905

7= -0.350386958 7'= -0.003473202 7= -4.402056132 7'= -0.00109439

Марс Плутон

х= -1.175550092 х'= -0.009235795 х= -29.09700146 х'= -0.00082459

У= 1.049657621 У'= -0.00809684 У= 13.24623843 У'= -0.00291887

7= 0.515432311 7'= -0.003443543 7= 12.99169925 7'= -0.00066926

Юпитер

х= 0.650552679 х'= 0.007407106

У= -4.752673184 У'= 0.00126012

7= -2.055191307 7'= 0.000359188

Т а б л и ц а 4

Локальная погрешность дискретизации на конечный момент времени 2073280.5

15 порядок, вложенный метод 15 порядок, метод экстраполяции 19 порядок, метод экстраполяции

координат скоростей координат скоростей Координат скоростей

Меркурий 0.0001459 1.29526Е-06 0.000145369 1.29049Е-06 3.7766Е-07 3.335Е-09

8.7496Е-06 7.86318Е-06 8.71776Е-06 7.8343Е-06 2.2751Е-08 2.032Е-08

2.0072Е-05 4.06121Е-06 1.99982Е-05 4.0463Е-06 5.2017Е-08 1.05Е-08

Венера 5.4286Е-08 2.48524Е-09 1.0384Е-07 4.75003Е-09 6.6425Е-09 3.031Е-10

8.0078Е-08 1.43699Е-09 1.52734Е-07 2.74232Е-09 9.7136Е-09 1.744Е-10

3.9302Е-08 4.83984Е-10 7.49635Е-08 9.23276Е-10 4.7689Е-09 5.871Е-11

Земля 3.9694Е-08 3.94973Е-10 4.86998Е-08 4.84438Е-10 6.0537Е-09 6.009Е-11

2.1649Е-08 6.22589Е-10 2.6576Е-08 7.63852Е-10 3.3077Е-09 9.48Е-11

9.4439Е-09 2.71768Е-10 1.15928Е-08 3.33423Е-10 1.4424Е-09 4.136Е-11

Марс 6.9797Е-09 5.72752Е-11 8.14672Е-09 6.69166Е-11 8.3321Е-10 6.873Е-12

6.0941Е-09 5.26367Е-11 7.1154Е-09 6.12184Е-11 7.3018Е-10 6.143Е-12

2.587Е-09 2.58867Е-11 3.02196Е-09 3.01006Е-11 3.1135Е-10 3.013Е-12

Юпитер 1.2705Е-10 2.3847Е-13 2.58129Е-10 2.57489Е-13 3.7685Е-11 8.4Е-15

2.697Е-11 1.10692Е-12 4.79501Е-11 9.2435Е-13 6.1995Е-12 4.699Е-14

1.067Е-11 5.89134Е-13 1.641Е-11 5.10072Е-13 1.64Е-12 2.004Е-14

Сатурн 2.0691Е-11 2.29083Е-13 2.71907Е-11 1.60636Е-13 8.5798Е-12 8.428Е-14

3.1539Е-11 1.30758Е-12 6.99529Е-11 1.30644Е-12 1.3688Е-10 1.264Е-14

1.5862Е-11 6.7514Е-13 2.6075Е-11 6.7686Е-13 5.641Е-11 9.18Е-15

Уран 1.2601Е-11 2.0579Е-13 2.05986Е-11 2.0346Е-13 1.2601Е-11 4.94Е-15

1.913Е-11 1.30668Е-12 4.91962Е-12 1.3003Е-12 2.855Е-11 9.298Е-16

6.1204Е-12 6.7582Е-13 1.39888Е-13 6.7297Е-13 1.26Е-11 6.199Е-16

Нептун 2.3302Е-11 2.1074Е-13 2.16005Е-11 2.0444Е-13 2.64Е-11 7.02Е-15

2.9008Е-12 1.30858Е-12 5.30065Е-12 1.30057Е-12 5.39Е-11 6.27Е-15

1.7595Е-12 6.7635Е-13 1.39444Е-13 6.7289Е-13 2.286Е-11 2.78Е-15

Плутон 1.6698Е-11 2.1295Е-13 3.47988Е-11 2.02977Е-13 1.4399Е-11 2.174Е-15

2.95Е-11 1.30768Е-12 2.57003Е-11 1.30164Е-12 4.6199Е-11 5.01Е-15

1.08Е-11 6.75554Е-13 1.69997Е-12 6.73791Е-13 1.1399Е-11 3.129Е-15

Результаты вычислений показывают, что локальные погрешности, полученные двумя различными методами, отличаются друг от друга незначительно^ Однако порядок точности координат и скоростей планет существенно зависит от порядка метода Эверхарта^ Как видно из таблицы, интегрирование уравнений движения методом Эверхарта 15 порядка дает погрешность в координатах и скоростях Меркурия на два порядка больше по сравнению с методом 19 поряд-ка^ Следовательно, метод Эверхарта 19 порядка предпочтительнее при получении координат планет с высокой степенью точности

Как видно из табл^ 3 и 4 наибольшая локальная погрешность получается в координатах Меркурия • Полная погрешность дискретизации часто оценивается с помощью локальной погрешности Однако для внутренних планет - Меркурия и Венеры - такое допущение не может быть оправданно, так как релятивистские эффекты, которые не учитывались в данном исследовании оказывают более существенное влияние на смещение координат, чем погрешность метода-

Хорошо известно, что основной вклад в эволюционные процессы малых тел солнечной системы вносят внешние планеты: Юпитер, Сатурн, Уран, Нептун Как показывают расчеты, значения максимальной локальной погрешности в координатах этих планет не превосходят 0•5' 10-10 а^ Для внешних планет локальная погрешность сопоставима с полной погрешностью дискретизации

Учитывая характер развития погрешности, можно утверждать, что метод Эверхарта позволяет создать банк данных координат больших планет с достаточно большой степенью точности на интервале времени порядка 10-20 тысяч лет Банк данных можно эффективно использовать при исследовании эволюции орбит малых тел Солнечной системы: астероидов, комет и метеоритных роев^

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

1. EverhartE. Implicit single methods for integrating orbits // Oelestial mechanics. 1974. V.10. Р. 35-55.

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

3. Справочное руководство по небесной механике и астродинамике // Под ред. Г.Н. Дубошина. М.: Наука, 1976. 862с.

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