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

 Одна из основных проблем, связанных с применением методов Рунге-Кутты, состоит в получении надежных и, в то же время простых и эффективных оценок погрешности решения задачи Коши. В решении данной проблемы значительные успехи достигнуты за счет введения вложенных методов Рунге-Кутты (ВМРК) или методов вложенных форм. Основные соображения, положенные в основу этих методов, таковы: нужно получить формулы из семейства Рунге-Кутты, которые использовали бы одни и те же значения функций – правой части уравнения – и определяли бы разные конкретные методы одного или смежных порядков.

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

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

Обыкновенные дифференциальные уравнения (ОДУ) широко используются для  математического моделирования  процессов и явлений в различных областях науки и техники. Переходные процессы в радиотехнических цепях, кинетика химических реакций, динамика биологических популяций, движение космических объектов, модели экономического развития исследуются с помощью ОДУ.

В дифференциальное уравнение  n-го порядка в качестве неизвестных величин входят функция у(х) и ее первые п производных по аргументу х

φ(х, y, y', …, у⁽ⁿ⁾)=0                                       (1)

Из теории ОДУ известно, что уравнение ф. 1 эквивалентно системе n уравнений первого порядка

φ(х, y₁, y'₁, y₂, y'₂, …, у_n, y'_n)=0                               (2)

 

где к = 1,..., n.

Уравнение (1) и эквивалентная ему система (2) имеют бесконечное множество решений. Единственные решения выделяют с помощью дополнительных условий, которым должны удовлетворять искомые решения. В зависимости от вида таких условий рассматривают три типа задач, для которых доказано существование и единственность решений.

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

у(х₀)= у₀ ,       у'(х₀)= у₁₀ ,…, у^(n-1)(x₀)=y_n-1,0

Для системы ОДУ типа ф. 2 начальные условия задаются в виде

у₁(х₀)= у₁₀ ,       у₂(х₀)= у₂₀ ,…, у_n(x₀)=y_n0

Ко второму типу задач  относятся так называемые граничные, или краевые задачи, в которых  дополнительные условия задаются в  виде функциональных соотношений между  искомыми решениями. Количество условий  должно совпадать с порядком п  уравнения или системы. Если решение  задачи определяется в интервале  х є [х₀, х_k], то такие условия могут быть заданы как на границах, так и внутри интервала. Минимальный порядок ОДУ, для которых может быть сформулирована граничная задача, равен двум.

Третий тип задач для  ОДУ - это задачи на собственные значения. Такие задачи отличаются тем, что  кроме искомых функций у(х) и  их производных в уравнения входят дополнительно m неизвестных параметров λ₁, λ₂,…, λ_m, которые называются собственными значениями. Для единственности решения на интервале [х₀, х_k] необходимо задать n + m граничных условий. В качестве примера можно назвать задачи определения собственных частот, коэффициентов диссипации, структуры электромагнитных полей и механических напряжений в колебательных системах, задачи нахождения фазовых коэффициентов, коэффициентов затухания, распределения напряженностей полей волновых процессов и т. д.

К численному решению ОДУ  приходится обращаться, когда не удается построить аналитическое решение задачи через известные функции. Хотя для некоторых задач численные методы оказываются более эффективными даже при наличии аналитических решений.

Большинство методов решения  ОДУ основано на задаче Коши.

Метод Рунге-Кутта 4 порядка

Пусть дано дифференциальное уравнение  первого порядка:

 с начальным условием y(x0)=y0.

Выберем шаг h и введем обозначения:

xi = x0 + ih  и yi = y(xi),   где   i = 0, 1, 2,....

Аналогично описанному выше методу производится решение дифференциального  уравнения. Отличие состоит в  делении шага на 4 части.

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

yi+1 = yi +∆yi, где i = 0, 1, 2...

∆y=(k1+2*k2+2*k3+k4)/6

a числа k1, k2 ,k3, k4 на каждом шаге вычисляются по формулам:

k1 = h*f(xi, yi )

k2 = f (xi +h/2, yi +k1 /2)*h

k3 = F(xi +h/2, yi +k2 /2)*h

k4 = F(xi +h, yi +k3 )*h

Это явный четырехэтапный метод 4 порядка  точности.

F(x, у) - заданная функция - должна  быть описана отдельно.

Входные параметры:

Х0, XК - начальное и конечное значения независимой переменной;

Y0 – значение y0 из начального условия y(x0)=y0;

N - количество отрезков разбиения;

Выходные параметры:

Y - массив значений искомого решения в узлах сетки.

Этот метод является наиболее распространенным методом решения  систем при шаге h=const. Его достоинством является высокая точность-погрешность - и меньшая склонность к возникновению неустойчивости решения. Алгоритм реализации метода заключается в циклических вычислениях Y_j(i+1) на каждом i+1 шаге по следующим формулам:

                                    (3)

                                        (4)

                                         (5)

                                 (6)

                         (7)

При переходе от одной формулы  к другой задаются или вычисляются  соответствующие значения x и Y_j и находятся по подпрограмме значения функции F_j(x,Y_j).

Метод Рунге-Кутта-Мерсона

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

Автоматическое изменение  шага в ходе решения систем дифференциальных уравнений необходимо, если решение  требуется получить с заданной точностью. При высокой точности (погрешность  ) и решении в виде кривых с сильно различающейся крутизной автоматическое изменение шага обеспечивает уменьшение общего числа шагов в несколько раз, резко уменьшается вероятность числовой неустойчивости, даёт более равномерное расположение точек графика кривых (решений) при их выводе на печать. Данный метод обеспечивает приближённую оценку погрешностей на каждом шаге интегрирования. Погрешность интегрирования имеет порядок h⁵. Этот метод реализуется следующим алгоритмом: Задаём число уравнений N, погрешность ε=E, начальный шаг интегрирования h=H и начальное значение y₁₀,…,y_N0. С помощью пяти циклов с управляющей переменной J=1,2,..,N вычисляем коэффициенты:

                  (8)

                                (9)

                       (10)

                       (11)

                (12)

Находим (в последнем цикле) значение (13)

                       (13)

И погрешность

                 (14)

Проверяем выполнения условий

                                              (15)

                                           (16)

Если условие (15) не выполняется, то делим шаг h на 2 и повторяем вычисления. Если это условие выполняется и выполняется условие (16), значение x_i+1=x_i+h и Y_j(i+1), то считаем, что решение системы дифференциальных уравнений найдено с заданной точностью. Если условие (16) не выполняется , шаг h увеличивается вдвое и вычисления повторяются.

Идея вложенных форм для методов  Рунге-Кутты реализована в различных  вычислительных схемах той или иной эффективности. Рассмотрим один из первых методов такого типа: пятиэтапный  метод Кутты-Мерсона четвертого порядка.

Метод Кутты-Мерсона (18) обладает той особенностью, что имеет порядок только для линейных уравнений с постоянными коэффициентами. Рассмотрим различные способы реализации метода вложенных форм Кутты-Мерсона при решении системы однородных линейных обыкновенных дифференциальных уравнений с постоянными коэффициентами (СЛОДУ), в силу того, что лучшие результаты схема дает именно в этом случае.

Общий вид задачи Коши для  однородных СЛОДУ с постоянными  коэффициентами:

 

 

                                                   (17) 

где – вектор неизвестных; – вектор начальных условий; A – матрица оэффициентов линейной системы.

Численное решение системы (17) последовательно можно получить по следующей схеме Кутты-Мерсона:

 

 

 

 

                       (18) 

Выполним элементарные преобразования системы (18) для оценки времени последовательной реализации одного шага интегрирования и оценим время вычисления с учетом:

 

 

 

                                          (19) 

Время, затраченное на выполнение последовательного вычисления шаговых  коэффициентов по схеме (19):

                                 (20)

где – время умножения матрицы на вектор; – время умножения вектора на скаляр; – время выполнения сложения векторов.

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

                             (21)

Предложенная схема имеет достаточно высокую вычислительную сложность. Для сокращения вычислений представим величины и на каждом шаге, как функции матрицы коэффициентов системы и шага интегрирования .

После упрощающих преобразований получим:

 

 

                        (22)

 

 

Вычислительная схема (22) позволяет уменьшить время каждого  шага интегрирования за счет создания подготовительного этапа, который  будет выполняться до начала основного  счета и содержать операции, не связанные непосредственно с  номером шага. Для вложенных методов  величина шага интегрирования h – не есть постоянная величина, поэтому основной операцией подготовительного этапа будет: вычисление степеней матрицы коэффициентов системы . Заметим, что эта операция являлась наиболее трудоемкой в схеме (22). Время подготовительного этапа вычисляется, как:

где – время вычисления скалярного произведения матриц; – время умножения матрицу на скаляр. Определим время выполнения одного шага интегрирования по схеме (23) с учетом того факта, что нет необходимости в вычислении обоих приближений, достаточно знать и :

                        (23)

Одной из модификаций метода Рунге-Кутта является метод Кутта-Мерсона (или пятиэтапный метод Рунге-Кутта четвертого порядка), который состоит в следующем.

1. На i-м шаге рассчитываются коэффициенты:

2. Вычисляем приближенное  значение  по формуле:

                             (24)

3. Вычисляем приближенное  значение  по формуле:

                                 (25)

4. Вычисляем оценочный коэффициент по формуле:

                                       (26)

5. Сравниваем R с точностью вычислений ε. Если R≥ ε , то уменьшаем шаг вдвое и возвращаемся к п.1. Если R< ε , то значение, вычисленное по формуле (25), и будет вычисленным значением (с точностью ε).

6. Перед переходом к  вычислению следующего значения  x, сравниваем R с ε/64 . Если R<ε/64, то дальнейшие вычисления можно проводить с удвоенным шагом h=2 h .