Материал: Понятие о численных методах решения обыкновенных дифференциальных уравнений

Полагая k=0,1,2,3 и т. д., мы будем получать последовательные значения угла и соответствующие ему значения угловой скорости u(t) в последовательные моменты времени

Совершенно аналогичным образом систему двух ОДУ второго порядка, которая описывает движения планет вокруг Солнца (31), можно свести к системе из четырёх уравнений первого порядка, если ввести две дополнительные переменные, представляющие собой скорости движения тела m вдоль координатных осей х и у соответственно:

Тогда исследуемая система принимает вид:

                                                                              (12)

Заменяя в системе (12) каждую производную её допредельным образом, получаем явные формулы для решения этой системы методом Эйлера.

. Недостатки метода Эйлера

Существует простая геометрическая интерпретация метода Эйлера. Рассмотрим снова задачу Коши (3) для одного ОДУ первой степени и соответствующее ему в методе Эйлера рекуррентное соотношение (8)

уравнение дифференциальный коши задача

Рис. 1

На плоскости (х,у) каждому частному решению задачи Коши, которая выделяется начальным условием , отвечает некоторая кривая, которая называется интегральной кривой. Изменяя , мы переходим от одной к другой интегральной кривой. Более того, можно доказать, что через каждую точку плоскости (х,у) проходит одна и только одна интегральная кривая.

Пусть решению рассматриваемой нами задачи Коши на рис. 1 отвечает жирная интегральная кривая. Рекуррентное уравнение метода Эйлера (8) можно записать в форме

                                                                                   (13)

Поскольку, согласно решаемому нами ОДУ, . Если считать параметр h непрерывной переменной, то легко понять, что уравнение (13) представляет собой уравнение касательной, проведённой в точке  к интегральной кривой у(х), проходящей через эту точку.

Таким образом, на каждом шаге метода Эйлера мы заменяем истинную интегральную кривую отрезком касательной, проведённой к этой кривой в начале микроинтервала [,]. Тем самым мы отклоняемся от искомой интегральной кривой на некоторую маленькую величину, причём её малость определяется малостью шага h.

В результате, на втором шаге метода Эйлера  мы строим касательную в точке  не к истинной, а к другой интегральной кривой - той, которая проходит через эту точку и, таким образом, ещё больше отклоняемся от искомой интегральной кривой, удовлетворяющей начальному условию .

Таким образом, в методе Эйлера мы заменяем искомую интегральную кривую некоторой ломаной линией, которая, по мере отдаления  от начальной точки , всё более и более отклоняется от этой интегральной кривой. Иными словами, абсолютная погрешность метода Эйлера имеет тенденцию к нарастанию по мере увеличения числа шагов этого метода.

Может показаться, что если выбрать шаг интегрирования h более маленьким, то можно существенным образом уменьшить эту погрешность. В общем случае это, однако, не так. Действительно, если мы должны получить решение исходной задачи Коши (3) на некотором заданном макроинтервале [,], то уменьшение шага h влечёт за собой увеличение числа шагов интегрирования N, поскольку . «Локальная погрешность» (погрешность на одном шаге) уменьшается при уменьшении h, но увеличение числа шагов может привести к росту «глобальной погрешности» на заданном интервале [,]. Вышеуказанная ситуация характерна для неустойчивых вычислительных процессов, использование которых на практике может привести к катастрофическим последствиям. В силу этого, необходимо рассматривать такие методы численного решения ОДУ, которые порождают достаточно устойчивые численные алгоритмы.

Весьма распространёнными и хорошо зарекомендовавшими себя на практике для решения ОДУ являются методы Рунге-Кутты. Это целый класс методов и мы, в качестве примера, рассмотрим так называемый четырёхточечный метод Рунге-Кутты

. Четырёхточечный метод Рунге-Кутты

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

Заметим, прежде всего, что решение дифференциального уравнения (3а) фактически определяет зависимость первой производной от двух независимых переменных - х и у. Это очень хорошо видно из рис. 1: фиксируя х, мы имеем бесконечное множество значение производной , поскольку каждая из пересекаемых вертикальной линией х=const интегральных кривых имеет своё направление касательной.

