Лекция 15.
Метод Рунге-Кутты

Метод Рунге-Кутты используют для расчета стандартных моделей достаточно часто, так как при небольшом объеме вычислений он обладает точностью метода Ο⁴(h).

Для построения разностной схемы интегрирования воспользуемся разложением функции

в ряд Тейлора:

Заменим вторую производную в этом разложении выражением

где

Причем Δx подбирается из условия достижения наибольшей точности записанного выражения. Для дальнейших выкладок произведем замену величины «y с тильдой» разложением в ряд Тейлора:

Для исходного уравнения (1) построим вычислительную схему:

которую преобразуем к виду:

Введем следующие обозначения:

Эти обозначения позволяют записать предыдущее выражение в форме:

Очевидно, что все введенные коэффициенты зависят от величины Δx и могут быть определены через коэффициент α, который в этом случае играет роль параметра:

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

Та же схема в форме разностного аналога уравнения (1):

При α = 0 получаем как частный случай уже известную схему Эйлера:

При α = 1:

При α = 1 проведение расчетов на очередном шаге интегрирования можно рассматривать как последовательность нижеследующих операций.

1. Вычисляется выражение, представляющее собой полушаг интегрирования по схеме Эйлера, то есть определяется приближенное значение искомой функции в точке x_k + h/2:

   

2. Для той же промежуточной точки находится приближенное значение производной:

   

3. Определяется уточненное значение функции в конечной точке всего шага, причем по схеме Эйлера с вычисленным на предыдущем шаге значением производной:

   

Геометрические построения (см. рис. 15.1) показывают, что получаемое в такой последовательности решение лежит «ближе» к истинному, чем вычисляемое по схеме Эйлера, то есть следует ожидать более высокой точности решения, получаемого методом Рунге-Кутты. Ранее мы назвали эту схему «модифицированным методом Эйлера».

Рис. 15.1. Иллюстрация расчета на шаге методом Рунге-Кутты при значении параметра α = 1

Теперь рассмотрим схему при α = 0.5 (геометрическая интерпретация результата приведена на рис. 15.2).

1. Выполняется полный шаг метода Эйлера с целью определения приближенного значения искомой функции на конце отрезка интегрирования:

   

2. Для этой же точки вычисляется приближенное значение производной:

   

3. Находится среднее значение двух производных, определенных на концах отрезка:

   

4. Вычисляется значение искомой функции в конечной точке всего шага по схеме Эйлера с усредненным значением производной:

   

Рис. 15.2. Иллюстрация расчета на шаге методом Рунге-Кутты при значении параметра α = 0.5

Иногда получающееся выражение называют схемой (методом) Эйлера-Коши. Геометрически понятно, что получаемый указанным способом результат также должен быть «ближе» к истинному решению, чем получаемый по схеме Эйлера.

Пример. Решить уравнение dy/dx = –y, y(0) = 1 методом Рунге-Кутты.

Поскольку правая часть дифференциального уравнения имеет вид: f(x, y) = –y, схема метода при α = 0.5 представляется следующим образом:

Построим последовательность значений искомой функции:

…

Pезультаты получаемого численного решения для значения аргумента x = 10 при различных шагах интегрирования приведены в табл. 15.1. Три верные значащие цифры получены для шага h = 0.01.

Оценим погрешность аппроксимации уравнения (1) разностной схемой метода Рунге-Кутты. Подставляем точное решение в разностный аналог исходного дифференциального уравнения и вычисляем невязку:

Подставим разложения функций

в полученное выражение:

Учитывая уравнение (1), а также выражение для производной

окончательно получаем, что ψ_k = Ο(h²), то есть метод Рунге-Кутты, независимо от значения параметра α, имеет второй порядок аппроксимации.

Методы Рунге-Кутты третьего и четвертого порядков

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

И две схемы Рунге-Кутты, имеющие четвертый порядок аппроксимации:

Пример. Решить методом Рунге-Кутты четвертого порядка уравнение dy/dx = –y, y(0) = 1.

В соответствии с приведенными выше соотношениями определяем коэффициенты:

Построим последовательность значений искомой функции:

…

Результаты получаемого численного решения для значения аргумента x = 10 при различных шагах интегрирования приведены в табл. 15.2. Три верные значащие цифры получены для шага h = 0.25.

Сравнение таблиц 15.1 и 15.2 с решениями одной и той же задачи позволяет сделать вывод, что более высокая степень аппроксимации дифференциального уравнения разностным аналогом позволяет получать более точное решение при более крупном шаге и, следовательно, меньшем числе шагов, то есть приводит к снижению требуемых ресурсов ЭВМ.

На сегодняшний день для грубого расчета вычисления производятся методом Эйлера, для точного расчета — методом Рунге-Кутты.