5.3. Методы Рунге-Кутта

Для уменьшения погрешности метода интегрирования ОДУ, исполь­зующего разложение искомого решения вряд Тейлора (5.6), необходи­мо учитывать большее количество членов ряда. Однако при этом воз­никает необходимость аппроксимации производных от правых частей ОДУ. Оснавная идея методов Рунге-Кутта заключается в том, что про­изводные аппроксимируются через значения функции f(x, у) в точках на интервале [x₀, x_0+h], которые выбираются из условия наибольшей близости алгоритма к ряду Тейлора. В зависимости от старшей сте­пени h, которой учитываются члены ряда, построены вычислительные схемы Рунге-Кутта разных порядков точности.

Так, например, для второго порядка получено однопараметрическое семействосхем вида [1,2]:

, (5.12)

где 0 <£ 1 — свободный параметр,

.

Локальная погрешность схемы (5.12) имеет третий порядок, глобальная — второй, т.е. решение ОДУ, полученное по этой схеме, рав­номерно сходится к точному решению с погрешностью O(h²).

Для параметра £ чаще всего используют значение £ = 0,5 и £ = 1. В первом случае формула (5.12) приобретает вид

, (5.13)

геометрическая интерпретация которой представлена на рис.5.2,а [1,2]. Вначале вычисляется приближенное решение ОДУ в точке х₀+h по формуле Эйлераy_э=y₀+hf₀ . Затем определяется наклон интеграль­ной кривой в найденной точке f(x₀+h,y_э) и после нахождения сред­него наклона на шаге h находится уточненное значение y_PK = y(x₀+h). Схемы подобного типа называют "прогноз — коррекция", что подра­зумевает грубое вычисление решения по формуле низкого порядка, а затем уточнение с учетом полученной информации о поведении инте­гральной кривой.

Рис. 5.2. Метод Рунге-Кутта второго порядка: а-£ = О,5; б-£ = 1

С целью экономии памяти при программировании алгоритма (5.13), обобщенного для решения системы ОДУ, изменим его запись с учетомтого, что :

, (5.14)

где k— номер решения для системы ОДУ. Теперь не придется держать в памяти ЭВМ массив начальных значений , его можно не сохра­нять после вычисления значений эйлеровских приближений , xотя по сравнению с методом Эйлера схема (5.14) требует дополнительного массива для запоминания значений . Во втором случае при £ = 1 от формулы (5.12) переходим к схеме

(5.15)

геометрический смысл которой отражает рис.5.2,б. Здесь при прогнозе определяется методом Эйлера решение в точке

,

а после вычисления наклона касательной к интегральной кривой в сред­ней точке решение корректируется по этому наклону.

Формула (5.15) обобщается для решения системы ОДУ аналогично схеме (5.14). По сравнению с программой метода Эйлера для сохране­ния начальных значений придется ввести дополнительный массив.

Схему (5.14) можно получить из метода Эйлера с помощью первой и второй формул Рунге без использования общего соотношения (5.12) для методов второго порядка.

Для построения вычислительных схем методов Рунге-Кутта дру­гих порядков в тейлоровском разложении искомого решения у(х) учитываются члены, содержащие степени шага hдо m включительно, где m — порядок схемы. После аппроксимации производных правой ча­сти ОДУ f(x, у)можно получить семейство схем Рунге-Кутта различ­ных порядков [1, 20, 21], наиболее используемые из которых третьего и четвертого порядка приведены ниже. Они обобщаются для решения систем ОДУ, записанных в форме Коши.

Одна из схем метода Рунге-Кутта третьего порядка имеет вид

, (5.16)

где

Также используют другую схему

(5.17)

где

Схемы (5.16) и (5.17) на каждом шаге вычислений требуют нахо­ждения правой части ОДУ в трех точках. Локальная погрешность схем имеет четвертый порядок, глобальной — третий.

Схема метода Рунге-Кутта четвертого порядка имеет вид

, (5.18)

где

Известна также схема метода Рунге-Кутта четвертого порядка:

, (5.19)

где

.

Наиболее часто в вычислительной практике используют схему ме­тода Рунге-Кутта четвертого порядка, поскольку она имеет рад пре­имуществ при программной реализации

, (5.20)

где

Схемы (5.18) — (5.20) на каждом шаге требуют вычисления правой части ОДУ в четырех точках. Локальная погрешность этих схем имеет пятый порядок, глобальная — четвертый. Для удобства программной реализации, особенно в случае системы ОДУ, формулы (5.20) можно преобразовать к виду

(5.21)

где

i = 1,2,...,n— номер уравнения в системе ОДУ из п уравнений.

Вычислительная схема метода Рунге-Кутта пятого порядка имеет вид

, (5.22)

где

Схема (5.22) на каждом шаге требует вычисления правой части ОДУ в пяти точках. Локальная погрешность этих схем имеет шестой поря­док, глобальная — пятый.