6.3. Методы Рунге-Кутты второго порядка

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

Так, например, для второго порядка получено однопараметрическое семейство схем вида [1]

где свободный параметр,

Локальная погрешность схем (6.14) имеет третий порядок, глобальная - второй; т.е. решение ОДУ полученное по этой схеме, равномерно сходится к точному решению с погрешностью

Для параметра наиболее часто используют значения В первом случае формула (6.14) приобретает вид

геометрическая интерпретация которой представлена на рис. 6.3 [1]. Вначале вычисляется приближенное решение ОДУ в точке по формуле Эйлера Затем определяется наклон интегральной кривой в найденной точке и после нахождения среднего наклона на шаге находится уточненное значение Схемы подобного типа называют "прогноз-коррекция”, что подразумевает грубое вычисление решения по формуле низкого порядка, а затем уточнение с учетом полученной информации о поведении интегральной кривой.

Рис. 6.3. Метод второго порядка (

Рис. 6.4. Метод Рунге-Кутты второго порядка

С целью экономии памяти при программировании алгоритма (6.15), обобщенного на системы ОДУ, изменим его запись с учетом того, что

где k — номер решения для системы ОДУ. Теперь не придется держать в памяти ЭВМ массив начальных значений его можно «забыть» после

вычисления значений эйлеровских приближений размещаемых на месте массива Хотя по сравнению с методом Эйлера схема (6.16) требует дополнительного массива для запоминания значений

Во втором случае при от формулы (6.14) переходим к схеме

геометрический смысл которой отражает рис. 6.4. Здесь при прогнозе определяется методом Эйлера решение в точке

а после вычисления наклона касательной к интегральной кривой в средней точке решение корректируется по этому наклону.

Формула (6.17) обобщается на системы ОДУ аналогично схеме (6.16). По сравнению с программой метода Эйлера для сохранения начальных значений придется ввести дополнительный массив.

Схему (6.17) можно получить из метода Эйлера с помощью первой и второй формул Рунге без использования общего соотношения (6.14) для методов второго порядка.

Программы 6.2 и 6.3 предназначены для решения задачи Коши для системы ОДУ с помощью первого и второго вариантов метода Рунге-Кутты второго порядка

Основной блок программы отличается от соответствующего блока программы наличием описания дополнительного массива для размещения правых частей системы ОДУ, вычисленных в начальной точке В строках 100-110 подпрограммы метода Рунге-Кутты вычисляется эйлеровское приближение в точке и запоминаются начальные производные в массиве С помощью операторов в строках 120-130 реализована правая часть формулы (6.16). Величина вычисляется вне цикла для экономии времени.

В программах 6.2 и 6.3 решается для примера та же система, полученная из уравнения Бесселя, что и в программах 6.1. В соответствующих программах на языках Фортран и Паскаль схемы Рунге-Кутты второго порядка реализованы в виде подпрограмм с именами с формальными параметрами, имеющими тот же смысл, что и в подпрограммах (программы 6.1). Обратите внимание на использование разных фактических параметров при обращениях к подпрограмме вычисления правых частей ОДУ с именем

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

Примеры расчетов по программам 6.2 и 6.3 сведены в табл. 6.1. Задача Коши для уравнения Бесселя при значении параметра на интервале

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

Таблица 5.1

(см. скан)

(см. скан)

(см. скан)

(см. скан)