Справочное руководство по небесной механике и астродинамике -
Скачать (прямая ссылка):
небесных тел получили очень большое распространение в связи с появлением
электронных вычислительных машин (ЭВМ). С помощью этих методов можно
получить таблицы численных значений координат небесных тел (или
оскулирующих элементов их орбит) на различные моменты времени.
В 1951 г. в США были опубликованы построенные таким путем таблицы
координат пяти внешних планет (Юпитера, Сатурна, Урана, Нептуна) на
период с 1653 по 2060 г. [17]. В 1950 и 1962 гг. в США были опубликованы
таблицы координат малых планет Цереры, Паллады, Юноны, Весты [18].
В Институте теоретической астрономии АН СССР (ИТА) регулярно публикуются
эфемериды малых планет, получаемые при помощи численного интегрирования.
Очень широко применяется численное интегрирование при изучении движения
комет и особенно в астродинамике для решения различных задач, относящихся
к движению искусственных небесных тел. При этом важное практическое
значение имеют так называемые краевые задачи.
Чаще всего рассматривают уравнения движения в прямоугольных координатах
ввиду простоты их правых частей. Тогда обычно имеют дело с системой
уравнений вида
'dF~==Fs(*i' Хл' ^ (s = 1....п), (7.3.01)
где Fs(x 1, ..., xn,t) - известные функции переменных хи ..., хп и
времени t. Значительно реже рассматривают уравнения относительно
оскулирующих элементов орбит. В этом случае рассматривается система
уравнений первого порядка-
= Fs (хи xn,t) (s = 1....п), (7.3.02)
однако правые части оказываются значительно сложнее.
668
Ч. VII. ЧИСЛЕННЫЕ МЕТОДЫ
[§ 3.01
Большинство методов численного интегрирования обыкновенных
дифференциальных уравнений рассчитано на системы вида
(7.3.02). При их применении к системам вида (7.3.01) последние приводят к
виду (7.3.02) заменой переменных
Имеются также методы, разработанные именно для систем вида (7.3.01).
§ 3.01. Метод Рунге- Кутта
Метод Рунге - Кутта непосредственно рассчитан на интегрирование систем
вида (7.3.02); он наиболее удобен при применении ЭВМ. Для того чтобы
начать вычисления, достаточно знать лишь начальные значения *s0 = *8(М
искомых функций.
Приведем сначала формулы для случая одного уравнения
Пусть дано начальное значение х0 = *(f0) и требуется вычислить значения
x(t) при fo = t0 + kh (k = 1, 2, 3, . ..). Число h называется шагом
интегрирования.
Если обозначить
то формулы, определяющие Х\ по заданному *0 с точностью до членов порядка
Л4, следующие:
Формулы для Адсо, имеющие погрешность порядка Л5 (эти формулы наиболее
употребительны при вычислениях на ЭВМ), следующие:
По формулам, аналогичным приведенным, вычисляются последовательно х2,
x3t... и т. д.
(7.3.03)
X (fo) = Xk, hF [х, t) = f (x, t),
(7.3.04)
A*0= -g- (&i + 2k2 + 2k3 + kt) ki=f{x0, t0),
(7.3.05)
§ 3.01] ГЛ. 3. ЧИСЛЕННОЕ РЕШЕНИЕ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЯ 669
Для контроля точности применяют вычисления с половинным шагом. А именно,
вычисляют Xi по приведенным формулам с шагом, равным А, а затем с шагом,
равным А/2. Разность между полученными двумя значениями Х\ принимают за
погрешность значения Х\, вычисленного с шагом А. Аналогичным образом
контролируют далее точность значения х2 и т. д.
Имеются программы вычислений на ЭВМ с автоматическим выбором шага при
заданной точности. В этом случае задается погрешность б (например, 1-10-
6) и некоторый первоначальный шаг интегрирования А0. ЭВМ вычисляет по
такой программе Х\ с шагом А0 и с шагом hi = А0/2 и сопоставляет
полученные значения (xi)/io, (xi)hi- Если разность между ними не
превышает по абсолютной величине б, то ЭВМ переходит к вычислению х2 с
шагом Ао- Если эта разность больше б, то ЭВМ выполняет вычисления
значения Х\ с шагом h2 = Ai/2. Если | (*i)m -(Х2К2] < б, то ЭВМ переходит
к вычислению х2 с шагом А. Если это условие не выполнено, то ЭВМ
вычисляет с шагом А3 = А2/2 и т. д.
Для системы уравнений вида (7.3.02) вычисления производятся параллельно
для каждого уравнения по формулам, аналогичным (7.3.04) или (7.3.05). .
Пусть дана система двух уравнений
Формулы, определяющие Х\, уi, например, с точностью до членов порядка А\
следующие:
^ = F(x, у, t), -g- = Ф (х, у, t) (7.3.06)
и начальные условия Хо = x(i0), уо = у (to) - Обозначим tk = to + kh,
hF(x, у, f) = f(x, у, t), АФ {х, у, t) = g(x, у, i), x(tk) = xk, у (tk) -
Уь (k = 1, 2, 3, ...).
*i = *о + А*о. У1 = Уо + &Уо,
Л*о=4(А. + ЗА3), At/0 = |(/1 + 3/3),
ki=f(x0, Уо, to), ll = g{x0, у0, to),
(7.3.07)
670
Ч. VII. ЧИСЛЕННЫЕ МЕТОДЫ
[§ 3.02
По таким же формулам вычисляются далее *2, У2, *з, Уз, и т. д.
Системы уравнений (7.3.01) расписываются при применения указанных выше
формул в виде систем уравнений первого порядка:
Чг=у*' .....*"• *) .....")¦ (7-3-08)
Ввиду простоты правых частей первой группы уравнений (для х'б) и
независимости функций Fa от уи уп формулы упрощаются.
При вычислениях по методу Рунге-Кутта значений искомых функций при каждом
последующем значении аргумента требуется вычислять несколько значений
правых частей уравнений в некоторых промежуточных точках. Поэтому объем
