Метод Рунге — Кутты

Отпатрулированная версия этой страницы, проверенная 14 мая 2012, была основана на этой версии.

Ме́тоды Ру́нге — Ку́тты (распространено неправильное название Ме́тоды Ру́нге — Ку́тта или даже Ме́тоды Ру́нге — Кутта́) — важное семейство численных алгоритмов решения обыкновенных дифференциальных уравнений и их систем. Данные итеративные методы явного и неявного приближённого вычисления были разработаны около 1900 года немецкими математиками К. Рунге и М. В. Куттой.

Формально, методом Рунге — Кутты является модифицированный и исправленный метод Эйлера, они представляют собой схемы второго порядка точности. Существуют стандартные схемы третьего порядка, не получившие широкого распространения. Наиболее часто используется и реализована в различных математических пакетах (Maple, MathCAD, Maxima) стандартная схема четвёртого порядка. Иногда при выполнении расчётов с повышенной точностью применяются схемы пятого и шестого порядков^[1]^[2]. Построение схем более высокого порядка сопряжено с большими вычислительными трудностями^[3]. Методы седьмого порядка должны иметь по меньшей мере девять стадий, в схему восьмого порядка входит 11 стадий. Хотя схемы девятого порядка не имеют большой практической значимости, неизвестно, сколько стадий необходимо для достижения этого порядка точности. Аналогичная задача существует для схем десятого и более высоких порядков^[3].

Классический метод Рунге — Кутты 4 порядка

Метод Рунге — Кутты 4 порядка столь широко распространён, что его часто называют просто методом Рунге — Кутты.

Рассмотрим задачу Коши

{\textbf {y}}'={\textbf {f}}(x,{\textbf {y}}),{\textbf {y}}(x_{0})={\textbf {y}}_{0}.

Тогда приближенное значение в последующих точках вычисляется по итерационной формуле:

{\textbf {y}}_{n+1}={\textbf {y}}_{n}+{1 \over 6}({\textbf {k}}_{1}+2{\textbf {k}}_{2}+2{\textbf {k}}_{3}+{\textbf {k}}_{4})

Вычисление нового значения проходит в четыре стадии:

{\textbf {k}}_{1}=h\ {\textbf {f}}\left(x_{n},{\textbf {y}}_{n}\right),

{\textbf {k}}_{2}=h\ {\textbf {f}}\left(x_{n}+{h \over 2},{\textbf {y}}_{n}+{1 \over 2}{\textbf {k}}_{1}\right),

{\textbf {k}}_{3}=h\ {\textbf {f}}\left(x_{n}+{h \over 2},{\textbf {y}}_{n}+{1 \over 2}{\textbf {k}}_{2}\right),

{\textbf {k}}_{4}=h\ {\textbf {f}}\left(x_{n}+h,{\textbf {y}}_{n}+{\textbf {k}}_{3}\right).

где $h$ — величина шага сетки по $x$

Этот метод имеет четвёртый порядок точности, то есть суммарная ошибка на конечном интервале интегрирования имеет порядок $O(h^{4})$ (ошибка на каждом шаге порядка $O(h^{5})$ ).

Прямые методы Рунге — Кутты

Семейство прямых методов Рунге — Кутты является обобщением метода Рунге — Кутты 4 порядка. Оно задаётся формулами

{\textbf {y}}_{n+1}={\textbf {y}}_{n}+\sum _{i=1}^{s}b_{i}{\textbf {k}}_{i},

где $h$ — величина шага сетки по $x$ и вычисление нового значения проходит в $s$ этапов:

{\begin{array}{ll}{\textbf {k}}_{1}=&h{\textbf {f}}(x_{n},{\textbf {y}}_{n}),\\{\textbf {k}}_{2}=&h{\textbf {f}}(x_{n}+c_{2}h,{\textbf {y}}_{n}+a_{21}{\textbf {k}}_{1}),\\\cdots &\\{\textbf {k}}_{s}=&h{\textbf {f}}(x_{n}+c_{s}h,{\textbf {y}}_{n}+a_{s1}{\textbf {k}}_{1}+a_{s2}{\textbf {k}}_{2}+\cdots +a_{s,s-1}{\textbf {k}}_{s-1})\end{array}}

Конкретный метод определяется числом $s$ и коэффициентами $b_{i},a_{ij}$ и $c_{i}$ . Эти коэффициенты часто упорядочивают в таблицу (называемую таблицей Бутчера)

{\begin{array}{c|ccccc}0&&&&&\\c_{2}&a_{21}&&&&\\c_{3}&a_{31}&a_{32}&&&\\\vdots &\vdots &\vdots &\ddots &&\\c_{s}&a_{s1}&a_{s2}&\dots &a_{ss-1}&\\\hline &b_{1}&b_{2}&\dots &b_{s-1}&b_{s}\end{array}}