Описываемый метод Рунге-Кутты, как и метод Эйлера, состоит из последовательности шагов величиной h, но, в отличие от последнего метода, на каждом шаге h находится не одно значение производной (в методе Эйлера находилось лишь ), а несколько. Они соответствуют разным значениям аргументов функции f(х,у). В четырёхточечном методе Рунге-Кутты (откуда и происходит его название) находится четыре различных значения и делается некоторое специфическое их усреднение (то есть берётся не простое среднее арифметическое значение этих производных!). После этого делается перемещение из точки в точку  по прямой в направлении тангенса угла наклона, которое определяется этим усреднённым значением производной. Заметим, что разные варианты методов Рунге-Кутты отличаются друг от друга стратегией выбора точек, в которых находятся производные на микроинтервале h и формулой усреднения значения этих производных.

Рис. 2

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

Этап I. Находим производную к интегральной кривой в точке 1 с координатами , подставляя эти координаты в правую часть дифференциального уравнения, т.е. вычисляя значение .

Из точки 1:  перемещаемся на полшага вперёд по прямой, направление которой задаётся этим значением производной (то есть по прямой с тангенсом угла наклона к оси абсцисс, равным значению  (этот этап метода Рунге-Кутты полностью аналогичен шагу метода Эйлера с шагом h/2). В результате, в плоскости (х,у) мы переходим в точку

.

Этап II. Через найденную таким образом точку 2 проходит своя интегральная кривая и мы находим направление касательной к ней, то есть вычисляем значение производной

Далее делается полшага вперёд с найденным значением производной (), но снова из начальной точки микроинтервала [,]. Таким образом, мы переходим в плоскости (х,у) в точку

.

Этап III. В этой точке находим значение производной , подставляя её координаты в правую часть дифференциального уравнения:

Эта производная определяет направление касательной к интегральной кривой, проходящей через точку 3.

Этап IV. Из начальной точки 1:  делаем на сей раз полный шаг вперёд (на величину h по оси х) по прямой, в направлении, которое определяется значением производной . В результате мы переходим в точку

Находим производную  в этой точке подстановкой её координат в правую часть дифференциального уравнения:

В результате четырёх описанных выше этапов мы нашли четыре значения производных. Производим их усреднение по формуле

                                                                    (14)

Таким образом,  является некоторым средневзвешенным значением найдённых четырёх производных: двум «внутренним» значениям производном соответствуют весовые множители 2, а двум крайним - множители 1 [деление в формуле (14) производится на сумму этих четырёх весовых множителей: 6=1+2+2+1].

Далее мы перемещаемся по прямой из начальной точки 1: в направлении, тангенс угла наклона которого к оси абсцисс определяется средним значением производной из формулы (14). Таким образом, из начальной точки с координатами переходим в точку плоскости (х,у) с координатами

.

Иными словами, полный шаг метода Рунге-Кутты определяется формулами

                                                                                  (15)

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

Более того, порядок точности рассматриваемого метода Рунге-Кутты на одном шаге величины h оценивается формулой

,                                                                                    (16)

где . Здесь  есть пятая производная от искомого решения дифференциального уравнения (3а) в некоторой точке на микроинтервале . Таким образом, локальная погрешность метода Рунге-Кутты (то есть погрешность на одном шаге h) пропорциональна пятой степени шага h и пятой производной искомого решения дифференциального уравнения.

Оценки точности типа (16) позволяют грубо оценить величину шага интегрирования шага h, необходимого для достижения требуемой точности решения исходного дифференциального уравнения. Для метода Эйлера аналогичная погрешность на одном шаге определяется формулой

Таким образом, четырёхточечный метод Рунге-Кутты на три порядка по шагу точнее метода Эйлера (например, при h=0.01 точность метода Рунге-Кутты в миллион раз выше точности метода Эйлера).

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

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

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