- Метод Рунге
-
Ме́тоды Ру́нге — Ку́тты (распространено неправильное название Ме́тоды Ру́нге — Ку́тта или даже Ме́тоды Ру́нге — Кутта́) — важное семейство численных алгоритмов решения обыкновенных дифференциальных уравнений и их систем. Данные итеративные методы явного и неявного приближённого вычисления были разработаны около 1900 года немецкими математиками К. Рунге и М. В. Куттой.
Формально, методом Рунге — Кутты является модифицированный и исправленный метод Эйлера, они представляют собой схемы второго порядка точности. Существуют стандартные схемы третьего порядка, не получившие широкого распространения. Наиболее часто используется и реализована в различных математических пакетах (Maple, MathCAD, Maxima) стандартная схема четвёртого порядка. Иногда при выполнении расчётов с повышенной точностью применяются схемы пятого и шестого порядков[1][2]. Построение схем более высокого порядка сопряжено с большими вычислительными трудностями[3]. Методы седьмого порядка должны иметь по меньшей мере девять стадий, в схему восьмого порядка входит 11 стадий. Хотя схемы девятого порядка не имеют большой практической значимости, неизвестно, сколько стадий необходимо для достижения этого порядка точности. Аналогичная задача существует для схем десятого и более высоких порядков[3].
Содержание
Классический метод Рунге — Кутты четвёртого порядка
Метод Рунге — Кутты четвёртого порядка столь широко распространён, что его часто называют просто методом Рунге — Кутты.
Рассмотрим задачу Коши
Тогда приближенное значение в последующих точках вычисляется по итерационной формуле:
Вычисление нового значения проходит в четыре стадии:
где — величина шага сетки по .
Этот метод имеет четвёртый порядок точности, то есть суммарная ошибка на конечном интервале интегрирования имеет порядок (ошибка на каждом шаге порядка ).
Прямые методы Рунге — Кутты
Семейство прямых методов Рунге — Кутты является обобщением метода Рунге — Кутты четвёртого порядка. Оно задаётся формулами
где — величина шага сетки по и вычисление нового значения проходит в этапов:
Конкретный метод определяется числом и коэффициентами и . Эти коэффициенты часто упорядочивают в таблицу (называемую таблицей Бутчера)
Для коэффициентов метода Рунге — Кутты должны быть выполнены условия для . Если требуется, чтобы метод имел порядок , то следует также обеспечить условие
где — приближение, полученное по методу Рунге — Кутты. После многократного дифференцирования это условие преобразуется в систему полиномиальных уравнений относительно коэффициентов метода.
Произношение
Согласно грамматическим нормам русского языка, фамилия Ку́тта склоняется, поэтому говорят: «Метод Ру́нге — Ку́тты». Правила русской грамматики предписывают склонять все мужские и женские фамилии, оканчивающиеся на -а, -я, которым предшествует согласный. Единственное исключение — фамилии французского происхождения с ударением на последнем слоге типа Дюма́, Золя́. Однако, иногда встречается несклоняемый вариант «Метод Ру́нге — Ку́тта» (например, в книге [4]).
Решение систем ОДУ
Метод Ру́нге — Ку́тты непосредственно обобщается на случай систем обыкновенных дифференциальных уравнений путём записи системы и метода в векторной форме.
Пример решения на алгоритмических языках программирования
производя замену и перенося в правую часть, получаем систему:
код на Java для решения системы дифференциальных уравнений методом Рунге-Кутты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); } }
Код на языке C#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) { Init(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):
MATLAB, runge.mfunction Dy = runge(x, y) Dy = y(:); Dy(1) = y(2); Dy(2) = cos(3*x) - 4*y(1);
Имя файла и имя функции должно совпадать, но оно может быть любым неиспользуемым ранее.
Затем необходимо создать главный файл c именем, например, main.m, который будет выполнять основные вычисления. Этот главный файл будет содержать следующий текст:
MATLAB, main.mclear; clc; % Очистка памяти и экрана x = 0.1; % Значение аргумента 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 даст искомое решение:
MATLAB, решениеy1 = y(find(x == 0.1))
Ответ: y1 = 0.98768
См. также
Ссылки
- ↑ Бахвалов Н. С., Жидков Н. П., Кобельков Г. М. Численные методы. — М.: Бином, 2001 — с. 363—375.
- ↑ Ильина В. А., Силаев П. К. Численные методы для физиков-теоретиков. т. 2. — Москва-Ижевск: Институт компьютерных исследований, 2004. — с. 16-30.
- ↑ 1 2 J. C. Butcher. Numerical Methods for Ordinary Differential Equations. The University of Auckland, New Zealand.
- ↑ Б. П. Демидович, И. А. Марон, Э. З. Шувалова. Численные методы анализа, 3-е изд. — М.: Наука, 1967.
Категории:- Численное интегрирование
- Численные методы
Wikimedia Foundation. 2010.