CC BY
40
Cloud of Science. 2017. T. 4. № 1 http:/ / cloudofscience.ru
Сравнительный анализ методов интегрирования обыкновенных дифференциальных уравнений
С. В. Корягин, И. С. Егорушкин
Московский технологический университет (МИРЭА) 119571, Москва, Вернадского пр-т, 78
e-mail: dongenealog2003@mail. ru
Аннотация. Данная статья посвящена реализации проблемно-ориентированного транслирующего средства в рамках систем непрерывного моделирования. Построение переходных процессов осуществляется по формальному описанию объекта, предлагаемому в виде набора обыкновенных дифференциальных и алгебраических уравнений. Предусмотрена возможность изменения методов и условий интегрирования.
Ключевые слова: непрерывное моделирование, методы трансляции, проблемно-ориентированные языки программирования.
1. Введение
При изучении сложных объектов и систем иногда невозможно провести эксперимент в лабораторных или натуральных условиях. С целью изучения таких объектов и систем возникло понятие имитационного моделирования [1].
В процессе моделирования исследуемого объекта на компьютере его модель воспроизводится таким образом, что имеется возможность управлять процессом имитации. На основе анализа полученных результатов делаются выводы о поведении модели и ее свойствах.
Процесс функционирования модели обычно формализуется в виде математической модели — описания объекта математическими средствами, позволяющего делать выводы об интересующих нас свойствах объекта при помощи формальных процедур. В большинстве случаев математические модели описаны в виде обыкновенных дифференциальных и/или алгебраических уравнений.
В данной статье представлено два варианта формального описания исследуемых объектов в виде наборов дифференциальных и алгебраических уравнений. При помощи разработанного транслирующего средства осуществляется получение переходных процессов с возможностью варьирования методов и условий интегрирования [2-5, 10].
2. Метод Эйлера
Метод является самым простым и популярным методом интегрирования, но не самым точным [11]. Согласно методу Эйлера, какждый следующий элемент вычисляется по формуле
У+1 = У + ¥ (X; Уi )• Более точное значение можно найти, разложив ряд Тейлора:
^ кк (к) , ' к2 . У,+1 = ^Т7У() = У, + кУг +V У" + ••• к =0 к ! 2
Графическое представление метода Эйлера представлено на рис 1. С целью увеличения точности вычисления необходимо уменьшить значение шага к.
птближенное .з начелие\
} /
4— юсгг ь
\
>ефу т реальное значег, иеу
»'Л VI
1
XI X 2 9
Рисунок 1. Метод Эйлера
3. Метод Рунге-Кутта 2-го порядка
Метод точнее, чем метод Эйлера, но вычисляется дольше из-за необходимости подсчета функции дважды [11]. Точность обеспечивается тем, что при вычислении следующего значения функции берется половина шага к/2 :
к1 = к х / ( х,; у, ); У+1 = У, + кх/{хг + к/2; уг + к/2 ). (3)
Графическое представление метода выглядит следующим образом (рис. 2).
4. Метод Рунге-Кутта 4-го порядка
С целью увеличения точности вычисления необходимо каждый раз уменьшать размер шага, разбивая его на более мелкие части [6]. Но с увеличением точности уве-
личивается и количество вычислений. Метод Рунге-Кутта 4-го порядка заключается в разбиении шага на 4 части:
к1 = к х / ( х^; у); к2 = к х / (хг + к/2; У + ^/2 ); кз= к х / (хг + к/2; yi + к^2 ); к4 = к х / (хг + к; у + к3);
кл кл кл к
Уг+1 = Уг + + + + • '+1 г 6 3 3 6
На графике приближение выглядит следующим образом (рис. 3).
гзр иб/1 'жен ное !нач ?ние V
/ 1 по Эйяе РУ
/ У пог оешъ ост
/1 )йле ? а
шчу
■ч
ре альъ оез юче шеу
нач ние фуш ции
(л)-
к
л Г1 £ X 2 Г
Рисунок 2. Метод Рунге-Кутта 2-го порядка
Рисунок 3. Метод Рунге-Кутта 4-го порядка
5. Метод Рунге-Кутта-Фельберга 4-5 порядка
При вычислении переходных процессов, важен момент перехода системы из 1 состояния в другое. Т. к. он может быть непредсказуем, то существуют методы интегрирования с динамическим шагом. Изначально шаг интегрирования мал для увеличения точности, а впоследствии увеличивается [7]. Нахождение следующего значения функции происходит следующим образом:
к/ = к х / (х; у);
к2 = к х / ^ х+к; Уi+к/ ^;
. . / 3к 3к/ 9к2 ]
к3 = кх/I х1 +—;у1 + + —^ ; 3 Ч 8 32 32 /
, У 12к 1932к 7200к 7296к ^
к4 = кх/I х,. +-;у1 +-1--2 +-3
4 ^ ' 13 2197 2197 2197
439к „, 3680к3 845к4
к = к х /1 х + к; у +-1 - 8к +
1 ' 216 513 4104
, , / к 8к 3544к 1859к, 11к
к = к х /I х +—;у.--- + 2к--3 +-4--5
^ ' 2 27 2565 4104 40
16к 6656к 28561к 9к 2к
У,-+1 = У, + 7ТТ + ■
135 12825 56430 50 55 6. Проблемно-ориентированный язык
Формальное описание проблемно-ориентированного языка имеет вид [8, 10]: Программа <имя программы>; Задано
<имя /-го дифференциального уравнения>=<дифференциальное уравнение>; Коэффициенты
<имя коэффициента>=<коэффициент>; Начальные данные
Range=[<начало интегрирования>, <конец интегрирования>]; й=<шаг интегрирования^ х=[(<начальные условия>)]; Метод
<метод интегрирования^ Конец
Структура, описанная выше, имеет следующие ключевые слова:
Программа — начало программы, где описывается имя программы;
Задано — блок, где описывается система дифференциальных уравнений;
Коэффициенты — блок, где описываются коэффициенты;
Начальные данные — блок, где описываются промежуток интегрирования, шаг и задача Коши (начальные условия);
Метод — блок в котором отмечаются числовые методы, по которым будет производиться интегрирование;
Конец — конец программы.
Внутри блоков каждая строка отделяется символом «;». Левая и правая части отделяются знаком «=». Для символов используется латинский алфавит. Для чисел используются как целые, так и вещественные числа с разделителем в виде точки «.».
Разработано транслирующее средство, позволяющее решать математические модели следующими методами: Эйлера, Рунге-Кутта 2-го порядка, Рунге-Кутта 4-го порядка и методом Рунге-Кутта-Фельберга. Воспользуемся разработанным транслятором для решения системы дифференциальных уравнений, применяемых в задачах расчета переходных процессов.
Уравнение Колмогорова, описывающее вероятность состояния системы для марковского процесса через промежуток времени используем в качестве эксперимента. Рассмотрим систему с 2 уязвимостями и 4 состояниями: Ц — начальное состояние системы; — выявлена первая уязвимость, но не устранена; Б3 — выявлена, но не устранена вторая уязвимость; — обе уязвимости выявлены, но не устранены. Коэффициенты I и 12 отображают вероятность нахождения 1-й и 2-й уязвимости соответственно. т1 и т2 показывают вероятность устранения уязвимости.
Система уравнения имеет следующий вид:
&Р(Ц)
щР ( ^2 )+ Щ2 Р ( )-(1 + 12 ) Р ( Р1 ) = 1Р ( Ц) + т2 Р ( Б4 )-(/2 + щ) Р ( ?2 ) = 12Р ( Ц ) + т,Р ( Б4 )-(/1 + т2) Р ( Б3 ) = 12Р ( ?2) + 1Р ( Б3)-( т1 + т2) Р ( Б4 ) =
&
&Р (Ц?)
& '
ЖР (¥3)
ж '
&Р(?А)
& '
Пример 1.
Для решения имеем следующие данные:
Ц = 1;
F2 = 0.3; F3 = 0.4; F4 = 0.2; ^ = 0.07; /2 = 0.09; m = 0.1;
m = 0.2.
Шаг интегрирования = 0.1, интервал интегрирования: от 0 до 30. Решим задачу методом Эйлера. программа Вероятность_Состояний; задано
dxdt (1) = m1 x x(2) + m2 x x(3) - (/1 + / 2) x x(1); dxdt (2) = /1 x x(1) + m2 x x(4) - (/ 2 + m1) x x(2); dxdt (3) = / 2 x x(1) + m1 x x(4) - (/1 + m2) x x(3); dxdt (4) = / 2 x x(2) + /1 x x(3) - (m1 + m2) x x(4);
коэффициенты /1 = 0.07; / 2 = 0.09; m1 = 0.1; m2 = 0.2;
начальные данные range = [0,30]; h = 0.1;
x = [1,0.3,0.4,0.2];
метод
euler;
конец
В результате отработки транслирующего средства получаем следующие графики, представленные на рис. 4. Из графика видно, что при низких значениях интенсивности обнаружения уязвимостей и при высокой интенсивности их обнаружения, вероятность, что система не изменит своего изначального состояния, 78%. При этом шанс обнаружения хотя бы одной уязвимости 54%, а второй или обеих — менее 50%.
Рисунок 4. Решение задачи определения состояний системы методом Эйлера (пример 1)
Пример 2.
Внесем изменения:
Б2= 0.5;
Бз= 0.2;
Р4= 0.3; 1 = 0.17; /2 = 0.39; т = 0.4; т2 = 0.5.
Шаг интегрирования = 0.2, интервал интегрирования: от 0 до 10. Используем метод Рунге-Кутта 4-го порядка.
В результате отработки транслирующего средства получаем следующие графики, представленные на рис. 5.
При увеличении интенсивности обнаружения почти в 3 раза шанс системы остаться неизменной составляет 68%, что всего лишь на 10% меньше, чем в предыдущем примере. Но вероятность обнаружения уязвимостей теперь не менее 35%.
0.6
I
0,5
0,4
0.3
1,5 3
7,5
Рисунок 5. Решение задачи определения состояний системы методом Рунге-Кутта 4-го порядка (пример 2)
В качестве второго эксперимента рассмотрим систему дифференциальных уравнений, описывающую конкуренцию двух популяций. а и в — коэффициенты
Пример 3.
Для решения имеем следующие данные: N(0) = 10;
М(0) = 50; а = 0.003; Ь = 0.004; п = 0.5; т = 0.3; пх = 0.07; тх = 0.05.
Шаг интегрирования = 0.05, интервал интегрирования: от 0 до 10. Используем метод Рунге-Кутта-Фельберга.
программа Конкуренция_популяций;
задано
dxdt(l) = x(l) + a x x(l) x(l - x(l)/nz - m x (x(2)/mz)); dxdt(2) = x(2) + b x x(2) x (l - x(2)/mz - n x (x(l)/nz)).
коэффициенты a = 0.003; b = 0.004; n = 0.5; m = 0.3; nz = 0.07; mz = 0.05; начальные данные range = [0,10];
h = 0.05; x = [10,50]; метод rkf4; конец
В результате отработки транслирующего средства получаем графики, представленные на рис. 6. При заданных начальных условиях на графике видно, что при низкой рождаемости, средней терпимости, но высокой емкостью среды — первая популяция обгонит вторую по численности, несмотря на то, что изначально была меньше в 5 раз.
Рисунок 6. Решение системы дифференциальных уравнений конкуренции популяций методом Рунге-Кутта-Фельберга (пример 3)
Пример 4.
N(0) = 100; М(0) = 70; а = 0.006; Ь = 0.005; п = 0.7; т = 0.2; пх = 0.03; тх = 0.06;
Шаг интегрирования равен 0.01, интервал интегрирования: от 0 до 5. Используем метод Рунге-Кутта 2-го порядка.
В результате отработки транслирующего средства получаем следующие графики, представленные на рис. 7. Можно заметить рост второй популяции, в первое время обусловленный терпимостью первой, но из-за низкой рождаемости и емкости среды дальнейший рост невозможен.
Рисунок 7. Решение системы дифференциальных уравнений конкуренции популяций методом Рунге-Кутта 2-го порядка (пример 4)
7. Заключение
При увеличении точности метода интегрирования увеличивается и скорость выполнения алгоритма [9]. Поэтому при моделировании процесса необходимо выбирать что важнее: точность или время. Метод Рунге-Кутта 4-го порядка позволяет найти компромисс, т. к. он не занимает много времени как методы более высокого порядка и более точен, чем методы низкого порядка.
Литература
[1] Коробов П. Н. Математическое программирование и моделирование экономических процессов. — СПб. : Изд-во ДНК, 2006.
[2] Корягин С. В., Ильмухин С. В., Цыпленков С. В. Система непрерывного моделирования // ВестникМГУПИ. 2013. № 47. С. 66-77.
[3] Ахо А., Сети Р., Ульман Д. Компиляторы: принципы, технология и инструменты. М., 2011.
[4] Ахо А. В., Хопкрофт Дж. Э., Ульман Дж. Д. Структуры данных и алгоритмы. — М. : Издательский дом «Вильямс», 2008.
[5] Свердлов С. З. Языки программирования и методы трансляции. — СПб. : Питер, 2009.
[6] 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.
[7] Балдин К. В., Башлыков В. Н., Рукосуев А. В. Основы теории вероятностей и математической статистики: учебник. — М. : Флинта: НОУ ВПО «МПСИ», 2010.
[8] Cellier F. E. Continuous System Simulation. — Springer, 2006.
[9] Birta Louis G., Gilbert A. Modelling and Simulation: Exploring Dynamic System Behaviour. — Ottawa : School of information technology and engineering, 2007.
[10] Корягин С. В., Яковлев А. А. Проблемно-ориентированный язык системы непрерывного моделирования // Вестник МГУПИ. 2014. № 6. С. 36-40.
[11] Корягин С. В., Яковлев А. А. Сравнительный анализ методов интегрирования с плавающим шагом // Cloud of Science. 2016. Т. 3. № 1. С. 95-105.
[12] Gunz M. Differential Equations and Numerical Approximation using C++, 2007.
Авторы:
Сергей Викторович Корягин — кандидат технических наук, доцент, доцент кафедры управления и моделирования систем, Московский технологический университет (МИРЭА)
Илья Сергеевич Егорушкин — студент 1-го курса магистратуры, Московский технологический университет (МИРЭА)
Comparative analysis of integration methods
S. V. Koriagin, I. S. Egorushkin
Moscow technological university MIREA 78, Prospect Vernadskogo, Moscow, Russia, 119571
e-mail: [email protected]
Abstract. This article discusses the options for implementing problem-oriented broadcasting media in the framework of the continuous simulation. Building transition process follow according formal description of the object, which is a set of ordinary differential and algebraic equations. Provide for possibility to change the integration methods and conditions for integration. Key words: Continuous simulation, translation methods, problem-oriented programming languages.
References
[1] Korobov P. N. (2006) Mathematical programming and simulation of economic processes. DNK [In Rus]
[2] Koriagin S. V., Il'mukhin S. V., Tsyplenkov S. V. (2013) VestnikMGUPI, 47: 66-77. [In Rus]
[3] Aho A., Sethi R., Ulman J. (2011) Compilyatory: principy, tekhnologiya i instrumenty. Williams [In Rus]
[4] Aho A., Hopcroft J., Ulman J. (2008) Struktury dannyh i algoritmy. Williams [In Rus]
[5] Sverdlov S. Z. (2009) Yazyki programmirovaniya i metody translyacii. Piter [In Rus]
[6] Harder D. W. (2012) 4th-order Runge Kutta and the Dormand-Prince Methods. University of Waterloo.
[7] Baldin K. V., Bashlykov V.N., Rukosuev A. V. (2010) Fundamentals of the theory of probability and mathematical statistics. Flinta [In Rus]
[8] Cellier F. E. (2006) Continuous system simulation. Springer.
[9] Birta Louis G., Gilbert A. (2007) Modelling and Simulation: Exploring Dynamic System Behaviour. School of information technology and engineering.
[10] Koriagin S. V., Yakovlev A. A. (2014) Vestnik MGUPI, 4:36-40. [In Rus]
[11] Koriagin S. V., Yakovlev A. A. (2016) Cloud of Science, 3(1):95-105. [In Rus]
[12] Gunz. M. (2007) Differential Equations and Numerical Approximation using C++.
