Материал: Математическое моделирование в экологии

11.3.5. Метод Рунге – Кутта

Пусть требуется найти решение дифференциального уравнения

y’ = f (t, y),

удовлетворяющее начальному условию

y’(t₀) = y₀.

Принцип, на котором основан метод Рунге – Кутта, можно пояснить, как и принцип, на котором основан метод Эйлера, с помощью разложения функции в ряд Тейлора

Чтобы удержать в ряде Тейлора член n-го порядка, необходимо вычислить n-ю производную зависимой переменной. При использовании модифицированного метода Эйлера для получения второй производной в конечно-разност-ной форме достаточно было знать наклон кривой на концах рассматриваемого интервала. Чтобы вычислить третью производную в конечно-разностном виде, необходимо иметь значения второй производной, по меньшей мере, в двух точках. Для этого необходимо дополнительно определить наклон кривой в некоторой промежуточной точке интервала h, т. е. между tn и n 1 t + . Очевидно, чем выше порядок вычисляяемой производной, тем больше дополнительных точек потребуется вычислить внутри интервала. Так как существует несколько способов расположения внутренних точек и выбора относительных весов для найденных производных, то метод Рунге – Кутта, в сущности, объединяет целое семейство методов решения дифференциальных уравнений.

Наиболее распространенным из них является метод четвёртого порядка точности, при котором удерживаются все члены ряда Тейлора, включая h⁴. Расчеты при использовании этого классического метода производятся по формулам:

где

Метод Эйлера и его модификация по сути дела являются методами Рунге – Кутта первого и второго порядка соответственно. Более высокая точность метода Рунге – Кутта позволяет увеличить шаг интегрирования h. Допустимая погрешность на шаге определяет его максимальную величину. В прикладных пакетах программ выбор шага часто осуществляется автоматически. Для этого проводят вычисления сначала с шагом h, а затем – с шагом h/2.

За оценку погрешности вычислений с шагом h/2 можно принять приближенную формулу

где - вычисленное значение с шагом h/2; y_n – вычисленное значение с шагом h. Пример: y’ = xy.

Решение:

При реализации методов Рунге – Кутта на ЭВМ для каждой точки проводят двойной счет. Если полученные при этом значения удовлетворяют выражению (5.4), то для точки t n+1 шаг удваивают, в противном случае уменьшают вдвое. Однако необходимо помнить, что выражение (5.4) приближенное и при неблагоприятных условиях можно получить совершенно ошибочные результаты, хотя в большинстве случаев дело обстоит благополучно.

11.3.6. Метод Рунге – Кутта для систем дифференциальных уравнений

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

можно принять

тогда

и получаем два уравнения первого порядка:

Задача Коши в этом случае содержит два начальных условия

y(t₀) = y₀ и z(t₀) = z₀.

Формулы Рунге – Кутта для рассматриваемого случая имеют

вид:

где

11.3.7. Общая характеристика одношаговых методов

Всем одношаговым методам присущи определенные общие черты:

1) чтобы получить информацию в новой точке, надо иметь данные лишь в одной предыдущей точке. Поэтому одношаговые методы называют «самостартующимися»;

2) в основе всех одношаговых методов лежит разложение функции в ряд Тейлора, в котором сохраняются члены, содержащие h в степени до k включительно. Целое число k называется порядком метода;

3) все одношаговые методы не требуют действительного вычисления производных, а вычисляется лишь сама функция в одной или нескольких промежуточных точках;

4) свойство «самостартования» позволяет легко менять величину шага.

3.8. Многошаговые методы

Вернемся к задаче Коши

y’ = f(t, y); y(0) = y₀. (3.13)

В предыдущих методах значение y_{n + 1} зависело только от информации в пре-дыдущей точке t_n. Кажется вполне вероятным, что можно добиться большей точности, если использовать информацию о нескольких предыдущих точках t_n, t_{n - 1},…. Именно так и поступают в многошаговых методах [4].

