CC BY
276
Cloud of Science. 2016. T. 3. № 1 http:/ / cloudofscience.ru ISSN 2409-031X
Сравнительный анализ методов интегрирования с плавающим шагом
С. В. Корягин, А. А. Яковлев
Московский технологический университет МИРЭА 119571, Москва, пр-кт Вернадского, 78
e-mail: dongenealog2003@mail.ru
Аннотация. В данной работе проводится анализ точности и сравнение двух популярных парных методов численного интегрирования 4-5 порядка: метод Рунге-Кутта-Фельберга и метод Рунге-Кутта-Дормана-Принса 4-5 порядка с автоматическим изменением шага интегрирования, реализованным в популярной среде для решения задач технических вычислений MATLAB с применением транслятора, автоматически формирующего необходимую статистическую информацию. Ключевые слова: численное интегрирование, метод Рунге-Кутта-Фельберга 4-5 порядка, модификация Дормана-Принса.
1. Введение
Актуальность темы обоснована широкой популярностью численного интегрирования, которое используется в случае отсутствия аналитического решения, либо громоздкостью найденного аналитического решения. Помимо того, преимущества численных методов по сравнению с аналитическими заключаются в относительной простоте их реализации в виде алгоритмов и вычислительных программ. У данных методов безусловно есть и недостатки и прежде всего это большие затраты времени вычислений и, в зависимости от параметров задачи и выбранного метода, высокие требования к объему памяти [1].
Проанализировав возможности современных методов численного интегрирования расчета, следует сделать вывод, что интегрирование нежестких систем дифференциальных уравнений, а также систем малой жесткости можно достаточно эффективно проводить явными парными методами с изменяющейся сеткой шага [2]. Проведем сравнительный анализ метода Рунге-Кутта-Фельберга 4-5 порядка, а также метода Рунге-Кутта в модификации Дормана-Принса 4-5 порядка. Оба метода предусматривают автоматическое изменение шага, а также возможность контроля погрешности интегрирования. Для методов 4-5 порядка погрешность метода имеет порядок И6. Широкое применение эти методы получили в системах непрерывного моделирования [3-6]. Метод Рунге-Кутта-Фельберга обеспечивает погрешность интегрирования того же порядка, что и методы Рунге-Кутта, но позволяет повысить устойчивость решения.
2. Метод Рунге-Кутта-Фельберга 4-5 порядка (РКФ45)
При интегрировании данным методом на каждом шаге рассчитываются 6 промежуточных точек [7]. Реализация метода предполагает построение неравенства для контроля точности вычислений, которое при правильном выборе погрешности не приводит к дополнительным вычислительным затратам. Экономия вычислительных ресурсов обеспечивается тем, что в алгоритме реализована оценка погрешности более точного решения 5-го порядка, позволяющая избежать двойного расчета дифференциальных уравнений. Это обеспечивается вычислением разности двух решений 4-го и 5-го порядка. Найденная оценка может использоваться для корректировки величины шага приращения аргумента.
Для нахождения нового значения неизвестной функции уи+1 последовательно вычисляются величины:
( г-1 \
к = К/
К +«Л , Уп
V з=1 У
г = 1,6. (1)
Значение пятого порядка точности вычисляется как взвешенная комбинация величин к.:
Уп
6
= Уп + £у, • к,, (2)
Значения всех коэффициентов даны в табл. 1. Если в формуле (2) коэффициенты у1 заменить на коэффициенты у*, то получим решение четвертого порядка.
Фактически на практике вычисляют решение пятого порядка и оценку погрешности:
6
Ап+1 (УГ-У,-) • к,. (3)
¿=1
При выборе допустимой погрешности метода следует учесть вероятность непредсказуемого возрастания накопленной ошибки. Дело в том, что все вычисления выполняются всегда с фиксированной точностью простых арифметических операций. Точность можно увеличить с помощью специальных способов отображения чисел в памяти компьютера и применением соответствующих алгоритмов арифметических действий. Это приведет к еще большим вычислительным затратам, но главное то, что в современных компьютерах максимальная точность представления чисел и точность выполнения арифметических операций фиксирована. При каждой операции совершается некоторая ошибка, так называемая ошибка округления. При последовательном выполнении операций с плавающей точкой ошибки округления
¿=1
накапливаются [8]. Следовательно, уменьшение шага интегрирования с целью увеличить точность решения имеет свой предел.
Таблица 1. Значения коэффициентов в формулах (1)-(3)
а, Ру У, * У*
0 16 135 25 216
1 1 0 0
4 4
3 3 9 6 656 1408
8 32 32 12 825 2 565
12 1932 7 200 7 296 28 561 2197
13 2197 2197 2197 56430 4104
1 439 -8 3 680 845 9 1
216 513 4104 50 " 5
1 8 2 3 544 1859 11 2 0
2 27 2 565 4104 40 55
Следующее, на что следует обратить внимание, это максимальная погрешность при использовании метода. Ошибочный выбор погрешности интегрирования может приводить к непомерному увеличению времени вычислений. Следует заметить, что данный алгоритм считается одним из лучших среди методов типа Рунге-Кутта пятого порядка точности.
3. Метод Рунге-Кутта, модификация Дормана-Принса 4-5 порядка (РК-ДП45)
Модификация Дормана-Принса 4-5 порядка метода Метод Рунге-Кутта (РК-ДП45) имеет таблицу бутчера со следующими коэффициентами [9] (рис. 1).
Таким образом, при расчете очередного значения решения метод РК-ДП45 опирается на 7 промежуточных точек. В целом метод схож с методом РКФ45. Необходимо отметить, что выбор шага в методе РК-ДП45 в системе МЛТЬЛБ реализован не тривиально, качество которого для нас представляет интерес.
Проверку расчета методом РКФ45 выполним в одной из систем обыкновенных дифференциальных уравнений из практикума Р.З. Даутова [10]. Рассмотрим полет снаряда, выпущенного с начальной скоростью V под углом 9 к горизонту. Будем придерживаться допущения, что земля плоская и вся траектория снаряда лежит в одной плоскости хОу.
0
1 5 1
3 10 1 4 3 4
4 5 11 9 14 3 40 9
8 9 4843 1458 3170 243 8056 729 53 162
1 9017 355 46732 49 5103
3168 33 5247 176 18656
1 35 384 0 500 113 125 192 2187 6784 11 84
5179 0 7571 393 92097 187 1
57600 16695 640 339200 2100 40
35 384 0 500 1113 125 192 2187 6784 11 84 0
4-й порядок
5-й порядок
Рисунок 1. Коэффициенты для метода Дормана-Принса
Уравнения движения центра масс снаряда в проекциях на оси координат запишутся в виде:
CpSv2cos(0)
2m ' (4)
CpSv2 sin (0)
y" = --
x = —
2т
где т — масса снаряда; v= (х '2 + у|2)12 — скорость движения; 9= агС§( у'/х') — угол между касательной к траектории и осью Ох; g — ускорение силы тяжести;
— площадь поперечного сечения снаряда; р — плотность воздуха; С — коэффициент лобового сопротивления снаряда.
Сведем приведенную систему уравнений (1) к 4-м уравнениям 1-го порядка:
х' = V • со^(9); у' = V • sin(9);
, Ср8^ (5)
V =-----Я • sln(v);
2т
9' = -(я • ос8(9))/V.
Теперь решим данную систему уравнений с помощью методов РК-ДП45 и РКФ45, используя следующие коэффициенты и начальные данные: масса снаряда т = 43.51; коэффициент лобового сопротивления с = 0.15; плотность воздуха р = 1.29; площадь поперечного сечения снаряда 5 = 0.35; ускорение свободного падения: Я=9.81; х0 = 0; у0 = 0; ^ = 655; 90 = 1.2.
Максимально возможная ошибка для метода РКФ45 Етах = 0.001, будем использовать переменный шаг. Интервал интегрирования: от 0 до 50. Метод РК-ДП45 возвращает результат, приведенный на рис. 2. Количество успешных шагов: 24. Количество возвратов: 0. Количество расчетов системы: 145. Значения на конце интервала: 1664.9, -945.4, 111.7, -1.5. Время расчета: 0.041 с.
У -if
1 v x
^_____ i
t
8
10 15 20 25 30 35 40 45 50
Рисунок 2. Решение задачи определения параметров полета снаряда методом Рунге-Кутта в модификации Дормана-Принса 4-5 порядка
Результаты интегрирования той же системы уравнений, но методом РКФ45 приведены на рис. 3. Количество успешных шагов: 46. Количество возвратов: 4. Количество расчетов системы: 294. Значения на конце интервала: 1664.1, -946, 111.7, -1.5. Время расчета: 0.014 с.
У
% X
t \
8 \
0 5 10 15 20 25 30 35 40 45 50
Рисунок 3. Решение задачи определения параметров полета снаряда методом Рунге-Кутта-Фельберга 4-5 порядка
Данный расчет показал, что метод РКФ45 при расчете данной системы дифференциальных уравнений эффективнее метода РК-ДП45 с точки зрения времени расчета. Это наблюдается наряду с тем, что количество расчетов системы дифференциальных уравнений метода РКФ45 существенно больше. Причина может быть в более сложной реализации контроля устойчивости метода РК-ДП45.
Далее протестируем методы интегрирования с точки зрения правильности расчета. За пример возьмем алгебраическое уравнение, имеющее аналитическое решение:
у' = ео8(х)-х2 +sin(x)■2x. (6)
Решением данного уравнения будет являться уравнение у=8т(х)-х2. Решим данное уравнение численно теми же двумя способами и сравним их результаты с аналитическим решением. Будем использовать следующие исходные данные: у0 = 0; х от 0 до 50. Шаг интегрирования: 0.1. Максимальная погрешность Етах = 0.0001.
Результат, полученный методом РК-ДП45, приведен на рис. 4. Количество успешных шагов: 28. Количество возвратов: 3. Количество расчетов системы 217. Значения на конце интервала: -655.6. Время расчета: 0.102 с.
2500
2000 -
1500 -
1000 -
500 -
-500 -■1000 -■1500 -■2000 --2500
Рисунок 4. Решение уравнения (6) методом Рунге-Кутта модификации Дормана-Принса 4-5 порядка
Метод РКФ45 показал результаты, приведенные на рис. 5. Количество успешных шагов: 164. Количество возвратов: 13. Количество расчетов системы: 1056. Значения на конце интервала: -655.9. Время расчета: 0.029.
В данном примере время расчета методом РКФ45 также существенно меньше, чем методом РК-ДП45. Но также интересует точность расчетов. Аналитическое решение со значение аргумента, равного 50, равно у = -655.9371.
Рисунок 5. Решение уравнения (6) методом Рунге-Кутта Фельберга 4-5 порядка
Таким образом, мы смогли убедиться не только в высокой производительности метода РКФ45, но также в его высокой точности.
Следующим примером продемонстрируем эффективность интегрирования жесткой системы уравнений. Покажем это на примере уравнения осциллятора Ван-Дер-Поля. Уравнение Ван-Дер-Поля имеет следующий вид [11]:
(2х ¡л 2\(х --ц,(1 - х ) — + х = 0 .
(2
(7)
(8)
Уравнение (7) преобразуется в систему уравнений
| х' = у;
| у = ц-(1 - х2 )• у - х-
В уравнении (8) ц — коэффициент, характеризующий нелинейность и силу затухания колебаний, х — координата точки, зависящая от времени (.
Решим систему уравнений (8) с коэффициентом ц = 5, интервал интегрирования: от 0 до 30; И = 0.0001; х0 = 0.5; у0 = -0.5.
Решение данной системы методом РК-ДП45 покано на рис. 6. Статистические показатели этого решения следующие. Количество успешных шагов: 168. Количество возвратов: 32. Количество расчетов системы: 1201. Значения на конце интервала: 1.9, -0.1. Время расчета: 0.77 с.
Для метода РКФ45 будем использовать следующую максимальную ошибку: Ещи = 0.000001. Оценим то, что нам покажет метод РКФ45 (рис. 7). Количество успешных шагов: 548. Количество возвратов: 68. Общее количество расчетов системы дифуров: 3690. Значения на конце интервала: 1.9-0.1. Время расчета: 0.127.
Таким образом, система рассчитана с теми же значениями на конце интервала, что в данном случае позволяет говорить о приблизительно такой же точности рас-
чета. В то же время скорость расчет значительно выше, чем у метода РК-ПД45. Разница во времени расчете составляет 0.643 мс.
Рунге-Кутта модификации Дормана-Принса 4-5 порядка
Рисунок 7. Решение системы дифференциальных уравнений Ван-Дер-Поля (8) методом Рунге-Кутта-Фелъберга 4-5 порядка
Проведенные исследования позволяют говорить о возможности качественно проводить интегрирование систем дифференциальных уравнений, в том числе с ограниченной жесткостью. Заметим, что в процессе интегрирования активно используется ограничение максимальной ошибки с целью варьирования точно-
сти/времени интегрирования. Данный «рычаг» позволяет нам выбирать наиболее подходящее соотношение данных показателей.
Литература
[1] Самарский А. А. Введение в численные методы : учебник. — СПб. : Лань, 2009.
[2] Гайдышев И. П. Решение научных и инженерных задач средствами Excel, VBA и C/C++. — СПб. : БХВ-Петербург, 2004.
[3] Корягин С. В. Перекодировщик языков непрерывного моделирования // Промышленные АСУ и контроллеры. 2013. № 10. С. 52-56.
[4] Cellier F. E., Kofman E. Continuous system simulation. — Springer, 2006.
[5] Birta Louis G., Gilbert A. Modelling and Simulation: Exploring Dynamic System Behaviour. — Ottawa : School of information technology and engineering, 2007.
[6] Корягин С. В., Яковлев А. А. Проблемно-ориентированный язык системы непрерывного моделирования // Промышленные АСУ и Контроллеры. 2014. № 6. С. 36-40.
[7] Галкин А. В., Дятчина Д. В. Численное решение математических моделей объектов, заданных составными системами дифференциальных уравнений // Современные проблемы науки и образования. 2011. № 6. (http://www.science-education.ru/100-5196)
[8] Емельянов Н. В. Практическая небесная механика. Спецкурс ГАИШ МГУ [Электронный ресурс]. Режим доступа: URL: http://www.sai.msu.ru/neb/pcm/pcm.htm, свободный.
[9] Harder D. W. 4th-order Runge Kutta and the Dormand-Prince Methods. — Ontario : M.Math. LEL Department of Electrical and Computer Engineering University of Waterloo, 2012.
[10] Даутов Р. З. Практикум по методам решения задачи Коши для систем ОДУ. — Казань : Казанский государственный университет, 2010.
[11] Эдвардс Ч. Г., Пенни Д. Э. Дифференциальные уравнения и краевые задачи: моделирование и вычисление с помощью Mathematica, Maple и MATLAB. — 3-е изд. — М. : Ви-льямс, 2008.
Авторы:
Сергей Викторович Корягин — кандидат технических наук, доцент, доцент кафедры управления и моделирования систем, Московский технологический университет МИРЭА Андрей Александрович Яковлев — аспирант кафедры управления и моделирования систем, Московский технологический университет МИРЭА
Comparative Analysis of the Integration Methods with a Floating Step
S. V. Koriagin, A. A. Yakovlev
Moscow Technological University MIREA 78, Prospect Vernadskogo, Moscow, Russia, 119571
e-mail: dongenealog2003@mail.ru
Abstract. In this article we will analyze the accuracy and compare two popular pair methods of numerical integration Runge-Kutta-Fehlberg of the 4-5 order and Runge-Kutta Dormand-Prince 4-5 order with automatic change of the integration step, implemented in Matlab, a popular environment for solving problems of technical computing. We will use external compiler, that automatically generates the necessary statistics.
Key words: numerical integration, Runge-Kutta-Fehlberg of the 4-5 order method, Dormand-Prince modification.
Referenses
[1] SamarskijA. A. (2009) Vvedenie v chislennye metody. Lan'. [In Rus]
[2] Gajdyshev I. P. (2004) Reshenie nauchnyh i inzhenemyh zadach sredstvami Excel, VBA i C/C++. BHV-Peterburg. [In Rus]
[3] Koriagin S. V. (2013) Promyshlennye ASUi kontrollery, 10:52-56. [In Rus]
[4] Cellier F. E., Kofman E. (2006) Continuous system simulation. Springer.
[5] Birta Louis G., Gilbert A. (2007) Modelling and Simulation: Exploring Dynamic System Behaviour. School of information technology and engineering.
[6] Koriagin S. V., Yakovlev A. A. (2014) Promyshlennye ASU i Kontrollery, 6:36-40.
[7] Galkin A. V., Djatchina D. V.(2011) Sovremennye problemy nauki i obrazovanija, 6. (http://www.science-education.ru/100-5196)
[8] Emelyanov N. V. Prakticheskaja nebesnaja mehanika. Speckurs GAISh MGU http://www.sai.msu.ru/neb/pcm/pcm.htm
[9] HarderD. W. (2012) 4th-order Runge Kutta and the Dormand-Prince Methods. University of Waterloo.
[10] Dautov R. Z. (2010) Praktikum po metodam reshenija zadachi koshi dlja sistem ODU. Kazanskij gosu-darstvennyj universitet.
[11] Edwards C. H., Penney D. E. (2008) Elementary Differential Equations: 7th eds. Pearson Education, Inc
