CC BY
45
УДК 521.1
А.Ф. Заусаев, А.Г. Ольхин
ОЦЕНКА ЛОКАЛЬНЫХ ОШИБОК ДИСКРЕТИЗАЦИИ В ЧИСЛЕННОМ ИНТЕГРИРОВАНИИ УРАВНЕНИЙ ДВИЖЕНИЯ БОЛЬШИХ ПЛАНЕТ (МЕРКУРИЙ-ПЛУТОН) МЕТОДОМ ЭВЕРХАРТА
Проведено исследование локальной погрешности в численном интегрировании уравнений больших планет (Меркурий-Плутон) методом Эверхарта. Показано, что для внешних планет (Юпитер, Сатурн, Уран, Нептун) координаты и скорость планет могут быть вычислены численным методом Эверхарта с высокой степенью точности на интервале времени порядка 10-20 тысяч лет.
Для повышения эффективности численного интегрирования уравнений движения небесных объектов: астероидов, комет, метеоритных роев, естественных и искусственных спутников и т.д., создаются банки данных координат больших планет на определенном интервале времени. В зависимости от применяемого метода, начальных данных, интервала интегрирования, точность координат больших планет может быть различной. Целью данной работы являлось нахождение интервала интегрирования, на котором ошибка в координатах больших планет не превосходит точности телескопических наблюдений.
Для решения этой задачи на интервале времени с 1950 по 965 г.г. проведено численное интегрирование уравнений движения больших планет (Меркурий-Плутон) методом Эверхарта [1] с переменным шагом и порядком. На основании анализа полученных результатов оценены локальные ошибки дискретизации на исследуемом интервале времени для каждой планеты в отдельности.
Метод Эверхарта является одним из самых эффективных по точности и быстродействию численным методом для решения задач небесной механики [2].
Рассмотрим основную идею построения алгоритмов Эверхарта на примере решения дифференциальных уравнений вида
^2 х ч /1Ч
—т = Р (х, г). (1)
ёг2
Представим правую часть (1) в виде ряда по степеням X в окрестности 11=0.
х = Р (х, г) = Р1 + А1 г + А2 г2 +...+Апгп (2)
Интегрируя уравнения (2) дважды, получим
х = Х1 + Р^ + А^2/2 + А2 г 3/з +... + Апгп+1 /п +1, (3)
х = х1 + х^ + 2/2 + А^ 3/б +... + Апгп+2 /(п + 1)(п + 2). (4)
Важно отметить, что полином, стоящий в правой части (2) не является отрезком ряда Тейлора. Коэффициенты А! вычисляются из условий наилучшего приближения х и х в момент времени
Т, соответствующий значению решения на конце шага Ь по формулам (3) и (4). Для связи А-
значений с Б-значениями используется выражение
Р = Р1 + а^ + а 2 г(г — г 2) + а 3 г(г — г 2)(г — /3) +... . (5)
В каждый момент времени 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 = ((((Р5 — Р1 ) / г5 — «1 ) ! ^52 — « 2 ) / ^53 — « 3 ) / г54 , (6)
Связь А-значений с а-значениями находится из следующих соотношений:
А1 = а1 + (—г 2 ^)а 2 + (г 2 ^)а 3 = Сца1 + с 21а1а 2 + с31а 3 +... ,
Aз —а з + ••• — Cззa з + •••
Коэффициенты с^ определяются из рекуррентных соотношений:
Сц — 1 , г — ],
Сц —-^С,-1,1 г > 1> (8)
Сц — С,-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)
Выражения (9) являются предсказывающими уравнениями для определения коэффициентов ЛР С помощью следующих двух исправляющих уравнений находятся значения решения на конце шага Ь:
х(Т) — Х1 + Х1Т + Р1Т2 /2 + Л1Т3 /2 * з + Л2Т4/з * 4 + Л3Т5 /4 * 5,
Х(Т) — Х1 + Р1Т + Л1Т2/2 + Л2Т 3/з + ЛзТ 4/4^ (10)
Значения 12 , 1з , 14 определяются из условий, которые позволили бы с помощью формул (10) получить решения с точностью до седьмого порядка, тх^ разности решений пятого и седьмого порядков должны обратиться в нуль:
Ах — 0 , Ах — 0 (11)
Вводя безразмерные величины к, — ^, можно получить систему алгебраических уравнений
для соответствующих им коэффициентов вида
с41/2 + с42/з + с 4з/4 +1/5 — 0, с41 /з + с42 /4 + с4з /5 + 1/6 — 0,
С41 /4 + с^ /5 + с4з / 6 + 1/7 — 0^ (12)
Решение этой системы
с441 — -4/з5 — -к2Нзк4 ; с^2 — 6/7 — Н2Нз + Нзк4 + Н2к4 ; с^3 — -12/7 — -Н2 -Нз - к4 (1з)
показывает, что значения Ь2, Ьз, Ь4 являются корнями следующего полинома третьей степени:
кз + (-12/7)к2 + (6/7)к - 4/з5 — 0^ (14)
Из (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(1 + т) — + > к2 т, (-А--------------3-), (15)
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 - 2 - р)
дискретизации; р-порядок метода; 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с
