Для уменьшения погрешности метода интегрирования ОДУ, использующего разложение искомого решения вряд Тейлора (5.6), необходимо учитывать большее количество членов ряда. Однако при этом возникает необходимость аппроксимации производных от правых частей ОДУ. Оснавная идея методов Рунге-Кутта заключается в том, что производные аппроксимируются через значения функции f(x, у) в точках на интервале [x0, x0+h], которые выбираются из условия наибольшей близости алгоритма к ряду Тейлора. В зависимости от старшей степени h, которой учитываются члены ряда, построены вычислительные схемы Рунге-Кутта разных порядков точности.
Так, например, для второго порядка получено однопараметрическое семействосхем вида [1,2]:
, (5.12)
где 0 <£ 1 — свободный параметр,
.
Локальная погрешность схемы (5.12) имеет третий порядок, глобальная — второй, т.е. решение ОДУ, полученное по этой схеме, равномерно сходится к точному решению с погрешностью O(h2).
Для параметра £ чаще всего используют значение £ = 0,5 и £ = 1. В первом случае формула (5.12) приобретает вид
, (5.13)
геометрическая интерпретация которой представлена на рис.5.2,а [1,2]. Вначале вычисляется приближенное решение ОДУ в точке х0+h по формуле Эйлераyэ=y0+hf0 . Затем определяется наклон интегральной кривой в найденной точке f(x0+h,yэ) и после нахождения среднего наклона на шаге h находится уточненное значение yPK = y(x0+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) на каждом шаге требует вычисления правой части ОДУ в пяти точках. Локальная погрешность этих схем имеет шестой порядок, глобальная — пятый.