Для коэффициентов метода Рунге — Кутты должны быть выполнены условия $\sum _{j=1}^{i-1}a_{ij}=c_{i}$ для $i=2,\ldots ,s$ . Если требуется, чтобы метод имел порядок $p$ , то следует так же обеспечить условие

{\bar {\textbf {y}}}(h+x_{0})-{\textbf {y}}(h+x_{0})=O(h^{p+1}),

где ${\bar {\textbf {y}}}(h+x_{0})$ — приближение, полученное по методу Рунге — Кутты. После многократного дифференцирования это условие преобразуется в систему полиномиальных уравнений на коэффициенты метода.

Произношение

Согласно грамматическим нормам русского языка, фамилия Ку́тта склоняется, поэтому говорят: «Метод Ру́нге — Ку́тты». Правила русской грамматики предписывают склонять все мужские и женские фамилии, оканчивающиеся на -а, -я, которым предшествует согласный. Единственное исключение — фамилии французского происхождения с ударением на последнем слоге типа Дюма́, Золя́. Однако, иногда встречается несклоняемый вариант «Метод Ру́нге — Ку́тта» (например, в книге ^[4]).

Решение систем ОДУ

Метод Ру́нге — Ку́тты непосредственно обобщается на случай систем обыкновенных дифференциальных уравнений путём записи системы и метода в векторной форме.

Пример программы

y" + 4y = cos(3x) | y(0) = 0.8, y'(0) = 2, x = [0,1], h = 0.1

y" = cos(3x) - 4y

замена y'=z

получаем систему

y' = z = g(x,y,z)

z' = cos(3x) - 4y = f(x,y,z)

код на 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 = f(Xo, Yo, Zo);
			q1 = g(Xo, Yo, Zo);
			
			k2 = f(Xo + h/2.0, Yo + (h*k1)/2.0, Zo + (h*q1)/2.0);
			q2 = g(Xo + h/2.0, Yo + (h*k1)/2.0, Zo + (h*q1)/2.0);
			
			k3 = f(Xo + h/2.0, Yo + (h*k2)/2.0, Zo + (h*q2)/2.0);
			q3 = g(Xo + h/2.0, Yo + (h*k2)/2.0, Zo + (h*q2)/2.0);
			
			k4 = f(Xo + h, Yo + h*k3, Zo + h*q3);
			q4 = g(Xo + h, Yo + h*k3, Zo + h*q3);
			
			Z1 = Zo + (h/6.0)*(k1 + 2.0*k2 + 2.0*k3 + k4);
			Y1 = Yo + (h/6.0)*(q1 + 2.0*q2 + 2.0*q3 + q4);
			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);
            // решаем до 10 секунд
            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 (версии 5.3):

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

В среде MATLAB (довольно распространённый и удобный язык для технических вычислений) для решения системы уравнений необходимо сначала записать функцию, вычисляющую производные, т.е. функции y = g(x,y,z) и z = cos(3x) - 4y = f(x,y,z), о чём сказано выше. Для этого в одной из папок, к которой имеется доступ из системы MATLAB нужно создать текстовый файл runge.m со следующим содержимым:

MATLAB, runge.m

function Dy = runge(x, y)
Dy = y(:);
Dy(1) = y(2);
Dy(2) = cos(3*x) - 4*y(1);

Имя файла и имя функции должно совпадать, но оно может быть любым неиспользуемым ранее.

Затем необходимо создать главный файл c именем, например, main.m, который который будет выполнять основные вычисления. Этот главный файл будет содержать следующий текст:

MATLAB, main.m

clear; 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.
↑ ¹ ² J. C. Butcher. Numerical Methods for Ordinary Differential Equations. The University of Auckland, New Zealand.
↑ Б. П. Демидович, И. А. Марон, Э. З. Шувалова. Численные методы анализа, 3-е изд. — М.: Наука, 1967.

[1] Бахвалов Н. С., Жидков Н. П., Кобельков Г. М. Численные методы. — М.: Бином, 2001 — с. 363—375.

[2] Ильина В. А., Силаев П. К. Численные методы для физиков-теоретиков. т. 2. — Москва-Ижевск: Институт компьютерных исследований, 2004. — с. 16-30.

[ReferenceA-3] ¹ ² J. C. Butcher. Numerical Methods for Ordinary Differential Equations. The University of Auckland, New Zealand.

[4] Б. П. Демидович, И. А. Марон, Э. З. Шувалова. Численные методы анализа, 3-е изд. — М.: Наука, 1967.

[1]

[2]

[3]

[4]

Метод Рунге — Кутты

Содержание

Классический метод Рунге — Кутты 4 порядка

Прямые методы Рунге — Кутты

Произношение

Решение систем ОДУ

Пример программы

См. также

Ссылки

Навигация

Метод Рунге — Кутты

Классический метод Рунге — Кутты 4 порядка

Прямые методы Рунге — Кутты

Произношение

Решение систем ОДУ

Пример программы

См. также

Ссылки

Навигация

Поиск