5.4. Метод Рунгe-Кутта-Мерсона

Мерсон предложил модификацию метода Рунге-Кутта четвертого порядка, позволяющую оценивать погрешность на каждом шаге и при­нимать решение об изменении шага. Схему Мерсона [22] с помощью эквивалентных преобразований приведем к виду, удобному для про­граммирования:

, (5.23)

где

Схема Мерсона требует на каждом шаге вычислять правую часть ОДУ в пяти точках, но за счет только одного дополнительного коэффи­циента k_iпо сравнению с классической схемой Рунге-Кутта на каждом шаге можно определить погрешность решения Rпо формуле

. (5.24)

Для автоматического выбора шага интегрирования рекомендуется следующий критерий. Если абсолютное значение величины R, вычис­ленное по формуле (5.24), больше допустимой погрешности ε, то шаг hуменьшается в два раза и вычисления по схеме (5.23) повторяются с точки . Шаг hможно удвоить при выполнении условия:

32|R| < . (5.25)

Автоматический выбор шага позволяет значительно сократить вре­мя решения ОДУ. Схема (5.23) обобщается для решения системы ОДУ аналогично классической схеме Рунге-Кутта.

5.5. Методы Адамса-Башфорта и Адамса-Маултона

При решении задачи Коши методами Рунге-Кутта необходимо вы­числять правые части ОДУ в нескольких точках на каждом шаге. Ко­личество точек зависит от порядка используемого метода. После того, как искомая интегральная кривая у (х) определена в нескольких точках , можно применить алгоритмы интерполяции и сократить количество правых частей ОДУ для получения решения в очередной точке x_n+1.Подобные методы называют многоточечными или много­шаговыми. Известно несколько типов таких методов [1,23, 24].

Алгоритмы многоточечных методов основываются на аппроксима­ции интерполяционными полиномами либо правых частей ОДУ, либо интегральных кривых .

Рассмотрим четырехточечный вариант одного из методов перво­го типа для задачи Коши, сформулированной в виде (5.5). С помо­щью любой из схем, рассмотренных в предыдущих разделах, вычис­лим решения заданного дифференциального уравнения в точ­ках . Правая часть уравнения f(x, у) на интегральной кривой, соответствующей начальному условию, будет функцией только одного аргумента х

f(x) = f(x,y(x)),

значения которой в рассматриваемых точках обозначим

В окрестности узлов функцию f(x) приближенно заменим интерполяционным полиномом Ньютона (2.6):

(5.26)

где — разделенные разности (2.10) — (2.12).

Представим искомое решение в точке в виде тейлоров­ского разложения около точки :

, (5.27)

где — производныепо х отправой части дифференциально­го уравнения в точке .

Дифференцируя полином (5.26), получим выражения для произ­водных

Последнее соотношение при в случае равноотстоящих узлов после подстановки в них разделенных разностей (2.10) — (2.12) при­нимает вид

(5.28)

Подставляя производные (5.28) в разложение (5.27), получим экс-траполяционную формулу Адамса-Башфорта [24]

(5.29)

имеющую пятый порядок локальной погрешности и четвертый — гло­бальной. Остаточный член формулы (5.29) равен [2]. Зна­чительная величина коэффициента в остаточном члене обусловлена тем,что точка лежит вне интервала расположения узлов , по зна­чениям функции f(x), в которых построен интерполяционный поли­ном. Таким образом, имеем дело с экстраполяцией, погрешность которой в соответствии с оценкой (2.15) всегда больше, чем при ин­терполяции.

Изменяя количество членов, учитываемых в ряде (5.27), можно по­лучить схемы Адамса-Башфорта различных порядков. Формула Адамса для переменного шага приведена в [2]. В табл. 5.1 приведены фор­мулы методов Адамса-Башфорта от первого до шестого порядка.

Таблица 5.1. Методы Адамса-Башфорта

Порядок	Алгоритмы	Локальная погрешность
1	2	3
Первый
Второй
Третий
Четвертый
Пятый	+ }
Шестой

Анализ табл. 5.1 показывает, что метод Адамса-Башфорта порядка k требует k начальных значений: Поэтому его называют k-шаговым. С целью уменьшения погрешности способом, аналогичным получению формулы (5.29), по узлам строится интерполяционная формула Адамса–Маултона:

(5.30)

Последняя формула является неявной, так как искомая величина необходима для вычисления значений функции , входя­щей в правую часть. Выражение(5.30) можно рассматривать как нели­нейное уравнение относительно неизвестной величины и решать его одним из методов, описанных выше. Наиболее часто здесь использу­ется метод простых итераций, хотя в некоторых случаях оказывается более предпочтительным метод Ньютона [24]. Следует иметь в виду, что каждая итерация потребует нового вычисления правой части диф­ференциального уравнения f(x, у). Решение, определенное по экстра-поляционной формуле (5.29), обычно выбирается в качестве начально­го приближения для итерационных методов, поэтому выражение (5.29) рассматривается как формула прогноза, тогда выражение (5.30) явля­ется формулой коррекции.

Таким образом, вычисления на каждом шаге интегрирования диф­ференциального уравнения осуществляются по схеме , где этапы вычислительного процесса обозначены буквами Р — прогноз, Е — вычисление функции f(x, у), С — коррекция, т— количество итераций коррекции [9]. В [2] приводится эмпирическое правило, со­гласно которому погрешность решения убывает только до тех пор, пока ,где p — порядок используемого неявного метода. Следователь­но, для метода четвертого порядка не следует выполнять более четы­рех итераций коррекции. С другой стороны, в [9] отмечается, что схе­ма является более устойчивой в смысле накопления вы­числительной погрешности по сравнению со схемой , следо­вательно, наиболее выгодной будет схема РЕСЕ. При реализации по­следней схемы на каждом шаге интегрирования осуществляется толь­ко одна коррекция.

Формулы (5.29) и (5.30) без изменения переносятся на системы ОДУ первого порядка, записанные в форме Коши.

В табл. 5.2 приведены формулы методов Адамса-Маултона от пер­вого до шестого порядка.

Таблица 5.2. Методы Адамса-Маултона

Порядок	Алгоритмы	Локальная погрешность
1	2	3
Первый
Второй
Третий
Четвертый
Пятый
Шестой