Метод Рунге — Кутты (Bymk; Jruiy — Trmmd)
Ме́тоды Ру́нге — Ку́тты (в литературе встречается название методы Рунге — Кутта) — большой класс численных методов решения задачи Коши для обыкновенных дифференциальных уравнений и их систем. Первые методы данного класса были предложены около 1900 года немецкими математиками К. Рунге и М. В. Куттой.
К классу методов Рунге — Кутты относятся явный метод Эйлера и модифицированный метод Эйлера с пересчётом, которые представляют собой соответственно методы первого и второго порядка точности. Существуют стандартные явные методы третьего порядка точности, не получившие широкого распространения. Наиболее часто используется и реализован в различных математических пакетах (Maple, MathCAD, Maxima) классический метод Рунге — Кутты, имеющий четвёртый порядок точности. При выполнении расчётов с повышенной точностью всё чаще применяются методы пятого и шестого порядков точности[1][2]. Построение схем более высокого порядка сопряжено с большими вычислительными трудностями[3].
Методы седьмого порядка должны иметь по меньшей мере девять стадий, а методы восьмого порядка — не менее 11 стадий. Для методов девятого и более высоких порядков (не имеющих, впрочем, большой практической значимости) неизвестно, сколько стадий необходимо для достижения соответствующего порядка точности[3].
Классический метод Рунге — Кутты четвёртого порядка
[править | править код]Метод Рунге — Кутты четвёртого порядка при вычислениях с постоянным шагом интегрирования столь широко распространён, что его часто называют просто методом Рунге — Кутты.
Рассмотрим задачу Коши для системы обыкновенных дифференциальных уравнений первого порядка (далее , а ).
Тогда приближенное значение в последующих точках вычисляется по итерационной формуле:
Вычисление нового значения проходит в четыре стадии:
где — величина шага сетки по .
Этот метод имеет четвёртый порядок точности. Это значит, что ошибка на одном шаге имеет порядок , а суммарная ошибка на конечном интервале интегрирования имеет порядок .
Явные методы Рунге — Кутты
[править | править код]Семейство явных методов Рунге — Кутты является обобщением как явного метода Эйлера, так и классического метода Рунге — Кутты четвёртого порядка. Оно задаётся формулами
где — величина шага сетки по , а вычисление нового значения проходит в этапов:
Конкретный метод определяется числом и коэффициентами и . Эти коэффициенты часто упорядочивают в таблицу (называемую таблицей Бутчера):
Для коэффициентов метода Рунге — Кутты должны быть выполнены условия для . Если требуется, чтобы метод имел порядок , то следует также обеспечить условие
где — приближение, полученное по методу Рунге — Кутты. После многократного дифференцирования это условие преобразуется в систему полиномиальных уравнений относительно коэффициентов метода.
Неявные методы Рунге — Кутты
[править | править код]Все до сих пор упомянутые методы Рунге — Кутты являются явными методами[англ.]. К сожалению, явные методы Рунге — Кутты, как правило, непригодны для решения жёстких систем уравнений из-за малой области их абсолютной устойчивости[4]. Неустойчивость явных методов Рунге — Кутты создаёт весьма серьёзные проблемы и при численном решении дифференциальных уравнений в частных производных[англ.].
Неустойчивость явных методов Рунге — Кутты мотивировала развитие неявных методов. Неявный метод Рунге — Кутты имеет вид[5][6]:
где
Явный метод характерен тем, что матрица коэффициентов для него имеет нижний треугольный вид (включая и нулевую главную диагональ) — в отличие от неявного метода, где матрица имеет произвольный вид. Это также видно по таблице Батчера[англ.].
Следствием этого различия является необходимость на каждом шагу решать систему уравнений для , где — число стадий. Это увеличивает вычислительные затраты, однако при достаточно малом можно применить принцип сжимающих отображений и решать данную систему методом простой итерации[7]. В случае одной итерации это увеличивает вычислительные затраты всего лишь в два раза.
С другой стороны, Жан Кунцма́н[фр.] (1961) и Джон Батчер[англ.] (1964) показали, что при любом количестве стадий существует неявный метод Рунге — Кутты с порядком точности . Это значит, например, что для описанного выше явного четырёхстадийного метода четвёртого порядка существует неявный аналог с вдвое большим порядком точности.
Неявный метод Рунге — Кутты второго порядка
[править | править код]Простейшим неявным методом Рунге — Кутты является модифицированный метод Эйлера «с пересчётом». Он задаётся формулой:
.
Для его реализации на каждом шаге необходимы как минимум две итерации (и два вычисления функции).
Прогноз:
- .
Коррекция:
- .
Вторая формула — это простая итерация решения системы уравнений относительно , записанной в форме сжимающего отображения. Для повышения точности итерацию-коррекцию можно сделать несколько раз, подставляя . Модифицированный метод Эйлера «с пересчётом» имеет второй порядок точности.
Устойчивость
[править | править код]Преимуществом неявных методов Рунге — Кутты в сравнении с явными является их бо́льшая устойчивость, что особенно важно при решении жестких уравнений. Рассмотрим в качестве примера линейное уравнение y' = λy. Обычный метод Рунге — Кутты, применённый к этому уравнению, сведётся к итерации , с r, равным
где обозначает вектор-столбец из единиц[8]. Функция называется функцией устойчивости[9]. Из формулы видно, что является отношением двух полиномов степени , если метод имеет стадий. Явные методы имеют строго нижнюю треугольную матрицу откуда следует, что и что функция устойчивости является многочленом[10].
Численное решение данного примера сходится к нулю при условии с . Множество таких называется областью абсолютной устойчивости. В частности, метод называется A-устойчивым, если все с находятся в области абсолютной устойчивости. Функция устойчивости явного метода Рунге — Кутты является многочленом, поэтому явные методы Рунге — Кутты в принципе не могут быть A-устойчивыми[10].
Если метод имеет порядок , то функция устойчивости удовлетворяет условию при . Таким образом, представляет интерес отношение многочленов данной степени, приближающее экспоненциальную функцию наилучшим образом. Эти отношения известны как аппроксимации Паде. Аппроксимация Паде с числителем степени и знаменателем степени А-устойчива тогда и только тогда, когда [11]
-стадийный метод Гаусса — Лежандра имеет порядок , поэтому его функция устойчивости является приближением Паде . Отсюда следует, что метод является A-устойчивым[12]. Это показывает, что A-устойчивые методы Рунге — Кутты могут иметь сколь угодно высокий порядок. В отличие от этого, порядок А-устойчивости методов Адамса не может превышать два.
Произношение
[править | править код]Согласно грамматическим нормам русского языка, фамилия Ку́тта склоняется, поэтому говорят: «Метод Ру́нге — Ку́тты». Правила русской грамматики предписывают склонять все фамилии (в том числе и мужские), оканчивающиеся на -а, -я, которым предшествует согласный. Единственное исключение — фамилии французского происхождения с ударением на последнем слоге типа Дюма́, Золя́[13]. Однако иногда встречается несклоняемый вариант «Метод Ру́нге — Ку́тта» (например, в книге[14]).
Пример решения на алгоритмических языках программирования
[править | править код]производя замену и перенося в правую часть, получаем систему:
public class MainClass {
public static void main(String[] args) {
int k = 2;
double Xo, Yo, Y1, Zo, Z1;
double k1, k2, k4, k3, h;
double q1, q2, q4, q3;
/*
*Начальные условия
*/
Xo = 0;
Yo = 0.8;
Zo = 2;
h = 0.1; // шаг
System.out.println("\tX\t\tY\t\tZ");
for(; r(Xo,2)<1.0; Xo += h){
k1 = h * f(Xo, Yo, Zo);
q1 = h * g(Xo, Yo, Zo);
k2 = h * f(Xo + h/2.0, Yo + q1/2.0, Zo + k1/2.0);
q2 = h * g(Xo + h/2.0, Yo + q1/2.0, Zo + k1/2.0);
k3 = h * f(Xo + h/2.0, Yo + q2/2.0, Zo + k2/2.0);
q3 = h * g(Xo + h/2.0, Yo + q2/2.0, Zo + k2/2.0);
k4 = h * f(Xo + h, Yo + q3, Zo + k3);
q4 = h * g(Xo + h, Yo + q3, Zo + k3);
Z1 = Zo + (k1 + 2.0*k2 + 2.0*k3 + k4)/6.0;
Y1 = Yo + (q1 + 2.0*q2 + 2.0*q3 + q4)/6.0;
System.out.println("\t" + r(Xo + h, k) + "\t\t" + r(Y1 ,k) + "\t\t" + r(Z1 ,k));
Yo = Y1;
Zo = Z1;
}
}
/**
* функция для округления и отбрасывания "хвоста"
*/
public static double r(double value, int k){
return (double)Math.round((Math.pow(10, k)*value))/Math.pow(10, k);
}
/**
* функции, которые получаются из системы
*/
public static double f(double x, double y, double z){
return (Math.cos(3*x) - 4*y);
}
public static double g(double x, double y, double z){
return (z);
}
}
using System;
using System.Collections.Generic;
namespace PRJ_RungeKutta
{
/// <summary>
/// Реализация метода Ру́нге — Ку́тты для обыкновенного дифференциального уравнения
/// </summary>
public abstract class RungeKutta
{
/// <summary>
/// Текущее время
/// </summary>
public double t;
/// <summary>
/// Искомое решение Y[0] — само решение, Y[i] — i-я производная решения
/// </summary>
public double[] Y;
/// <summary>
/// Внутренние переменные
/// </summary>
double[] YY, Y1, Y2, Y3, Y4;
protected double[] FY;
/// <summary>
/// Конструктор
/// </summary>
/// <param name="N">размерность системы</param>
public RungeKutta(uint N)
{
Init(N);
}
/// <summary>
/// Конструктор
/// </summary>
public RungeKutta(){}
/// <summary>
/// Выделение памяти под рабочие массивы
/// </summary>
/// <param name="N">Размерность массивов</param>
protected void Init(uint N)
{
Y = new double[N];
YY = new double[N];
Y1 = new double[N];
Y2 = new double[N];
Y3 = new double[N];
Y4 = new double[N];
FY = new double[N];
}
/// <summary>
/// Установка начальных условий
/// </summary>
/// <param name="t0">Начальное время</param>
/// <param name="Y0">Начальное условие</param>
public void SetInit(double t0, double[] Y0)
{
t = t0;
if (Y == null)
Init((uint)Y0.Length);
for (int i = 0; i < Y.Length; i++)
Y[i] = Y0[i];
}
/// <summary>
/// Расчет правых частей системы
/// </summary>
/// <param name="t">текущее время</param>
/// <param name="Y">вектор решения</param>
/// <returns>правая часть</returns>
abstract public double[] F(double t, double[] Y);
/// <summary>
/// Следующий шаг метода Рунге-Кутта
/// </summary>
/// <param name="dt">текущий шаг по времени (может быть переменным)</param>
public void NextStep(double dt)
{
int i;
if (dt < 0) return;
// рассчитать Y1
Y1 = F(t, Y);
for (i = 0; i < Y.Length; i++)
YY[i] = Y[i] + Y1[i] * (dt / 2.0);
// рассчитать Y2
Y2 = F(t + dt / 2.0, YY);
for (i = 0; i < Y.Length; i++)
YY[i] = Y[i] + Y2[i] * (dt / 2.0);
// рассчитать Y3
Y3 = F(t + dt / 2.0, YY);
for (i = 0; i < Y.Length; i++)
YY[i] = Y[i] + Y3[i] * dt;
// рассчитать Y4
Y4 = F(t + dt, YY);
// рассчитать решение на новом шаге
for (i = 0; i < Y.Length; i++)
Y[i] = Y[i] + dt / 6.0 * (Y1[i] + 2.0 * Y2[i] + 2.0 * Y3[i] + Y4[i]);
// рассчитать текущее время
t = t + dt;
}
}
class TMyRK : RungeKutta
{
public TMyRK(uint N) : base(N) { }
/// <summary>
/// пример математический маятник
/// y''(t)+y(t)=0
/// </summary>
/// <param name="t">Время</param>
/// <param name="Y">Решение</param>
/// <returns>Правая часть</returns>
public override double[] F(double t, double[] Y)
{
FY[0] = Y[1];
FY[1] = -Y[0];
return FY;
}
/// <summary>
/// Пример использования
/// </summary>
static public void Test()
{
// Шаг по времени
double dt = 0.001;
// Объект метода
TMyRK task = new TMyRK(2);
// Определим начальные условия y(0)=0, y'(0)=1 задачи
double[] Y0 = { 0, 1 };
// Установим начальные условия задачи
task.SetInit(0, Y0);
// решаем до 15 секунд
while (task.t <= 15)
{
Console.WriteLine("Time = {0:F5}; Func = {1:F8}; d Func / d x = {2:F8}", task.t, task.Y[0], task.Y[1]); // вывести t, y, y'
// рассчитать на следующем шаге, шаг интегрирования
task.NextStep(dt);
}
Console.ReadLine();
}
}
class Program
{
static void Main(string[] args)
{
TMyRK.Test();
}
}
}
В программе на С# используется абстрактный класс RungeKutta, в котором следует переопределить абстрактный метод F, задающий правые части уравнений.
Пример решения в среде MATLAB
[править | править код]Решение систем дифференциальных уравнений методом Рунге — Кутты является одним из самых распространённых численных методов решений в технике. В среде MATLAB реализована его одна из разновидностей — метод Дормана — Принса[англ.]. Для решения системы уравнений необходимо сначала записать функцию, вычисляющую производные, то есть функции y = g(x, y, z) и z = cos(3x) − 4y = f(x, y, z), о чём сказано выше. В одном из каталогов, к которому имеется доступ из системы MATLAB, нужно создать текстовый файл с именем (например) runge.m со следующим содержимым (для MATLAB версии 5.3):
function Dy = runge(x, y)
Dy = y(:);
Dy(1) = y(2);
Dy(2) = cos(3*x) - 4*y(1);
Имя файла и имя функции должно совпадать, но оно может быть любым неиспользуемым ранее.
Затем необходимо создать главный файл c именем, например, main.m, который будет выполнять основные вычисления. Этот главный файл будет содержать следующий текст:
clear; clc; % Очистка памяти и экрана
h = 0.1; % Шаг интегрирования
x_fin = 8; % Конечное время интегрирования
y0 = 0.8; % Начальное значение функции
Dy0 = 2; % Начальное значение производной функции
[x, y] = ode45('runge', [0:h:x_fin], [y0 Dy0]); % Метод Рунге — Кутты
plot(x, y, 'LineWidth', 2); grid; % Построение графика и сетки
legend('y(x)', 'y''(x)', 0); % Легенда на графике
Так как MATLAB ориентирован на работу с матрицами, решение по методу Рунге — Кутты очень легко выполняется для целого ряда x, как, например, в приведённом примере программы. Здесь решение — график функции в пределах времён от 0 до x_fin.
Переменные x и y, полученные в результате работы функции ODE45, есть векторы значений. Очевидно, что решение конкретно заданного выше примера — второй элемент x, так как первое значение — 0, шаг интегрирования h = 0,1, а интересующее значение x = 0,1. Следующая запись в командном окне MATLAB даст искомое решение:
y1 = y(find(x == 0.1))
Ответ: y1 = 0,98768
Примечания
[править | править код]- ↑ Бахвалов Н. С., Жидков Н. П., Кобельков Г. М. . Численные методы. — М.: Лаборатория Базовых Знаний, 2001. — 630 с. — ISBN 5-93208-043-4. — С. 363—375.
- ↑ Ильина В. А., Силаев П. К. . Численные методы для физиков-теоретиков. Ч. 2. — Москва-Ижевск: Институт компьютерных исследований, 2004. — 118 с. — ISBN 5-93972-320-9. — С. 16—30.
- ↑ 1 2 Butcher J. C.[англ.] . Numerical Methods for Ordinary Differential Equations. — New York: John Wiley & Sons, 2008. — ISBN 978-0-470-72335-7.
- ↑ Süli & Mayers, 2003, p. 349—351.
- ↑ Iserles, 1996, p. 41.
- ↑ Süli & Mayers, 2003, p. 351—352.
- ↑ Неявные методы Рунге — Кутты Архивная копия от 6 марта 2019 на Wayback Machine.
- ↑ Hairer & Wanner, 1996, p. 40—41.
- ↑ Hairer & Wanner, 1996, p. 40.
- ↑ 1 2 Iserles, 1996, p. 60.
- ↑ Iserles, 1996, p. 62—63.
- ↑ Iserles, 1996, p. 63.
- ↑ Как склонять фамилии (трудные случаи) — «Грамота.ру» — справочно-информационный Интернет-портал «Русский язык» . Дата обращения: 5 июля 2016. Архивировано 23 мая 2018 года.
- ↑ Демидович Б. П., Марон И. А., Шувалова Э. З. . Численные методы анализа. 3-е изд. — М.: Наука, 1967.
Литература
[править | править код]- Hairer E., Wanner G. . Solving ordinary differential equations II: Stiff and differential-algebraic problems. 2nd ed. — Berlin, New York: Springer-Verlag, 1996. — ISBN 978-3-540-60452-5.
- Iserles A. . A First Course in the Numerical Analysis of Differential Equations. — Cambridge: Cambridge University Press, 1996. — ISBN 978-0-521-55655-2.
- Süli E., Mayers D. . An Introduction to Numerical Analysis. — Cambridge: Cambridge University Press, 2003. — ISBN 0-521-00794-1.