Если проинтегрировать уравнение (3.13) на отрезке [t_n, t_n+1], то получим

или

где P(t) – полином, аппроксимирующий f(t, y).

Чтобы построить полином степени N, используем предыдущие решения y_n, y_{n – 1}, …. в точках t_n,…, t_n-1, …, t_{n - N}, ... . Мы по-прежнему считаем, что узлы t расположены равномерно с шагом h. В принципе, можно проинтегрировать этот полином явно, что ведет к следующему методу:

(3.14)

В простейшем случае, когда N = 0, полином P – есть константа, равная f_n, и (3.14) превращается в обычный метод Эйлера.

Если N = 1, то P – есть линейная функция, проходящая через

точки (t_n-1, f_{n -1}) и (t_n, f_n) , т. е.

Интегрируя этот полином от t_n до t_{n - 1}, получаем следующий метод:

(3.15)

который является двухшаговым, поскольку использует информацию в двух точках - t_n и t_{n - 1}. Аналогично, если N = 2 , то P есть квадратичный полином, интерполирующий данные (t_n-2 , f_n-2 ), (t_n-1 , f_n-1), (t_n , f_n), а соответствующий метод имеет вид

(3.16)

Если N = 3, то интерполяционный полином является кубическим, а соответствующий метод определяется формулой

(3.17)

Отметим, что метод (3.16) является трехшаговым, а (3.17) – четырехшаговым. Формулы (3.15) – (3.17) известны как методы Адамса –Башфорта. Метод (3.15) имеет второй порядок точности, поэтому его называют методом Адамса – Башфорта второго порядка. Аналогично, методы (3.16) и (3.17) называют соответственно методами Адамса – Башфорта третьего и четвертого порядков.

Этот процесс, в принципе, можно бы продолжить, используя все большее число предыдущих точек, а следовательно, и интерполяционный полином P более высокой степени, и получить Адамса – Башфорта сколь угодно высокого порядка. Однако точность вычислений с увеличением порядка возрастает нелинейно. Чем дальше отстоит предыдущая точка от текущей точки, тем слабее она влияет на точность. Многошаговые методы порождают проблему, которая не возникала при использовании одношаговых методов. Так как в

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

11.3.9. Методы прогноза и коррекции

Методы Адамса – Башфорта используют уже сосчитанное значение в точке n, t и в предыдущих точках. В принципе, при построении интерполяционного полинома мы можем использовать и точки t_n+1, t_n+2 и т.д. Простейший случай при этом состоит в использовании точек t_n+1, t_n ,…, t_n-N и построении интерполяционного полинома степени N +1. При этом возникает класс методов, известный как методы Адамса - Моултона. Если N = 0, то P - линейная функция, проходящая через точки (t_n , f_n) и (t_n+1 , f_n+1), и соответствующий метод

(3.18)

является методом Адамса - Моултона второго порядка.

Если, N = 2 , то P - кубический полином, построенный по точкам (t_n+1, f_n+1), (t_n, f_n), (t_n-1, f_n-1) и (t_n-2, f_n-2) и соответствующий метод

(3.19)

является методом Адамса-Моултона четвёртого порядка.

Заметим, что в формулах (3.18) и (3.19) значение f_n+1 неизвестно. Дело в том, что для вычисления f(t_n+1, y_n+1)=f_n+1 нужно знать три значения y_n+1, которое само пока является неизвестным. Например, соотношение (3.18) является уравнением

(3.20)

относительно неизвестного значения y_n+1. То же самое справедливо и относительно (3.19). Следовательно, методы Адамса-Моултона определяют y_n+1 неявно и в силу этого называются неявными. В то же время, методы Адамса – Башфорта называются явными, поскольку они для нахождения значения y_n+1 не требуют решения никаких уравнений. На практике обычно не решают уравнение (3.20), а используют совместно явную и неявную формулы, что приводит к методу прогноза и коррекции. Одним из широко используемых методов прогноза и коррекции является объединение методов Адамса четвёртого порядка (3.17) и (3.19):