CC BY
40
К ВЫБОРУ ЧИСЛЕННЫХ МЕТОДОВ ИНТЕГРИРОВАНИЯ УРАВНЕНИЙ ДВИЖЕНИЯ НАВИГАЦИОННЫХ СПУТНИКОВ
Артем Андреевич Карауш
Новосибирский государственный технический университет, 630092, г. Новосибирск, пр. К. Маркса, 20, аспирант кафедры автоматики, e-mail: [email protected]
Александр Сергеевич Толстиков
Сибирская государственная геодезическая академия, 630108, г. Новосибирск, ул. Плахотного, 10, доцент кафедры метрологии, стандартизации и сертификации, ФГУП «СНИИМ», 630004, г. Новосибирск, пр -т Димитрова 4, начальник службы времени и частоты, e-mail: [email protected]
В статье проведен сравнительный анализ численных методов. Приведены обоснования для использования метода Эверхарта в качестве наиболее эффективного метода численного интегрирования орбит спутников.
Ключевые слова: численное интегрирование, метод Эверхарта, восстановление орбит
ИСЗ.
TO THE CHOICE OF NUMERICAL METHODS OF INTEGRATION OF THE EQUATIONS OF MOVEMENT OF NAVIGATING COMPANIONS
Artem A. Karaush
Novosibirsk State Technical University, 630092, Novosibirsk, K. Marksa, 20, post-graduate student, automation department, e-mail: [email protected]
Aleksandr S. Tolstikov
The Siberian state geodetic academy, 630108, Novosibirsk, Street Plahotnogo, 10, the senior lecturer of chair of metrology, standardization and certification, FGUP «SNIIM»", 630004, Novosibirsk, Dimitrova 4, head of National time and frequency service, e-mail: [email protected]
Comparative analysis of digital computational methods is presented in the article. It is shown, that Evehart method is the most effective method of digital integration of satellite orbits.
Key words: digital integration, Everhart method, recovering satellite orbits.
Центральной задачей эфемеридно-временного обеспечения спутниковых навигационных систем (ГЛОНАСС, GPS, GALILEO) является прогнозирование движения навигационных спутников (НС). Точность такого прогнозирования зависит от степени адекватности математических моделей, применяемых для описания движения НС, от правильности выбора стартовой точки прогноза и, в значительной степени, от характеристик выбранного численного метода интегрирования уравнений движения НС [1].
Целью настоящей работы является сравнительный анализ численных методов, применяемых для интегрирования уравнений движения НС по
околокруговым орбитам, и обоснованный выбор на этой основе наиболее эффективных численных схем.
Для указанного сравнительного анализа численных методов выбраны математические модели движения в классе дифференциальных уравнений
Ц 5
= —з77Ги(0 + Хд(0> u(0 = uo>u(0 = V С1)
Р (О м
uT (t) = [ x(t), y(t), z(t)] - вектор текущих координат НС в
квазиинерциальной системе координат (ИСК), p(t) = yjx2(t) + y\t) + z2(t) -текущий радиус орбиты НС, л - геоцентрическая постоянная гравитационного поля Земли, 8г. (t) - вектор действующих на НС возмущений: i = 1 - от несферичности гравитационного поля Земли, i = 2,3 - от гравитационного воздействия на НС Луны и Солнца, i = 4 - от радиационного давления на НС солнечного излучения, i = 5 - немоделируемые возмущения случайной природы для которых в лучшем случае могут быть получены статистические характеристики.
Возмущения s.(i),z = 1,...,4образуют группу моделируемых возмущений, которые с той или иной точностью могут быть представлены математическими моделями и учтены при расчете орбитального движения НС.
Возмущения si (t) приводятся к центру масс НС в объектоцентрической системе координат и преобразуются в ИСК с помощью известных матричных преобразований [2].
Наряду с общими требованиями точности к численным схемам интегрирования предъявляются дополнительные требования, порожденные особенностями уравнения движения НС вида (1).
1. Решение уравнения (1), не содержащего демпфирующих членов, находится на границе устойчивости. Выбранная численная схема не должна нарушать это условие и не накапливать погрешностей, при интегрировании на больших интервалах времени.
2. Численная схема не должна терять точность интегрирования при скачкообразных возмущениях в правой части уравнения (1). Такое изменение
возмущения s4(t) возможно при вхождении НС в тень Земли, когда скачком
снимается радиационное давление на НС солнечного излучения.
Анализ точности численных схем принято строить на основе сравнения аналитического решения уравнения движения и решений этого уравнения, полученных с помощью сравниваемых численных схем. Однако для общего
случая, при условии действия на НС всех возмущений (t), уравнение (1) в квадратурах не интегрируется [2, 3]. По этой причине для сравнительного анализа используется упрощенное уравнение
'¿(О = -со1 ■ г(() + 5!(О + 54(0, г(/0) = г0,г(г0) = ¿0 , ®2 = (2)
решением которого является проекция на ось Ъ ИСК движения НС в плоскости (ХОЪ), для которого ^ (?) имитирует возмущение от второй зональной гармоники амплитудой а, (?) - возмущение от радиационного
давления на НС солнечного излучения в виде прямоугольного импульса 6(£ — Т) амплитудой Ь с моментом появления Т и моментом снятия Т.
Метод $ конечное $ среднее арифм $ max time feval
Рунге-Кутты 1,2000E-02 6,7600E-02 7,5570E-01 0,049 4000
Адамса-Штёрмера 2,6659E-06 -4,271Ш-08 2,6934E-06 0,053 6006
Эверхарта 2,7940E-09 1Д19Ш-09 3,1665E-08 0,506 22051
Адамса-Мультона-Коуэлла 4,4694E-06 8,8574E-07 8,3372E-06 0,204 1991
Аналитическое решение уравнения (2) представляется выражением [4].
Таблица 1. Шаг = 100 с
а-Бт22а>С) а-5УП(<гС)-57П(3юС) эсоб2{оС} a■cos(3cot')■cos(a>t')
г(С) =-----^ 2-------1-----------7—2-----------------^Г~2-----1-----------7—2----------
2о 6о 2о 6о
Ъ-в(С — Т) Ъ-вС — Т) -Б1п(Тгю) -Б1п(а>С) Ъ-вС — Т) соБ(Тгю) -соб(фС)
+
2 2 2 СО О со
Ъ -в(С — Т2) - БШ2(а>С) Ъ -в(С — Т2)-бш(Т2™') б1П(ю£) Ъ -в(С — Т2) ■ соб2(оС)
2 2 2 2 ОЭ 22 Ф
Ъ -вС — Т) -собТо) -собОС) . , , , ,
2--- -----2-----2 2 2---- --- + с - Б1П(о>С) + с - СОБ(о>С).
со
Параметры С'1 и °2 определяются из начальных условий 2^{) ) ~ 2{]’> _ .
Сравнение численного интегрирования уравнения (2) с его аналитическим решением на интервале времени 105 с. проводилось для схем:
— Рунге-Кутты 4-го порядка [4.5],
— Адамса-Штёрмера 5-го порядка [4],
— Эверхарта 15- го порядка [6],
— Адамса-Мультона-Коуэлла 16- го порядка [5,7].
Численные схемы 2, 3, 4 не требуют сведения уравнения (2) к системе уравнений 1-го порядка. Схемы 1, 3 одношаговые, Для схем 2, 4 разгонная сетка решений получена на основе аналитического решения уравнения (2).
При моделировании выбирались: радиус орбиты Р = 25-10 м-,
Л = 398600 -109 м3 / с2 а = и -10—6 Ь = и -10—7. Начальные условия
? 5
2(^о) ~ 2а--¿(*о) “ ¿овыбирались так, чтобы движение НС начиналось от экватора.
Таблица 2. Шаг = 1000 с
Метод 8 конечное 8 среднее арифм 8 max time feval
Рунге-Кутты 1,8273E+03 1,3073E+01 1,8273E+03 0,005 400
Адамса-Штёрмера 1,8959E+00 -1,4830E-01 2,5660E+00 0,005 606
Эверхарта 2,6077E-08 -5,1868E-10 4,8429E-08 0,050 2251
Адамса-Мультона- Коуэлла 6,2514E+02 7,1383E+01 7,1766E+02 0,043 191
Таблица 3. Со ступенчатым воздействием. Шаг = 100 с
Метод 8 конечное 8 среднее арифм 8 max time feval
Рунге-Кутты 3,1900E-02 7,0400E-02 8,0340E-01 0,054 4000
Адамса-Штёрмера 1,0100E-02 1,1000E-03 2,0500E-02 0,052 6006
Эверхарта 1,3970E-08 -7,2660E-10 5,2154E-08 0,492 22051
Адамса-Мультона- Коуэлла 8,4000E-03 1,2000E-03 2,0500E-02 0,276 1991
Результаты сравнительного анализа численных схем приведены в таблицах 1, 2, 3. Для каждого из сравниваемых методов оценивалась погрешность интегрирования в метрах в виде:
- Абсолютного значения погрешности на конец интегрирования 8 -конечное,
- Среднего арифметического на интервале интегрирования 8 -среднее арифм.,
- Максимального значения на интервале интегрирования 8 -max.
Основные результаты и выводы
1. Схема Рунге-Кутты показала свою работоспособность при малом шаге интегрирования 100 с. и неприемлемую точность при шаге 1 000 с. Как правило одношаговые схемы применяются в сочетании с многошаговыми схемами для заполнения с малым шагом разгонной сетки решений для запуска много -шаговых схем.
2. Многошаговые методы Адамса-Штёрмера и Адамса-Мультона-Коуэлла обеспечивают приемлемую точность интегрирования при шаге 100 с. и шаге 1 000 с. и требуют знания решения в узлах разгонной сетки. Количество узлов зависит от порядка применяемой схемы.
Это создает определенные неудобства, связанные с необходимостью сгущения сетки решений из за риска потери точности интегрирования на участках вхождения НС в теневые зоны Земли, где возникают скачкообразные изменения радиационного давления солнечного излучения на НС.
3. Наилучшие по точности результаты показала одношаговая схема Эверхарта в условиях интегрирования с шагом 100 с. и с шагом 1000 с. гладкой функции и при скачкообразных возмущениях в правой части уравнения (2).
Дополнительный положительный эффект применения схемы Эверхарта заключается в том, что погрешность интегрирования оказывается знакопеременной функцией времени и обеспечивает тенденцию не накапливаться на интервале интегрирования. Об этом свидетельствует малость среднего арифметического погрешности 8 по сравнению с выборочными значениями этой погрешности.
Определенным неудобством применения схемы Эверхарта, по сравнению с рассмотренными схемами, является большое число обращений к вычислению правой части уравнения (2), характеризующееся числом feval в таблицах 1, 2, 3.
БИБЛИОГРАФИЧЕСКИЙ СПИСОК
1. Сурнин Ю.В., Кужелев С.В. Модели движения ИСЗ и точность численного прогнозирования орбит // Геодезия и картография. - 1982. - № 10. - С. 8-13.
2. Абалакин В.К., Аксенов Е.П., Гребенников Е.А., Демин В.Г., Рябов Ю.А. Справочное руководство по небесной механике и астродинамике. - М., 1976.
3. Чеботарев, Г. А. (1965). Аналитические и численные методы небесной механики. -М., 1965
4. Камке Э. Справочник по обыкновенным дифференциальным уравнениям. - М.,
1961.
5. Бордовицына, Т. В. Современные численные методы в задачах небесной механики. - Томск, 1984.
6. Everhart, E. Dept. of Physics and Astronomy, University of Denver, Colo., U.S.A. Celestial Mechanics, 1974, 10, p. 35-55.
7. Бордовицына, Т. В., Авдюшев, В. А. (2007). Теория движения искусственных спутников Земли. - Томск, 2007.
© А.А. Карауш, А.С. Толстиков, 2012
