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

Предположим, что функция х от t затабулирована через равные промежутки времени (интервал равен ), так что задается формулой , где — эпоха, в которую х принимает значение

Тогда можно составить таблицу, подобную табл. 7.1.

Таблица 7.1

Первая разность получается вычитанием из вторая разность — вычитанием из и т. д.

Первая сумма находится по формуле

вторая сумма — по формуле

Полуразности в промежутках между нечетными разностями и столбцами сумм и между четными разностями и столбцами сумм вводятся по формулам

Эти полуразности отличаются от предыдущих разностей буквой .

Если t лежит внутри диапазона таблицы, то с ее помощью можно проводить интерполяцию (т. е. получать значения х для любого значения независимой переменной t, даже если это значение t и не отвечает целому ).

Для этой цели существуют различные формулы, использующие затабулированные величины, например формула Бесселя

где В — интерполяционные коэффициенты Бесселя. Они являются функциями и приведены во многих работах.

Существует также формула Эверетта

Здесь Е — коэффициенты Эверетта (функции ), которые также можно найти во многих работах, например в [16].

Разности более высоких порядков в табл. вида 7.1 связаны с производными соответствующих порядков х по формулами, при помощи которых можно выполнять численное дифференцирование. Так, формула численного дифференцирования Бесселя имеет вид

(значения В затабулированы).

Во многих задачах нужны значения производных только в узловых и половинных точках. Если это так, то уравнение (7.78) принимает вид

Рассмотрим теперь численное интегрирование дифференциального уравнения, которое нельзя проинтегрировать аналитически. Пусть уравнение имеет вид

где F — некоторая функция

Предположим, что мы подставили в уравнение ряд

и получили несколько первых постоянных выразив их через начальное условие при . Тогда при помощи этого ряда

для малого диапазона значений мы сможем вычислить значения х и F. Можно составить таблицу значений х и F (типа табл. 7.1), в которой принадлежит этому диапазону и принимает значения положительное или отрицательное целое число, шаг таблицы).

Обычно вычисляют и заносят в таблицу функцию, равную не F, a F, умноженную на Н, т. е.

Тогда имеем

и можно показать, что

В эпоху, соответствующую следующему шагу таблицы, имеем

По мере надобности разности можно оценить, основываясь на известной процедуре их формирования в таблице.

Кроме того, существуют экстраполяционные формулы, например

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

представляет собой одно из уравнений системы, то мы можем написать

Процедуру решения, как и в случае уравнения первого порядка, можно начать с того, что при помощи решения в виде ряда

системы уравнений типа (7.83), справедливого на малом интервале времени, формируется таблица сумм и разностей. На практике при составлении таблицы для заданного интервала времени, с которого начинают решение, обычно используют невозмущенную кеплеровскую орбиту.

В случае когда уравнения имеют второй порядок

где — значения при следующем узле таблицы имеем

Как и раньше, можно оценить разности, причем на практике оценку можно провести настолько точно, что после определения х, при помощи которого (совместно с у и ) из уравнения (7.84) вычисляется значение X, часто оказывается, что дальнейшие итерации на этом шаге уже не нужны.

Уравнение (7.88) можно использовать для получения экстраполированного значения х, если оценить значения величины X и линейных разностей. Можно воспользоваться также другим соотношением:

    (7.89)

Иногда при приближении тела (например, кометы) к перигелию возникает необходимость деления шага таблицы пополам. После прохождения перигелия шаг снова может быть удвоен.

Чтобы проиллюстрировать некоторые из обсуждавшихся выше методов, рассмотрим численное интегрирование уравнения второго порядка

с начальными условиями при

Подставляя в (7.90) ряд

и вычисляя коэффициенты при степенях получаем

где

Положим шаг таблицы равным 0,1. Ряд (7.91) можно использовать для вычисления значений х во втором столбце табл. 7.2.

Таблица 7.2

Значения функции X можно записать в третьем столбце таблицы.

Теперь мы можем составить таблицу для функции X, включая в нее, когда это возможно, разности а также полуразности . При помощи уравнений (7.85) и (7.86) вычисляются и заносятся в таблицу значения и

Последовательные значения первой и второй сумм получаются при помощи (7.74) и (7.75). В результате на данном этапе мы получим новую таблицу (см. табл. 7.3 выше ступенчатой линии).

Чтобы получить значение х в момент мы можем оценить разности, входящие в (7.88). Считая, что равно нулю, выпишем значения Они равны соответственно . Если к тому же предположить, что равно нулю, то можно также выписать значения для . Это соответственно 0,0000048 и 0,0000362. Значение нам известно, так что при помощи (7.81) можно вычислить первое приближение которое оказывается равным 0,35145. Исходя из этого значения, можно записать новые значения разностей. Оказалось, что эти новые разности отличаются своих прежних значений только в пятом знаке после запятой, так что новое значение не надо вычислять из (7.88). Новые разности вписаны в табл. 7.3 ниже ступенчатой линии. Теперь можно записать первую и вторую суммы и и начать новый шаг (вычисление ).

Другой способ состоит в том, что для получения первого приближения используется уравнение (7.89) в виде

При этом нужно оценить только

Если требуется найти к, то для того, чтобы можно было воспользоваться (7.87), надо заносить в таблицу еще и последовательные полуразности.

Численным методам посвящена обширная литература; многие математики, такие, как Ньютон, Гаусс, Лагранж, Бессель, Стирлинг и др., разработали изящные методы интерполяции, численного дифференцирования и интегрирования, решения дифференциальных уравнений, аппроксимации данных и т. д.

Ранее уже упоминалось, что одним из лучших методов численного интегрирования уравнений второго порядка, чаще всего встречающихся в задачах орбитального движения, является метод Гаусса—Джексона. Если придерживаться введенных выше обозначений, то для двукратного интегрирования используется формула

Таблица 7.3 (см. скан)

а для однократного интегрирования

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